text:
"""Manage survey planning and schedule configuration data.

The normal usage is::

    >>> config = Configuration()
    >>> config.programs.BRIGHT.max_sun_altitude()
    <Quantity -13.0 deg>

Use dot notation to specify nodes in the configuration hieararchy and
function call notation to access terminal node values.

Terminal node values are first converted according to YAML rules. Strings
containing a number followed by valid astropy units are subsequently converted
to astropy quantities.  Strings of the form YYYY-MM-DD are converted to
datetime.date objects.

To change a value after the configuration has been loaded into memory use,
for example::

    >>> config.full_moon_nights.set_value(5)

Assigned values must have the appropriate converted types, for example::

    >>> import datetime
    >>> config.last_day.set_value(datetime.date(2024, 1, 1))
    >>> import astropy.units as u
    >>> config.location.temperature.set_value(-5 * u.deg_C)

The configuration is implemented as a singleton so the YAML file is only
loaded and parsed the first time a Configuration() is built.  Subsequent
calls to Configuration() always return the same object.
"""
from __future__ import print_function, division

import os
import re

import yaml

import astropy.units
import astropy.utils.data

# Extract a number from a string with optional leading and
# trailing whitespace.
_float_pattern = re.compile(
    r'\s*([-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)?)\s*')

class Node(object):
    """A single node of a configuration data structure.

    The purpose of this class is to allow nested dictionaries to be
    accessed using attribute dot notation, and to implement automatic
    conversion of strings of the form "<value> <units>" into corresponding
    astropy quantities.
    """
    def __init__(self, value, path=[]):
        self._path = path
        if isinstance(value, dict):
            # Remember our keys.
            self._keys = value.keys()
            # Recursively add sub-dictionaries as new child attributes.
            for name in self._keys:
                child_path = path + [name]
                self.__dict__[name] = Node(value[name], child_path)
        else:
            # Define the value of a leaf node.
            try:
                # Try to interpret as an astropy quantity with units.
                found_number = _float_pattern.match(value)
                if found_number:
                    number = float(found_number.group(1))
                    unit = value[found_number.end():]
                    try:
                        self._value = astropy.units.Quantity(number, unit)
                    except ValueError:
                        raise ValueError(
                            'Invalid unit for {0}: {1}'
                            .format('.'.join(self._path), unit))
                else:
                    self._value = value
            except (TypeError, ValueError):
                self._value = value

    @property
    def path(self):
        """Return the full path to this node using dot notation.
        """
        return '.'.join(self._path)

    @property
    def keys(self):
        """Return the list of keys for a non-leaf node or raise a RuntimeError
        for a terminal node.
        """
        try:
            return self._keys
        except AttributeError:
            raise RuntimeError(
                '{0} is a terminal config node.'.format(self.path))

    def __call__(self):
        """Return a terminal node's value or raise a RuntimeError for
        a non-terminal node.
        """
        try:
            return self._value
        except AttributeError:
            raise RuntimeError(
                '{0} is a non-terminal config node.'.format(self.path))

    def set_value(self, new_value):
        """Set a terminal node's value or raise a RuntimeError for
        a non-terminal node.
        """
        try:
            old_value = self._value
            if not isinstance(new_value, type(old_value)):
                raise RuntimeError(
                'new type ({}) does not match old type ({}).'
                .format(type(new_value), type(old_value)))
            self._value = new_value
        except AttributeError:
            raise RuntimeError(
                '{0} is a non-terminal config node.'.format(self.path))

class Configuration(Node):
    """Top-level configuration data node.
    """
    __instance = None

    @staticmethod
    def reset():
        """Forget our singleton instance.  Mainly intended for unit tests."""
        Configuration.__instance = None

    @staticmethod
    def _get_full_path(file_name):
        # Locate the config file in our pkg data/ directory if no path is given.
        if os.path.split(file_name)[0] == '':
            full_path = astropy.utils.data._find_pkg_data_path(
                os.path.join('data', file_name))
        else:
            full_path = file_name
        return full_path

    def __new__(cls, file_name=None):
        """Implement a singleton access pattern.
        """
        if Configuration.__instance is None:
            Configuration.__instance = object.__new__(cls)
            Configuration.__instance._initialize(file_name)
        elif file_name is not None and file_name != Configuration.__instance.file_name:
            raise RuntimeError('Configuration already loaded from {0}'
                               .format(Configuration.__instance.file_name))
        return Configuration.__instance

    def __init__(self, file_name='config.yaml'):
        """Return the unique configuration object for this session.

        The configuration will be loaded from the specified file when this
        constructor is called for the first time.  Subsequent calls with
        a different file name will result in a RuntimeError.

        Parameters
        ----------
        file_name : string
            Name of a YAML file including a valid YAML extension.  The file
            is assumed to be under this package's data/ directory unless
            a path is included (either relative or absolute).
        """
        pass

    def _initialize(self, file_name=None):
        """Initialize a configuration data structure from a YAML file.
        """
        if file_name is None:
            file_name = 'config.yaml'

        # Remember the file name since it is not allowed to change.
        self.file_name = file_name

        full_path = self._get_full_path(file_name)

        # Validate that all mapping keys are valid python identifiers
        # and that there are no embedded sequences.
        valid_key = re.compile(r'^[a-zA-Z_][a-zA-Z0-9_]*\Z')
        with open(full_path) as f:
            next_value_is_key = False
            for token in yaml.scan(f):
                if isinstance(
                    token,
                    (yaml.BlockSequenceStartToken,
                     yaml.FlowSequenceStartToken)):
                    raise RuntimeError('Config sequences not supported.')
                if next_value_is_key:
                    if not isinstance(token, yaml.ScalarToken):
                        raise RuntimeError(
                            'Invalid config key type: {0}'.format(token))
                    if not valid_key.match(token.value):
                        raise RuntimeError(
                            'Invalid config key name: {0}'.format(token.value))
                next_value_is_key = isinstance(token, yaml.KeyToken)

        # Load the config data into memory.
        with open(full_path) as f:
            Node.__init__(self, yaml.safe_load(f))

        # Output path is not set until it is first used.
        self._output_path = None

    def set_output_path(self, output_path):
        """Set the output directory for relative paths.

        The path must exist when this method is called. Called by
        :meth:`get_path` for a non-absolute path. This method updates the
        configuration output_path value.

        Parameters
        ----------
        output_path : str
            A path possibly including environment variables enclosed in {...}
            that will be substituted from the current environment.

        Raises
        ------
        ValueError
            Path uses undefined environment variable or does not exist.
        """
        try:
            self._output_path = output_path.format(**os.environ)
        except KeyError as e:
            raise ValueError(
                'Environment variable not set for output_path: {0}'.format(e))
        if not os.path.isdir(self._output_path):
            raise ValueError(
                'Non-existent output_path: {0}'.format(self._output_path))
        # Update our config node.
        self.output_path._value = output_path

    def get_path(self, name):
        """Prepend this configuration's output_path to non-absolute paths.

        Configured by the ``output_path`` node and :meth:`set_output_path`.

        An absolute path is returned immediately so an environment variable
        used in output_path only needs to be defined if relative paths
        are used.

        Parameters
        ----------
        name : str
            Absolute or relative path name, which does not need to exist yet.

        Returns
        -------
        str
            Path name to use. Relative path names will have our output_path
            prepended.  Absolute path names will be unchanged.
        """
        if os.path.isabs(name):
            return name
        if self._output_path is None:
            self.set_output_path(self.output_path())
        return os.path.join(self._output_path, name)

max_stars_repo_name: desi-bgs/desisurvey

max_issues_repo_name: desi-bgs/desisurvey

max_forks_repo_name: desi-bgs/desisurvey

reason: import astropy


text:
from itertools import cycle
from json import load

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import rc

with open('bgm_anime_dataset.json', 'r', encoding='utf8') as f:
    data = load(f)

scores = np.array(
    [bangumi['rating']['score'] for bangumi in data],
    dtype=np.float64
)

count = scores.size
mean = np.mean(scores, dtype=np.float64)
median = np.median(scores)
min_score = np.min(scores)
max_score = np.max(scores)

n, bins = np.histogram(scores, bins=np.arange(-0.05, 10.05, 0.1))
hist = np.vstack((bins[:-1] + 0.05, n))

for i in range(hist.shape[1]):
    print('%.1f %3d' % (hist[0, i], hist[1, i]))

rc('font', family='FZZhunYuan-M02', size=14)

fig, ax = plt.subplots()
for i, c in zip(range(20), cycle(('#fca2ae', '#f6c2d0'))):
    ax.bar(hist[0, i*5:i*5+5], hist[1, i*5:i*5+5], width=0.05, color=c)
ax.annotate(
    '最高：%d @ %.1f' % (hist[1, 66], hist[0, 66]),
    (hist[0, 66], hist[1, 66]),
    xycoords='data', xytext=(0.25, 0.56), textcoords='axes fraction',
    arrowprops={'arrowstyle': '->'}
)
ax.annotate(
    '局部最低：%d @ %.1f' % (hist[1, 69], hist[0, 69]),
    (hist[0, 69], hist[1, 69]),
    xycoords='data', xytext=(0.25, 0.43), textcoords='axes fraction',
    arrowprops={'arrowstyle': '->'}
)
ax.annotate(
    '平均：%.3f' % (mean),
    (mean, 0),
    xycoords='data', xytext=(0.58, -0.12), textcoords='axes fraction',
    arrowprops={'arrowstyle': '->'}
)
ax.annotate(
    '中位：%.1f' % (median),
    (median, 0),
    xycoords='data', xytext=(0.68, -0.12), textcoords='axes fraction',
    arrowprops={'arrowstyle': '->'}
)
ax.annotate(
    '最低：%.1f' % (min_score),
    (min_score, 0),
    xycoords='data', xytext=(0.125, -0.12), textcoords='axes fraction',
    arrowprops={'arrowstyle': '->'}
)
ax.annotate(
    '最高：%.1f' % (max_score),
    (max_score, 0),
    xycoords='data', xytext=(0.87, -0.12), textcoords='axes fraction',
    arrowprops={'arrowstyle': '->'}
)
ax.grid(True, axis='y')
ax.set_xticks(np.arange(0, 10 + 0.1, 0.5))
ax.set_xlabel('评分')
ax.set_ylabel('动画数')
ax.set_title(f'Bangumi 动画评分分布    总计{count}部')
fig.tight_layout()
plt.show()

ext: py

reason: import numpy


text:
#!/usr/bin/env python
# -*- coding: utf-8 -*-
# File: DCGAN.py
# Author: Yuxin Wu

import argparse
import glob
import numpy as np
import os
import tensorflow as tf

from tensorpack import *
from tensorpack.tfutils.scope_utils import auto_reuse_variable_scope
from tensorpack.utils.viz import stack_patches

from GAN import GANModelDesc, GANTrainer, RandomZData

"""
1. Download the 'aligned&cropped' version of CelebA dataset
   from http://mmlab.ie.cuhk.edu.hk/projects/CelebA.html

2. Start training:
    ./DCGAN-CelebA.py --data /path/to/img_align_celeba/ --crop-size 140
    Generated samples will be available through tensorboard

3. Visualize samples with an existing model:
    ./DCGAN-CelebA.py --load path/to/model --sample

You can also train on other images (just use any directory of jpg files in
`--data`). But you may need to change the preprocessing.

A pretrained model on CelebA is at http://models.tensorpack.com/GAN/
"""

class Model(GANModelDesc):
    def __init__(self, shape, batch, z_dim):
        self.shape = shape
        self.batch = batch
        self.zdim = z_dim

    def inputs(self):
        return [tf.TensorSpec((None, self.shape, self.shape, 3), tf.float32, 'input')]

    def generator(self, z):
        """ return an image generated from z"""
        nf = 64
        l = FullyConnected('fc0', z, nf * 8 * 4 * 4, activation=tf.identity)
        l = tf.reshape(l, [-1, 4, 4, nf * 8])
        l = BNReLU(l)
        with argscope(Conv2DTranspose, activation=BNReLU, kernel_size=4, strides=2):
            l = Conv2DTranspose('deconv1', l, nf * 4)
            l = Conv2DTranspose('deconv2', l, nf * 2)
            l = Conv2DTranspose('deconv3', l, nf)
            l = Conv2DTranspose('deconv4', l, 3, activation=tf.identity)
            l = tf.tanh(l, name='gen')
        return l

    @auto_reuse_variable_scope
    def discriminator(self, imgs):
        """ return a (b, 1) logits"""
        nf = 64
        with argscope(Conv2D, kernel_size=4, strides=2):
            l = (LinearWrap(imgs)
                 .Conv2D('conv0', nf, activation=tf.nn.leaky_relu)
                 .Conv2D('conv1', nf * 2)
                 .BatchNorm('bn1')
                 .tf.nn.leaky_relu()
                 .Conv2D('conv2', nf * 4)
                 .BatchNorm('bn2')
                 .tf.nn.leaky_relu()
                 .Conv2D('conv3', nf * 8)
                 .BatchNorm('bn3')
                 .tf.nn.leaky_relu()
                 .FullyConnected('fct', 1)())
        return l

    def build_graph(self, image_pos):
        image_pos = image_pos / 128.0 - 1

        z = tf.random_uniform([self.batch, self.zdim], -1, 1, name='z_train')
        z = tf.placeholder_with_default(z, [None, self.zdim], name='z')

        with argscope([Conv2D, Conv2DTranspose, FullyConnected],
                      kernel_initializer=tf.truncated_normal_initializer(stddev=0.02)):
            with tf.variable_scope('gen'):
                image_gen = self.generator(z)
            tf.summary.image('generated-samples', image_gen, max_outputs=30)
            with tf.variable_scope('discrim'):
                vecpos = self.discriminator(image_pos)
                vecneg = self.discriminator(image_gen)

        self.build_losses(vecpos, vecneg)
        self.collect_variables()

    def optimizer(self):
        lr = tf.get_variable('learning_rate', initializer=2e-4, trainable=False)
        return tf.train.AdamOptimizer(lr, beta1=0.5, epsilon=1e-3)

def get_augmentors():
    augs = []
    if args.load_size:
        augs.append(imgaug.Resize(args.load_size))
    if args.crop_size:
        augs.append(imgaug.CenterCrop(args.crop_size))
    augs.append(imgaug.Resize(args.final_size))
    return augs

def get_data():
    assert args.data
    imgs = glob.glob(args.data + '/*.jpg')
    ds = ImageFromFile(imgs, channel=3, shuffle=True)
    ds = AugmentImageComponent(ds, get_augmentors())
    ds = BatchData(ds, args.batch)
    ds = PrefetchDataZMQ(ds, 5)
    return ds

def sample(model, model_path, output_name='gen/gen'):
    pred = PredictConfig(
        session_init=get_model_loader(model_path),
        model=model,
        input_names=['z'],
        output_names=[output_name, 'z'])
    pred = SimpleDatasetPredictor(pred, RandomZData((100, args.z_dim)))
    for o in pred.get_result():
        o = o[0] + 1
        o = o * 128.0
        o = np.clip(o, 0, 255)
        o = o[:, :, :, ::-1]
        stack_patches(o, nr_row=10, nr_col=10, viz=True)

def get_args(default_batch=128, default_z_dim=100):
    parser = argparse.ArgumentParser()
    parser.add_argument('--gpu', help='comma separated list of GPU(s) to use.')
    parser.add_argument('--load', help='load model')
    parser.add_argument('--sample', action='store_true', help='view generated examples')
    parser.add_argument('--data', help='a jpeg directory')
    parser.add_argument('--load-size', help='size to load the original images', type=int)
    parser.add_argument('--crop-size', help='crop the original images', type=int)
    parser.add_argument(
        '--final-size', default=64, type=int,
        help='resize to this shape as inputs to network')
    parser.add_argument('--z-dim', help='hidden dimension', type=int, default=default_z_dim)
    parser.add_argument('--batch', help='batch size', type=int, default=default_batch)
    global args
    args = parser.parse_args()
    if args.gpu:
        os.environ['CUDA_VISIBLE_DEVICES'] = args.gpu
    return args

if __name__ == '__main__':
    args = get_args()
    M = Model(shape=args.final_size, batch=args.batch, z_dim=args.z_dim)
    if args.sample:
        sample(M, args.load)
    else:
        logger.auto_set_dir()
        GANTrainer(
            input=QueueInput(get_data()),
            model=M).train_with_defaults(
            callbacks=[ModelSaver()],
            steps_per_epoch=300,
            max_epoch=200,
            session_init=SaverRestore(args.load) if args.load else None
        )

max_stars_repo_name: bzamecnik/tensorpack


text:
import torch
import numpy as np
import torch.backends.cudnn as cudnn
import os
from tqdm import tqdm
from net.models import deeplabv3plus
from sklearn.metrics import accuracy_score
from net import loss
# import matplotlib as mpl
# mpl.use('Agg')
import matplotlib.pyplot as plt
from apex import amp
from tensorboardX import SummaryWriter
from dataset.my_datasets import MyDataSet_seg, MyValDataSet_seg
from torch.utils import data

model_urls = {'deeplabv3plus_xception': 'models/deeplabv3plus_xception_VOC2012_epoch46_all.pth'}

INPUT_SIZE = '224, 224'
w, h = map(int, INPUT_SIZE.split(','))
LEARNING_RATE = 0.0001
MOMENTUM = 0.9
POWER = 0.9
WEIGHT_DECAY = 0.0005
NUM_CLASSES = 2
TRAIN_NUM = 2000
BATCH_SIZE = 16
EPOCH = 500
STEPS = (TRAIN_NUM/BATCH_SIZE)*EPOCH
FP16 = True
NAME = 'DR_CoarseSN/'

def lr_poly(base_lr, iter, max_iter, power):
    return base_lr * ((1 - float(iter) / max_iter) ** (power))

def adjust_learning_rate(optimizer, i_iter):
    lr = lr_poly(LEARNING_RATE, i_iter, STEPS, POWER)
    optimizer.param_groups[0]['lr'] = lr
    return lr

def val_mode_seg(valloader, model, path, epoch):
    dice = []
    sen = []
    spe = []
    acc = []
    jac_score = []
    for index, batch in enumerate(valloader):

        data, mask, name = batch
        data = data.cuda()
        mask = mask[0].data.numpy()
        val_mask = np.int64(mask > 0)
        # print(name)

        model.eval()
        with torch.no_grad():
            pred = model(data)

        pred = torch.softmax(pred, dim=1).cpu().data.numpy()
        pred_arg = np.argmax(pred[0], axis=0)

        #y_pred
        y_true_f = val_mask.reshape(val_mask.shape[0]*val_mask.shape[1], order='F')
        y_pred_f = pred_arg.reshape(pred_arg.shape[0]*pred_arg.shape[1], order='F')

        intersection = np.float(np.sum(y_true_f * y_pred_f))
        dice.append((2. * intersection) / (np.sum(y_true_f) + np.sum(y_pred_f)))
        sen.append(intersection / np.sum(y_true_f))
        intersection0 = np.float(np.sum((1 - y_true_f) * (1 - y_pred_f)))
        spe.append(intersection0 / np.sum(1 - y_true_f))
        acc.append(accuracy_score(y_true_f, y_pred_f))
        jac_score.append(intersection / (np.sum(y_true_f) + np.sum(y_pred_f) - intersection))

        if index in [100]:
            fig = plt.figure()
            ax = fig.add_subplot(131)
            ax.imshow(data[0].cpu().data.numpy().transpose(1, 2, 0))
            ax.axis('off')
            ax = fig.add_subplot(132)
            ax.imshow(mask)
            ax.axis('off')
            ax = fig.add_subplot(133)
            ax.imshow(pred_arg)
            ax.axis('off')
            fig.suptitle('RGB image,ground truth mask, predicted mask',fontsize=6)
            fig.savefig(path + name[0][:-4] + '_e' + str(epoch) + '.png', dpi=200, bbox_inches='tight')
            ax.cla()
            fig.clf()
            plt.close()

    return np.array(acc), np.array(dice), np.array(sen), np.array(spe), np.array(jac_score)

def Jaccard(pred_arg, mask):
    pred_arg = np.argmax(pred_arg.cpu().data.numpy(), axis=1)
    mask = mask.cpu().data.numpy()

    y_true_f = mask.reshape(mask.shape[0] * mask.shape[1] * mask.shape[2], order='F')
    y_pred_f = pred_arg.reshape(pred_arg.shape[0] * pred_arg.shape[1] * pred_arg.shape[2], order='F')

    intersection = np.float(np.sum(y_true_f * y_pred_f))
    jac_score = intersection / (np.sum(y_true_f) + np.sum(y_pred_f) - intersection)

    return jac_score

def main():
    """Create the network and start the training."""
    writer = SummaryWriter('models/' + NAME)

    cudnn.enabled = True

    ############# Create coarse segmentation network
    model = deeplabv3plus(num_classes=NUM_CLASSES)
    optimizer = torch.optim.Adam(model.parameters(), lr=LEARNING_RATE)
    model.cuda()
    if FP16 is True:
        model, optimizer = amp.initialize(model, optimizer, opt_level="O1")
    model = torch.nn.DataParallel(model)

    ############# Load pretrained weights
    pretrained_dict = torch.load(model_urls['deeplabv3plus_xception'])
    net_dict = model.state_dict()
    pretrained_dict = {k: v for k, v in pretrained_dict.items() if (k in net_dict) and (v.shape == net_dict[k].shape)}
    net_dict.update(pretrained_dict)
    model.load_state_dict(net_dict)

    print(len(net_dict))
    print(len(pretrained_dict))

    model.train()
    model.float()

    DR_loss = loss.Fusin_Dice_rank()

    cudnn.benchmark = True

    ############# Load training and validation data
    data_train_root = 'dataset/seg_data/Training_resize_seg/'
    data_train_list = 'dataset/ISIC/Training_seg.txt'
    trainloader = data.DataLoader(MyDataSet_seg(data_train_root, data_train_list, crop_size=(w, h)),
                                  batch_size=BATCH_SIZE, shuffle=True, num_workers=8, pin_memory=True)

    data_val_root = 'dataset/seg_data/ISIC-2017_Validation_Data/'
    data_val_list = 'dataset/ISIC/Validation_seg.txt'
    valloader = data.DataLoader(MyValDataSet_seg(data_val_root, data_val_list), batch_size=1, shuffle=False, num_workers=8,
                                pin_memory=True)

    path = 'models/' + NAME
    if not os.path.isdir(path):
        os.mkdir(path)
    f_path = path + 'outputxx.txt'

    val_jac = []

    ############# Start the training
    for epoch in range(EPOCH):

        train_loss_D = []
        train_loss_R = []
        train_loss_total = []
        train_jac = []

        for i_iter, batch in tqdm(enumerate(trainloader)):

            step = (TRAIN_NUM/BATCH_SIZE)*epoch+i_iter

            images, labels, name = batch
            images = images.cuda()
            labels = labels.cuda().squeeze(1)

            optimizer.zero_grad()
            lr = adjust_learning_rate(optimizer, step)

            model.train()
            preds = model(images)

            loss_D, loss_R = DR_loss(preds, labels)
            term = loss_D + 0.05 * loss_R

            if FP16 is True:
                with amp.scale_loss(term, optimizer) as scaled_loss:
                    scaled_loss.backward()
            else:
                term.backward()

            optimizer.step()

            writer.add_scalar('learning_rate', lr, step)
            writer.add_scalar('loss', term.cpu().data.numpy(), step)

            train_loss_D.append(loss_D.cpu().data.numpy())
            train_loss_R.append(loss_R.cpu().data.numpy())
            train_loss_total.append(term.cpu().data.numpy())
            train_jac.append(Jaccard(preds, labels))

        print("train_epoch%d: lossTotal=%f, lossDice=%f, lossRank=%f, Jaccard=%f \n" % (epoch, np.nanmean(train_loss_total), np.nanmean(train_loss_D), np.nanmean(train_loss_R), np.nanmean(train_jac)))

        ############# Start the validation
        [vacc, vdice, vsen, vspe, vjac_score] = val_mode_seg(valloader, model, path, epoch)
        line_val = "val%d: vacc=%f, vdice=%f, vsensitivity=%f, vspecifity=%f, vjac=%f \n" % \
                (epoch, np.nanmean(vacc), np.nanmean(vdice), np.nanmean(vsen), np.nanmean(vspe),
                    np.nanmean(vjac_score))

        print(line_val)
        f = open(f_path, "a")
        f.write(line_val)

        ############# Plot val curve
        val_jac.append(np.nanmean(vjac_score))
        plt.figure()
        plt.plot(val_jac, label='val jaccard', color='blue', linestyle='--')
        plt.legend(loc='best')

        plt.savefig(os.path.join(path, 'jaccard.png'))
        plt.clf()
        plt.close()
        plt.show()

        plt.close('all')

        writer.add_scalar('val_Jaccard', np.nanmean(vjac_score), epoch)

        ############# Save network
        torch.save(model.state_dict(), path + 'CoarseSN_e' + str(epoch) + '.pth')

if __name__ == '__main__':
    main()

max_stars_repo_name: jianpengz/MB-DCNN

max_issues_repo_name: DLWK/MB-DCNN

max_forks_repo_name: DLWK/MB-DCNN


text:
import math
from ..lgssm import LGSSM
import jax.numpy as np
from jax import vmap
from jax.lax import scan, cond
from jax.numpy import concatenate
from jax.scipy.linalg import cho_factor, cho_solve

LOG2PI = math.log(2 * math.pi)
INV2PI = (2 * math.pi) ** -1

def first_filtering_element(m0, P0, F, Q, H, R, y):
    """
    Equation 11 in "GPR in Logarithmic Time"
    """

    def _res_nan():
        A = np.zeros_like(F)
        b = m0
        C = P0
        eta = np.zeros_like(m0)
        J = np.zeros_like(F)

        return A, b, C, J, eta

    def _res_not_nan():
        S1 = H @ (P0 @ H.T) + R
        S1_chol = cho_factor(S1)
        K1t = cho_solve(S1_chol, H @ P0)

        A = np.zeros_like(F)
        b = m0 - K1t.T @ (y - H @ m0)
        C = P0 - (K1t.T @ S1) @ K1t

        S = H @ (Q @ H.T) + R
        chol = cho_factor(S)

        HF = H @ F
        eta = HF.T @ np.squeeze(cho_solve(chol, np.expand_dims(y, 1)))
        J = HF.T @ cho_solve(chol, HF)
        return A, b, C, J, eta

    res = cond(np.isnan(y), _res_nan, _res_not_nan)
    return res

def _generic_filtering_element_nan(F, Q):
    A = F
    b = np.zeros(F.shape[:2], dtype=F.dtype)
    C= Q
    eta = np.zeros(F.shape[:2], dtype=F.dtype)
    J = np.zeros_like(F)

    return A, b, C, J, eta

def _generic_filtering_element(F, Q, H, R, y):
    """
    Equation 10 in "GPR in Logarithmic Time"
    """

    S = H @ (Q @ H.T) + R

    chol = cho_factor(S)
    Kt = cho_solve(chol, H @ Q)
    A = F - (Kt.T @ H) @ F

    b= Kt.T @ y
    C = Q - (Kt.T @ H) @ Q

    HF = H @ F
    eta = HF.T @ np.squeeze(cho_solve(chol, np.expand_dims(y, 1)))
    J = HF.T @ cho_solve(chol, HF)

    return A, b, C, J, eta

def make_associative_scan_filtering_elements(m0, P0, Fs,  Qs, H, R, observations):
    init_res = first_filtering_element(m0, P0, Fs[0], Qs[0], H, R observations[0])

max_stars_repo_name: AntixK/parallel-SSVI

max_issues_repo_name: AntixK/parallel-SSVI

max_forks_repo_name: AntixK/parallel-SSVI

reason: import jax,from jax


text:
from __future__ import absolute_import
from __future__ import division
from __future__ import print_function
from __future__ import unicode_literals

import os
import numpy as np
from six.moves import xrange

import argparse
import logging
import tensorflow as tf
from tensorflow.python.platform import flags

from cleverhans.utils_mnist import data_mnist
from cleverhans.utils import to_categorical
from cleverhans.utils import set_log_level
from cleverhans.utils import build_model_save_path, parse_model_settings
from cleverhans.utils_tf import model_train, model_eval, batch_eval, tf_model_load
from cleverhans.attacks import FastGradientMethod
from cleverhans.attacks_tf import jacobian_graph, jacobian_augmentation

from cleverhans_tutorials.tutorial_models import make_basic_cnn, make_basic_binary_cnn, MLP
from cleverhans_tutorials.tutorial_models import Flatten, Linear, SimpleLinear, ReLU, SoftmaxT1

FLAGS = flags.FLAGS
MAX_BATCH_SIZE = 100
TEST_SET_SIZE = 1000

INIT_T = 1.0
# enum adversarial training types
ADVERSARIAL_TRAINING_MADRYETAL = 1
ADVERSARIAL_TRAINING_FGSM = 2
MAX_EPS = 0.3

def setup_tutorial():
    """
    Helper function to check correct configuration of tf for tutorial
    :return: True if setup checks completed
    """

    # Set TF random seed to improve reproducibility
    tf.set_random_seed(1234)

    return True

def prep_bbox(sess, logits_scalar, x, y, X_train, Y_train, X_test, Y_test,
              img_rows, img_cols, channels, nb_epochs, batch_size, learning_rate,
              rng, phase=None, binary=False, scale=False, nb_filters=64,
              model_path=None, adv=0, delay=0, eps=0.3):
    """
    Define and train a model that simulates the "remote"
    black-box oracle described in the original paper.
    :param sess: the TF session
    :param x: the input placeholder for MNIST
    :param y: the ouput placeholder for MNIST
    :param X_train: the training data for the oracle
    :param Y_train: the training labels for the oracle
    :param X_test: the testing data for the oracle
    :param Y_test: the testing labels for the oracle
    :param nb_epochs: number of epochs to train model
    :param batch_size: size of training batches
    :param learning_rate: learning rate for training
    :param rng: numpy.random.RandomState
    :return:
    """

    # Define TF model graph (for the black-box model)
    save = False
    train_from_scratch = False

    if model_path is not None:
        if os.path.exists(model_path):
            # check for existing model in immediate subfolder
            if any(f.endswith('.meta') for f in os.listdir(model_path)):
                binary, scale, nb_filters, batch_size, learning_rate, nb_epochs, adv = parse_model_settings(
                    model_path)
                train_from_scratch = False
            else:
                model_path = build_model_save_path(
                    model_path, binary, batch_size, nb_filters, learning_rate, nb_epochs, adv, delay, scale)
                print(model_path)
                save = True
                train_from_scratch = True
    else:
        train_from_scratch = True  # train from scratch, but don't save since no path given

    if binary:
        if scale:
            #from cleverhans_tutorials.tutorial_models import make_scaled_binary_cnn
            # model = make_scaled_binary_cnn(phase, 'bb_binsc_', input_shape=(
            from cleverhans_tutorials.tutorial_models import make_scaled_binary_rand_cnn
            model = make_scaled_binary_rand_cnn(phase, logits_scalar, 'bb_binsc_', input_shape=(
                None, img_rows, img_cols, channels), nb_filters=nb_filters)
        else:
            from cleverhans_tutorials.tutorial_models import make_basic_binary_cnn
            model = make_basic_binary_cnn(phase, logits_scalar, 'bb_bin_', input_shape=(
                None, img_rows, img_cols, channels), nb_filters=nb_filters)
    else:
        from cleverhans_tutorials.tutorial_models import make_basic_cnn
        model = make_basic_cnn(phase, logits_scalar, 'bb_fp_', input_shape=(
            None, img_rows, img_cols, channels), nb_filters=nb_filters)

    preds = model(x, reuse=False)
    print("Defined TensorFlow model graph.")

    def evaluate():
        # Print out the accuracy on legitimate data
        eval_params = {'batch_size': batch_size}
        acc = model_eval(
            sess, x, y, preds, X_test, Y_test, phase=phase, args=eval_params)
        print('Test accuracy of black-box on legitimate test '
              'examples: %.4f' % acc)

    # Train an MNIST model
    train_params = {
        'binary': binary,
        'nb_epochs': nb_epochs,
        'batch_size': batch_size,
        'learning_rate': learning_rate,
        'loss_name': 'bb train loss',
        'filename': 'bb_model',
        'train_scope': 'bb_model',
        'reuse_global_step': False,
        'is_training': True
    }

    if adv != 0:
        if adv == ADVERSARIAL_TRAINING_MADRYETAL:
            from cleverhans.attacks import MadryEtAl
            nb_iter = 20
            train_attack_params = {'eps': MAX_EPS, 'eps_iter': 0.01,
                                   'nb_iter': nb_iter}
            train_attacker = MadryEtAl(model, sess=sess)

        if adv == ADVERSARIAL_TRAINING_FGSM:
            from cleverhans.attacks import FastGradientMethod
            train_attacker = FastGradientMethod(model, back='tf', sess=sess)

        # create the adversarial trainer
        train_attack_params.update({'clip_min': 0., 'clip_max': 1.})
        adv_x_train = train_attacker.generate(x, phase, **train_attack_params)
        preds_adv = model.get_probs(adv_x_train)

    if train_from_scratch:
        if save:
            train_params.update({'log_dir': model_path})
            if adv and delay > 0:
                train_params.update({'nb_epochs': delay})

        # do clean training for 'nb_epochs' or 'delay' epochs
        model_train(sess, x, y, preds, X_train, Y_train, phase=phase,
                    evaluate=evaluate, args=train_params, save=save, rng=rng)

        # optionally do additional adversarial training
        if adv:
            print("Adversarial training for %d epochs" % (nb_epochs - delay))
            train_params.update({'nb_epochs': nb_epochs - delay})
            train_params.update({'reuse_global_step': True})
            model_train(sess, x, y, preds, X_train, Y_train, phase=phase,
                        predictions_adv=preds_adv, evaluate=evaluate, args=train_params,
                        save=save, rng=rng)
    else:
        tf_model_load(sess, model_path)
        print('Restored model from %s' % model_path)

    accuracy = evaluate()

    return model, preds, accuracy, model_path

def substitute_model(img_rows=28, img_cols=28, nb_classes=10):
    """
    Defines the model architecture to be used by the substitute. Use
    the example model interface.
    :param img_rows: number of rows in input
    :param img_cols: number of columns in input
    :param nb_classes: number of classes in output
    :return: tensorflow model
    """
    input_shape = (None, img_rows, img_cols, 1)

    # Define a fully connected model (it's different than the black-box)
    layers = [Flatten(),
              SimpleLinear(200),
              ReLU(),
              SimpleLinear(200),
              ReLU(),
              SimpleLinear(nb_classes),
              SoftmaxT1()]

    return MLP(layers, input_shape)

def train_sub(sess, logits_scalar, x, y, bbox_preds, X_sub, Y_sub, nb_classes,
              nb_epochs_s, batch_size, learning_rate, data_aug, lmbda,
              rng, binary=False, phase=None, model_path=None):
    """
    This function creates the substitute by alternatively
    augmenting the training data and training the substitute.
    :param sess: TF session
    :param x: input TF placeholder
    :param y: output TF placeholder
    :param bbox_preds: output of black-box model predictions
    :param X_sub: initial substitute training data
    :param Y_sub: initial substitute training labels
    :param nb_classes: number of output classes
    :param nb_epochs_s: number of epochs to train substitute model
    :param batch_size: size of training batches
    :param learning_rate: learning rate for training
    :param data_aug: number of times substitute training data is augmented
    :param lmbda: lambda from arxiv.org/abs/1602.02697
    :param rng: numpy.random.RandomState instance
    :param phase: placeholder for batch_norm phase (training or testing)
    :param phase_val: True if training, False if testing
    :return:
    """
    # Define TF model graph (for the black-box model)
    model_sub = substitute_model()
    preds_sub = model_sub(x)
    print("Defined TensorFlow model graph for the substitute.")

    # Define the Jacobian symbolically using TensorFlow
    grads = jacobian_graph(preds_sub, x, nb_classes)

    train_params = {
        'binary': False,
        'nb_epochs': nb_epochs_s,
        'batch_size': batch_size,
        'learning_rate': learning_rate,
        'filename': 'sub_model',
        'train_scope': 'sub_model',
        'reuse_global_step': False,
        'is_training': True
    }

    # Train the substitute and augment dataset alternatively
    for rho in xrange(data_aug):
        print("Substitute training epoch #" + str(rho))

        if rho > 0:
            train_params.update({'reuse_global_step': True})
        if model_path is not None:
            train_params.update({'log_dir': model_path})
            model_train(sess, x, y, preds_sub, X_sub, to_categorical(Y_sub),
                        phase=phase, save=True, init_all=False, args=train_params,
                        rng=rng)
        else:
            model_train(sess, x, y, preds_sub, X_sub, to_categorical(Y_sub),
                        phase=phase, init_all=False, args=train_params,
                        rng=rng)

        # If we are not at last substitute training iteration, augment dataset
        if rho < data_aug - 1:
            print("Augmenting substitute training data.")
            # Perform the Jacobian augmentation
            X_sub = jacobian_augmentation(sess, x, X_sub, Y_sub, grads, lmbda)

            print("Labeling substitute training data.")
            # Label the newly generated synthetic points using the black-box
            Y_sub = np.hstack([Y_sub, Y_sub])
            X_sub_prev = X_sub[int(len(X_sub) / 2):]
            eval_params = {'batch_size': batch_size}
            bbox_val = batch_eval(sess, [x], [bbox_preds], [X_sub_prev],
                                  feed={phase: False}, args=eval_params)[0]
            # Note here that we take the argmax because the adversary
            # only has access to the label (not the probabilities) output
            # by the black-box model
            Y_sub[int(len(X_sub) / 2):] = np.argmax(bbox_val, axis=1)

    return model_sub, preds_sub

def mnist_blackbox(train_start=0, train_end=60000, test_start=0,
                   test_end=10000, nb_classes=10, batch_size=128, nb_samples=10, nb_filters=64,
                   eps=0.3, learning_rate=0.001, nb_epochs=10, holdout=150, data_aug=6,
                   nb_epochs_s=10, lmbda=0.1, binary=False, scale=False, model_path=None,
                   targeted=False, data_dir=None, adv=False, delay=0):
    """
    MNIST tutorial for the black-box attack from arxiv.org/abs/1602.02697
    :param train_start: index of first training set example
    :param train_end: index of last training set example
    :param test_start: index of first test set example
    :param test_end: index of last test set example
    :return: a dictionary with:
             * black-box model accuracy on test set
             * substitute model accuracy on test set
             * black-box model accuracy on adversarial examples transferred
               from the substitute model
    """

    # Set logging level to see debug information
    set_log_level(logging.DEBUG)

    # Dictionary used to keep track and return key accuracies
    accuracies = {}

    # Perform tutorial setup
    assert setup_tutorial()

    # Create TF session
    sess = tf.Session()

    # Get MNIST data
    X_train, Y_train, X_test, Y_test = data_mnist(datadir=data_dir,
                                                  train_start=train_start,
                                                  train_end=train_end,
                                                  test_start=test_start,
                                                  test_end=test_end)

    # Initialize substitute training set reserved for adversary
    X_sub = X_test[:holdout]
    Y_sub = np.argmax(Y_test[:holdout], axis=1)
    print("Y_sub shape=")
    print(Y_sub.shape)

    # Redefine test set as remaining samples unavailable to adversaries
    X_test = X_test[holdout:]
    Y_test = Y_test[holdout:]

    # MNIST-specific dimensions
    img_rows = 28
    img_cols = 28
    channels = 1
    nb_classes = 10

    # Define input and output TF placeholders
    x = tf.placeholder(tf.float32, shape=(None, 28, 28, 1))
    y = tf.placeholder(tf.float32, shape=(None, 10))
    phase = tf.placeholder(tf.bool, name='phase')

    logits_scalar = tf.placeholder_with_default(
        INIT_T, shape=(), name="logits_temperature")

    # Seed random number generator so tutorial is reproducible
    rng = np.random.RandomState([2017, 8, 30])

    # Simulate the black-box model locally
    # You could replace this by a remote labeling API for instance
    print("Preparing the black-box model.")
    prep_bbox_out = prep_bbox(sess, logits_scalar, x, y, X_train, Y_train, X_test, Y_test,
                              img_rows, img_cols, channels, nb_epochs, batch_size, learning_rate,
                              rng=rng, phase=phase, binary=binary, scale=scale,
                              nb_filters=nb_filters, model_path=model_path,
                              adv=adv, delay=delay, eps=eps)

    model, bbox_preds, accuracies['bbox'], model_path = prep_bbox_out

    # Train substitute using method from https://arxiv.org/abs/1602.02697
    print("Training the substitute model.")
    train_sub_out = train_sub(sess, logits_scalar, x, y, bbox_preds, X_sub, Y_sub,
                              nb_classes, nb_epochs_s, batch_size,
                              learning_rate, data_aug, lmbda, rng=rng,
                              phase=phase, model_path=model_path)
    model_sub, preds_sub = train_sub_out

    # Evaluate the substitute model on clean test examples
    eval_params = {'batch_size': batch_size}
    acc = model_eval(sess, x, y, preds_sub, X_test,
                     Y_test, phase=phase, args=eval_params)
    accuracies['sub'] = acc

    # Initialize the Fast Gradient Sign Method (FGSM) attack object.
    fgsm_par = {'eps': eps, 'ord': np.inf, 'clip_min': 0., 'clip_max': 1.}
    fgsm = FastGradientMethod(model_sub, sess=sess)

    if targeted:
        from cleverhans.utils import build_targeted_dataset
        adv_inputs, true_labels, adv_ys = build_targeted_dataset(
            X_test, Y_test, np.arange(nb_samples), nb_classes, img_rows, img_cols, channels)
        att_batch_size = np.clip(
            nb_samples * (nb_classes - 1), a_max=MAX_BATCH_SIZE, a_min=1)
        nb_adv_per_sample = nb_classes - 1
        yname = "y_target"

    else:
        att_batch_size = np.minimum(nb_samples, MAX_BATCH_SIZE)
        nb_adv_per_sample = 1
        adv_ys = None
        yname = "y"

    # Craft adversarial examples using the substitute
    eval_params = {'batch_size': att_batch_size}

    if targeted:
        fgsm_par.update({yname: adv_ys})
        x_adv_sub = fgsm.generate_np(adv_inputs, phase, **fgsm_par)
        accuracy = model_eval(sess, x, y, model(x, reuse=True), x_adv_sub, true_labels,
                              phase=phase, args=eval_params)
    else:
        x_adv_sub = fgsm.generate(x, phase, **fgsm_par)

        # Evaluate the accuracy of the "black-box" model on adversarial
        # examples
        accuracy = model_eval(sess, x, y, model(x_adv_sub), X_test, Y_test,
                              phase=phase, args=eval_params)
    print('Test accuracy of oracle on adversarial examples generated '
          'using the substitute: ' + str(accuracy))
    accuracies['bbox_on_sub_adv_ex'] = accuracy

    return accuracies

def main(argv=None):
    mnist_blackbox(nb_classes=FLAGS.nb_classes, batch_size=FLAGS.batch_size,
                   nb_samples=FLAGS.nb_samples, nb_filters=FLAGS.nb_filters,
                   eps=FLAGS.eps, learning_rate=FLAGS.learning_rate,
                   nb_epochs=FLAGS.nb_epochs, holdout=FLAGS.holdout,
                   data_aug=FLAGS.data_aug, nb_epochs_s=FLAGS.nb_epochs_s,
                   lmbda=FLAGS.lmbda, binary=FLAGS.binary, scale=FLAGS.scale, model_path=FLAGS.model_path,
                   targeted=FLAGS.targeted, data_dir=FLAGS.data_dir,
                   adv=FLAGS.adv, delay=FLAGS.delay)

if __name__ == '__main__':
    parser = argparse.ArgumentParser()

    # General flags
    parser.add_argument('--gpu', help='id of GPU to use')
    parser.add_argument("--model_path", help='Path to save or load model')
    parser.add_argument(
        '--data_dir', help='Path to training data', default='/scratch/gallowaa/mnist')
    parser.add_argument("--nb_classes", type=int, default=10)
    parser.add_argument("--batch_size", type=int, default=128)
    parser.add_argument("--learning_rate", type=float, default=0.001)

    # For oracle
    parser.add_argument("--nb_epochs", type=int, default=15)
    parser.add_argument(
        '--binary', help='Use a binary model?', action="store_true")
    parser.add_argument('--scale', help='Scale activations of the binary model?',
                        action="store_true")
    parser.add_argument('--nb_filters', type=int, default=64,
                        help='Number of filters in first layer')

    # For substitute
    parser.add_argument("--holdout", type=int, default=150)
    parser.add_argument("--data_aug", type=int, default=6)
    parser.add_argument("--nb_epochs_s", type=int, default=10)
    parser.add_argument("--lmbda", type=float, default=0.1)

    # Attack
    parser.add_argument("--eps", type=float, default=0.3)
    parser.add_argument(
        '--targeted', help='Run a targeted attack?', action="store_true")
    parser.add_argument('--nb_samples', type=int,
                        default=10, help='Nb of inputs to attack')

    # Adversarial training flags
    parser.add_argument(
        '--adv', help='Adversarial training type?', type=int, default=0)
    parser.add_argument('--delay', type=int,
                        default=10, help='Nb of epochs to delay adv training by')

    FLAGS, _ = parser.parse_known_args()

    if FLAGS.gpu:
        os.environ['CUDA_VISIBLE_DEVICES'] = FLAGS.gpu

    tf.app.run()

max_stars_repo_name: AngusG/cleverhans-attacking-bnns

max_issues_repo_name: AngusG/cleverhans-attacking-bnns

max_forks_repo_name: AngusG/cleverhans-attacking-bnns


text:
# -*- coding: utf-8 -*-
"""
Created on Sat Mar 17 19:37:14 2018

@author: Marina Mautner
"""
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import seaborn as sns
import os
from pathlib import Path

# Set current working directory to script folder
os.chdir(os.path.dirname(__file__))

# Set cwd to repository folder
path = Path(os.getcwd()).parent.parent
os.chdir(path)

df = pd.read_csv(Path.cwd()/'data_processed'/'wells'/'POPULATION.csv')
df = df.set_index(pd.to_datetime(df.Year.apply(lambda x: str(x)+'/01/01')))

t1 = np.array([1950,1960,1970,1980,1990,1995,2000,2005,2010])
p1 = df.PopZMVM_Census.dropna()['1950/01/01':]
t2 = np.array([2010,2011,2012,2013,2014,2015,2016,2017,2018,2019,2020,2021,2022,
               2023,2024,2025,2026,2027,2028,2029,2030])
p2 = df.PopZMVM_Proj.dropna()

# Fit population curve
# Census
g1 = np.polyfit(t1, p1, 2)
f1 = np.poly1d(g1)

# CONAPO
g2 = np.polyfit(t2, p2, 2)
f2 = np.poly1d(g2)

x = np.arange(1950,2030,1)
Pt1 = np.zeros(len(x))
Pt2 = np.zeros(len(x))
for i in range(0,len(x)):
    Pt1[i] = f1(x[i])
    Pt2[i] = f2(x[i])

plt.plot(np.array(df.Year),np.array(df.PopZMVM_Proj)/1000000,'.',color='black')
plt.plot(t1,p1/1000000,'s',color='grey',zorder=3)
plt.plot(x,Pt1/1000000,'--',color='red')
plt.plot(x,Pt2/1000000,'--',color='blue')
plt.grid(True)
plt.xlabel('Year')
plt.ylabel('Population (millions)')
plt.legend(['CONAPO 2014 Projections','1895 to 2010 Census Data', 'Census Model','CONAPO Projection Model'])

#%%
PM = pd.DataFrame(np.transpose(np.array([x,Pt1])),columns=['Year','P'])
PM = PM.set_index(pd.to_datetime(PM.Year.apply(lambda x: str(int(x))+'/01/01')))
PM.to_csv('C:/Users/MM/Google Drive/ValleMexico/Data/Pumping/POP_CM.csv')

Pop = PM.P['2005-01-01':'2016-01-01']
Pop['2011-01-01'] = np.nan
Pop = Pop.dropna().values
Q = df.REPDA_m3d.dropna().values

WEL = pd.DataFrame(np.transpose(np.array([Q,Pop])/10**6),index=df.REPDA_m3d.dropna().index,columns=['Q','P'])
#
#Population
g = np.polyfit(Pop, Q, 1)
f = np.poly1d(g)

##Time
#year = np.arange(2005,2017)
#year = np.delete(year,6)
#g = np.polyfit(year, Q, 1)
#f = np.poly1d(g)

residuals = f(Pop) - Q
#residuals = f(year) - Q
ss_res = np.sum(residuals**2)
ss_tot = np.sum((Q-np.mean(Q))**2)
r_squared = 1 - (ss_res / ss_tot)
rmse = np.sqrt(ss_res/len(Pop))
#rmse = np.sqrt(ss_res/len(year))

plt.grid(True)
plt.xlim(1.4e7,2.2e7)
#plt.xlim(1980,2020)
sns.regplot(x=Pop, y='Q', data=WEL,marker='.',order=1)
#sns.regplot(x=year, y='Q', data=WEL,marker='.',order=1)
plt.xlabel('Population (millions)')
#plt.xlabel('Year')
plt.ylabel('Groundwater Concessions (' r'$hm^3/d$' ')')
print(rmse)
print(r_squared)
print(f)

max_stars_repo_name: mrlmautner/UrbanGW

max_issues_repo_name: mrlmautner/UrbanGW

max_forks_repo_name: mrlmautner/UrbanGW


text:
import logging
import numpy as np
from sklearn.neighbors import KDTree
from sklearn.neighbors import KNeighborsClassifier
from typing import Tuple, Any

logger = logging.getLogger(__name__)

class TrustScore(object):

    def __init__(self, k_filter: int = 10, alpha: float = 0., filter_type: str = None,
                 leaf_size: int = 40, metric: str = 'euclidean', dist_filter_type: str = 'point') -> None:
        """
        Initialize trust scores.

        Parameters
        ----------
        k_filter
            Number of neighbors used during either kNN distance or probability filtering.
        alpha
            Fraction of instances to filter out to reduce impact of outliers.
        filter_type
            Filter method; either 'distance_knn' or 'probability_knn'
        leaf_size
            Number of points at which to switch to brute-force. Affects speed and memory required to build trees.
            Memory to store the tree scales with n_samples / leaf_size.
        metric
            Distance metric used for the tree. See sklearn's DistanceMetric class for a list of available metrics.
        dist_filter_type
            Use either the distance to the k-nearest point (dist_filter_type = 'point') or
            the average distance from the first to the k-nearest point in the data (dist_filter_type = 'mean').
        """
        self.k_filter = k_filter
        self.alpha = alpha
        self.filter = filter_type
        self.eps = 1e-12
        self.leaf_size = leaf_size
        self.metric = metric
        self.dist_filter_type = dist_filter_type

    def filter_by_distance_knn(self, X: np.ndarray) -> np.ndarray:
        """
        Filter out instances with low kNN density. Calculate distance to k-nearest point in the data for each
        instance and remove instances above a cutoff distance.

        Parameters
        ----------
        X
            Data

        Returns
        -------
        Filtered data.
        """
        kdtree = KDTree(X, leaf_size=self.leaf_size, metric=self.metric)
        knn_r = kdtree.query(X, k=self.k_filter + 1)[0]  # distances from 0 to k-nearest points
        if self.dist_filter_type == 'point':
            knn_r = knn_r[:, -1]
        elif self.dist_filter_type == 'mean':
            knn_r = np.mean(knn_r[:, 1:], axis=1)  # exclude distance of instance to itself
        cutoff_r = np.percentile(knn_r, (1 - self.alpha) * 100)  # cutoff distance
        X_keep = X[np.where(knn_r <= cutoff_r)[0], :]  # define instances to keep
        return X_keep

    def filter_by_probability_knn(self, X: np.ndarray, Y: np.ndarray) -> Tuple[np.ndarray, np.ndarray]:
        """
        Filter out instances with high label disagreement amongst its k nearest neighbors.

        Parameters
        ----------
        X
            Data
        Y
            Predicted class labels

        Returns
        -------
        Filtered data and labels.
        """
        if self.k_filter == 1:
            logger.warning('Number of nearest neighbors used for probability density filtering should '
                           'be >1, otherwise the prediction probabilities are either 0 or 1 making '
                           'probability filtering useless.')
        # fit kNN classifier and make predictions on X
        clf = KNeighborsClassifier(n_neighbors=self.k_filter, leaf_size=self.leaf_size, metric=self.metric)
        clf.fit(X, Y)
        preds_proba = clf.predict_proba(X)
        # define cutoff and instances to keep
        preds_max = np.max(preds_proba, axis=1)
        cutoff_proba = np.percentile(preds_max, self.alpha * 100)  # cutoff probability
        keep_id = np.where(preds_max >= cutoff_proba)[0]  # define id's of instances to keep
        X_keep, Y_keep = X[keep_id, :], Y[keep_id]
        return X_keep, Y_keep

    def fit(self, X: np.ndarray, Y: np.ndarray, classes: int = None) -> None:
        """
        Build KDTrees for each prediction class.

        Parameters
        ----------
        X
            Data
        Y
            Target labels, either one-hot encoded or the actual class label.
        classes
            Number of prediction classes, needs to be provided if Y equals the predicted class.
        """
        self.classes = classes if classes is not None else Y.shape[1]
        self.kdtrees = [None] * self.classes  # type: Any

        # KDTree and kNeighborsClassifier need 2D data
        if len(X.shape) > 2:
            logger.warning('Reshaping data from {0} to {1} so k-d trees can '
                           'be built.'.format(X.shape, X.reshape(X.shape[0], -1).shape))
            X = X.reshape(X.shape[0], -1)

        # make sure Y represents predicted classes, not one-hot encodings
        if len(Y.shape) > 1:
            Y = np.argmax(Y, axis=1)

        if self.filter == 'probability_knn':
            X_filter, Y_filter = self.filter_by_probability_knn(X, Y)

        for c in range(self.classes):

            if self.filter is None:
                X_fit = X[np.where(Y == c)[0]]
            elif self.filter == 'distance_knn':
                X_fit = self.filter_by_distance_knn(X[np.where(Y == c)[0]])
            elif self.filter == 'probability_knn':
                X_fit = X_filter[np.where(Y_filter == c)[0]]

            no_x_fit = len(X_fit) == 0
            if no_x_fit and len(X[np.where(Y == c)[0]]) == 0:
                logger.warning('No instances available for class %s', c)
            elif no_x_fit:
                logger.warning('Filtered all the instances for class %s. Lower alpha or check data.', c)

            self.kdtrees[c] = KDTree(X_fit, leaf_size=self.leaf_size, metric=self.metric)  # build KDTree for class c

    def score(self, X: np.ndarray, Y: np.ndarray, k: int = 2, dist_type: str = 'point') -> np.ndarray:
        """
        Calculate trust scores = ratio of distance to closest class other than the
        predicted class to distance to predicted class.

        Parameters
        ----------
        X
            Instances to calculate trust score for.
        Y
            Either prediction probabilities for each class or the predicted class.
        k
            Number of nearest neighbors used for distance calculation.
        dist_type
            Use either the distance to the k-nearest point (dist_type = 'point') or
            the average distance from the first to the k-nearest point in the data (dist_type = 'mean').

        Returns
        -------
        Batch with trust scores.
        """
        # make sure Y represents predicted classes, not probabilities
        if len(Y.shape) > 1:
            Y = np.argmax(Y, axis=1)

        # KDTree needs 2D data
        if len(X.shape) > 2:
            logger.warning('Reshaping data from {0} to {1} so k-d trees can '
                           'be queried.'.format(X.shape, X.reshape(X.shape[0], -1).shape))
            X = X.reshape(X.shape[0], -1)

        d = np.tile(None, (X.shape[0], self.classes))  # init distance matrix: [nb instances, nb classes]

        for c in range(self.classes):
            d_tmp = self.kdtrees[c].query(X, k=k)[0]  # get k nearest neighbors for each class
            if dist_type == 'point':
                d[:, c] = d_tmp[:, -1]
            elif dist_type == 'mean':
                d[:, c] = np.mean(d_tmp, axis=1)

        sorted_d = np.sort(d, axis=1)  # sort distance each instance in batch over classes
        # get distance to predicted and closest other class and calculate trust score
        d_to_pred = d[range(d.shape[0]), Y]
        d_to_closest_not_pred = np.where(sorted_d[:, 0] != d_to_pred, sorted_d[:, 0], sorted_d[:, 1])
        trust_score = d_to_closest_not_pred / (d_to_pred + self.eps)
        return trust_score

max_stars_repo_name: NunoEdgarGFlowHub/alibi

max_issues_repo_name: NunoEdgarGFlowHub/alibi

max_forks_repo_name: NunoEdgarGFlowHub/alibi


text:
#!/usr/bin/env python

__author__ = 'jessada'

import sys
import argparse
import tempfile
from os.path import join as join_path
import subprocess

# for creating temp file
from random import choice
from string import ascii_letters

# for parsing into ExAC format
import numpy as np
import numpy.lib.recfunctions as rfn

np.set_printoptions(threshold=np.inf)

def count_depth(bam_list, region, out_file):
    cmd = "samtools depth"
    cmd += " -f " + bam_list
    cmd += " -r " + region
    cmd += " > " + out_file

    print >> sys.stderr, "##INFO executing: " + cmd

    # error handling part might be incomplete
    p = subprocess.Popen(cmd,
                         shell=True,
                         stdout=subprocess.PIPE,
                         stderr=subprocess.PIPE,
                         )
    stdout_data, stderr_data = p.communicate()
    print stdout_data
    print >> sys.stderr, stderr_data

def samtools_depth_to_panel_coverage(depth_file, out_file):
    chrom = np.genfromtxt(depth_file,
                          delimiter="\t",
                          dtype="S",
                          )[:,0:1]
    pos = np.genfromtxt(depth_file,
                        delimiter="\t",
                        dtype="S",
                        )[:,1:2]
    depths = np.genfromtxt(depth_file,
                           delimiter="\t",
                           )[:,2:]
    n_samples = float(depths.shape[1])
    print >> sys.stderr, "##INFO number of samples: " + str(int(n_samples))
    merged_arrays = [chrom,
                     pos,
                     np.mean(depths, axis=1),
                     np.median(depths, axis=1),
                     (depths > 1).sum(axis=1)/n_samples,
                     (depths > 5).sum(axis=1)/n_samples,
                     (depths > 10).sum(axis=1)/n_samples,
                     (depths > 15).sum(axis=1)/n_samples,
                     (depths > 20).sum(axis=1)/n_samples,
                     (depths > 25).sum(axis=1)/n_samples,
                     (depths > 30).sum(axis=1)/n_samples,
                     (depths > 50).sum(axis=1)/n_samples,
                     (depths > 100).sum(axis=1)/n_samples,
                     ]
    out_exac = rfn.merge_arrays(merged_arrays,
                                flatten = True,
                                usemask = False,
                                )

    # write calculated coverage to file
    header="#chrom"
    header+="\tpos"
    header+="\tmean"
    header+="\tmedian"
    header+="\t1"
    header+="\t5"
    header+="\t10"
    header+="\t15"
    header+="\t20"
    header+="\t25"
    header+="\t30"
    header+="\t50"
    header+="\t100"
    fmt="%s\t%s\t%.2f\t%.2f\t%.4f\t%.4f\t%.4f\t%.4f\t%.4f\t%.4f\t%.4f\t%.4f\t%.4f"
    np.savetxt(out_file,
               out_exac,
               fmt=fmt,
               delimiter="\t",
               header=header,
               comments="",
               )

def main(args):
    # tmp_depth_file is a temp file to keep depth count
    # please feel free to modify its name and location as you like
    tmp_depth_file = join_path(tempfile.mkdtemp(),
                               "tmp_" + ''.join([choice(ascii_letters) for i in range(6)]))
    count_depth(args.bam_list, args.region, tmp_depth_file)

    panel_coverage_out = "Panel."
    panel_coverage_out += args.region.replace(":", "_").replace("-", "_")
    panel_coverage_out += ".coverage.txt"
    samtools_depth_to_panel_coverage(tmp_depth_file,
                                     join_path(args.out_dir, panel_coverage_out)
                                     )

if __name__ == '__main__':
    parser = argparse.ArgumentParser()

    parser.add_argument('--region', '-r', metavar="CHR:FROM-TO", help='Only report depth in specified region.', required=True)
    parser.add_argument('--bam_list', '-b', metavar="FILE", help='Compute depth at list of positions or regions in specified BAM FILE.', required=True)
    parser.add_argument('--out_dir', '-o', metavar="DIR", help='Output directory.', required=True)
    args = parser.parse_args()
    main(args)