text:
"""SunPy utility functions"""
from __future__ import absolute_import
from . util import *

__all__ = []
__all__ += util.__all__

def system_info():
    """Prints system information.

    Prints information about the runtime environment that SunPy lives in.
    Information about the OS, architecture, Python, and all major dependencies
    are included.

    The goal of this function is to provide enough information for someone
    running SunPy code or replicating a bug to setup a comparible environment
    to that which was originally used.

    Author: `Keith Hughitt <keith.hughitt@nasa.gov>`
    """
    import platform
    import datetime
    import sunpy

    print("==========================================================")
    print(" SunPy Installation Information\n")
    print(" " + datetime.datetime.utcnow().strftime("%A, %d. %B %Y %I:%M%p UT"))
    print("==========================================================\n")

    system = platform.system()
    proc = platform.processor()

    print("###########")
    print(" General")
    print("###########")

    # OS and architecture information
    if system == "Linux":
        distro = " ".join(platform.linux_distribution())
        print("OS: %s (Linux %s %s)" %  (distro, platform.release(), proc))
    elif system == "Darwin":
        print("OS: Mac OS X %s (%s)" %  (platform.mac_ver()[0], proc))
    elif system == "Windows":
        print("OS: Windows %s %s (%s)" %  (platform.release(),
                                        platform.version(), proc))
    else:
        print ("Unknown OS (%s)" % proc)

    # Python version
    arch = platform.architecture()[0]
    print("Python: %s (%s)\n" % (platform.python_version(), arch))

    # Dependencies
    try:
        from numpy import __version__ as numpy_version
    except ImportError:
        numpy_version = "NOT INSTALLED"

    try:
        from scipy import __version__ as scipy_version
    except ImportError:
        scipy_version = "NOT INSTALLED"

    try:
        from matplotlib import __version__ as matplotlib_version
    except ImportError:
        matplotlib_version = "NOT INSTALLED"

    try:
        from pyfits import __version__ as pyfits_version
    except ImportError:
        pyfits_version = "NOT INSTALLED"

    try:
        from suds import __version__ as suds_version
    except ImportError:
        suds_version = "NOT INSTALLED"

    try:
        from pandas import __version__ as pandas_version
    except ImportError:
        pandas_version = "NOT INSTALLED"

    try:
        from PyQt4.QtCore import PYQT_VERSION_STR as pyqt_version
    except ImportError:
        pyqt_version = "NOT INSTALLED"

    print("###########")
    print(" Libraries")
    print("###########")

    print("SunPy: %s" % sunpy.__version__)
    print("NumPy: %s" % numpy_version)
    print("SciPy: %s" % scipy_version)
    print("Matplotlib: %s" % matplotlib_version)
    print("PyFITS: %s" % pyfits_version)
    print("pandas: %s" % pandas_version)
    print("SUDS: %s" % suds_version)
    print("PyQt: %s\n" % pyqt_version)

max_stars_repo_name: astrofrog/sunpy

max_issues_repo_name: astrofrog/sunpy

max_forks_repo_name: astrofrog/sunpy

reason: from numpy,from scipy


text:
from __future__ import print_function, division

import torch
import torch.optim as optim
from torch.optim import lr_scheduler
from torchvision import transforms
from torch.utils.data import sampler

import time
import copy
import torch.nn.functional as F  # useful stateless functions
from torch.utils.data import DataLoader

from ssdd.dataset import MaskDataset
from ssdd.model import DecisionNet, SegmentationNet

torch.backends.cudnn.benchmark = True

img_transform = transforms.Compose([
    transforms.ToTensor(),
])

msk_transform = transforms.Compose([
    transforms.ToTensor(),
])

class SequentialSampler(sampler.Sampler):
    r"""Samples elements sequentially, always in the same order.

    Arguments:
        data_source (Dataset): dataset to sample from
    """

    def __init__(self, data_source):
        self.data_source = data_source

    def __iter__(self):
        return iter(self.data_source)

    def __len__(self):
        return len(self.data_source)

print("scale: ", 2)
#   img_resize=(128 * 2, 128 * 2),
#   msk_resize=(16 * 2, 16 * 2),
dataset = MaskDataset("./pre_data",
                      img_resize=(128 * 2, 128 * 2),
                      msk_resize=(128, 128),
                      img_transform=img_transform,
                      msk_transform=msk_transform)

ALL_TOTAL = 400
num_train = 250
num_val = 75
num_test = 75

dataloaders = dict()
dataloaders['train'] = DataLoader(dataset, batch_size=1, shuffle=False, num_workers=5,
                                  sampler=SequentialSampler(range(num_train)))

dataloaders['val'] = DataLoader(dataset, batch_size=1, shuffle=False, num_workers=5,
                                sampler=SequentialSampler(range(num_train, num_train + num_val)))

dataloaders['test'] = DataLoader(dataset, batch_size=1, shuffle=False, num_workers=5,
                                sampler=SequentialSampler(range(num_train + num_val, num_train + num_val + num_test)))
dataset_sizes = {
    'train': num_train,
    'val': num_val,
    'test': num_test
}

print(dataset_sizes)
device = torch.device("cuda:0" if torch.cuda.is_available() else "cpu")
print(device)

def train_model_seg(seg_model, optimizer, scheduler, num_epochs=25):
    since = time.time()
    train_loss_list = []
    val_loss_list = []
    for epoch in range(num_epochs):
        print('Epoch {}/{}'.format(epoch, num_epochs - 1))
        print('-' * 10)

        # Each epoch has a training and validation phase
        for phase in ['train', 'val']:
            if phase == 'train':
                seg_model.train()  # Set model to training mode
            else:
                pass
                # model.eval()  # Set model to evaluate mode

            running_loss = 0.0

            # Iterate over data.
            filename = 0
            dataset_iter = dataloaders[phase].__iter__()
            for _ in range(len(dataloaders[phase])):
                inputs, masks, labels = dataset_iter.__next__()

                inputs = inputs.cuda()
                masks = masks.cuda()
                # labels = labels.cuda()

                # zero the parameter gradients
                optimizer.zero_grad()
                # forward
                # track history if only in train
                with torch.set_grad_enabled(phase == 'train'):
                    feat, pred_mask = seg_model.forward(inputs)
                    # pred_mask = torch.sigmoid(pred_mask)
                    loss = F.binary_cross_entropy(pred_mask.reshape(1, -1), masks.reshape(1, -1))
                    # loss = F.mse_loss(pred_mask.reshape(1, -1), masks.reshape(1, -1))
                    if phase == 'train':
                        loss.backward()
                        optimizer.step()
                    if epoch % 3 == 0:
                        with torch.no_grad():
                            save_dir = "./visualization/{}_epoch-{}".format(phase, epoch)
                            visualization(save_dir, str(filename) + ".jpg", inputs[0][0], masks[0][0], pred_mask[0][0])
                            filename += 1
                # statistics
                running_loss += loss.item() * inputs.size(0)

            if phase == 'train':
                scheduler.step()

            epoch_loss = running_loss / dataset_sizes[phase]
            if phase == 'train':
                train_loss_list.append(epoch_loss)
            else:
                val_loss_list.append(epoch_loss)

            print('{} Loss: {:.4f}'.format(phase, epoch_loss))

        print()

    time_elapsed = time.time() - since
    print('Training complete in {:.0f}m {:.0f}s\n'.format(
        time_elapsed // 60, time_elapsed % 60))

    print("train_loss_list =", train_loss_list)
    print("val_loss_list =", val_loss_list)
    # load best model weights
    return seg_model

def visualization(save_dir, filename, img, mask, pred_mask):
    import os
    import numpy as np
    if not os.path.exists(save_dir):
        os.makedirs(save_dir)
    img = img.cpu()
    mask = mask.cpu()
    pred_mask = pred_mask.cpu()
    image = np.array(img) * 255
    mask = np.array(mask) * 255
    pred_mask = np.array(pred_mask.detach().numpy()) * 255
    img_visual = concat_image([image, mask, pred_mask])
    visualization_path = os.path.join(save_dir, filename)
    img_visual.save(visualization_path)

def concat_image(images, mode="L"):
    from PIL import Image
    if not isinstance(images, list):
        raise Exception('images must be a list  ')
    count = len(images)
    size = Image.fromarray(images[0]).size
    target = Image.new(mode, (size[0] * count, size[1] * 1))
    for i in range(count):
        image = Image.fromarray(images[i]).resize(size, Image.BILINEAR)
        target.paste(image, (i*size[0], 0, (i+1)*size[0], size[1]))
    return target

def run_seg():
    seg_model = SegmentationNet(1)
    # seg_model = torch.nn.DataParallel(seg_model.cuda())
    seg_model = seg_model.cuda()

    # Observe that all parameters are being optimized
    optimizer = optim.Adam(seg_model.parameters(), lr=0.001)

    # Decay LR by a factor of 0.1 every 7 epochs
    exp_lr_scheduler = lr_scheduler.StepLR(optimizer, step_size=7, gamma=0.1)

    seg_model = train_model_seg(seg_model, optimizer, exp_lr_scheduler, num_epochs=7)
    import os
    save_dir = "./saved/"
    if not os.path.exists(save_dir):
        os.makedirs(save_dir)
    torch.save(seg_model, save_dir + 'segmentation')

def train_model_cls(seg_model, cls_model, optimizer, scheduler, num_epochs=25):
    train_acc_list = []
    val_acc_list = []
    train_loss_list = []
    val_loss_list = []
    since = time.time()
    best_model_wts = copy.deepcopy(cls_model.state_dict())
    best_acc = 0.0
    for epoch in range(num_epochs):
        print('Epoch {}/{}'.format(epoch, num_epochs - 1))
        print('-' * 10)

        # Each epoch has a training and validation phase
        for phase in ['train', 'val']:
            if phase == 'train':
                cls_model.train()  # Set model to training mode
            else:
                cls_model.train()  # Set model to training mode
                pass
                # model.eval()  # Set model to evaluate mode

            running_loss = 0.0
            running_corrects = torch.tensor(0)
            tp = 0
            tn = 0
            fp = 0
            fn = 0

            dataset_iter = dataloaders[phase].__iter__()
            for _ in range(len(dataset_iter)):
                inputs, masks, labels = dataset_iter.__next__()

                inputs = inputs.cuda()
                # masks = masks.cuda()
                labels = labels.cuda()
                # zero the parameter gradients
                optimizer.zero_grad()
                # forward
                with torch.no_grad():
                    feats, pred_masks = seg_model.forward(inputs)

                with torch.set_grad_enabled(phase == 'train'):
                    scores = cls_model.forward(feats, pred_masks)
                    cls_loss = F.cross_entropy(scores, labels)

                    _, preds = torch.max(scores, 1)

                    # loss = cls_loss + seg_loss
                    if phase == 'train':
                        cls_loss.backward()
                        optimizer.step()

                # statistics
                if preds[0] == labels[0]:
                    if preds[0] == 1:
                        tp += 1
                    else:
                        tn += 1
                else:
                    if preds[0] == 1:
                        fp += 1
                    else:
                        fn += 1

                running_loss += cls_loss.item() * inputs.size(0)
                running_corrects += torch.sum(preds == labels.data)
            if phase == 'train':
                scheduler.step()
            epoch_loss = running_loss / dataset_sizes[phase]
            epoch_acc = running_corrects.float() / dataset_sizes[phase]

            if phase == 'train':
                train_acc_list.append(epoch_acc.item())
                train_loss_list.append(epoch_loss)
            else:
                val_acc_list.append(epoch_acc.item())
                val_loss_list.append(epoch_loss)

            print('{} Loss: {:.4f} Acc: {:.4f}'.format(
                phase, epoch_loss, epoch_acc))

            print('tp: {} fn: {} fp: {} tn: {}'.format(tp, fn, fp, tn))
            # deep copy the model
            if phase == 'val' and epoch_acc > best_acc:
                best_acc = epoch_acc
                best_model_wts = copy.deepcopy(cls_model.state_dict())

        print()

    time_elapsed = time.time() - since
    print('Training complete in {:.0f}m {:.0f}s'.format(
        time_elapsed // 60, time_elapsed % 60))
    print('Best val Acc: {:4f}\n'.format(best_acc))

    print("train_loss_list =", train_loss_list)
    print("val_loss_list =", val_loss_list)
    print("train_acc_list =", train_acc_list)
    print("val_acc_list =", val_acc_list)
    print("time_elapsed =", time_elapsed)

    # load best model weights
    cls_model.load_state_dict(best_model_wts)
    return cls_model

def test_model(seg_model, cls_model):
    since = time.time()
    for phase in ['test']:
        running_corrects = torch.tensor(0)
        tp = 0
        tn = 0
        fp = 0
        fn = 0
        filename = 0

        dataset_iter = dataloaders[phase].__iter__()
        for _ in range(len(dataset_iter)):
            inputs, masks, labels = dataset_iter.__next__()

            inputs = inputs.cuda()
            # inputs = inputs.float()
            labels = labels.cuda()
            with torch.no_grad():
                feats, pred_masks = seg_model.forward(inputs)
                scores = cls_model.forward(feats, pred_masks)
                _, preds = torch.max(scores, 1)

                save_dir = "./visualization/{}".format(phase)
                visualization(save_dir, str(filename) + ".jpg", inputs[0][0], masks[0][0], pred_masks[0][0])
                filename += 1

            if preds[0] == labels[0]:
                if preds[0] == 1:
                    tp += 1
                else:
                    tn += 1
            else:
                if preds[0] == 1:
                    fp += 1
                else:
                    fn += 1

            running_corrects += torch.sum(preds == labels.data)
        epoch_acc = running_corrects.float() / dataset_sizes[phase]

        print('{} Acc: {:.4f}'.format(phase, epoch_acc))
        print('tp: {} fn: {} fp: {} tn: {}'.format(tp, fn, fp, tn))
        print()

    time_elapsed = time.time() - since
    print("time_elapsed =", time_elapsed, "\n")

def run_cls(save_dir):
    seg_model = torch.load(save_dir)
    seg_model = seg_model.cuda()
    cls_model = DecisionNet()
    cls_model = cls_model.cuda()

    # Observe that all parameters are being optimized
    # with regularization
    # optimizer = optim.Adam(cls_model.parameters(), lr=0.001, weight_decay=0.1)
    # without regularization
    optimizer = optim.Adam(cls_model.parameters(), lr=0.001)

    # Decay LR by a factor of 0.1 every 7 epochs
    exp_lr_scheduler = lr_scheduler.StepLR(optimizer, step_size=7, gamma=0.1)

    cls_model = train_model_cls(seg_model, cls_model, optimizer, exp_lr_scheduler, num_epochs=14)
    print("Test Network!\n")
    test_model(seg_model, cls_model)
    import os
    save_dir = "./saved/"
    if not os.path.exists(save_dir):
        os.makedirs(save_dir)
    torch.save(cls_model, save_dir + 'classification')

if __name__ == '__main__':
    import sys
    if len(sys.argv) == 1:
        print()
        print('Train Segmentation Network!\n')
        run_seg()
        print()
        print('Train Classification Network!\n')
        run_cls('./saved/segmentation')
    else:
        if sys.argv[1] == 'seg':
            run_seg()
        elif sys.argv[1] == 'cls':
            run_cls('./saved/segmentation')
        elif sys.argv[1] == 'all':
            run_seg()
            run_cls('./saved/segmentation')
        else:
            print('run `python train.py seg` or `python train.py cls`')

ext: py

max_stars_repo_name: OYMiss/ssdd-net

max_issues_repo_name: OYMiss/ssdd-net

max_forks_repo_name: OYMiss/ssdd-net

reason: import numpy


text:
"""
 Copyright (C) 2018-2021 Intel Corporation

 Licensed under the Apache License, Version 2.0 (the "License");
 you may not use this file except in compliance with the License.
 You may obtain a copy of the License at

      http://www.apache.org/licenses/LICENSE-2.0

 Unless required by applicable law or agreed to in writing, software
 distributed under the License is distributed on an "AS IS" BASIS,
 WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 See the License for the specific language governing permissions and
 limitations under the License.
"""

import unittest

import numpy as np
from generator import generator, generate

from mo.graph.graph import Node
from mo.ops.expand_dims import ExpandDims
from mo.utils.unittest.graph import build_graph

nodes_attributes = {
    'data_1': {
        'kind': 'data',
        'shape': np.array([2, 3, 224, 224]),
        'value': None,
    },
    'expand_dims': {
        'type': 'None',
        'kind': 'op',
    },
    'data_2': {
        'kind': 'data',
        'shape': None,
        'value': None,
    }
}

@generator
class ExpandDimsOp(unittest.TestCase):
    @generate(*[(0, [1, 2, 3, 224, 224]),
                (1, [2, 1, 3, 224, 224]),
                (2, [2, 3, 1, 224, 224]),
                (3, [2, 3, 224, 1, 224]),
                (4, [2, 3, 224, 224, 1]),
                ])
    def test_expand_dims_infer(self, axis, ref_out_shape):
        graph = build_graph(nodes_attributes,
                            [('data_1', 'expand_dims'),
                             ('expand_dims', 'data_2')],
                            {'expand_dims': {'expand_axis': axis}})
        expand_dims_node = Node(graph, 'expand_dims')

        ExpandDims.infer(expand_dims_node)

        self.assertTrue(np.array_equal(expand_dims_node.out_node().shape, np.array(ref_out_shape)))

@generator
class ExpandDimsOpValueInfer(unittest.TestCase):
    @generate(*[(0, [2, 3, 224, 224], [1, 2, 3, 224, 224]),
                (1, [2, 3, 224, 224], [2, 1, 3, 224, 224]),
                (2, [2, 3, 224, 224], [2, 3, 1, 224, 224]),
                (3, [2, 3, 224, 224], [2, 3, 224, 1, 224]),
                (4, [2, 3, 224, 224], [2, 3, 224, 224, 1]),
                ])
    def test_expand_dims_infer_value(self, axis, in_shape, ref_out_shape):
        in_value = np.random.rand(*in_shape)
        graph = build_graph(nodes_attributes,
                            [('data_1', 'expand_dims'),
                             ('expand_dims', 'data_2')],
                            {'data_1': {'value': in_value},
                             'expand_dims': {'expand_axis': axis}})
        expand_dims_node = Node(graph, 'expand_dims')

        ExpandDims.infer(expand_dims_node)

        self.assertTrue(np.array_equal(expand_dims_node.out_node().shape, np.array(ref_out_shape)))
        self.assertTrue(np.array_equal(expand_dims_node.out_node().value, np.array(in_value.reshape(ref_out_shape))))

max_stars_repo_name: calvinfeng/openvino

max_issues_repo_name: calvinfeng/openvino

max_forks_repo_name: calvinfeng/openvino


text:
# ------------------------------------------------------------------------------
# Loads Cityscapes semantic dataset.
# Written by Bowen Cheng (bcheng9@illinois.edu)
# ------------------------------------------------------------------------------

import glob
import os

import numpy as np

from .base_dataset import BaseDataset
from .utils import DatasetDescriptor
from ..transforms import build_transforms

_CITYSCAPES_INFORMATION = DatasetDescriptor(
    splits_to_sizes={'train': 2975,
                     'trainval': 3475,
                     'val': 500,
                     'test': 1525},
    num_classes=19,
    ignore_label=255,
)

_CITYSCAPES_TRAIN_ID_TO_EVAL_ID = [7, 8, 11, 12, 13, 17, 19, 20, 21, 22,
                                   23, 24, 25, 26, 27, 28, 31, 32, 33]

# A map from data type to folder name that saves the data.
_FOLDERS_MAP = {
    'image': 'leftImg8bit',
    'label': 'gtFine',
}

# A map from data type to filename postfix.
_POSTFIX_MAP = {
    'image': '_leftImg8bit',
    'label': '_gtFine_labelTrainIds',
}

# A map from data type to data format.
_DATA_FORMAT_MAP = {
    'image': 'png',
    'label': 'png',
}

class Cityscapes(BaseDataset):
    """
    Cityscapes semantic segmentation dataset.
    Arguments:
        root: Str, root directory.
        split: Str, data split, e.g. train/val/test.
        is_train: Bool, for training or testing.
        crop_size: Tuple, crop size.
        mirror: Bool, whether to apply random horizontal flip.
        min_scale: Float, min scale in scale augmentation.
        max_scale: Float, max scale in scale augmentation.
        scale_step_size: Float, step size to select random scale.
        mean: Tuple, image mean.
        std: Tuple, image std.
    """
    def __init__(self,
                 root,
                 split,
                 is_train=True,
                 crop_size=(513, 1025),
                 mirror=True,
                 min_scale=0.5,
                 max_scale=2.,
                 scale_step_size=0.25,
                 mean=(0.485, 0.456, 0.406),
                 std=(0.229, 0.224, 0.225),
                 **kwargs):
        super(Cityscapes, self).__init__(root, split, is_train, crop_size, mirror, min_scale, max_scale,
                                         scale_step_size, mean, std)

        self.num_classes = _CITYSCAPES_INFORMATION.num_classes
        self.ignore_label = _CITYSCAPES_INFORMATION.ignore_label
        self.label_pad_value = (self.ignore_label, )

        # Get image and annotation list.
        self.img_list = self._get_files('image', self.split)
        self.ann_list = self._get_files('label', self.split)

        assert len(self) == _CITYSCAPES_INFORMATION.splits_to_sizes[self.split]

        self.transform = build_transforms(self, is_train)

    def _get_files(self, data, dataset_split):
        """Gets files for the specified data type and dataset split.
        Args:
            data: String, desired data ('image' or 'label').
            dataset_split: String, dataset split ('train', 'val', 'test')
        Returns:
            A list of sorted file names or None when getting label for test set.
        """
        if data == 'label' and dataset_split == 'test':
            return None
        pattern = '*%s.%s' % (_POSTFIX_MAP[data], _DATA_FORMAT_MAP[data])
        search_files = os.path.join(
            self.root, _FOLDERS_MAP[data], dataset_split, '*', pattern)
        filenames = glob.glob(search_files)
        return sorted(filenames)

    @staticmethod
    def train_id_to_eval_id():
        return _CITYSCAPES_TRAIN_ID_TO_EVAL_ID

    def _convert_train_id_to_eval_id(self, prediction):
        """Converts the predicted label for evaluation.
        There are cases where the training labels are not equal to the evaluation
        labels. This function is used to perform the conversion so that we could
        evaluate the results on the evaluation server.
        Args:
            prediction: Semantic segmentation prediction.
        Returns:
            Semantic segmentation prediction whose labels have been changed.
        """
        converted_prediction = prediction.copy()
        for train_id, eval_id in enumerate(self.train_id_to_eval_id()):
            converted_prediction[prediction == train_id] = eval_id

        return converted_prediction

    @staticmethod
    def create_label_colormap():
        """Creates a label colormap used in CITYSCAPES segmentation benchmark.
        Returns:
            A colormap for visualizing segmentation results.
        """
        colormap = np.zeros((256, 3), dtype=np.uint8)
        colormap[0] = [128, 64, 128]
        colormap[1] = [244, 35, 232]
        colormap[2] = [70, 70, 70]
        colormap[3] = [102, 102, 156]
        colormap[4] = [190, 153, 153]
        colormap[5] = [153, 153, 153]
        colormap[6] = [250, 170, 30]
        colormap[7] = [220, 220, 0]
        colormap[8] = [107, 142, 35]
        colormap[9] = [152, 251, 152]
        colormap[10] = [70, 130, 180]
        colormap[11] = [220, 20, 60]
        colormap[12] = [255, 0, 0]
        colormap[13] = [0, 0, 142]
        colormap[14] = [0, 0, 70]
        colormap[15] = [0, 60, 100]
        colormap[16] = [0, 80, 100]
        colormap[17] = [0, 0, 230]
        colormap[18] = [119, 11, 32]
        return colormap


text:
#######################################################################
# Copyright (C) 2017 Shangtong Zhang(zhangshangtong.cpp@gmail.com)    #
# Permission given to modify the code as long as you keep this        #
# declaration at the top                                              #
#######################################################################
from typing import List, Optional

import numpy as np
import torch
import torch.nn as nn
import torch.nn.functional as F
from regym.networks.utils import BaseNet, layer_init, tensor, entropy
from regym.networks.bodies import DummyBody

class CategoricalDQNet(nn.Module, BaseNet):

    def __init__(self,
                 body: nn.Module,
                 action_dim: int,
                 actfn=F.leaky_relu,
                 use_cuda=False):
        BaseNet.__init__(self)
        super(CategoricalDQNet, self).__init__()

        self.body = body
        self.action_dim = action_dim
        self.use_cuda = use_cuda

        self.actfn = actfn

        body_output_features = self.body.feature_dim
        self.qsa = nn.Linear(body_output_features, self.action_dim)

        if self.use_cuda:
            self = self.cuda()

    def forward(self, x: torch.Tensor, action: torch.Tensor = None,
                legal_actions: List[int] = None):
        # Forward pass till last layer
        x = self.body(x)
        # Q values for all actions
        q_values = self.qsa(x)

        # TODO: check if this works
        if legal_actions is not None:
            q_values = self._mask_ilegal_action_logits(q_values,
                                 legal_actions, self.action_dim)
        if action is None:
            q_value, action = q_values.max(dim=1)

        probs = F.softmax(q_values, dim=-1)
        return {'a': action,
                'Q': q_values,
                'entropy': entropy(probs)}

class CategoricalDuelingDQNet(nn.Module, BaseNet):

    def __init__(self,
                 body: nn.Module,
                 action_dim: int,
                 actfn=F.leaky_relu):
        BaseNet.__init__(self)
        super(CategoricalDuelingDQNet, self).__init__()
        self.action_dim = action_dim

        self.body = body

        self.value = layer_init(nn.Linear(body.feature_dim, 1))
        self.advantage = layer_init(nn.Linear(body.feature_dim, action_dim))

    def forward(self, x, action: torch.Tensor = None,
                legal_actions: List[int] = None):
        x = self.body(tensor(x))
        V = self.value(x)
        A = self.advantage(x)

        Q = V.expand_as(A) + (A - A.mean(1, keepdim=True))

        probs = F.softmax(Q, dim=-1)
        return {'V': V, 'A': A, 'Q': Q, 'a': Q.max(dim=1)[1],
                'entropy': entropy(probs)}

class CategoricalHead(nn.Module, BaseNet):
    def __init__(self, input_dim, output_dim, body: nn.Module = None):
        super().__init__()
        self.body = body
        self.fc_categorical = layer_init(nn.Linear(input_dim, output_dim))
        self.input_dim = input_dim
        self.output_dim = output_dim

    def forward(self, x: torch.Tensor, legal_actions: List[int] = None):
        body_embedding = self.body(x) if self.body else x
        logits = self.fc_categorical(body_embedding)

        if legal_actions:
            logits = self._mask_ilegal_action_logits(logits, legal_actions, self.output_dim)

        probs = F.softmax(logits, dim=-1)
        log_probs = probs.log()
        action = torch.distributions.Categorical(logits=logits).sample()
        action = action.view(-1, 1)
        return {'probs': probs, 'log_probs': log_probs, 'a': action,
                'entropy': entropy(probs)}

class QuantileNet(nn.Module, BaseNet):
    def __init__(self, action_dim, num_quantiles, body):
        super(QuantileNet, self).__init__()
        self.fc_quantiles = layer_init(nn.Linear(body.feature_dim, action_dim * num_quantiles))
        self.action_dim = action_dim
        self.num_quantiles = num_quantiles
        self.body = body
        self.to(Config.DEVICE)

    def forward(self, x):
        phi = self.body(tensor(x))
        quantiles = self.fc_quantiles(phi)
        quantiles = quantiles.view((-1, self.action_dim, self.num_quantiles))
        return quantiles

class ActorCriticNet(nn.Module):
    def __init__(self, state_dim, action_dim, body, actor_body, critic_body):
        super(ActorCriticNet, self).__init__()
        if body is None: body = DummyBody(state_dim)
        if actor_body is None: actor_body = DummyBody(body.feature_dim)
        if critic_body is None: critic_body = DummyBody(body.feature_dim)
        self.body = body
        self.actor_body = actor_body
        self.critic_body = critic_body
        self.fc_action = layer_init(nn.Linear(actor_body.feature_dim, action_dim), 1e-3)
        self.fc_critic = layer_init(nn.Linear(critic_body.feature_dim, 1), 1e-3)

        self.actor_params = list(self.actor_body.parameters()) + list(self.fc_action.parameters())
        self.critic_params = list(self.critic_body.parameters()) + list(self.fc_critic.parameters())
        self.phi_params = list(self.body.parameters())

class GaussianActorCriticNet(nn.Module, BaseNet):
    def __init__(self,
                 state_dim,
                 action_dim,
                 body=None,
                 actor_body=None,
                 critic_body=None):
        super(GaussianActorCriticNet, self).__init__()
        self.network = ActorCriticNet(state_dim, action_dim, body, actor_body, critic_body)
        self.std = nn.Parameter(torch.zeros(action_dim))

    def forward(self, obs, action=None, rnn_states=None):
        obs = tensor(obs)

        if rnn_states is not None and 'phi_arch' in rnn_states:
            phi, rnn_states['phi_arch'] = self.network.body( (obs, rnn_states['phi_arch']) )
        else:
            phi = self.network.body(obs)

        if rnn_states is not None and 'actor_arch' in rnn_states:
            phi_a, rnn_states['actor_arch'] = self.network.actor_body( (phi, rnn_states['actor_arch']) )
        else:
            phi_a = self.network.actor_body(phi)

        if rnn_states is not None and 'critic_arch' in rnn_states:
            phi_v, rnn_states['critic_arch'] = self.network.critic_body( (phi, rnn_states['critic_arch']) )
        else:
            phi_v = self.network.critic_body(phi)

        mean = F.tanh(self.network.fc_action(phi_a))
        v = self.network.fc_critic(phi_v)
        dist = torch.distributions.Normal(mean, F.softplus(self.std))
        if action is None:
            action = dist.sample()
        log_prob = dist.log_prob(action).sum(-1).unsqueeze(-1)
        entropy = dist.entropy().sum(-1).unsqueeze(-1)

        if rnn_states is not None:
            return {'a': action,
                    'log_pi_a': log_prob,
                    'ent': entropy,
                    'V': v,
                    'rnn_states': rnn_states}
        else:
            return {'a': action,
                    'log_pi_a': log_prob,
                    'ent': entropy,
                    'V': v}

class CategoricalActorCriticNet(nn.Module, BaseNet):
    def __init__(self,
                 state_dim: int,
                 action_dim: int,
                 critic_gate_fn: Optional[str] = None,
                 body: nn.Module = None,
                 actor_body: nn.Module = None,
                 critic_body: nn.Module = None):
        BaseNet.__init__(self)
        super(CategoricalActorCriticNet, self).__init__()
        self.action_dim = action_dim

        self.critic_gate_fn = critic_gate_fn

        self.network = ActorCriticNet(state_dim, action_dim, body, actor_body, critic_body)

    # TODO: type hint rnn_states
    def forward(self, obs: torch.Tensor, action: int = None, rnn_states=None,
                legal_actions: List[int] = None):
        if rnn_states is not None:
            next_rnn_states = {k: None for k in rnn_states}

        if rnn_states is not None and 'phi_arch' in rnn_states:
            phi, next_rnn_states['phi_arch'] = self.network.body( (obs, rnn_states['phi_arch']) )
        else:
            phi = self.network.body(obs)

        if rnn_states is not None and 'actor_arch' in rnn_states:
            phi_a, next_rnn_states['actor_arch'] = self.network.actor_body( (phi, rnn_states['actor_arch']) )
        else:
            phi_a = self.network.actor_body(phi)

        if rnn_states is not None and 'critic_arch' in rnn_states:
            phi_v, next_rnn_states['critic_arch'] = self.network.critic_body( (phi, rnn_states['critic_arch']) )
        else:
            phi_v = self.network.critic_body(phi)

        logits = self.network.fc_action(phi_a)
        if legal_actions:
            logits = self._mask_ilegal_action_logits(logits, legal_actions, self.action_dim)
        # Size: batch x action_dim
        v = self.network.fc_critic(phi_v)
        if self.critic_gate_fn: v = self.critic_gate_fn(v)
        # Size: batch x 1

        dist = torch.distributions.Categorical(logits=logits)
        if action is None:
            # Size: batch x 1
            action = dist.sample()
            action = action.view(-1, 1)
        log_prob = dist.log_prob(action).unsqueeze(-1)
        # estimates the log likelihood of each action against each batched distributions... : batch x batch x 1
        log_prob = torch.cat([log_prob[idx][idx].view((1,-1)) for idx in range(log_prob.size(0))], dim=0)
        # retrieve the log likelihood of each batched action against its relevant batched distribution: batch x 1
        entropy = dist.entropy().unsqueeze(-1)
        # retrieve the the entropy of each batched distribution: batch x 1

        prediction = {'a': action,
                      'log_pi_a': log_prob,
                      'entropy': entropy,
                      'V': v,
                      'probs': dist.probs}
        if rnn_states is not None: prediction.update({'rnn_states': rnn_states, 'next_rnn_states': next_rnn_states})
        return prediction

class PolicyInferenceActorCriticNet(nn.Module, BaseNet):
    '''
    Inspired from A3C-AMF architecture:
        'Agent Modeling as Auxiliary Task for Deep Reinforcement Learning'
        Pablo Hernandez-Leal et al. (arXiv:1907.09597)
    Originally suggested as DPIQN architecture in:
        'A Deep Policy Inference Q-Network for Multi-Agent Systems'
        Zhang-Wei Hong et al (ariv:1712.07893)
    '''
    def __init__(self,
                 num_policies: int,
                 num_actions: int,
                 feature_extractor: BaseNet,
                 policy_inference_body: BaseNet,
                 actor_critic_body: BaseNet):
        '''
        TODO
        '''
        BaseNet.__init__(self)
        super().__init__()

        self.num_policies = num_policies
        # NOTE: all three
        # (policy_inference_body feature_dim, actor_critic_body feature_dim, actor_critic_head input_dim)
        # should be the same!

        self.feature_extractor = feature_extractor

        self.policy_inference_body = policy_inference_body
        self.actor_critic_body = actor_critic_body

        self.policy_inference_heads = nn.ModuleList([
                CategoricalHead(
                    policy_inference_body.feature_dim, num_actions)
                for _ in range(num_policies)])

        self.actor_critic_head = CategoricalActorCriticNet(
                state_dim=self.actor_critic_body.feature_dim,
                action_dim=num_actions)

    def forward(self, x: torch.Tensor, legal_actions: List[int] = None):
        feature_embeddings = self.feature_extractor(x)

        # Opponent modelling branch
        policy_type_embedding = self.policy_inference_body(feature_embeddings)
        # How about adding a non-linearity here?
        policy_inference_predictions = map(lambda head: head(policy_type_embedding),
                                           self.policy_inference_heads)

        # Actor critic branch
        actor_critic_body_embedding = self.actor_critic_body(feature_embeddings)
        # element wise vector multiplication
        actor_critic_head_input = policy_type_embedding * actor_critic_body_embedding

        actor_critic_prediction = self.actor_critic_head(
                actor_critic_head_input, legal_actions=legal_actions)

        # TODO: think how to regain this.
        # Current non-commented out approach only uses
        # policy_inference_predictions_dict = {
        #         f'policy_{i}': prediction
        #         for i, prediction
        #         in zip(range(self.num_policies), policy_inference_predictions)}
        policy_inference_predictions_dict = {
                 f'policy_{i}': prediction['probs']
                 for i, prediction
                 in zip(range(self.num_policies), policy_inference_predictions)}

        # Merging 2 dicts together
        aggregated_predictions = {
                **policy_inference_predictions_dict, **actor_critic_prediction}

        return aggregated_predictions


text:
import pandas as pd
import numpy as np
import pytest
from helper_functions_test import Wrangle, Randomize

### Instantiate Classes
wrangle_methods = Wrangle()
rand_methods = Randomize()

### Create DF for testing
df = pd.DataFrame(
    np.random.randint(0, 100, size=(15, 4)),
    columns=list('ABCD')
)

def test_nc():
    """
    Test Null_Count Method from Wrangle Class
    """
    assert wrangle_methods.null_count(df) == 0

    ### Compare null count to ensure proper accuracy of function.
    assert wrangle_methods.null_count(df) == df.isnull().sum().sum()

def test_rand():
    """
    Test Randomize Method from Randomize Class
    """
    assert len(rand_methods.randomize(df, seed=42)) == len(df)


text:
import pytest
import numpy as np

@pytest.fixture(scope="module")
def get_rotation_matrix():
    def rmat(alpha, beta, gamma):
        """
        alpha : rotation around z axis - tait-bryan angle
        beta : rotation around y axis - tait-bryan
        gamma : rotation around x axis - tair-bryan
        """
        ca = np.cos(alpha)
        cb = np.cos(beta)
        cy = np.cos(gamma)

        sa = np.sin(alpha)
        sb = np.sin(beta)
        sy = np.sin(gamma)

        r = np.zeros((3, 3))
        r[0, 0] = ca * cb
        r[0, 1] = ca * sb * sy - sa * cy
        r[0, 2] = ca * sb * cy + sa * sy
        r[1, 0] = sa * cb
        r[1, 1] = sa * sb * sy + ca * cy
        r[1, 2] = sa * sb * cy - ca * sy
        r[2, 0] = -sb
        r[2, 1] = cb * sy
        r[2, 2] = cb * cy

        return r

    return rmat

@pytest.fixture(scope="class")
def dummy_rotated_accel(np_rng, get_rotation_matrix):
    x = np_rng.standard_normal((5000, 3))
    x *= 0.1  # make the noise small
    x[:, 2] += 1.0  # make one axis gravity

    r = get_rotation_matrix(0.1, 0.13, -0.07)

    x_rot = (r @ x.T).T

    return x, x_rot

@pytest.fixture(scope="module")
def day_ends():
    return 2000, 4000

@pytest.fixture(scope="module")
def sleep_ends():
    # treat sleep as exclusionary
    sleep_starts = {
        1: np.array([200, 1200, 2200, 4200]),
        2: np.array([200, 1800, 3800]),
        3: np.array([200, 1500, 4200]),
    }
    sleep_stops = {
        1: np.array([800, 1800, 2800, 5000]),
        2: np.array([800, 2500, 4400]),
        3: np.array([200, 1900, 5000]),
    }
    return sleep_starts, sleep_stops

@pytest.fixture(scope="module")
def wear_ends():
    wear_starts = np.array([200, 2300, 3000])
    wear_stops = np.array([1800, 2900, 3900])
    return wear_starts, wear_stops

@pytest.fixture(scope="module")
def true_intersect_ends():
    starts = {
        1: np.array([2800, 3000]),
        2: np.array([2500, 3000]),
        3: np.array([2300, 3000]),
    }
    stops = {
        1: np.array([2900, 3900]),
        2: np.array([2900, 3800]),
        3: np.array([2900, 3900]),
    }
    return starts, stops

@pytest.fixture(scope="module")
def true_sleep_only_ends():
    starts = {1: np.array([2000, 2800]), 2: np.array([2500]), 3: np.array([2000])}
    stops = {1: np.array([2200, 4000]), 2: np.array([3800]), 3: np.array([4000])}
    return starts, stops

@pytest.fixture(scope="class")
def rle_arr():
    return [0] * 5 + [1] * 3 + [0] * 4 + [1] * 7 + [0] * 2 + [1] * 6 + [0] * 1

@pytest.fixture(scope="class")
def rle_truth():
    lens = np.asarray([5, 3, 4, 7, 2, 6, 1])
    indices = np.asarray([0, 5, 8, 12, 19, 21, 27])
    vals = np.asarray([0, 1, 0, 1, 0, 1, 0])
    return lens, indices, vals

@pytest.fixture(scope="class")
def dummy_time():
    return np.arange(0, 10, 0.02)  # 50 hz, 500 samples

@pytest.fixture(scope="class")
def dummy_acc(np_rng):
    return np_rng.random((500, 3))

@pytest.fixture(scope="class")
def dummy_idx_1d():
    # original and truth value
    return np.array([0, 20, 113, 265, 481]), np.array([0, 4, 23, 53, 96])

@pytest.fixture(scope="class")
def dummy_idx_2d():
    # original and truth value
    orig = np.array([[0, 20], [113, 265], [481, 499]])
    truth = np.array([[0, 4], [23, 53], [96, 99]])
    return orig, truth

@pytest.fixture(scope="module")
def dummy_frag_predictions():
    return np.array(
        [
            0,
            0,
            0,
            1,
            1,
            1,
            1,
            1,
            1,
            1,
            0,
            0,
            1,
            1,
            1,
            0,
            0,
            0,
            0,
            0,
            1,
            1,
            1,
            0,
            0,
            0,
            0,
            0,
            0,
            0,
        ]
    )

max_stars_repo_name: PfizerRD/scikit-digital-health

max_issues_repo_name: PfizerRD/scikit-digital-health

max_forks_repo_name: PfizerRD/scikit-digital-health


text:
from abc import abstractmethod, abstractproperty
from contextlib import contextmanager
from multiprocessing import BoundedSemaphore
from pathlib import Path
from typing import Callable, Iterable, Union

import matplotlib
import matplotlib.pyplot as plt
import seaborn as sns
import nibabel as nib
import numpy as np
import pandas as pd
from neurolang.frontend import NeurolangDL, NeurolangPDL
from neurolang.frontend.neurosynth_utils import StudyID
from neurolang.regions import ExplicitVBR, ExplicitVBROverlay, region_union
from nilearn import datasets, image

class NeurolangEngineSet:
    """
    Utility class to hold a set of Neurolang engines (either NeurolangDL
    or NeurolangPDL). It keeps track of the total number of engines added
    to the set, as well as a semaphore to manage acquiring and releasing the
    engines in a safe concurrent manner.
    """

    def __init__(self, engine: Union[NeurolangDL, NeurolangPDL]) -> None:
        self.engines = set((engine,))
        self.counter = 1
        self.sema = BoundedSemaphore(value=1)

    def add_engine(self, engine: Union[NeurolangDL, NeurolangPDL]) -> None:
        """
        Add an engine to the set and increase the total number of engines
        in the set.

        Parameters
        ----------
        engine : Union[NeurolangDL, NeurolangPDL]
            the engine to add.
        """
        self.engines.add(engine)
        self.counter += 1
        self.sema = BoundedSemaphore(value=self.counter)

    @contextmanager
    def engine(self, timeout: int = None):
        """
        ContextManager to safely acquire an engine from the set.
        This method will block on the semaphore until an engine can be
        safely used from this set, and then yield it.
        At the end, the engine is put back in the set and the semaphore
        value is released.

        Parameters
        ----------
        timout : int
            When invoked with a timeout other than None, it will block for
            at most timeout seconds.

        Yields
        -------
        Union[NeurolangDL, NeurolangPDL]
            an engine.
        """
        lock = self.sema.acquire(timeout=timeout)
        if lock:
            try:
                engine = self.engines.pop()
                yield engine
            finally:
                self.engines.add(engine)
                self.sema.release()
        else:
            yield None

class NeurolangEngineConfiguration:
    """
    A NeurolangEngineConfiguration is a combination of an id key and
    a method which returns a Neurolang instance.
    """

    @abstractproperty
    def key(self):
        pass

    @abstractproperty
    def atlas(self):
        pass

    @abstractmethod
    def create(self) -> Union[NeurolangDL, NeurolangPDL]:
        pass

    def __eq__(self, other: object) -> bool:
        if not isinstance(other, NeurolangEngineConfiguration):
            return super().__eq__(other)
        return self.key == other.key

    def __hash__(self) -> int:
        return hash(self.key)

class NeurosynthEngineConf(NeurolangEngineConfiguration):
    def __init__(self, data_dir: Path, resolution=None) -> None:
        super().__init__()
        self.data_dir = data_dir
        self.resolution = resolution
        self._mni_atlas = None
        self._mni_brain_mask = None

    @property
    def key(self):
        return "neurosynth"

    @property
    def atlas(self):
        if self._mni_atlas is None:
            self._mni_atlas = nib.load(
                datasets.fetch_icbm152_2009(
                    data_dir=str(self.data_dir / "icbm")
                )["t1"]
            )
        return self._mni_atlas

    @property
    def brain_mask(self):
        if self._mni_brain_mask is None:
            self._mni_brain_mask = nib.load(
                datasets.fetch_icbm152_2009(
                    data_dir=str(self.data_dir / "icbm")
                )["mask"]
            )
        return self._mni_brain_mask

    def create(self) -> NeurolangPDL:
        mask = self.brain_mask
        if self.resolution is not None:
            mask = image.resample_img(mask, np.eye(3) * self.resolution)
        nl = init_frontend(mask)
        add_ploting_functions(nl)
        load_neurosynth_data(self.data_dir, nl, mask)
        return nl

def load_neurosynth_data(data_dir: Path, nl, mni_mask: nib.Nifti1Image):
    ns_database_fn, ns_features_fn = datasets.utils._fetch_files(
        data_dir / "neurosynth",
        [
            (
                "database.txt",
                "https://github.com/neurosynth/neurosynth-data/raw/e8f27c4a9a44dbfbc0750366166ad2ba34ac72d6/current_data.tar.gz",
                {"uncompress": True},
            ),
            (
                "features.txt",
                "https://github.com/neurosynth/neurosynth-data/raw/e8f27c4a9a44dbfbc0750366166ad2ba34ac72d6/current_data.tar.gz",
                {"uncompress": True},
            ),
        ],
    )

    activations = pd.read_csv(ns_database_fn, sep="\t")
    activations["id"] = activations["id"].apply(StudyID)
    mni_peaks = activations.loc[activations.space == "MNI"][
        ["x", "y", "z", "id"]
    ]
    non_mni_peaks = activations.loc[activations.space == "TAL"][
        ["x", "y", "z", "id"]
    ]
    proj_mat = np.linalg.pinv(
        np.array(
            [
                [0.9254, 0.0024, -0.0118, -1.0207],
                [-0.0048, 0.9316, -0.0871, -1.7667],
                [0.0152, 0.0883, 0.8924, 4.0926],
                [0.0, 0.0, 0.0, 1.0],
            ]
        ).T
    )
    projected = np.round(
        np.dot(
            np.hstack(
                (
                    non_mni_peaks[["x", "y", "z"]].values,
                    np.ones((len(non_mni_peaks), 1)),
                )
            ),
            proj_mat,
        )[:, 0:3]
    )
    projected_df = pd.DataFrame(
        np.hstack([projected, non_mni_peaks[["id"]].values]),
        columns=["x", "y", "z", "id"],
    )
    peak_data = pd.concat([projected_df, mni_peaks]).astype(
        {"x": int, "y": int, "z": int}
    )
    study_ids = peak_data[["id"]].drop_duplicates()

    features = pd.read_csv(ns_features_fn, sep="\t")
    features.rename(columns={"pmid": "id"}, inplace=True)

    term_data = pd.melt(
        features,
        var_name="term",
        id_vars="id",
        value_name="tfidf",
    ).query("tfidf > 0")[["term", "tfidf", "id"]]
    term_data["id"] = term_data["id"].apply(StudyID)

    nl.add_tuple_set(peak_data, name="PeakReported")
    nl.add_tuple_set(study_ids, name="Study")
    nl.add_tuple_set(term_data, name="TermInStudyTFIDF")
    nl.add_uniform_probabilistic_choice_over_set(
        study_ids, name="SelectedStudy"
    )
    nl.add_tuple_set(
        np.round(
            nib.affines.apply_affine(
                mni_mask.affine,
                np.transpose(mni_mask.get_fdata().astype(int).nonzero()),
            )
        ).astype(int),
        name="Voxel",
    )

def init_frontend(mni_mask):
    """
    Create a Neurolang Probabilistic engine and add some aggregation methods.

    Returns
    -------
    NeurolangPDL
        the Neurolang engine
    """
    nl = NeurolangPDL()

    @nl.add_symbol
    def agg_create_region_ijk(
        i: Iterable, j: Iterable, k: Iterable
    ) -> ExplicitVBR:
        voxels = np.c_[i, j, k]
        return ExplicitVBR(voxels, mni_mask.affine, image_dim=mni_mask.shape)

    @nl.add_symbol
    def agg_create_region_overlay_ijk(
        i: Iterable, j: Iterable, k: Iterable, p: Iterable
    ) -> ExplicitVBR:
        voxels = np.c_[i, j, k]
        return ExplicitVBROverlay(
            voxels, mni_mask.affine, p, image_dim=mni_mask.shape
        )

    @nl.add_symbol
    def agg_create_region(
        x: Iterable, y: Iterable, z: Iterable
    ) -> ExplicitVBR:
        mni_coords = np.c_[x, y, z]
        voxels = np.round(
            nib.affines.apply_affine(
                np.linalg.inv(mni_mask.affine),
                mni_coords.astype(float),
            )
        ).astype(int)
        return ExplicitVBR(voxels, mni_mask.affine, image_dim=mni_mask.shape)

    @nl.add_symbol
    def agg_create_region_overlay(
        x: Iterable, y: Iterable, z: Iterable, p: Iterable
    ) -> ExplicitVBR:
        mni_coords = np.c_[x, y, z]
        voxels = np.round(
            nib.affines.apply_affine(
                np.linalg.inv(mni_mask.affine),
                mni_coords.astype(float),
            )
        ).astype(int)
        return ExplicitVBROverlay(
            voxels, mni_mask.affine, p, image_dim=mni_mask.shape
        )

    @nl.add_symbol
    def startswith(prefix: str, s: str) -> bool:
        """Describe the prefix of string `s`.

        Parameters
        ----------
        prefix : str
            prefix to query.
        s : str
            string to check whether its
            prefixed by `s`.

        Returns
        -------
        bool
            whether `s` is prefixed by
            `prefix`.
        """
        return s.startswith(prefix)

    @nl.add_symbol
    def principal_direction(s: ExplicitVBR, direction: str, eps=1e-6) -> bool:
        """Describe the principal direction of
        the extension of a volumetric region.

        Parameters
        ----------
        s : ExplicitVBR
            region to analyse the principal
            direction of its extension.
        direction : str
            principal directions, one of
            `LR`, `AP`, `SI`, for the directions
            left-right, anterio-posterior, and
            superior inferior respectively.
        eps : float, optional
            minimum difference on between
            directional standard deviations,
            by default 1e-6.

        Returns
        -------
        bool
            wether the principal variance of
            `s` is `direction`.
        """
        # Assuming RAS coding os the xyz space.
        c = ["LR", "AP", "SI"]

        s_xyz = s.to_xyz()
        cov = np.cov(s_xyz.T)
        evals, evecs = np.linalg.eig(cov)
        i = np.argmax(np.abs(evals))
        abs_max_evec = np.abs(evecs[:, i].squeeze())
        sort_dir = np.argsort(abs_max_evec)
        if (
            np.abs(abs_max_evec[sort_dir[-1]] - abs_max_evec[sort_dir[-2]])
            < eps
        ):
            return False
        else:
            main_dir = c[sort_dir[-1]]
        return (direction == main_dir) or (direction[::-1] == main_dir)

    return nl

class DestrieuxEngineConf(NeurolangEngineConfiguration):
    def __init__(self, data_dir: Path, resolution: int = None) -> None:
        super().__init__()
        self.data_dir = data_dir
        self.resolution = resolution
        self._mni_atlas = None

    @property
    def key(self):
        return "destrieux"

    @property
    def atlas(self):
        if self._mni_atlas is None:
            self._mni_atlas = nib.load(
                datasets.fetch_icbm152_2009(
                    data_dir=str(self.data_dir / "icbm")
                )["t1"]
            )
        return self._mni_atlas

    def create(self) -> NeurolangPDL:
        mask = self.atlas
        if self.resolution is not None:
            mask = image.resample_img(mask, np.eye(3) * self.resolution)
        nl = init_frontend(mask)
        nl.add_symbol(
            region_union,
            name="region_union",
            type_=Callable[[Iterable[ExplicitVBR]], ExplicitVBR],
        )

        load_destrieux_atlas(self.data_dir, nl)
        return nl

def load_destrieux_atlas(data_dir, nl):
    destrieux_atlas = datasets.fetch_atlas_destrieux_2009(
        data_dir=str(data_dir / "destrieux")
    )

    nl.new_symbol(name="destrieux")
    destrieux_atlas_image = nib.load(destrieux_atlas["maps"])
    destrieux_labels = dict(destrieux_atlas["labels"])
    destrieux_set = set()
    for k, v in destrieux_labels.items():
        if k == 0:
            continue
        destrieux_set.add(
            (
                v.decode("utf8").replace("-", " ").replace("_", " "),
                ExplicitVBR.from_spatial_image_label(destrieux_atlas_image, k),
            )
        )
    nl.add_tuple_set(destrieux_set, name="destrieux")

def add_ploting_functions(nl: Union[NeurolangDL, NeurolangPDL]):
    matplotlib.use('Agg')

    @nl.add_symbol
    def agg_kde(terms: Iterable, probs: Iterable) -> matplotlib.figure.Figure:
        """
        Create a kde plot showing prob distribution per term.

        Parameters
        ----------
        terms : Iterable[str]
            the terms
        probs: Iterable[float]
            the probs

        Returns
        -------
        matplotlib.figure.Figure
            a Figure
        """
        df = pd.DataFrame(
            {
                "terms": terms,
                "probs": probs,
            }
        )
        fig, ax = plt.subplots()
        fig.suptitle("Distribution of probs / term")
        sns.kdeplot(ax=ax, data=df[df.probs > 0.002], x="probs", hue="terms")
        return fig

max_stars_repo_name: gzanitti/NeuroLang

max_issues_repo_name: gzanitti/NeuroLang

max_forks_repo_name: jonasrenault/NeuroLang


text:
# -*- coding: utf-8 -*-
"""
Created on Tue May 17 12:46:20 2016

@author: Hossam Faris
"""

import numpy
import math

# define the function blocks
def prod(it):
    p = 1
    for n in it:
        p *= n
    return p

def Ufun(x, a, k, m):
    y = k * ((x - a) ** m) * (x > a) + k * ((-x - a) ** m) * (x < (-a))
    return y

def F1(x):
    s = numpy.sum(x ** 2)
    return s

def F2(x):
    o = sum(abs(x)) + prod(abs(x))
    return o

def F3(x):
    dim = len(x) + 1
    o = 0
    for i in range(1, dim):
        o = o + (numpy.sum(x[0:i])) ** 2
    return o

def F4(x):
    o = max(abs(x))
    return o

def F5(x):
    dim = len(x)
    o = numpy.sum(
        100 * (x[1:dim] - (x[0 : dim - 1] ** 2)) ** 2 + (x[0 : dim - 1] - 1) ** 2
    )
    return o

def F6(x):
    o = numpy.sum(abs((x + 0.5)) ** 2)
    return o

def F7(x):
    dim = len(x)

    w = [i for i in range(len(x))]
    for i in range(0, dim):
        w[i] = i + 1
    o = numpy.sum(w * (x ** 4)) + numpy.random.uniform(0, 1)
    return o

def F8(x):
    o = sum(-x * (numpy.sin(numpy.sqrt(abs(x)))))
    return o

def F9(x):
    dim = len(x)
    o = numpy.sum(x ** 2 - 10 * numpy.cos(2 * math.pi * x)) + 10 * dim
    return o

def F10(x):
    dim = len(x)
    o = (
        -20 * numpy.exp(-0.2 * numpy.sqrt(numpy.sum(x ** 2) / dim))
        - numpy.exp(numpy.sum(numpy.cos(2 * math.pi * x)) / dim)
        + 20
        + numpy.exp(1)
    )
    return o

def F11(x):
    dim = len(x)
    w = [i for i in range(len(x))]
    w = [i + 1 for i in w]
    o = numpy.sum(x ** 2) / 4000 - prod(numpy.cos(x / numpy.sqrt(w))) + 1
    return o

def F12(x):
    dim = len(x)
    o = (math.pi / dim) * (
        10 * ((numpy.sin(math.pi * (1 + (x[0] + 1) / 4))) ** 2)
        + numpy.sum(
            (((x[: dim - 1] + 1) / 4) ** 2)
            * (1 + 10 * ((numpy.sin(math.pi * (1 + (x[1 :] + 1) / 4)))) ** 2)
        )
        + ((x[dim - 1] + 1) / 4) ** 2
    ) + numpy.sum(Ufun(x, 10, 100, 4))
    return o

def F13(x):
    if x.ndim==1:
        x = x.reshape(1,-1)

    o = 0.1 * (
        (numpy.sin(3 * numpy.pi * x[:,0])) ** 2
        + numpy.sum(
            (x[:,:-1] - 1) ** 2
            * (1 + (numpy.sin(3 * numpy.pi * x[:,1:])) ** 2), axis=1
        )
        + ((x[:,-1] - 1) ** 2) * (1 + (numpy.sin(2 * numpy.pi * x[:,-1])) ** 2)
    ) + numpy.sum(Ufun(x, 5, 100, 4))
    return o

def F14(x):
    aS = [
        [
            -32, -16, 0, 16, 32,
            -32, -16, 0, 16, 32,
            -32, -16, 0, 16, 32,
            -32, -16, 0, 16, 32,
            -32, -16, 0, 16, 32,
        ],
        [
            -32, -32, -32, -32, -32,
            -16, -16, -16, -16, -16,
            0, 0, 0, 0, 0,
            16, 16, 16, 16, 16,
            32, 32, 32, 32, 32,
        ],
    ]
    aS = numpy.asarray(aS)
    bS = numpy.zeros(25)
    v = numpy.matrix(x)
    for i in range(0, 25):
        H = v - aS[:, i]
        bS[i] = numpy.sum((numpy.power(H, 6)))
    w = [i for i in range(25)]
    for i in range(0, 24):
        w[i] = i + 1
    o = ((1.0 / 500) + numpy.sum(1.0 / (w + bS))) ** (-1)
    return o

def F15(L):
    aK = [
        0.1957, 0.1947, 0.1735, 0.16, 0.0844, 0.0627,
        0.0456, 0.0342, 0.0323, 0.0235, 0.0246,
    ]
    bK = [0.25, 0.5, 1, 2, 4, 6, 8, 10, 12, 14, 16]
    aK = numpy.asarray(aK)
    bK = numpy.asarray(bK)
    bK = 1 / bK
    fit = numpy.sum(
        (aK - ((L[0] * (bK ** 2 + L[1] * bK)) / (bK ** 2 + L[2] * bK + L[3]))) ** 2
    )
    return fit

def F16(L):
    o = (
        4 * (L[0] ** 2)
        - 2.1 * (L[0] ** 4)
        + (L[0] ** 6) / 3
        + L[0] * L[1]
        - 4 * (L[1] ** 2)
        + 4 * (L[1] ** 4)
    )
    return o

def F17(L):
    o = (
        (L[1] - (L[0] ** 2) * 5.1 / (4 * (numpy.pi ** 2)) + 5 / numpy.pi * L[0] - 6)
        ** 2
        + 10 * (1 - 1 / (8 * numpy.pi)) * numpy.cos(L[0])
        + 10
    )
    return o

def F18(L):
    o = (
        1
        + (L[0] + L[1] + 1) ** 2
        * (
            19
            - 14 * L[0]
            + 3 * (L[0] ** 2)
            - 14 * L[1]
            + 6 * L[0] * L[1]
            + 3 * L[1] ** 2
        )
    ) * (
        30
        + (2 * L[0] - 3 * L[1]) ** 2
        * (
            18
            - 32 * L[0]
            + 12 * (L[0] ** 2)
            + 48 * L[1]
            - 36 * L[0] * L[1]
            + 27 * (L[1] ** 2)
        )
    )
    return o

# map the inputs to the function blocks
def F19(L):
    aH = [[3, 10, 30], [0.1, 10, 35], [3, 10, 30], [0.1, 10, 35]]
    aH = numpy.asarray(aH)
    cH = [1, 1.2, 3, 3.2]
    cH = numpy.asarray(cH)
    pH = [
        [0.3689, 0.117, 0.2673],
        [0.4699, 0.4387, 0.747],
        [0.1091, 0.8732, 0.5547],
        [0.03815, 0.5743, 0.8828],
    ]
    pH = numpy.asarray(pH)
    o = 0
    for i in range(0, 4):
        o = o - cH[i] * numpy.exp(-(numpy.sum(aH[i, :] * ((L - pH[i, :]) ** 2))))
    return o

def F20(L):
    aH = [
        [10, 3, 17, 3.5, 1.7, 8],
        [0.05, 10, 17, 0.1, 8, 14],
        [3, 3.5, 1.7, 10, 17, 8],
        [17, 8, 0.05, 10, 0.1, 14],
    ]
    aH = numpy.asarray(aH)
    cH = [1, 1.2, 3, 3.2]
    cH = numpy.asarray(cH)
    pH = [
        [0.1312, 0.1696, 0.5569, 0.0124, 0.8283, 0.5886],
        [0.2329, 0.4135, 0.8307, 0.3736, 0.1004, 0.9991],
        [0.2348, 0.1415, 0.3522, 0.2883, 0.3047, 0.6650],
        [0.4047, 0.8828, 0.8732, 0.5743, 0.1091, 0.0381],
    ]
    pH = numpy.asarray(pH)
    o = 0
    for i in range(0, 4):
        o = o - cH[i] * numpy.exp(-(numpy.sum(aH[i, :] * ((L - pH[i, :]) ** 2))))
    return o

def F21(L):
    aSH = [
        [4, 4, 4, 4],
        [1, 1, 1, 1],
        [8, 8, 8, 8],
        [6, 6, 6, 6],
        [3, 7, 3, 7],
        [2, 9, 2, 9],
        [5, 5, 3, 3],
        [8, 1, 8, 1],
        [6, 2, 6, 2],
        [7, 3.6, 7, 3.6],
    ]
    cSH = [0.1, 0.2, 0.2, 0.4, 0.4, 0.6, 0.3, 0.7, 0.5, 0.5]
    aSH = numpy.asarray(aSH)
    cSH = numpy.asarray(cSH)
    fit = 0
    for i in range(5):
        v = numpy.matrix(L - aSH[i, :])
        fit = fit - ((v) * (v.T) + cSH[i]) ** (-1)
    o = fit.item(0)
    return o

def F22(L):
    aSH = [
        [4, 4, 4, 4],
        [1, 1, 1, 1],
        [8, 8, 8, 8],
        [6, 6, 6, 6],
        [3, 7, 3, 7],
        [2, 9, 2, 9],
        [5, 5, 3, 3],
        [8, 1, 8, 1],
        [6, 2, 6, 2],
        [7, 3.6, 7, 3.6],
    ]
    cSH = [0.1, 0.2, 0.2, 0.4, 0.4, 0.6, 0.3, 0.7, 0.5, 0.5]
    aSH = numpy.asarray(aSH)
    cSH = numpy.asarray(cSH)
    fit = 0
    for i in range(7):
        v = numpy.matrix(L - aSH[i, :])
        fit = fit - ((v) * (v.T) + cSH[i]) ** (-1)
    o = fit.item(0)
    return o

def F23(L):
    aSH = [
        [4, 4, 4, 4],
        [1, 1, 1, 1],
        [8, 8, 8, 8],
        [6, 6, 6, 6],
        [3, 7, 3, 7],
        [2, 9, 2, 9],
        [5, 5, 3, 3],
        [8, 1, 8, 1],
        [6, 2, 6, 2],
        [7, 3.6, 7, 3.6],
    ]
    cSH = [0.1, 0.2, 0.2, 0.4, 0.4, 0.6, 0.3, 0.7, 0.5, 0.5]
    aSH = numpy.asarray(aSH)
    cSH = numpy.asarray(cSH)
    fit = 0
    for i in range(10):
        v = numpy.matrix(L - aSH[i, :])
        fit = fit - ((v) * (v.T) + cSH[i]) ** (-1)
    o = fit.item(0)
    return o

def getFunctionDetails(a):
    # [name, lb, ub, dim]
    param = {
        "F1": ["F1", -100, 100, 30],
        "F2": ["F2", -10, 10, 30],
        "F3": ["F3", -100, 100, 30],
        "F4": ["F4", -100, 100, 30],
        "F5": ["F5", -30, 30, 30],
        "F6": ["F6", -100, 100, 30],
        "F7": ["F7", -1.28, 1.28, 30],
        "F8": ["F8", -500, 500, 30],
        "F9": ["F9", -5.12, 5.12, 30],
        "F10": ["F10", -32, 32, 30],
        "F11": ["F11", -600, 600, 30],
        "F12": ["F12", -50, 50, 30],
        "F13": ["F13", -50, 50, 30],
        "F14": ["F14", -65.536, 65.536, 2],
        "F15": ["F15", -5, 5, 4],
        "F16": ["F16", -5, 5, 2],
        "F17": ["F17", -5, 15, 2],
        "F18": ["F18", -2, 2, 2],
        "F19": ["F19", 0, 1, 3],
        "F20": ["F20", 0, 1, 6],
        "F21": ["F21", 0, 10, 4],
        "F22": ["F22", 0, 10, 4],
        "F23": ["F23", 0, 10, 4],
    }
    return param.get(a, "nothing")

'''
# ESAs space mission design benchmarks https://www.esa.int/gsp/ACT/projects/gtop/
from fcmaes.astro import (
    MessFull,
    Messenger,
    Gtoc1,
    Cassini1,
    Cassini2,
    Rosetta,
    Tandem,
    Sagas,
)
def Ca1(x):
    return Cassini1().fun(x)
def Ca2(x):
    return Cassini2().fun(x)
def Ros(x):
    return Rosetta().fun(x)
def Tan(x):
    return Tandem(5).fun(x)
def Sag(x):
    return Sagas().fun(x)
def Mef(x):
    return MessFull().fun(x)
def Mes(x):
    return Messenger().fun(x)
def Gt1(x):
    return Gtoc1().fun(x)

def getFunctionDetails(a):
    # [name, lb, ub, dim]
    param = {
        "F1": ["F1", -100, 100, 30],
        "F2": ["F2", -10, 10, 30],
        "F3": ["F3", -100, 100, 30],
        "F4": ["F4", -100, 100, 30],
        "F5": ["F5", -30, 30, 30],
        "F6": ["F6", -100, 100, 30],
        "F7": ["F7", -1.28, 1.28, 30],
        "F8": ["F8", -500, 500, 30],
        "F9": ["F9", -5.12, 5.12, 30],
        "F10": ["F10", -32, 32, 30],
        "F11": ["F11", -600, 600, 30],
        "F12": ["F12", -50, 50, 30],
        "F13": ["F13", -50, 50, 30],
        "F14": ["F14", -65.536, 65.536, 2],
        "F15": ["F15", -5, 5, 4],
        "F16": ["F16", -5, 5, 2],
        "F17": ["F17", -5, 15, 2],
        "F18": ["F18", -2, 2, 2],
        "F19": ["F19", 0, 1, 3],
        "F20": ["F20", 0, 1, 6],
        "F21": ["F21", 0, 10, 4],
        "F22": ["F22", 0, 10, 4],
        "F23": ["F23", 0, 10, 4],
        "Ca1": [
            "Ca1",
            Cassini1().bounds.lb,
            Cassini1().bounds.ub,
            len(Cassini1().bounds.lb),
        ],
        "Ca2": [
            "Ca2",
            Cassini2().bounds.lb,
            Cassini2().bounds.ub,
            len(Cassini2().bounds.lb),
        ],
        "Gt1": ["Gt1", Gtoc1().bounds.lb, Gtoc1().bounds.ub, len(Gtoc1().bounds.lb)],
        "Mes": [
            "Mes",
            Messenger().bounds.lb,
            Messenger().bounds.ub,
            len(Messenger().bounds.lb),
        ],
        "Mef": [
            "Mef",
            MessFull().bounds.lb,
            MessFull().bounds.ub,
            len(MessFull().bounds.lb),
        ],
        "Sag": ["Sag", Sagas().bounds.lb, Sagas().bounds.ub, len(Sagas().bounds.lb)],
        "Tan": [
            "Tan",
            Tandem(5).bounds.lb,
            Tandem(5).bounds.ub,
            len(Tandem(5).bounds.lb),
        ],
        "Ros": [
            "Ros",
            Rosetta().bounds.lb,
            Rosetta().bounds.ub,
            len(Rosetta().bounds.lb),
        ],
    }
    return param.get(a, "nothing")
'''


text:
from Bio import SeqIO
from Bio.SeqUtils import GC,ProtParam
from Bio.Seq import Seq
import numpy
import os,re
from lib.findcds import FindCDS
from lib.fickett import Fickett
from lib.eiip import *
import pandas as pd
from sklearn.ensemble import RandomForestClassifier
import pickle
## (1) farward frame translations
## result: A translated sequence will be generated
def frame_translation(seq, shift, genetic_code=1):
    from Bio.Seq import translate
    length = len(seq)
    frame = {}
    fragment_length = 3 * ((length-shift) // 3)
    frame = translate(seq[shift:shift+fragment_length], genetic_code)
    return frame

def mRNA_translate(mRNA):
    return Seq(mRNA).translate()

def protein_param(putative_seqprot):
    return putative_seqprot.isoelectric_point()

def txcds(FAfile):
   # ### run txCdsPredict on the input FASTA file ####
    cmd = "txCdsPredict " + FAfile + " -anyStart tmp.cds"
    os.system(cmd)
    cds_len = dict()
    cds_score = dict()
    temp = open("tmp.cds")
    for line in temp:
        line_array = line.split()
        id_array = line_array[0].split('|')
        trans_id = id_array[0]
        pred_start = int(line_array[1])
        pred_end = int(line_array[2])
        pred_len = pred_end-pred_start
        cds_len[trans_id] = pred_len
        cds_score[trans_id] = float(line_array[5])
    temp.close()
    os.system("rm tmp.cds")
    return cds_len, cds_score
    ############ end of running txCdsPredict ################

def getStop(seq):
    # 5
    frame_0 = frame_translation(seq, 0, genetic_code=1)
    frame_1 = frame_translation(seq, 1, genetic_code=1)
    frame_2 = frame_translation(seq, 2, genetic_code=1)
    stop = (frame_0.count("*"), frame_1.count("*"), frame_2.count("*"))
    std_stop = numpy.std(stop)
    return std_stop

def getex(seqRNA,strinfoAmbiguous):
    '''seqCDS:longest ORF'''
    seqCDS, _ , _ , orf_fullness = FindCDS(seqRNA).longest_orf("+")
    seqprot = mRNA_translate(seqCDS)
    pep_len = len(seqprot)  # pep_len = len(seqprot.strip("*"))
    newseqprot = strinfoAmbiguous.sub("", str(seqprot))
    '''exclude ambiguous amio acid X, B, Z, J, Y in peptide sequence'''
    # 7 6
    protparam_obj = ProtParam.ProteinAnalysis(str(newseqprot.strip("*")))
    if pep_len > 0:
        isoelectric_point = protein_param(protparam_obj)
    else:
        orf_fullness = -1
        isoelectric_point = 0.0
    return isoelectric_point, orf_fullness

def getlenperc(cds_len,first,seq_len):
    cds = cds_len[first[0]]
    len_perc = float(cds)/seq_len
    return len_perc

def getcdss(cds_score,first):
    return cds_score[first[0]]

def write_files(seq,cds_len,cds_score,first,seqRNA,strinfoAmbiguous,fickett_obj, f_out,class_type=0):
    seq_len = len(seq)  # 1
    cdsscore = getcdss(cds_score,first) # 2
    fickett_score = fickett_obj.fickett_value(seqRNA) # 3
    gc = GC(seq) # 4
    std_stop = getStop(seq) # 5
    isoelectric_point,orf_fullness = getex(seqRNA,strinfoAmbiguous) # 6 7
    len_perc = getlenperc(cds_len,first,seq_len) # 8
    cwt = get_Ewtf(get_numseq(seq, 1), wavlet)  # 9

    # print features
    print("%d" % seq_len, end='\t', file=f_out)  # 1
    print("%d" % cdsscore, end='\t', file=f_out)  # 2
    print("%s" % str(fickett_score), end='\t', file=f_out)  # 3
    print("%0.2f" % gc, end='\t', file=f_out)  # 4
    print("%.6f" % std_stop, end='\t', file=f_out)  # 5
    print("%s" % str(orf_fullness), end='\t', file=f_out)  # 6
    print("%s" % str(isoelectric_point), end='\t', file=f_out)  # 7
    print("%0.2f" % len_perc, end='\t', file=f_out)  # 8
    for o in cwt:
        print("%0.8f" % o, end='\t', file=f_out)  # 9
    if class_type != -1:
        print("%d" % class_type, file=f_out)

def get_feat(FAfile,outFileLink,islncs = 0):
    strinfoAmbiguous = re.compile("X|B|Z|J|U", re.I)
    ptU = re.compile("U", re.I)
    fickett_obj = Fickett()
    # ------

    f_in = open(FAfile)
    cds_len, cds_score = txcds(FAfile) # 2 # ### run txCdsPredict on the input FASTA file ####
    print("Extra Informations Start:")
    proIndex = 0
    for record in SeqIO.parse(f_in, "fasta"):
        proIndex = proIndex + 1
        print("%s\r" % str(proIndex).zfill(7),end='')
        ID = record.id
        first = ID.split("|")
        seq = record.seq
        seqRNA = ptU.sub("T", str(seq).strip())
        seqRNA = seqRNA.upper()
        write_files(seq,cds_len,cds_score,first,seqRNA,strinfoAmbiguous,fickett_obj, outFileLink,islncs)
    print("".zfill(7))
    print("Extra Informations Finished")
    f_in.close()

help_docs = '''
Usage:
1 . \033[31;1mIf you want to train your model\033[0m
    \033[31;1mExample:\033[0m prelnc2.py \033[31;1m-l\033[0m ath_lnc.fa \033[31;1m-p\033[0m ath_pc.fa \033[31;1m-t\033[0m train \033[31;1m-o\033[0m ath_model.pkl \033[31;1m0 100 0.3 10\033[0m
    -l lncRNAfile    -p pcRNAfile    -train for trainning   -o outputfile
                        \033[31;1m 0 100 0.3 10 \033[0m
    Default max_depth=\033[31;1mNone\033[0m, min_samples_leaf=\033[31;1m100\033[0m, min_weight_fraction_leaf=\033[31;1m0.3\033[0m,n_estimators=\033[31;1m10\033[0m,
2 . \033[31;1mIf you want to predict your RNA\033[0m
    \033[31;1mExample:\033[0mprelnc2.py \033[31;1m-f\033[0m ath_lnc.fa \033[31;1m-m\033[0m ath_model.pkl \033[31;1m-t\033[0m predict \033[31;1m-o\033[0m result.txt
    -f  file_in for predict  -m modelfile   -predict   -o outputfile
'''

def model_create(model_path, train_f, r_para):
    # create model
    data = pd.read_table(train_f, header=None)
    x = numpy.array(data[:-1])
    y = numpy.array(data[-1])
    rf = RandomForestClassifier(n_jobs=-1,
                                max_depth=r_para[0],
                                min_samples_leaf=r_para[1],
                                min_weight_fraction_leaf=r_para[2],
                                n_estimators = r_para[3]
                                # max_depth=None, min_samples_leaf=100, min_weight_fraction_leaf=0.3,n_estimators=10,
                                )

    model = rf.fit(x, y)
    with open(model_path, 'wb') as f:
        pickle.dump(model, f)

def model_predict(model_path, test_f):
    with open(model_path, 'rb') as m:
        model = pickle.load(m)

    # predict model
    data = pd.read_table(test_f, header=None)
    x = numpy.array(data)
    y_pre = model.predict(x)
    with open('temp', 'w') as t:
        numpy.savetxt(t, y_pre, fmt='%0.8f')
    cl = "paste " + test_f + " temp > " + test_f + ".predict.txt"
    os.system(cl)
    os.system("rm temp")


text:
import numpy as np
import cv2
import torch
from .post_parser import remove_subjects

def padding_image_overlap(image, overlap_ratio=0.46):
    h, w = image.shape[:2]
    pad_length = int(h* overlap_ratio)
    pad_w = w+2*pad_length
    pad_image = np.zeros((h, pad_w, 3), dtype=np.uint8)
    top, left = 0, pad_length
    bottom, right = h, w+pad_length
    pad_image[top:bottom, left:right] = image

    # due to BEV takes square input, so we convert top, bottom to the state that assuming square padding
    pad_height = (w - h)//2
    top = pad_height
    bottom = w - top
    left = 0
    right = w
    image_pad_info = torch.Tensor([top, bottom, left, right, h, w])
    return pad_image, image_pad_info, pad_length

def get_image_split_plan(image, overlap_ratio=0.46):
    h, w = image.shape[:2]
    aspect_ratio = w / h
    slide_time = int(np.ceil((aspect_ratio - 1) / (1 - overlap_ratio))) + 1

    crop_box = [] # left, right, top, bottom
    move_step = (1 - overlap_ratio) * h
    for ind in range(slide_time):
        if ind == (slide_time-1):
            left = w-h
        else:
            left = move_step * ind
            right = left+h
        crop_box.append([left, right, 0, h])

    return np.array(crop_box).astype(np.int32)

def exclude_boudary_subjects(outputs, drop_boundary_ratio, ptype='left', torlerance=0.05):
    if ptype=='left':
        drop_mask = outputs['cam'][:, 2] > (1 - drop_boundary_ratio + torlerance)
    elif ptype=='right':
        drop_mask = outputs['cam'][:, 2] < (drop_boundary_ratio - 1 - torlerance)
    remove_subjects(outputs, torch.where(drop_mask)[0])

def convert_crop_cam_params2full_image(cam_params, crop_bbox, image_shape):
    h, w = image_shape
    # adjust scale, cam 3: depth, y, x
    scale_adjust = (crop_bbox[[1,3]]-crop_bbox[[0,2]]).max() / max(h, w)
    cam_params *= scale_adjust

    # adjust x
    # crop_bbox[:2] -= pad_length
    bbox_mean_x = crop_bbox[:2].mean()
    cam_params[:,2] += bbox_mean_x / (w /2) - 1
    return cam_params

def collect_outputs(outputs, all_outputs):
    keys = list(outputs.keys())
    for key in keys:
        if key not in all_outputs:
            all_outputs[key] = outputs[key]
        else:
            if key in ['smpl_face']:
                continue
            if key in ['center_map']:
                all_outputs[key] = torch.cat([all_outputs[key], outputs[key]],3)
                continue
            if key in ['center_map_3d']:
                all_outputs[key] = torch.cat([all_outputs[key], outputs[key]],2)
                continue
            all_outputs[key] = torch.cat([all_outputs[key], outputs[key]],0)

max_stars_repo_name: jjandnn/ROMP

max_issues_repo_name: jjandnn/ROMP

max_forks_repo_name: jjandnn/ROMP


text:
import argparse
import os
import pathlib
from tkinter import BOTTOM
import urllib.request
import numpy as np
import pandas as pd
from joblib import Parallel, delayed
import mne
from mne.io.brainvision.brainvision import _aux_vhdr_info

from mne_bids import write_raw_bids, print_dir_tree, make_report, BIDSPath

lemon_info = pd.read_csv(
    "./META_File_IDs_Age_Gender_Education_Drug_Smoke_SKID_LEMON.csv")
lemon_info = lemon_info.set_index("ID")
eeg_subjects = pd.read_csv('./lemon_eeg_subjects.csv')
lemon_info = lemon_info.loc[eeg_subjects.subject]
lemon_info['gender'] = lemon_info['Gender_ 1=female_2=male'].map({1: 2, 2: 1})
lemon_info['age_guess'] = np.array(
  lemon_info['Age'].str.split('-').tolist(), dtype=np.int).mean(1)
subjects = list(lemon_info.index)

def convert_lemon_to_bids(lemon_data_dir, bids_save_dir, n_jobs=1, DEBUG=False):
    """Convert TUAB dataset to BIDS format.

    Parameters
    ----------
    lemon_data_dir : str
        Directory where the original LEMON dataset is saved, e.g.
        `/storage/store3/data/LEMON_RAW`.
    bids_save_dir : str
        Directory where to save the BIDS version of the dataset.
    n_jobs : None | int
        Number of jobs for parallelization.
    """
    subjects_ = subjects
    if DEBUG:
        subjects_ = subjects[:1]

    good_subjects = Parallel(n_jobs=n_jobs)(
        delayed(_convert_subject)(subject, lemon_data_dir, bids_save_dir)
        for subject in subjects_)
    subjects_ = [sub for sub in good_subjects if not isinstance(sub, tuple)]
    _, bad_subjects, errs = zip(*[
        sub for sub in good_subjects if isinstance(sub, tuple)])
    bad_subjects = pd.DataFrame(
        dict(subjects= bad_subjects, error=errs))
    bad_subjects.to_csv(
        '/storage/store3/data/LEMON_EEG_BIDS/bids_conv_erros.csv')
    # update the participants file as LEMON has no official age data
    participants = pd.read_csv(
        "/storage/store3/data/LEMON_EEG_BIDS/participants.tsv", sep='\t')
    participants = participants.set_index("participant_id")
    participants.loc[subjects_, 'age'] = lemon_info.loc[subjects_, 'age_guess']
    participants.to_csv(
        "/storage/store3/data/LEMON_EEG_BIDS/participants.tsv", sep='\t')

def _convert_subject(subject, data_path, bids_save_dir):
    """Get the work done for one subject"""
    try:
        fname = pathlib.Path(data_path) / subject / "RSEEG" / f"{subject}.vhdr"
        raw = mne.io.read_raw_brainvision(fname)

        raw.set_channel_types({"VEOG": "eog"})
        montage = mne.channels.make_standard_montage('standard_1005')
        raw.set_montage(montage)
        sub_id = subject.strip("sub-")
        raw.info['subject_info'] = {
            'participant_id': sub_id,
            'sex': lemon_info.loc[subject, 'gender'],
            'age': lemon_info.loc[subject, 'age_guess'],
            # XXX LEMON shares no public age
            'hand': lemon_info.loc[subject, 'Handedness']
        }
        events, _ = mne.events_from_annotations(raw)

        events = events[(events[:, 2] == 200) | (events[:, 2] == 210)]
        event_id = {"eyes/open": 200, "eyes/closed": 210}
        bids_path = BIDSPath(
            subject=sub_id, session=None, task='RSEEG',
            run=None,
            root=bids_save_dir, datatype='eeg', check=True)

        write_raw_bids(
            raw,
            bids_path,
            events_data=events,
            event_id=event_id,
            overwrite=True
        )
    except Exception as err:
        print(err)
        return ("BAD", subject, err)
    return subject

if __name__ == '__main__':
    parser = argparse.ArgumentParser(description='Convert LEMON to BIDS.')
    parser.add_argument(
        '--lemon_data_dir', type=str,
        default='/storage/store3/data/LEMON_RAW',
        help='Path to the original data.')
    parser.add_argument(
        '--bids_data_dir', type=str,
        default=pathlib.Path("/storage/store3/data/LEMON_EEG_BIDS"),
        help='Path to where the converted data should be saved.')
    parser.add_argument(
        '--n_jobs', type=int, default=1,
        help='number of parallel processes to use (default: 1)')
    parser.add_argument(
        '--DEBUG', type=bool, default=False,
        help='activate debugging mode')
    args = parser.parse_args()

    convert_lemon_to_bids(
        args.lemon_data_dir, args.bids_data_dir, n_jobs=args.n_jobs,
        DEBUG=args.DEBUG)

    print_dir_tree(args.bids_data_dir)
    print(make_report(args.bids_data_dir))

max_stars_repo_name: agramfort/brain-age-benchmark-paper

max_issues_repo_name: agramfort/brain-age-benchmark-paper

max_forks_repo_name: agramfort/brain-age-benchmark-paper


text:
from collections import OrderedDict

import numpy as np
from ipyatom.bonds import compute_bonds
from ipyatom.transforms import apply_transforms
from jsonextended import edict
from scipy.spatial.qhull import ConvexHull

from ipyatom.repeat_density import cube_frac2cart

# don't make required just yet
try:
    import ipyvolume.pylab as p3
except ImportError:
    pass
try:
    import ipywidgets as widgets
except ImportError:
    pass

def _create_mesh(points, color="black", solid=False, lines=True, triangle_indices=None, line_indices=None):
    """ create an ipyvolume mesh from a number of points

    Parameters
    ----------
    points: list
        [[x, y, z], ...]
    color: str
    solid: bool
    triangle_indices: list or None
        [[i, j, k], ...], if None then computed from scipy.spatial.qhull.ConvexHull triangles
    line_indices: list or None
        [[i, j], ...], if None then computed from scipy.spatial.qhull.ConvexHull triangles

    Returns
    -------

    """
    x, y, z = np.array(points).T
    try:
        hull = ConvexHull(points)
    except:
        hull = ConvexHull(points, qhull_options='QJ')
    if triangle_indices is None:
        triangle_indices = hull.simplices.tolist()
    if line_indices is None:
        line_indices = []
        for i, j, k in hull.simplices.tolist():
            line_indices += [[i, j], [i, k], [j, k]]

    mesh = p3.plot_trisurf(x, y, z, triangles=triangle_indices, lines=line_indices, color=color)
    if not solid:
        mesh.visible_faces = False
    if not lines:
        mesh.visible_lines = False
    return mesh

def _get_cylinder_mesh(p0, p1, radius=1, segments=50):
    """

    Parameters
    ----------
    p0: list or numpy.array((3,))
    p1: list or numpy.array((3,))
    radius: float
    segments: int

    Returns
    -------

    """
    p0 = np.array(p0)
    p1 = np.array(p1)
    # vector in direction of axis
    v = p1 - p0
    # find magnitude of vector
    mag = np.linalg.norm(v)
    # unit vector in direction of axis
    v = v / mag
    # make some vector not in the same direction as v
    not_v = np.array([1, 0, 0])
    if (v == not_v).all() or (v == -not_v).all():
        not_v = np.array([0, 1, 0])
    # make vector perpendicular to v
    n1 = np.cross(v, not_v)
    # normalize n1
    n1 /= np.linalg.norm(n1)
    # make unit vector perpendicular to v and n1
    n2 = np.cross(v, n1)
    # surface ranges over t from 0 to length of axis and 0 to 2*pi
    t = np.linspace(0, mag, 2)
    theta = np.linspace(0, 2 * np.pi, segments)
    # use meshgrid to make 2d arrays
    t, theta = np.meshgrid(t, theta)
    # generate coordinates for surface
    X, Y, Z = [p0[i] + v[i] * t + radius * np.sin(theta) * n1[i] + radius * np.cos(theta) * n2[i] for i in [0, 1, 2]]
    return list(zip(X[:, 0], Y[:, 0], Z[:, 0])) + list(zip(X[:, 1], Y[:, 1], Z[:, 1]))

def _create_bond(p1, p2, c1, c2=None, radius=.1, segments=20):
    """ create a bond

    Parameters
    ----------
    p1: list or numpy.array((3,))
    p2: list or numpy.array((3,))
    c1: str or numpy.array((3,))
    c2: str or numpy.array((3,))
    radius: float
    segments: int

    """
    if c2 is None or np.all(c1 == c2):
        return [_create_mesh(_get_cylinder_mesh(p1, p2, radius=radius, segments=segments),
                             c1, solid=True, lines=False)]
    else:
        p1 = np.array(p1)
        p2 = np.array(p2)
        pmid = p1 + 0.5 * (p2 - p1)
        m1 = _create_mesh(_get_cylinder_mesh(p1, pmid, radius=radius, segments=segments),
                          c1, solid=True, lines=False)
        m2 = _create_mesh(_get_cylinder_mesh(pmid, p2, radius=radius, segments=segments),
                          c2, solid=True, lines=False)
        return [m1, m2]

def create_ivol(vstruct,
                width=500, height=400, ssize=5,
                min_voxels=None,  max_voxels=None, **volargs):
    """

    Parameters
    ----------
    vstruct: dict
    width: int
    height: int
    ssize: int
    min_voxels : int
        minimum number of voxels in density cube
    max_voxels : int
        maximum number of voxels in density cube
    volargs: dict

    Returns
    -------

    Examples
    --------

    >>> from jsonextended import edict

    >>> dstruct = {
    ...  'type': 'repeat_density',
    ...  'dtype': 'charge',
    ...  'name': '',
    ...  'dcube':np.ones((3,3,3)),
    ...  'centre':[0,0,0],
    ...  'cell_vectors':{
    ...      'a':[2.,0,0],
    ...      'b':[0,2.,0],
    ...      'c':[0,0,2.]},
    ...   'color_bbox': 'black',
    ...   'transforms': []
    ... }
    >>> cstruct = {
    ...  'type': 'repeat_cell',
    ...  'name': '',
    ...  'centre':[0,0,0],
    ...  'cell_vectors':{
    ...      'a':[2.,0,0],
    ...      'b':[0,2.,0],
    ...      'c':[0,0,2.]},
    ...   'color_bbox': 'black',
    ...   'sites': [{
    ...         'label': "Fe",
    ...         'ccoord': [1,1,1],
    ...         'color_fill': "red",
    ...         'color_outline': None,
    ...         'transparency': 1,
    ...         'radius': 1,
    ...   }],
    ...   'bonds': [],
    ...   'transforms': []
    ... }
    >>> vstruct = {"elements": [dstruct, cstruct], "transforms": []}
    >>> new_struct, fig, controls = create_ivol(vstruct)

    >>> print(edict.apply(edict.filter_keys(new_struct, ["ivol"], list_of_dicts=True),
    ...                   "ivol", lambda x: [v.__class__.__name__ for v in x], list_of_dicts=True))
    {'elements': [{'ivol': ['Figure', 'Mesh']}, {'ivol': ['Mesh', 'Scatter']}]}

    """
    new_struct = apply_transforms(vstruct)
    bonds = compute_bonds(new_struct)#edict.filter_keyvals(vstruct, {"type": "repeat_cell"}, keep_siblings=True))

    # ivolume currently only allows one volume rendering per plot
    # voltypes = edict.filter_keyvals(vstructs,[('type','repeat_density')])
    vol_index = [i for i, el in enumerate(vstruct['elements'])
                 if el['type'] == 'repeat_density']
    assert len(vol_index) <= 1, "ipyvolume only allows one volume rendering per scene"

    p3.clear()
    fig = p3.figure(width=width, height=height, controls=True)
    fig.screen_capture_enabled = True

    # the volume rendering must be created first,
    # for appropriately scaled axis
    if vol_index:
        volstruct = new_struct['elements'][vol_index[0]]
        a = volstruct['cell_vectors']['a']
        b = volstruct['cell_vectors']['b']
        c = volstruct['cell_vectors']['c']
        centre = volstruct['centre']
        #print(centre)

        # convert dcube to cartesian
        out = cube_frac2cart(volstruct['dcube'], a, b, c, centre,
                             max_voxels=max_voxels, min_voxels=min_voxels, make_cubic=True)
        new_density, (xmin, ymin, zmin), (xmax, ymax, zmax) = out

        vol = p3.volshow(new_density, **volargs)

        if volstruct["color_bbox"] is not None:
            a = np.asarray(a)
            b = np.asarray(b)
            c = np.asarray(c)
            o = np.asarray(centre) - 0.5*(a+b+c)
            mesh = _create_mesh([o, o+a, o+b, o+c, o+a+b, o+a+c, o+c+b, o+a+b+c], color=volstruct["color_bbox"],
                                line_indices=[[0, 1], [0, 2], [0, 3], [2, 4], [2, 6],
                                              [1, 4], [1, 5], [3, 5], [3, 6], [7, 6], [7, 4], [7, 5]])
            vol = [vol, mesh]

        # todo better way of storing ivol components?
        volstruct['ivol'] = vol

        # appropriately scale axis
        p3.xlim(xmin, xmax)
        p3.ylim(ymin, ymax)
        p3.zlim(zmin, zmax)

    for element in new_struct['elements']:
        if element['type'] == 'repeat_density':
            continue
        elif element['type'] == 'repeat_cell':
            scatters = []
            if element["color_bbox"] is not None:
                a = np.asarray(element['cell_vectors']['a'])
                b = np.asarray(element['cell_vectors']['b'])
                c = np.asarray(element['cell_vectors']['c'])
                centre = element['centre']
                o = np.asarray(centre) - 0.5 * (a + b + c)
                mesh = _create_mesh([o, o + a, o + b, o + c, o + a + b, o + a + c, o + c + b, o + a + b + c],
                                    color=element["color_bbox"],
                                    line_indices=[[0, 1], [0, 2], [0, 3], [2, 4], [2, 6],
                                                  [1, 4], [1, 5], [3, 5], [3, 6], [7, 6], [7, 4], [7, 5]])
                scatters.append(mesh)

            for color, radius in set([(s['color_fill'], s['radius']) for s in element["sites"]]):

                scatter = edict.filter_keyvals(element, {"color_fill": color, "radius": radius}, "AND",
                                               keep_siblings=True, list_of_dicts=True)
                scatter = edict.combine_lists(scatter, ["sites"], deepcopy=False)
                scatter = edict.remove_keys(scatter, ["sites"], deepcopy=False)

                x, y, z = np.array(scatter['ccoord']).T
                s = p3.scatter(x, y, z, marker='sphere', size=ssize * radius, color=color)
                scatters.append(s)

            element['ivol'] = scatters
        elif element['type'] == 'repeat_poly':
            polys = []
            for poly in element['polys']:
                mesh = _create_mesh(poly, element['color'], element['solid'])
                polys.append(mesh)
            element['ivol'] = polys
        else:
            raise ValueError("unknown element type: {}".format(element['type']))

    for bond in bonds:
        p1, p2, c1, c2, radius = bond
        meshes = _create_bond(p1, p2, c1, c2, radius)

    # split up controls
    if vol_index:
        (level_ctrls, figbox,
         extractrl1, extractrl2) = p3.gcc().children
        controls = OrderedDict([('transfer function', [level_ctrls]),
                                ('lighting', [extractrl1, extractrl2])])
    else:
        # figbox = p3.gcc().children
        controls = {}

    return new_struct, fig, controls

# def create_ivol_control(vstruct, vparam, ctype, cparam='value', **ckwargs):
#     """
#
#     Parameters
#     ----------
#     vstruct
#     vparam
#     ctype
#     cparam
#     ckwargs
#
#     Returns
#     -------
#
#     """
#     import ipywidgets as widgets
#     assert 'ivol' in vstruct
#     ctrl = getattr(widgets, ctype)(**ckwargs)
#     widgets.jslink((vstruct['ivol'], vparam), (ctrl, cparam))
#
#     return ctrl

def add_controls(fig, ctrl_layout,
                 top=False):
    """

    Parameters
    ----------
    fig
    ctrl_layout
    top

    Returns
    -------

    """
    if not ctrl_layout:
        return fig
    if isinstance(ctrl_layout, list):
        layout = []
        for cntrl in ctrl_layout:
            if isinstance(cntrl, list):
                cntrl = widgets.HBox(cntrl)
            layout.append(cntrl)
        options = widgets.VBox(layout)
    else:
        tabs = OrderedDict()
        for name, cntrls in ctrl_layout.items():
            layout = []
            for cntrl in cntrls:
                if isinstance(cntrl, list):
                    cntrl = widgets.HBox(cntrl)
                layout.append(cntrl)
            tabs[name] = widgets.VBox(layout)
        options = widgets.Tab(children=tuple(tabs.values()))
        for i, name in enumerate(tabs):
            options.set_title(i, name)

    if top:
        return widgets.VBox([options, fig])
    else:
        return widgets.VBox([fig, options])

def ipy_style(fig, rescale=True, init_angles=None, hide_axes=False, fov=1):
    """ style the ipyvolume figure

    Parameters
    ----------
    fig: ipyvolume.Figure
    rescale: bool
        rescale axes to same aspect ratio
    init_angles: list
        initial angles [x, y, z]
    hide_axes: bool
    fov: float
    """
    if hide_axes:
        fig.xlabel = ' '
        fig.ylabel = ' '
        fig.zlabel = ' '
    fig.camera_fov = fov
    fig.style = {'axes': {'color': 'grey',
                          'label': {'color': 'grey'},
                          'ticklabel': {'color': 'grey'},
                          'visible': not hide_axes},
                 'background-color': 'white',
                 'box': {'visible': False}}

    x0, x1 = fig.xlim
    y0, y1 = fig.ylim
    z0, z1 = fig.zlim
    # currently ipyvolume doesn't scale; https://github.com/maartenbreddels/ipyvolume/issues/35
    if rescale:
        cube_side = max([abs(x1 - x0), abs(y1 - y0), abs(z1 - z0)])
        fig.xlim = [x0 - 0.5 * (cube_side - abs(x1 - x0)), x1 + 0.5 * (cube_side - abs(x1 - x0))]
        fig.ylim = [y0 - 0.5 * (cube_side - abs(y1 - y0)), y1 + 0.5 * (cube_side - abs(y1 - y0))]
        fig.zlim = [z0 - 0.5 * (cube_side - abs(z1 - z0)), z1 + 0.5 * (cube_side - abs(z1 - z0))]
    if init_angles:
        x, y, z = init_angles
        fig.anglex = x
        fig.angley = y
        fig.anglez = z

    return fig

max_stars_repo_name: chrisjsewell/ipyatom

max_issues_repo_name: chrisjsewell/ipyatom

max_forks_repo_name: chrisjsewell/ipyatom

reason: import numpy,from scipy