text:
from __future__ import print_function
from __future__ import division
from __future__ import unicode_literals
from __future__ import absolute_import
from builtins import *  # NOQA
from future import standard_library
standard_library.install_aliases()  # NOQA
import argparse
import os

import chainer
import numpy as np

import chainerrl
from chokozainerrl import experiments
from chokozainerrl import tools
from chainerrl import explorers
from chainerrl import misc
from chainerrl import replay_buffer

from chainerrl.wrappers import atari_wrappers

def make_args(argstr):
    parser = argparse.ArgumentParser()
    parser.add_argument('--mode', type=str, default='check')
    parser.add_argument('--env', type=str, default='BreakoutNoFrameskip-v4')
    parser.add_argument('--outdir', type=str, default='results')
    parser.add_argument('--gpu', type=int, default=0)
    parser.add_argument('--load-agent', type=str, default=None)
    parser.add_argument('--log-type',type=str,default="full_stream")
    parser.add_argument('--save-mp4',type=str,default="test.mp4")

    parser.add_argument('--steps', type=int, default=5 * 10 ** 7)
    parser.add_argument('--step-offset', type=int, default=0)
    parser.add_argument('--checkpoint-frequency', type=int,default=None)
    parser.add_argument('--max-frames', type=int,default=30 * 60 * 60)  # 30 minutes with 60 fps
    parser.add_argument('--eval-interval', type=int, default=250000)
    parser.add_argument('--eval-n-runs', type=int, default=10)

    parser.add_argument('--seed', type=int, default=0)
    parser.add_argument('--use-sdl', action='store_true', default=False)
    parser.add_argument('--final-exploration-frames',type=int, default=10 ** 6)
    parser.add_argument('--final-epsilon', type=float, default=0.01)
    parser.add_argument('--eval-epsilon', type=float, default=0.001)
    parser.add_argument('--replay-start-size', type=int, default=5 * 10 ** 4)
    parser.add_argument('--target-update-interval',type=int, default=10 ** 4)
    parser.add_argument('--update-interval', type=int, default=4)
    parser.add_argument('--batch-size', type=int, default=32)
    parser.add_argument('--lr', type=float, default=2.5e-4)
    parser.add_argument('--render', action='store_true', default=False)
    parser.add_argument('--monitor', action='store_true', default=False)
    myargs = parser.parse_args(argstr)
    return myargs

def main(args):
    import logging
    logging.basicConfig(level=logging.INFO, filename='log')

    if(type(args) is list):
        args=make_args(args)
    if not os.path.exists(args.outdir):
        os.makedirs(args.outdir)

    # Set a random seed used in ChainerRL.
    misc.set_random_seed(args.seed, gpus=(args.gpu,))

    # Set different random seeds for train and test envs.
    train_seed = args.seed
    test_seed = 2 ** 31 - 1 - args.seed

    def make_env(test):
        # Use different random seeds for train and test envs
        env_seed = test_seed if test else train_seed
        env = atari_wrappers.wrap_deepmind(
            atari_wrappers.make_atari(args.env, max_frames=args.max_frames),
            episode_life=not test,
            clip_rewards=not test)
        env.seed(int(env_seed))
        if test:
            # Randomize actions like epsilon-greedy in evaluation as well
            env = chainerrl.wrappers.RandomizeAction(env, args.eval_epsilon)
        if args.monitor:
            env = chainerrl.wrappers.Monitor(
                env, args.outdir,
                mode='evaluation' if test else 'training')
        if args.render:
            env = chainerrl.wrappers.Render(env)
        return env

    env = make_env(test=False)
    eval_env = make_env(test=True)

    n_actions = env.action_space.n

    n_atoms = 51
    v_max = 10
    v_min = -10
    q_func = chainerrl.links.Sequence(
        chainerrl.links.NatureDQNHead(),
        chainerrl.q_functions.DistributionalFCStateQFunctionWithDiscreteAction(
            None, n_actions, n_atoms, v_min, v_max,
            n_hidden_channels=0, n_hidden_layers=0),
    )

    # Draw the computational graph and save it in the output directory.
    chainerrl.misc.draw_computational_graph(
        [q_func(np.zeros((4, 84, 84), dtype=np.float32)[None])],
        os.path.join(args.outdir, 'model'))

    # Use the same hyper parameters as https://arxiv.org/abs/1707.06887
    opt = chainer.optimizers.Adam(2.5e-4, eps=1e-2 / args.batch_size)
    opt.setup(q_func)

    rbuf = replay_buffer.ReplayBuffer(10 ** 6)

    explorer = explorers.LinearDecayEpsilonGreedy(
        1.0, args.final_epsilon,
        args.final_exploration_frames,
        lambda: np.random.randint(n_actions))

    def phi(x):
        # Feature extractor
        return np.asarray(x, dtype=np.float32) / 255

    agent = chainerrl.agents.CategoricalDQN(
        q_func, opt, rbuf, gpu=args.gpu, gamma=0.99,
        explorer=explorer, replay_start_size=args.replay_start_size,
        target_update_interval=args.target_update_interval,
        update_interval=args.update_interval,
        batch_accumulator='mean',
        phi=phi,
    )

    if args.load_agent:
        agent.load(args.load_agent)

    if (args.mode=='train'):
        experiments.train_agent_with_evaluation(
            agent=agent, env=env, steps=args.steps,
            checkpoint_freq=args.checkpoint_frequency,
            step_offset=args.step_offset,
            eval_n_steps=None,
            eval_n_episodes=args.eval_n_runs, eval_interval=args.eval_interval,
            outdir=args.outdir,
            save_best_so_far_agent=False,
            eval_env=eval_env,
            log_type=args.log_type
        )
    elif (args.mode=='check'):
        return tools.make_video.check(env=env,agent=agent,save_mp4=args.save_mp4)

    elif (args.mode=='growth'):
        return tools.make_video.growth(env=env,agent=agent,outdir=args.outdir,max_num=args.max_frames,save_mp4=args.save_mp4)

max_stars_repo_name: chokozainer/chokozainerrl

max_issues_repo_name: chokozainer/chokozainerrl

max_forks_repo_name: chokozainer/chokozainerrl


text:
#the purpose of this file/package is to be able to debug the camera feed to make sure it is outputting results as expected
import cv2
import numpy as np
import pudb
import rclpy
from cv_bridge import CvBridge, CvBridgeError
from rclpy.node import Node
from sensor_msgs.msg import Image
from std_msgs.msg import Header

class DebugCamera(Node):
    def __init__(self):
        super().__init__('DebugCamera')
        self.subscription = self.create_subscription(Image, 'camera', self.image_callback, 100)
        self.bridge = CvBridge()

    def image_callback(self, msg):
        self.get_logger().info('I heard {0}'.format(str(msg.header)))
        self.get_logger().info('This is an upscaled image, the original is 64x64')
        image = self.bridge.imgmsg_to_cv2(msg, "bgr8")
        #The above converts the image to RGB, but we want it to stay BGR
        image = cv2.cvtColor(image, cv2.COLOR_RGB2BGR)
        image = cv2.resize(image, dsize=(500,500), interpolation=cv2.INTER_CUBIC)
        image = np.fliplr(image)
        #The arrow is strictly for visualization, not used for navigation
        cv2.imshow('frame', image)
        key = cv2.waitKey(1)
        if key == ord('q'):
            self.__del__()

def main(args=None):
    rclpy.init(args=args)

    waypoint_publisher = DebugCamera()

    rclpy.spin(waypoint_publisher)

    # Destroy the node explicitly
    # (optional - otherwise it will be done automatically
    # when the garbage collector destroys the node object)
    waypoint_publisher.destroy_node()
    rclpy.shutdown()

ext: py

reason: import numpy


text:
import numpy as np
#import pylab as pl
import matplotlib.pyplot as pl

pl.style.use('seaborn-talk')
pl.rc('axes', linewidth=2)
ax = pl.gca()

ax.get_yaxis().set_tick_params(which='both', direction='in', width=2, right='on')
ax.get_xaxis().set_tick_params(which='both', direction='in', width=2, top='on')
ax.get_yaxis().set_tick_params(which='major', size=10)
ax.get_xaxis().set_tick_params(which='major', size=10)
ax.get_yaxis().set_tick_params(which='minor', size=5)
ax.get_xaxis().set_tick_params(which='minor', size=5)

data = np.loadtxt('LeptonicModulation.dat')

pl.title(r"Leptonic modulation at $Q = M_Z$ and $y = 0$", fontsize=20)

pl.xlabel(r'$p_{T,\ell}$ [GeV]', fontsize=22, labelpad=-2)
pl.ylabel(r'$\frac{d\sigma}{dQ dy dq_T dp_T} / \frac{d\sigma}{dQ dy dq_T}\;\times 100$', fontsize=22, labelpad=-2)
pl.xlim(41.5, 45.5)
pl.ylim(0, 20)
pl.xticks(size=20)
pl.yticks(size=20)

np = 10000

pl.plot(data[0*np:1*np,1], 100*data[0*np:1*np,2], 'r-', label=r'$q_T = 1$ GeV')
pl.plot(data[1*np:2*np,1], 100*data[1*np:2*np,2], 'b-', label=r'$q_T = 2$ GeV')
pl.plot(data[2*np:3*np,1], 100*data[2*np:3*np,2], 'm-', label=r'$q_T = 3$ GeV')
pl.plot(data[3*np:4*np,1], 100*data[3*np:4*np,2], 'g-', label=r'$q_T = 4$ GeV')
pl.plot(data[4*np:5*np,1], 100*data[4*np:5*np,2], 'c-', label=r'$q_T = 5$ GeV')

pl.legend(fontsize=22)

pl.savefig('LeptonicModulation.pdf')
pl.show()


text:
"""RefineNet-LightWeight

RefineNet-LigthWeight PyTorch for non-commercial purposes

Copyright (c) 2018, Vladimir Nekrasov (vladimir.nekrasov@adelaide.edu.au)
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

* Redistributions of source code must retain the above copyright notice, this
  list of conditions and the following disclaimer.

* Redistributions in binary form must reproduce the above copyright notice,
  this list of conditions and the following disclaimer in the documentation
  and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
"""

import torch.nn as nn
import torch.nn.functional as F
import torch

import numpy as np

from utils.helpers import maybe_download
from utils.layer_factory import conv1x1, conv3x3, CRPBlock

data_info = {
    7 : 'Person',
    21: 'VOC',
    40: 'NYU',
    60: 'Context'
    }

models_urls = {
    '50_person'  : 'https://cloudstor.aarnet.edu.au/plus/s/mLA7NxVSPjNL7Oo/download',
    '101_person' : 'https://cloudstor.aarnet.edu.au/plus/s/f1tGGpwdCnYS3xu/download',
    '152_person' : 'https://cloudstor.aarnet.edu.au/plus/s/Ql64rWqiTvWGAA0/download',

    '50_voc'     : 'https://cloudstor.aarnet.edu.au/plus/s/2E1KrdF2Rfc5khB/download',
    '101_voc'    : 'https://cloudstor.aarnet.edu.au/plus/s/CPRKWiaCIDRdOwF/download',
    '152_voc'    : 'https://cloudstor.aarnet.edu.au/plus/s/2w8bFOd45JtPqbD/download',

    '50_nyu'     : 'https://cloudstor.aarnet.edu.au/plus/s/gE8dnQmHr9svpfu/download',
    '101_nyu'    : 'https://cloudstor.aarnet.edu.au/plus/s/VnsaSUHNZkuIqeB/download',
    '152_nyu'    : 'https://cloudstor.aarnet.edu.au/plus/s/EkPQzB2KtrrDnKf/download',

    '101_context': 'https://cloudstor.aarnet.edu.au/plus/s/hqmplxWOBbOYYjN/download',
    '152_context': 'https://cloudstor.aarnet.edu.au/plus/s/O84NszlYlsu00fW/download',

    '50_imagenet' : 'https://download.pytorch.org/models/resnet50-19c8e357.pth',
    '101_imagenet': 'https://download.pytorch.org/models/resnet101-5d3b4d8f.pth',
    '152_imagenet': 'https://download.pytorch.org/models/resnet152-b121ed2d.pth',
    }

stages_suffixes = {0 : '_conv',
                   1 : '_conv_relu_varout_dimred'}

class BasicBlock(nn.Module):
    expansion = 1

    def __init__(self, inplanes, planes, stride=1, downsample=None):
        super(BasicBlock, self).__init__()
        self.conv1 = conv3x3(inplanes, planes, stride)
        self.bn1 = nn.BatchNorm2d(planes)
        self.relu = nn.ReLU(inplace=True)
        self.conv2 = conv3x3(planes, planes)
        self.bn2 = nn.BatchNorm2d(planes)
        self.downsample = downsample
        self.stride = stride

    def forward(self, x):
        residual = x

        out = self.conv1(x)
        out = self.bn1(out)
        out = self.relu(out)

        out = self.conv2(out)
        out = self.bn2(out)

        if self.downsample is not None:
            residual = self.downsample(x)

        out += residual
        out = self.relu(out)

        return out

class Bottleneck(nn.Module):
    expansion = 4

    def __init__(self, inplanes, planes, stride=1, downsample=None):
        super(Bottleneck, self).__init__()
        self.conv1 = nn.Conv2d(inplanes, planes, kernel_size=1, bias=False)
        self.bn1 = nn.BatchNorm2d(planes)
        self.conv2 = nn.Conv2d(planes, planes, kernel_size=3, stride=stride,
                               padding=1, bias=False)
        self.bn2 = nn.BatchNorm2d(planes)
        self.conv3 = nn.Conv2d(planes, planes * 4, kernel_size=1, bias=False)
        self.bn3 = nn.BatchNorm2d(planes * 4)
        self.relu = nn.ReLU(inplace=True)
        self.downsample = downsample
        self.stride = stride

    def forward(self, x):
        residual = x

        out = self.conv1(x)
        out = self.bn1(out)
        out = self.relu(out)

        out = self.conv2(out)
        out = self.bn2(out)
        out = self.relu(out)

        out = self.conv3(out)
        out = self.bn3(out)

        if self.downsample is not None:
            residual = self.downsample(x)

        out += residual
        out = self.relu(out)

        return out

### This class was made to abstract the RefineNet code from the original
### ResNet (originally ResNetLW) module into its own module. Further, I have
### made the following modifications to convert from LW-RefineNet to RefineNet:
###     - conv1x1 changed to conv3x3 in all lines ending in _varout_dimred,
###       which are part of the multi-resolution fusion blocks
###     - RCU blocks added
###     - heavy reorganization, renaming, etc
### I would like to rename the blocks in this module but I think the names are important for
###   other parts of the code
class RefineNet(nn.Module):

    # inplanes_hr/lr are the number of channels of the high resolution input and
    #   low resolution input respectively. if inplanes_hr is omitted, no fusion takes
    #   place: this is useful for layer 4 in our case, to prep for fusing it with layer 3.
    # the inputs will be passed through:
    #   - if HR is provided: 2 RCU blocks for HR only
    #     otherwise, 2 RCU blocks for LR only
    #   - fusion of HR and LR (if HR is provided)
    #   - CRP block + 1 RCU block for result of fusion (or LR if no HR provided)
    def __init__(self, inplanes_lr, inplanes_hr=None, fancy_upsample=False):
        super(RefineNet, self).__init__()

        self.fancy_upsample = fancy_upsample

        # first set of RCU blocks
        if inplanes_hr == None:
            self.inplanes = inplanes_lr
            self.rcu1_lr = self._make_layer(BasicBlock, planes=inplanes_lr, blocks=2)
        else:
            self.inplanes = inplanes_hr
            self.rcu1_hr = self._make_layer(BasicBlock, planes=inplanes_hr, blocks=2)

        # fusion
        if inplanes_hr != None:
            self.adapt_stage2_b2_joint_varout_dimred = conv3x3(inplanes_hr, inplanes_hr, bias=False)
            self.mflow_conv_g1_b3_joint_varout_dimred = conv3x3(inplanes_lr, inplanes_hr, bias=False)
            if fancy_upsample:
                # learnable upsample with a single (strided) transpose convolution (for now)
                self.upsample = nn.ConvTranspose2d(inplanes_hr, inplanes_hr, kernel_size=3, stride=2, padding=1)
                #self.upsample = nn.ConvTranspose2d(inplanes_hr, inplanes_hr, kernel_size=5, stride=2, padding=2)

        # CRP and RCU for fusion result
        outplanes = inplanes_hr if inplanes_hr != None else inplanes_lr
        self.mflow_conv_g1_pool = self._make_crp(outplanes, outplanes, 4)
        self.inplanes = outplanes
        self.rcu2 = self._make_layer(BasicBlock, planes=outplanes, blocks=1)

    # _make_layer and _make_crp copied from the ResNetLW class

    def _make_crp(self, in_planes, out_planes, stages):
        layers = [CRPBlock(in_planes, out_planes,stages)]
        return nn.Sequential(*layers)

    def _make_layer(self, block, planes, blocks, stride=1):
        downsample = None
        if stride != 1 or self.inplanes != planes * block.expansion:
            downsample = nn.Sequential(
                nn.Conv2d(self.inplanes, planes * block.expansion,
                          kernel_size=1, stride=stride, bias=False),
                nn.BatchNorm2d(planes * block.expansion),
            )

        layers = []
        layers.append(block(self.inplanes, planes, stride, downsample))
        self.inplanes = planes * block.expansion
        for i in range(1, blocks):
            layers.append(block(self.inplanes, planes))

        return nn.Sequential(*layers)

    def forward(self, x_lr, x_hr=None):

        # first set of RCU blocks
        if x_hr is None:
            x_lr = self.rcu1_lr(x_lr)
        else:
            x_hr = self.rcu1_hr(x_hr)

        # fusion (or not)
        x = x_lr
        if x_hr is not None:
            x_lr = self.mflow_conv_g1_b3_joint_varout_dimred(x_lr)
            if self.fancy_upsample:
                x_lr = self.upsample(x_lr, output_size=x_hr.size())  # might have to do [2:]
            else:
                x_lr = nn.Upsample(size=x_hr.size()[2:], mode='bilinear', align_corners=True)(x_lr)

            x_hr = self.adapt_stage2_b2_joint_varout_dimred(x_hr)

            x = x_lr + x_hr
            x = F.relu(x)

        # CRP and RCU for fusion result
        x = self.mflow_conv_g1_pool(x)  # CRP
        x = self.rcu2(x)

        out = x
        return out

### added for RDFNet implementation
class MMFNet(nn.Module):

    # in_channels will be 512 for MMFNet-4 and 256 for others, as per RDFNet paper
    def __init__(self, in_channels):
        super(MMFNet, self).__init__()

        self.do = nn.Dropout(p=0.5)

        # these next blocks are the same for both RGB and HHA because
        #  the input volumes are exactly the same dimensions (including # channels)
        # conv3 is also used for the convolution before and after fusion

        # pre-fusion RGB blocks
        self.conv1_rgb = conv1x1(in_planes=in_channels, out_planes=in_channels)
        self.RCUs_rgb = nn.Sequential(  # 2 RCU blocks
            BasicBlock(inplanes=in_channels, planes=in_channels),  # expansion=1, no downsampling
            BasicBlock(inplanes=in_channels, planes=in_channels)   # expansion=1, no downsampling
        )
        self.conv3_rgb = conv3x3(in_planes=in_channels, out_planes=in_channels)

        # pre-fusion HHA blocks
        self.conv1_hha = conv1x1(in_planes=in_channels, out_planes=in_channels)
        self.RCUs_hha = nn.Sequential(  # 2 RCU blocks
            BasicBlock(inplanes=in_channels, planes=in_channels),  # expansion=1, no downsampling
            BasicBlock(inplanes=in_channels, planes=in_channels)   # expansion=1, no downsampling
        )
        self.conv3_hha = conv3x3(in_planes=in_channels, out_planes=in_channels)

        # post-fusion block
        #self.maxpool = nn.MaxPool2d(kernel_size=5, stride=1)
        #self.conv3 = conv3x3(in_planes=in_channels, out_planes=in_channels)
        self.crp = CRPBlock(in_planes=in_channels, out_planes=in_channels, n_stages=1)

    def forward(self, x_rgb, x_depth):
        # pre-fusion RGB
        x_rgb = self.do(x_rgb)
        x_rgb = self.conv1_rgb(x_rgb)
        x_rgb = self.RCUs_rgb(x_rgb)
        x_rgb = self.conv3_rgb(x_rgb)

        # pre-fusion HHA
        x_depth = self.do(x_depth)
        x_depth = self.conv1_hha(x_depth)
        x_depth = self.RCUs_hha(x_depth)
        x_depth = self.conv3_hha(x_depth)

        # fusion
        x = x_rgb + x_depth

        # post-fusion
        x = F.relu(x)
        #residual = x
        #x = self.maxpool(x)
        #x = self.conv3(x)
        #out = x + residual
        out = self.crp(x)

        return out

### This module now represents the RDFNet
###
### changes and additions to the original version of this class:
###  - reorganization and commenting for clarity
###  - RefineNet component isolated into its own RefineNet module: see above
###  - to convert from RefineNet to RDFNet:
###     - depth track added (all blocks ending in "_depth")
###     - MMFNet added to fuse RGB and depth tracks
class ResNet(nn.Module):

    def __init__(self, block, layers, num_classes=40):#21):
        super(ResNet, self).__init__()

        self.do = nn.Dropout(p=0.5)
        self.conv1 = nn.Conv2d(3, 64, kernel_size=7, stride=2, padding=3, bias=False)
        self.conv1_depth = nn.Conv2d(3, 64, kernel_size=7, stride=2, padding=3, bias=False)
        self.bn1 = nn.BatchNorm2d(64)
        self.bn1_depth = nn.BatchNorm2d(64)
        #self.relu = nn.ReLU(inplace=True) -- switched to F.relu to avoid errors
        self.maxpool = nn.MaxPool2d(kernel_size=3, stride=2, padding=1)

        # backbone
        self.inplanes = 64
        self.layer1 = self._make_layer(block, 64, layers[0])
        self.layer2 = self._make_layer(block, 128, layers[1], stride=2)
        self.layer3 = self._make_layer(block, 256, layers[2], stride=2)
        self.layer4 = self._make_layer(block, 512, layers[3], stride=2)
        self.inplanes = 64
        self.layer1_depth = self._make_layer(block, 64, layers[0])
        self.layer2_depth = self._make_layer(block, 128, layers[1], stride=2)
        self.layer3_depth = self._make_layer(block, 256, layers[2], stride=2)
        self.layer4_depth = self._make_layer(block, 512, layers[3], stride=2)

        # dimensionality reduction right off the backbone
        # note that for some reason, x4 uses the outl1 block, x3 uses outl2, etc
        # I would like to rename these as well, but I think the names are used for some other stuff
        self.p_ims1d2_outl1_dimred = conv1x1(2048, 512, bias=False)  # really l4
        self.p_ims1d2_outl2_dimred = conv1x1(1024, 256, bias=False)  # really l3
        self.p_ims1d2_outl3_dimred = conv1x1(512, 256, bias=False)   # really l2
        self.p_ims1d2_outl4_dimred = conv1x1(256, 256, bias=False)   # really l1
        self.p_ims1d2_outl1_dimred_depth = conv1x1(2048, 512, bias=False)  # really l4
        self.p_ims1d2_outl2_dimred_depth = conv1x1(1024, 256, bias=False)  # really l3
        self.p_ims1d2_outl3_dimred_depth = conv1x1(512, 256, bias=False)   # really l2
        self.p_ims1d2_outl4_dimred_depth = conv1x1(256, 256, bias=False)   # really l1

        # MMFNets
        self.MMFNet_l4 = MMFNet(in_channels=512)
        self.MMFNet_l3 = MMFNet(in_channels=256)
        self.MMFNet_l2 = MMFNet(in_channels=256)
        self.MMFNet_l1 = MMFNet(in_channels=256)

        # RefineNets
        self.RefineNet_l4 = RefineNet(inplanes_lr=512)  # no fusion step
        self.RefineNet_l4_l3 = RefineNet(inplanes_lr=512, inplanes_hr=256, fancy_upsample=True)
        self.RefineNet_l3_l2 = RefineNet(inplanes_lr=256, inplanes_hr=256, fancy_upsample=True)
        self.RefineNet_l2_l1 = RefineNet(inplanes_lr=256, inplanes_hr=256, fancy_upsample=True)

        # CLF convolutional step applied to layer 1 output to get class predictions
        self.clf_conv = nn.Conv2d(256, num_classes, kernel_size=3, stride=1,
                                  padding=1, bias=True)

    def _make_crp(self, in_planes, out_planes, stages):
        layers = [CRPBlock(in_planes, out_planes,stages)]
        return nn.Sequential(*layers)

    def _make_layer(self, block, planes, blocks, stride=1):
        downsample = None
        if stride != 1 or self.inplanes != planes * block.expansion:
            downsample = nn.Sequential(
                nn.Conv2d(self.inplanes, planes * block.expansion,
                          kernel_size=1, stride=stride, bias=False),
                nn.BatchNorm2d(planes * block.expansion),
            )

        layers = []
        layers.append(block(self.inplanes, planes, stride, downsample))
        self.inplanes = planes * block.expansion
        for i in range(1, blocks):
            layers.append(block(self.inplanes, planes))

        return nn.Sequential(*layers)

    def forward(self, x, x_depth):
        x = self.conv1(x)
        x = self.bn1(x)
        x = F.relu(x)
        x = self.maxpool(x)
        x_depth = self.conv1_depth(x_depth)
        x_depth = self.bn1_depth(x_depth)
        x_depth = F.relu(x_depth)
        x_depth = self.maxpool(x_depth)

        # backbone
        l1 = self.layer1(x)
        l2 = self.layer2(l1)
        l3 = self.layer3(l2)
        l4 = self.layer4(l3)
        l1_depth = self.layer1(x_depth)
        l2_depth = self.layer2(l1_depth)
        l3_depth = self.layer3(l2_depth)
        l4_depth = self.layer4(l3_depth)

        # not sure why this is only for l3 and l4, but it's in the original code
        # so I added the same thing added for depth
        l4 = self.do(l4)
        l3 = self.do(l3)
        l4_depth = self.do(l4_depth)
        l3_depth = self.do(l3_depth)

        # dimensionality reduction right off the backbone
        x4 = F.relu(self.p_ims1d2_outl1_dimred(l4))
        x3 = F.relu(self.p_ims1d2_outl2_dimred(l3))
        x2 = F.relu(self.p_ims1d2_outl3_dimred(l2))
        x1 = F.relu(self.p_ims1d2_outl4_dimred(l1))
        x4_depth = F.relu(self.p_ims1d2_outl1_dimred_depth(l4_depth))
        x3_depth = F.relu(self.p_ims1d2_outl2_dimred_depth(l3_depth))
        x2_depth = F.relu(self.p_ims1d2_outl3_dimred_depth(l2_depth))
        x1_depth = F.relu(self.p_ims1d2_outl4_dimred_depth(l1_depth))

        # MMFNets
        x4 = self.MMFNet_l4(x4, x4_depth)
        x3 = self.MMFNet_l3(x3, x3_depth)
        x2 = self.MMFNet_l2(x2, x2_depth)
        x1 = self.MMFNet_l1(x1, x1_depth)

        # RefineNets
        x4 = self.RefineNet_l4(x4)
        x3 = self.RefineNet_l4_l3(x4, x3)
        x2 = self.RefineNet_l3_l2(x3, x2)
        x1 = self.RefineNet_l2_l1(x2, x1)

        # CLF convolutional step to x1 to get class predictions
        out = self.clf_conv(x1)
        return out

### new: use ImageNet pretrained weights for keys containing "_depth" as well
### might help marginally
def use_pretrained_depth_track(model_dict, pretrained_model_dict):
    for key in model_dict:
        if "_depth" in key:
            pretrained_key = "".join(key.split("_depth"))
            if pretrained_key in pretrained_model_dict:
                model_dict[key] = pretrained_model_dict[pretrained_key]

### the functions below have been modified to call use_pretrained_depth_track

def rf_lw50(num_classes, imagenet=False, pretrained=True, **kwargs):
    model = ResNet(Bottleneck, [3, 4, 6, 3], num_classes=num_classes, **kwargs)
    if imagenet:
        key = '50_imagenet'
        url = models_urls[key]
        pretrained_model_dict = maybe_download(key, url)
        model.load_state_dict(pretrained_model_dict, strict=False)
        use_pretrained_depth_track(model.state_dict(), pretrained_model_dict) ### new
    elif pretrained:
        dataset = data_info.get(num_classes, None)
        if dataset:
            bname = '50_' + dataset.lower()
            key = 'rf_lw' + bname
            url = models_urls[bname]
            model.load_state_dict(maybe_download(key, url), strict=False)
    return model

def rf_lw101(num_classes, imagenet=False, pretrained=True, **kwargs):
    model = ResNet(Bottleneck, [3, 4, 23, 3], num_classes=num_classes, **kwargs)
    if imagenet:
        key = '101_imagenet'
        url = models_urls[key]
        pretrained_model_dict = maybe_download(key, url)
        model.load_state_dict(pretrained_model_dict, strict=False)
        use_pretrained_depth_track(model.state_dict(), pretrained_model_dict) ### new
    elif pretrained:
        dataset = data_info.get(num_classes, None)
        if dataset:
            bname = '101_' + dataset.lower()
            key = 'rf_lw' + bname
            url = models_urls[bname]
            model.load_state_dict(maybe_download(key, url), strict=False)
    return model

def rf_lw152(num_classes, imagenet=False, pretrained=True, **kwargs):
    model = ResNet(Bottleneck, [3, 8, 36, 3], num_classes=num_classes, **kwargs)
    if imagenet:
        key = '152_imagenet'
        url = models_urls[key]
        pretrained_model_dict = maybe_download(key, url)
        model.load_state_dict(pretrained_model_dict, strict=False)
        use_pretrained_depth_track(model.state_dict(), pretrained_model_dict) ### new
    elif pretrained:
        dataset = data_info.get(num_classes, None)
        if dataset:
            bname = '152_' + dataset.lower()
            key = 'rf_lw' + bname
            url = models_urls[bname]
            model.load_state_dict(maybe_download(key, url), strict=False)
    return model

max_stars_repo_name: egvincent/rgbd-sem-seg

max_issues_repo_name: egvincent/rgbd-sem-seg

max_forks_repo_name: egvincent/rgbd-sem-seg


text:
import unittest
import json
import os
import tempfile
import shutil
import urllib2

import numpy as np

from geopandas import GeoDataFrame

class TestDataFrame(unittest.TestCase):

    def setUp(self):
        # Data from http://www.nyc.gov/html/dcp/download/bytes/nybb_13a.zip
        # saved as geopandas/examples/nybb_13a.zip.
        if not os.path.exists(os.path.join('examples', 'nybb_13a.zip')):
            with open(os.path.join('examples', 'nybb_13a.zip'), 'w') as f:
                response = urllib2.urlopen('http://www.nyc.gov/html/dcp/download/bytes/nybb_13a.zip')
                f.write(response.read())
        self.df = GeoDataFrame.from_file(
            '/nybb_13a/nybb.shp', vfs='zip://examples/nybb_13a.zip')
        self.tempdir = tempfile.mkdtemp()
        self.boros = np.array(['Staten Island', 'Queens', 'Brooklyn',
                               'Manhattan', 'Bronx'])

    def tearDown(self):
        shutil.rmtree(self.tempdir)

    def test_from_file_(self):
        self.assertTrue('geometry' in self.df)
        self.assertTrue(len(self.df) == 5)
        self.assertTrue(np.alltrue(self.df['BoroName'].values == self.boros))

    def test_to_json(self):
        text = self.df.to_json()
        data = json.loads(text)
        self.assertTrue(data['type'] == 'FeatureCollection')
        self.assertTrue(len(data['features']) == 5)

    def test_to_file(self):
        tempfilename = os.path.join(self.tempdir, 'boros.shp')
        self.df.to_file(tempfilename)
        # Read layer back in?
        df = GeoDataFrame.from_file(tempfilename)
        self.assertTrue('geometry' in df)
        self.assertTrue(len(df) == 5)
        self.assertTrue(np.alltrue(df['BoroName'].values == self.boros))

max_stars_repo_name: jwass/geopandas

max_issues_repo_name: jwass/geopandas

max_forks_repo_name: jwass/geopandas


text:
# -*- coding: utf-8 -*-
"""

"""
import os
from osgeo import gdal
import numpy as np
import sys
import subprocess
from sklearn.ensemble import RandomForestClassifier
from sklearn import metrics
import matplotlib.pyplot as plt
import raster_array_funcspy35 as ra

site = "Site4"
#global variarbles
no_data_flag = -9999
pos_val = 0 #pixel value for the positive class
neg_val = 1 #pixel vlaue for the negative class
pos_prop = 0.5 #percent of positive class to sample for training
neg_prop = 0.5 #percent of negative class to sample for trainig

#set filepaths

outdir = r"D:\3rdStudy_ONeil\Results\RF\{}".format(site)
data_tif = os.path.join(outdir, "{}_comp4.tif".format(site))
#test_labels = r"D:\2ndStudy_ONeil\Tool_testing\results\Site4\Run3\test_w0.850_nw0.920.tif"
test_labels = r"D:\3rdStudy_ONeil\Results\RF\Site4\test_labels_DL_test.tif"
train_labels = r"D:\2ndStudy_ONeil\Tool_testing\results\Site4\Run3\train_w0.150_nw0.080.tif"
f_path_in = os.path.dirname(data_tif)
f_in = os.path.basename(data_tif)
rf_results_name = "{}_rf_predict.tif".format(site)

wdir = r"D:\2ndStudy_ONeil\Tool_testing\data\{}\input_variables".format(site)
twi = os.path.join(wdir, "pmSLP50_0.9p_grass_at_twi.tif")
curv = os.path.join(wdir, "pmSLP50_0.9p_curv.tif")
dtw = os.path.join(wdir, "pmSLP50_0.9p_slp_grassrcost_dtw.tif")
ndvi = r"D:\3rdStudy_ONeil\Results\RF\Site4\NVDI.tif"
###flist1 = [twi, curv, dtw, ndvi]
#roi_lim = r"D:\2ndStudy_ONeil\Tool_testing\data\Site2\lims_prj.shp"
#arr, meta = ra.gtiff_to_arr(dtw, "float")
#print(np.min(arr))
#arr_ma = np.ma.masked_invalid(arr)
#print(np.min(arr_ma))
#arr_ma2 = np.ma.masked_values(arr, -9999.)
#print(np.min(arr_ma2))
#arr_out = np.ma.filled(arr_ma, -9999.)
#print(np.min(arr_out))
#arr_tif = ra.arr_to_gtiff(arr, meta, outdir, "dtw.tif")
#sys.exit(0)
##clip = ra.clip_geotif(dtw, outdir, roi_lim)
#sys.exit(0)
#flist=[]
#for f in flist1:
#    clip = ra.clip_geotif(f, outdir, roi_lim)
#    flist.append(clip)
#
#comp = comp_bands(outdir, flist, "Site1_comp.tif")
#
##clip = ra.clip_geotif(comp, outdir, roi_lim)
#sys.exit(0)

def create_tt_feats(input_feat, train_labels_2d, test_labels_2d, out_dir):

    feat_arr, feat_meta = ra.gtiff_to_arr(input_feat, "float")
#    for b in range(4):
#        print(np.min(feat_arr[:, :, b]))
#
#        band = np.ma.masked_values(feat_arr[:, :, b], -9999.)
#        print(np.min(band))
#
#        band1 = np.ma.filled(band, np.nan)
#        print(np.min(band1))
#        feat_arr[:, :, b] = band1

    train_fname = os.path.join(out_dir, "{}_train_feats.tif".format(site))
    test_fname = os.path.join(out_dir, "{}_test_feats.tif".format(site))

    #create arrays of input variables within training and testing limits
    train_features = np.copy(feat_arr)
    test_features = np.copy(feat_arr)
    print(np.min(train_labels_2d))
    print(np.min(np.ma.masked_invalid(train_labels_2d)))

    train_features[np.isnan(train_labels_2d), : ] = -9999.
    test_features[np.isnan(test_labels_2d), : ] = -9999.

    _ = ra.arr_to_gtiff(train_features, feat_meta, out_dir, train_fname)
    _ = ra.arr_to_gtiff(test_features, feat_meta, out_dir, test_fname)

    return train_features, test_features

train_labels_2d, meta = ra.gtiff_to_arr(train_labels, dtype="float")
test_labels_2d, meta = ra.gtiff_to_arr(test_labels, dtype="float")

#train_features, test_features = create_tt_feats(data_tif, train_labels_2d, test_labels_2d, outdir)
#train_features, t_meta = ra.gtiff_to_arr(os.path.join(outdir, "{}_train_feats.tif".format(site)), "float")
#test_features, te_meta = ra.gtiff_to_arr(os.path.join(outdir, "{}_test_feats.tif".format(site)), "float")

def classify(train_features, train_labels_2d, test_features, n_trees=300, tree_depth=None, max_f = 'auto'):

    train_features1 = np.copy(train_features)
    for b in range(4):
        print(np.min(train_features[:, :, b]))

        band = np.ma.masked_invalid(train_features[:, :, b])
        print(np.min(band))

        band1 = np.ma.filled(band, -9999.)
        print(np.min(band1))
        train_features1[:, :, b] = band1

#    #train RF model
    train_X = train_features1[~np.ma.masked_values(train_features1[:, :, 0], -9999).mask]
    print(np.min(train_X))
    train_Y = train_labels_2d[~np.ma.masked_values(train_features1[:, :, 0], -9999).mask]
    print(np.min(train_Y))
    print(np.shape(train_X), np.shape(train_Y))
    #train_X is an array of nonNaN values from train_features, 3D array --> 2D array ([X*Y, Z])
    #train_Y is an array of nonNaN values from train_labels, 2D array --> 1D array (X*Y)
    #both have NO NaN VALUES, which is required for sklearn rf

    #Initialize RF model
    print ("\nInitializing Random Forest model.....\n")

    rf_clf = RandomForestClassifier( n_estimators = n_trees, max_depth = tree_depth, \
                                    oob_score = True, n_jobs = -1, class_weight = 'balanced',\
                                    max_features = max_f, random_state = 21 )

    print ("Training model...\n")
    rf_fit = rf_clf.fit(train_X, train_Y)

    #save feature importance
    n_feats = len(rf_fit.feature_importances_)
    importance=[]
    bands = np.arange(1, n_feats+1)

    for b, imp in zip(bands, rf_fit.feature_importances_):
        print('Band {b} importance: {imp}'.format(b=b, imp=imp.round(3)))
        importance.append('Band {b} importance: {imp} \n'.format(b=b, imp=imp.round(3)))

    #execute RF classification and fuzzy classification

    test_feat_masked = np.ma.masked_invalid(test_features)
    test_feat_nonan = np.ma.filled(test_feat_masked, -9999.) #still 3D
    temp_shape = ( test_feat_nonan.shape[0] * test_feat_nonan.shape[1], test_feat_nonan.shape[2] )

    test_shaped = test_feat_nonan[:, :, :].reshape(temp_shape)

    print ("Executing prediction...\n")
    rf_predict = rf_fit.predict(test_shaped)
    fuzzy_predict = rf_fit.predict_proba(test_shaped)

    #reshape outputs to prepare for geotif export
    rf_predict = rf_predict.reshape(test_features[:, :, 0].shape)
    fuzzy_predict_pos = fuzzy_predict[:,0].reshape(test_features[:, :, 0].shape)

    return rf_predict, fuzzy_predict_pos, importance

#rf_predict, fuzzy_predict_pos, importance = classify(train_features, train_labels_2d, test_features)
#rf_tif = ra.arr_to_gtiff(rf_predict, meta, outdir, "{}_rf_predict.tif".format(site))
#rf_tif = os.path.join(outdir, "{}_rf_predict.tif".format(site))
rf_tif = r"D:\3rdStudy_ONeil\Results\RF\Site4\rf_predict_DL_test.tif"
rf_predict, _ = ra.gtiff_to_arr(rf_tif, "int")

def get_acc(test_labels_2d, rf_predict, out_dir, pos_val=0, fuzzy_predict_pos=None, importance=None):
    """
    :return: various accuracy metrics
    """
    #true labels
    tl_flat = test_labels_2d.ravel()
    tl = tl_flat[~np.isnan(tl_flat)]

    #predicted values
    pv_flat = rf_predict.ravel()
    pv = pv_flat[~np.isnan(tl_flat)]

    print ("Performing accuracy assessment... \n")

    #compute stats
    conf_matrix= metrics.confusion_matrix(tl, pv)
    acc_score= metrics.accuracy_score(tl, pv)
    class_report = metrics.classification_report(tl, pv)
    iou = metrics.jaccard_score(tl, pv, average=None)

    print("Confusion Matrix:\n")
    print(conf_matrix)
    print("----------------------")

    print("Accuracy Score:\n")
    print(acc_score)
    print("----------------------")

    print("Classification report:\n")
    print(class_report)
    print("----------------------")

    print("Jaccard Score (IoU):\n")
    print(iou)
    print("----------------------")

    return conf_matrix, class_report

conf_matrix, class_report = get_acc(test_labels_2d, rf_predict, outdir)

#*******************************run code**************************************#
#
#def run():
#
#    #open geotiffs
#    array, meta = gtiff_to_arr(f_path_in, f_in)
#    verif_arr, verif_meta = gtiff_to_arr(os.path.dirname(verification), os.path.basename(verification))
#
#    #use the verification geotiff to create training and testing files, input percentage from each class to sample for training
#
#    train_labels, test_labels = create_tt_labels(os.path.dirname(verification), os.path.basename(verification), \
#                     pos_prop=pos_prop, neg_prop=neg_prop, pos_val=pos_val, neg_val=neg_val, out_dir=outdir, no_data=no_data_flag)
#
#    #use train and test labels to index training and testing pixels from data raster
#    train_features, test_features = create_tt_feats(array, meta, train_labels, test_labels, outdir)
#
#    #execute classification
#    rf_predict, fuzzy_predict_pos, importance = classify(train_features, train_labels, test_features)
#
#    rf_out = arr_to_gtiff(rf_predict, verif_meta, outdir, rf_results_name, dtype='int')
#    fuz_out = arr_to_gtiff(fuzzy_predict_pos, verif_meta, outdir, fuz_results_name, dtype='int')
#
#    #get classifier accuracy
#    get_acc(test_labels, rf_predict, fuzzy_predict_pos, importance, outdir, meta, 0)
#
#    return rf_out, fuz_out
#
#rf_out, fuz_out = run()


text:
#!/mnt/data1/Matt/anaconda_install/anaconda2/envs/matt_EMAN2/bin/python
################################################################################
# imports
import argparse; import sys
parser = argparse.ArgumentParser('Generate noisy and noiseless 2D projections of randomly oriented and translated MRC files.')
parser.add_argument('--input_mrc_dir', type=str, help='input directory containing MRC files to be rotated, translated, CTF-convolved, and projected')
parser.add_argument('--output_noise_dir', type=str, help='output directory to store noisy 2D projections')
parser.add_argument('--output_noiseless_dir', type=str, help='output directory to store noiseless 2D projections')
parser.add_argument('--output_semMap_dir', type=str, help='output directory to store semantic segmentation targets')
parser.add_argument('--numProjs', type=int, help='total number of projections to make')
parser.add_argument('--nProcs', type=int, help='total number of parallel threads to launch')
args = parser.parse_args()
print('')
if(args.input_mrc_dir == None or args.output_noise_dir == None or args.output_noiseless_dir == None):
	print('Please enter an input_mrc_dir, AND an output_noise_dir, AND an output_noiseless_dir')
	sys.exit('The preferred input style may be found with ./projection_generator.py -h')

if(args.numProjs == None):
	sys.exit('Please enter the number of projection images you would like with the --numProjs flag')

if(args.nProcs == None):
	print('No process number specified, using one thread')
	nProcs = 1
else:
	nProcs = args.nProcs

if(args.numProjs % nProcs != 0):
	print('The numProjs that you specified was not a multiple of nProcs.')
	print('Instead of %d 2D projections, this program will generate %d 2D projections'%(args.numProjs, args.numProjs/nProcs*nProcs))

folder = args.input_mrc_dir
noNoise_outputDir = args.output_noiseless_dir
noise_outputDir = args.output_noise_dir
semMap_outputDir = args.output_semMap_dir
TOTAL_NUM_TO_MAKE = args.numProjs

if(folder[-1] != '/'): folder = folder + '/'
if(noNoise_outputDir[-1] != '/'): noNoise_outputDir = noNoise_outputDir + '/'
if(noise_outputDir[-1] != '/'): noise_outputDir = noise_outputDir + '/'
if(semMap_outputDir[-1] != '/'): semMap_outputDir = semMap_outputDir + '/'
print('The program will now generate %d 2D projections'%(args.numProjs/args.nProcs*args.nProcs))
################################################################################
# import of python packages
import numpy as np
from EMAN2 import *; from sparx import *; import mrcfile
import json; import glob
from multiprocessing import Pool
import os; from tqdm import tqdm
from skimage.morphology import erosion,dilation; from skimage.morphology import disk
################################################################################
################################################################################
# import data
actin_orig = []
file_names = sorted(os.listdir(folder))
for file_name in file_names:
	if(file_name[-4:] == '.mrc' and file_name[:10] == 'bent_actin'):
		actin_orig.append(EMData(folder+file_name))

################################################################################
# image sizes
box_len = 384; BL = 96; cropBox = BL*2
################################################################################
def launch_parallel_process(thread_idx):
	index=num_per_proc*thread_idx
	for i in tqdm(range(0,num_per_proc)):
		local_random_state = np.random.RandomState(None)
		# First: randomly pick one of the actin mrc files that were loaded into actin_orig
		r0 = local_random_state.randint(0,len(actin_orig))
		num_filaments = local_random_state.randint(0,4)
		if(num_filaments == 0):
			target = np.concatenate((np.ones((1,cropBox,cropBox)), np.zeros((2,cropBox,cropBox))), axis=0)
			with mrcfile.new(semMap_outputDir + 'actin_rotated%05d.mrcs'%(i+num_per_proc*thread_idx), overwrite=True) as mrc:
				mrc.set_data(target.astype('float32'))

			r7 = local_random_state.uniform(1.0, 4.0) #defocus
			r8 = max(local_random_state.normal(0.050, 0.010),0) # noise amplitude
			target_eman = EMNumPy.numpy2em(np.zeros((cropBox,cropBox)))
			target_eman.process_inplace('math.simulatectf',{'ampcont':10.0,'apix':4.12,'bfactor':0.0,'cs':2.7,'defocus':r7,'noiseamp':r8,'purectf':False,'voltage':300.0})
			target_noise = EMNumPy.em2numpy(target_eman)
			with mrcfile.new(noise_outputDir + 'actin_rotated%05d.mrc'%(i+num_per_proc*thread_idx), overwrite=True) as mrc:
				mrc.set_data(target_noise.astype('float32'))
			with mrcfile.new(noNoise_outputDir + 'actin_rotated%05d.mrc'%(i+num_per_proc*thread_idx), overwrite=True) as mrc:
				mrc.set_data(target[1].astype('float32'))

		else:
			target = np.zeros((5,cropBox,cropBox)); target_lp50 = np.zeros((5,cropBox,cropBox))
			bundle_idxs = []
			for j in range(1, num_filaments+1):
				r0_name = file_names[r0]
				rotated_actin = actin_orig[r0].copy()
				# handle bundles
				r9 = local_random_state.uniform()
				fil_is_bundle = r9 < 0.65
				if(fil_is_bundle): # i.e. if this filament-unit is a bundle
					r10 = int(local_random_state.normal(loc=90, scale=1.5) + int(local_random_state.randint(0,2)*180))-90 #
					r11 = int(local_random_state.random_sample()*360) #
					r12 = local_random_state.normal(loc=36.8, scale=3.2) # bundle distance, measured by Alfred
					r13 = local_random_state.uniform(-44,44) # should cover full crossover
					rotated_actin2 = actin_orig[r0].copy()
					t2 = Transform()
					t2.set_params({'type':'eman', 'az':0, 'alt':r10, 'phi':r11, 'tx':0, 'ty':r12, 'tz':r13}) #tz is shift up and down filament axis # y and x control distance, alt is (anti-)parallel
					rotated_actin2.transform(t2)
					rotated_actin.add(rotated_actin2)
					bundle_idxs.append(j)

				# Rotation angles: azimuth, alt, phi, then Translations: tx, ty,tz
				r1, r2, r3 = int(local_random_state.random_sample()*360),int(local_random_state.normal(loc=90, scale=2.5)),int(local_random_state.random_sample()*360)
				#r4, r5, r6 = local_random_state.normal(0, 25), local_random_state.normal(0, 25), local_random_state.normal(0, 25)
				r4, r5, r6 = local_random_state.uniform(-60, 60), local_random_state.uniform(-60, 60), local_random_state.uniform(-60, 60)
				t = Transform()
				t.set_params({'type':'eman','az':r1, 'alt':r2, 'phi':r3, 'tx':r4, 'ty':r5, 'tz':r6})
				rotated_actin.transform(t) # apply rotation and translation
				proj_eman = rotated_actin.project('standard',Transform()) # project
				proj_eman_lp50 = proj_eman.process('filter.lowpass.gauss', {'apix':4.12, 'cutoff_freq':1.0/40.0})
				proj_np_lp50 = EMNumPy.em2numpy(proj_eman_lp50)
				proj_np_lp50 = dilation(proj_np_lp50, selem=disk(8))
				proj_np_lp50 = erosion(proj_np_lp50, selem=disk(16))
				proj_np = EMNumPy.em2numpy(proj_eman)
				center = box_len/2
				# Save the target image
				target_filament = proj_np[center-BL:center+BL, center-BL:center+BL]
				#target_filament = (target_filament - np.mean(target_filament)) / np.std(target_filament) # normalize
				target[j] = target_filament

				target_filament_lp50 = proj_np_lp50[center-BL:center+BL, center-BL:center+BL]
				target_filament_lp50 = (target_filament_lp50 - np.mean(target_filament_lp50)) / np.std(target_filament_lp50) # normalize
				target_lp50[j] = target_filament_lp50

			# Generate noisy image
			target_proj = np.sum(target, axis=0)
			target_proj = (target_proj - np.mean(target_proj)) / np.std(target_proj)
			r7 = local_random_state.uniform(1.0, 4.0) #defocus
			r8 = max(local_random_state.normal(0.027, 0.005),0) # noise amplitude 0.050, 0.010 is default
			target_eman = EMNumPy.numpy2em(target_proj)
			target_eman.process_inplace('math.simulatectf',{'ampcont':10.0,'apix':4.12,'bfactor':0.0,'cs':2.7,'defocus':r7,'noiseamp':r8,'purectf':False,'voltage':300.0})
			target_noise = EMNumPy.em2numpy(target_eman)
			target_noise = (target_noise - np.mean(target_noise)) / np.std(target_noise) # normalize
			with mrcfile.new(noise_outputDir + 'actin_rotated%05d.mrc'%(i+num_per_proc*thread_idx), overwrite=True) as mrc:
				mrc.set_data(target_noise.astype('float32'))
			with mrcfile.new(noNoise_outputDir + 'actin_rotated%05d.mrc'%(i+num_per_proc*thread_idx), overwrite=True) as mrc:
				mrc.set_data(target_proj.astype('float32'))

			# make masks from each proj
			for j in range(1, num_filaments+1):
				target_lp50[j] = (target_lp50[j]>0)
			temp = np.sum(target_lp50, axis=0)

			semMap = np.array(np.zeros((3,cropBox,cropBox)))
			semMap[0] = (temp==0)>0
			#semMap[1] = -1*semMap[0] + 1
			semMap[2] = np.sum(target_lp50[bundle_idxs], axis=0) > 0
			semMap[1] = -1*((semMap[0] + semMap[2]) - 1)
			with mrcfile.new(semMap_outputDir + 'actin_rotated%05d.mrcs'%(i+num_per_proc*thread_idx), overwrite=True) as mrc:
				mrc.set_data(semMap.astype('float32'))

################################################################################
# run in parallel
num_per_proc = TOTAL_NUM_TO_MAKE / nProcs
if __name__ == '__main__':
	p=Pool(nProcs)
	p.map(launch_parallel_process, range(0, nProcs))

################################################################################
# Now all files are written, combine all json files into one master json file
read_files = glob.glob(noise_outputDir+'params_*.json')
output_list = []
for f in read_files:
	for line in open(f, 'r'):
		output_list.append(json.loads(line))

#sort the json dictionaries based on iteration number
output_list = sorted(output_list, key=lambda i: i['iteration'])
for line in output_list:
	with open(noise_outputDir+'master_params.json', 'a') as fp:
		data_to_write = json.dumps(line)
		fp.write(data_to_write + '\n')


text:
import numpy as np
import tensorflow as tf
from tensorflow.keras import layers

from nfp.preprocessing import zero, SmilesPreprocessor
from nfp.preprocessing.features import Tokenizer

class CifPreprocessor(SmilesPreprocessor):
    def __init__(self, radius=None, num_neighbors=12):
        self.site_tokenizer = Tokenizer()
        self.radius = radius
        self.num_neighbors = num_neighbors
        self.max_sites = 0
        self.max_distance = 0

    @property
    def site_classes(self):
        return self.site_tokenizer.num_classes + 1

    @staticmethod
    def site_features(site):
        species = site.as_dict()['species']
        assert len(species) == 1
        return species[0]['element']

    def construct_feature_matrices(self, crystal, train=False):

        self.site_tokenizer.train = train
        site_features = np.array([self.site_tokenizer(self.site_features(site))
                                  for site in crystal.sites])

        # Record for array sizing later
        if train & (crystal.num_sites > self.max_sites):
            self.max_sites = crystal.num_sites

        connectivity = []
        distances = []

        if self.radius is None:
            # Get the expected number of sites / volume, then find a radius
            # expected to yield 2x the desired number of neighbors
            desired_vol = (crystal.volume / crystal.num_sites) * self.num_neighbors
            radius = 2 * (desired_vol / (4 * np.pi / 3))**(1/3)
        else:
            radius = self.radius

        # Iterate over each site's neighbors
        for i, neighbors in enumerate(crystal.get_all_neighbors(radius)):
            assert len(neighbors) >= self.num_neighbors, \
                f"Only {len(neighbors)} neighbors for site {i}"

            sorted_neighbors = sorted(neighbors, key=lambda x: x[1])[:self.num_neighbors]
            for _, distance, j, _ in sorted_neighbors:
                connectivity += [(i, j)]
                distances += [distance]

                if train & (distance > self.max_distance):
                    self.max_distance = distance

        connectivity = np.array(connectivity, dtype='int')
        distances = np.array(distances, dtype='float')

        return {'site_features': site_features,
                'distance': distances,
                'connectivity': connectivity}

    tfrecord_features = {
        'site_features': tf.io.FixedLenFeature([], dtype=tf.string),
        'distance': tf.io.FixedLenFeature([], dtype=tf.string),
        'connectivity': tf.io.FixedLenFeature([], dtype=tf.string)
    }

    output_types = {'site_features': tf.int64,
                    'distance': tf.float64,
                    'connectivity': tf.int64}

    output_shapes = {'site_features': tf.TensorShape([None]),
                     'distance': tf.TensorShape([None]),
                     'connectivity': tf.TensorShape([None, 2])}

    @staticmethod
    def padded_shapes(max_sites=-1, max_bonds=-1):

        return {
            'site_features': [max_sites],
            'distance': [max_bonds],
            'connectivity': [max_bonds, 2]
        }

    padding_values = {
        'site_features': zero,
        'distance': tf.constant(np.nan, dtype=tf.float64),
        'connectivity': zero}

class RBFExpansion(layers.Layer):
    def __init__(self, dimension=128, init_gap=10, init_max_distance=7, trainable=False):
        super(RBFExpansion, self).__init__()
        self.init_gap = init_gap
        self.init_max_distance = init_max_distance
        self.dimension = dimension
        self.trainable = trainable

    def build(self, input_shape):
        self.centers = tf.Variable(
            name='centers',
            initial_value=tf.range(0, self.init_max_distance,
                                   delta=self.init_max_distance/self.dimension),
            trainable=self.trainable,
            dtype=tf.float32)

        self.gap = tf.Variable(
            name='gap',
            initial_value=tf.constant(self.init_gap, dtype=tf.float32),
            trainable=self.trainable,
            dtype=tf.float32)

    def call(self, inputs):
        distances = tf.where(
            tf.math.is_nan(inputs), tf.zeros_like(inputs, dtype=inputs.dtype), inputs)
        offset = tf.expand_dims(distances, -1) - tf.cast(self.centers, inputs.dtype)
        logits = -self.gap * (offset)**2
        return tf.exp(logits)

    def compute_mask(self, inputs, mask=None):
        return tf.logical_not(tf.math.is_nan(inputs))

    def get_config(self):
        return {'init_gap': self.init_gap,
                'init_max_distance': self.init_max_distance,
                'dimension': self.dimension,
                'trainable': self.trainable}

max_stars_repo_name: harrysorensennrel/rlmolecule

max_issues_repo_name: harrysorensennrel/rlmolecule

max_forks_repo_name: harrysorensennrel/rlmolecule


text:
import numpy as np
import _pickle as cPickle
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import cv2
import random
import tensorflow as tf

from utilities import label_img_to_color

from model import ENet_model

project_dir = "/home/pikl/segmentation/"

data_dir = project_dir + "data/"

# change this to not overwrite all log data when you train the model:
model_id = "1"

batch_size = 4
img_height = 512
img_width = 1024

model = ENet_model(model_id, img_height=img_height, img_width=img_width,
            batch_size=batch_size)

no_of_classes = model.no_of_classes

# load the mean color channels of the train imgs:
train_mean_channels = cPickle.load(open("data/mean_channels.pkl"))

# load the training data from disk:
train_img_paths = cPickle.load(open(data_dir + "train_img_paths.pkl"))
train_trainId_label_paths = cPickle.load(open(data_dir + "train_trainId_label_paths.pkl"))
train_data = zip(train_img_paths, train_trainId_label_paths)

# compute the number of batches needed to iterate through the training data:
no_of_train_imgs = len(train_img_paths)
no_of_batches = int(no_of_train_imgs/batch_size)

# load the validation data from disk:
val_img_paths = cPickle.load(open(data_dir + "val_img_paths.pkl"))
val_trainId_label_paths = cPickle.load(open(data_dir + "val_trainId_label_paths.pkl"))
val_data = zip(val_img_paths, val_trainId_label_paths)

# compute the number of batches needed to iterate through the val data:
no_of_val_imgs = len(val_img_paths)
no_of_val_batches = int(no_of_val_imgs/batch_size)

# define params needed for label to onehot label conversion:
layer_idx = np.arange(img_height).reshape(img_height, 1)
component_idx = np.tile(np.arange(img_width), (img_height, 1))

def evaluate_on_val():
    random.shuffle(val_data)
    val_img_paths, val_trainId_label_paths = zip(*val_data)

    val_batch_losses = []
    batch_pointer = 0
    for step in range(no_of_val_batches):
        batch_imgs = np.zeros((batch_size, img_height, img_width, 3), dtype=np.float32)
        batch_onehot_labels = np.zeros((batch_size, img_height, img_width,
                    no_of_classes), dtype=np.float32)

        for i in range(batch_size):
            # read the next img:
            img = cv2.imread(val_img_paths[batch_pointer + i], -1)
            img = img - train_mean_channels
            batch_imgs[i] = img

            # read the next label:
            trainId_label = cv2.imread(val_trainId_label_paths[batch_pointer + i], -1)

            # convert the label to onehot:
            onehot_label = np.zeros((img_height, img_width, no_of_classes), dtype=np.float32)
            onehot_label[layer_idx, component_idx, trainId_label] = 1
            batch_onehot_labels[i] = onehot_label

        batch_pointer += batch_size

        batch_feed_dict = model.create_feed_dict(imgs_batch=batch_imgs,
                    early_drop_prob=0.0, late_drop_prob=0.0,
                    onehot_labels_batch=batch_onehot_labels)

        # run a forward pass, get the batch loss and the logits:
        batch_loss, logits = sess.run([model.loss, model.logits],
                    feed_dict=batch_feed_dict)

        val_batch_losses.append(batch_loss)
        print ("epoch: %d/%d, val step: %d/%d, val batch loss: %g" % (epoch+1,
                    no_of_epochs, step+1, no_of_val_batches, batch_loss))

        if step < 4:
            # save the predicted label images to disk for debugging and
            # qualitative evaluation:
            predictions = np.argmax(logits, axis=3)
            for i in range(batch_size):
                pred_img = predictions[i]
                label_img_color = label_img_to_color(pred_img)
                cv2.imwrite((model.debug_imgs_dir + "val_" + str(epoch) + "_" +
                            str(step) + "_" + str(i) + ".png"), label_img_color)

    val_loss = np.mean(val_batch_losses)
    return val_loss

def train_data_iterator():
    random.shuffle(train_data)
    train_img_paths, train_trainId_label_paths = zip(*train_data)

    batch_pointer = 0
    for step in range(no_of_batches):
        # get and yield the next batch_size imgs and onehot labels from the train data:
        batch_imgs = np.zeros((batch_size, img_height, img_width, 3), dtype=np.float32)
        batch_onehot_labels = np.zeros((batch_size, img_height, img_width,
                    no_of_classes), dtype=np.float32)

        for i in range(batch_size):
            # read the next img:
            img = cv2.imread(train_img_paths[batch_pointer + i], -1)
            img = img - train_mean_channels
            batch_imgs[i] = img

            # read the next label:
            trainId_label = cv2.imread(train_trainId_label_paths[batch_pointer + i], -1)

            # convert the label to onehot:
            onehot_label = np.zeros((img_height, img_width, no_of_classes), dtype=np.float32)
            onehot_label[layer_idx, component_idx, trainId_label] = 1
            batch_onehot_labels[i] = onehot_label

        batch_pointer += batch_size

        yield (batch_imgs, batch_onehot_labels)

no_of_epochs = 100

# create a saver for saving all model variables/parameters:
saver = tf.train.Saver(tf.trainable_variables(), write_version=tf.train.SaverDef.V2)

# initialize all log data containers:
train_loss_per_epoch = []
val_loss_per_epoch = []

# initialize a list containing the 5 best val losses (is used to tell when to
# save a model checkpoint):
best_epoch_losses = [1000, 1000, 1000, 1000, 1000]
print("Num GPUs Available: ", len(tf.config.experimental.list_physical_devices('GPU')))
tf.debugging.set_log_device_placement(True)
with tf.Session() as sess:
    # initialize all variables/parameters:
    init = tf.global_variables_initializer()
    sess.run(init)

    for epoch in range(no_of_epochs):
        print ("###########################")
        print ("######## NEW EPOCH ########")
        print ("###########################")
        print ("epoch: %d/%d" % (epoch+1, no_of_epochs))

        # run an epoch and get all batch losses:
        batch_losses = []
        for step, (imgs, onehot_labels) in enumerate(train_data_iterator()):
            # create a feed dict containing the batch data:
            batch_feed_dict = model.create_feed_dict(imgs_batch=imgs,
                        early_drop_prob=0.01, late_drop_prob=0.1,
                        onehot_labels_batch=onehot_labels)

            # compute the batch loss and compute & apply all gradients w.r.t to
            # the batch loss (without model.train_op in the call, the network
            # would NOT train, we would only compute the batch loss):
            batch_loss, _ = sess.run([model.loss, model.train_op],
                        feed_dict=batch_feed_dict)
            batch_losses.append(batch_loss)

            print ("step: %d/%d, training batch loss: %g" % (step+1, no_of_batches, batch_loss))

        # compute the train epoch loss:
        train_epoch_loss = np.mean(batch_losses)
        # save the train epoch loss:
        train_loss_per_epoch.append(train_epoch_loss)
        # save the train epoch losses to disk:
        cPickle.dump(train_loss_per_epoch, open("%strain_loss_per_epoch.pkl"
                    % model.model_dir, "w"))
        print ("training loss: %g" % train_epoch_loss)

        # run the model on the validation data:
        val_loss = evaluate_on_val()

        # save the val epoch loss:
        val_loss_per_epoch.append(val_loss)
        # save the val epoch losses to disk:
        cPickle.dump(val_loss_per_epoch, open("%sval_loss_per_epoch.pkl"\
                    % model.model_dir, "w"))
        print ("validation loss: %g" % val_loss)

        if val_loss < max(best_epoch_losses): # (if top 5 performance on val:)
            # save the model weights to disk:
            checkpoint_path = (model.checkpoints_dir + "model_" +
                        model.model_id + "_epoch_" + str(epoch + 1) + ".ckpt")
            saver.save(sess, checkpoint_path)
            print ("checkpoint saved in file: %s" % checkpoint_path)

            # update the top 5 val losses:
            index = best_epoch_losses.index(max(best_epoch_losses))
            best_epoch_losses[index] = val_loss

        # plot the training loss vs epoch and save to disk:
        plt.figure(1)
        plt.plot(train_loss_per_epoch, "k^")
        plt.plot(train_loss_per_epoch, "k")
        plt.ylabel("loss")
        plt.xlabel("epoch")
        plt.title("training loss per epoch")
        plt.savefig("%strain_loss_per_epoch.png" % model.model_dir)
        plt.close(1)

        # plot the val loss vs epoch and save to disk:
        plt.figure(1)
        plt.plot(val_loss_per_epoch, "k^")
        plt.plot(val_loss_per_epoch, "k")
        plt.ylabel("loss")
        plt.xlabel("epoch")
        plt.title("validation loss per epoch")
        plt.savefig("%sval_loss_per_epoch.png" % model.model_dir)
        plt.close(1)

max_stars_repo_name: alekpikl/segmentation

max_issues_repo_name: alekpikl/segmentation

max_forks_repo_name: alekpikl/segmentation


text:
import sys
import numpy as np

data = np.loadtxt(sys.argv[1])
print(np.average(data))
print(np.median(data))
print(np.min(data))
print(np.max(data))
print(np.percentile(data, 95))
print(np.percentile(data, 99))
print(np.percentile(data, 99.9))

max_stars_repo_name: stefan-ctrl/paxi

max_issues_repo_name: stefan-ctrl/paxi

max_forks_repo_name: stefan-ctrl/paxi


text:
#%%
"""
Created on July 12 2021
Shifted GBM and pricing of caplets/floorlets

This code is purely educational and comes from "Financial Engineering" course by L.A. Grzelak
The course is based on the book “Mathematical Modeling and Computation
in Finance: With Exercises and Python and MATLAB Computer Codes”,
by C.W. Oosterlee and L.A. Grzelak, World Scientific Publishing Europe Ltd, 2019.
@author: Lech A. Grzelak
"""

import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as st
import enum
import scipy.optimize as optimize

# This class defines puts and calls
class OptionType(enum.Enum):
    CALL = 1.0
    PUT = -1.0

def GeneratePathsGBMShifted(NoOfPaths,NoOfSteps,T,r,sigma,S_0,shift):
    S_0shift = S_0 + shift
    if (S_0shift < 0.0):
        raise('Shift is too small!')

    paths =GeneratePathsGBM(NoOfPaths,NoOfSteps,T,r,sigma,S_0shift)
    Sshifted = paths["S"] - shift
    time = paths["time"]
    return {"time":time,"S":Sshifted}

def GeneratePathsGBM(NoOfPaths,NoOfSteps,T,r,sigma,S_0):
    Z = np.random.normal(0.0,1.0,[NoOfPaths,NoOfSteps])
    X = np.zeros([NoOfPaths, NoOfSteps+1])
    W = np.zeros([NoOfPaths, NoOfSteps+1])
    time = np.zeros([NoOfSteps+1])

    X[:,0] = np.log(S_0)

    dt = T / float(NoOfSteps)
    for i in range(0,NoOfSteps):
        # making sure that samples from normal have mean 0 and variance 1
        if NoOfPaths > 1:
            Z[:,i] = (Z[:,i] - np.mean(Z[:,i])) / np.std(Z[:,i])
        W[:,i+1] = W[:,i] + np.power(dt, 0.5)*Z[:,i]
        X[:,i+1] = X[:,i] + (r - 0.5 * sigma * sigma) * dt + sigma * (W[:,i+1]-W[:,i])
        time[i+1] = time[i] +dt

    #Compute exponent of ABM
    S = np.exp(X)
    paths = {"time":time,"S":S}
    return paths

# Shifted Black-Scholes Call option price
def BS_Call_Put_Option_Price_Shifted(CP,S_0,K,sigma,tau,r,shift):
    K_new = K + shift
    S_0_new = S_0 + shift
    return BS_Call_Put_Option_Price(CP,S_0_new,K_new,sigma,tau,r)

# Black-Scholes Call option price
def BS_Call_Put_Option_Price(CP,S_0,K,sigma,tau,r):
    if K is list:
        K = np.array(K).reshape([len(K),1])
    d1    = (np.log(S_0 / K) + (r + 0.5 * np.power(sigma,2.0)) * tau) / (sigma * np.sqrt(tau))
    d2    = d1 - sigma * np.sqrt(tau)
    if CP == OptionType.CALL:
        value = st.norm.cdf(d1) * S_0 - st.norm.cdf(d2) * K * np.exp(-r * tau)
    elif CP == OptionType.PUT:
        value = st.norm.cdf(-d2) * K * np.exp(-r * tau) - st.norm.cdf(-d1)*S_0
    return value

# Implied volatility method
def ImpliedVolatilityBlack76(CP,marketPrice,K,T,S_0):
    # To determine initial volatility we interpolate define a grid for sigma
    # and interpolate on the inverse
    sigmaGrid = np.linspace(0.0,2.0,5000)
    optPriceGrid = BS_Call_Put_Option_Price(CP,S_0,K,sigmaGrid,T,0.0)
    sigmaInitial = np.interp(marketPrice,optPriceGrid,sigmaGrid)
    print("Initial volatility = {0}".format(sigmaInitial))

    # Use determined input for the local-search (final tuning)
    func = lambda sigma: np.power(BS_Call_Put_Option_Price(CP,S_0,K,sigma,T,0.0) - marketPrice, 1.0)
    impliedVol = optimize.newton(func, sigmaInitial, tol=1e-15)
    print("Final volatility = {0}".format(impliedVol))
    return impliedVol

# Implied volatility method
def ImpliedVolatilityBlack76Shifted(CP,marketPrice,K,T,S_0,shift):
    # To determine initial volatility we interpolate define a grid for sigma
    # and interpolate on the inverse
    sigmaGrid = np.linspace(0.0,2.0,5000)
    optPriceGrid = BS_Call_Put_Option_Price_Shifted(CP,S_0,K,sigmaGrid,T,0.0,shift)
    sigmaInitial = np.interp(marketPrice,optPriceGrid,sigmaGrid)
    print("Initial volatility = {0}".format(sigmaInitial))

    # Use determined input for the local-search (final tuning)
    func = lambda sigma: np.power(BS_Call_Put_Option_Price_Shifted(CP,S_0,K,sigma,T,0.0,shift) - marketPrice, 1.0)
    impliedVol = optimize.newton(func, sigmaInitial, tol=1e-15)
    print("Final volatility = {0}".format(impliedVol))
    return impliedVol

def mainCalculation():
    NoOfPaths = 10000
    NoOfSteps = 500
    T         = 3.0
    sigma     = 0.2
    L0        = -0.05
    shift     = 0.1
    K         = [0.95]
    CP        = OptionType.CALL

    P0T =lambda T: np.exp(-0.1*T)

    np.random.seed(4)
    Paths = GeneratePathsGBMShifted(NoOfPaths,NoOfSteps,T,0.0,sigma,L0,shift)
    time  = Paths["time"]
    L     = Paths["S"]

    print(np.mean(L[:,-1]))

    # Plot first few paths
    plt.figure(1)
    plt.plot(time,np.transpose(L[0:20,:]))
    plt.grid()

    # Shifted lognormal for different shift parameters
    plt.figure(2)
    shiftV = [1.0,2.0,3.0,4.0,5.0]
    legend = []
    for shiftTemp in shiftV:
        x = np.linspace(-shiftTemp,10,1000)
        lognnormPDF = lambda x,t :  st.lognorm.pdf(x+shiftTemp, scale = np.exp(np.log(L0+shiftTemp) + (- 0.5 * sigma * sigma)*t), s= np.sqrt(t) * sigma)
        pdf_x= lognnormPDF(x,T)
        plt.plot(x,pdf_x)
        legend.append('shift={0}'.format(shiftTemp))
    plt.legend(legend)
    plt.xlabel('x')
    plt.ylabel('pdf')
    plt.title('shifted lognormal density')
    plt.grid()

    # Call/Put option prices, MC vs. Analytical
    plt.figure(3)
    K = np.linspace(-shift,np.abs(L0)*3,25)
    optPriceMCV=np.zeros([len(K),1])
    for idx in range(0,len(K)):
        optPriceMCV[idx] =0.0
        if CP == OptionType.CALL:
            optPriceMCV[idx] = P0T(T)*np.mean(np.maximum(L[:,-1]-K[idx],0.0))
        elif CP == OptionType.PUT:
            optPriceMCV[idx] = P0T(T)*np.mean(np.maximum(K[idx]-L[:,-1],0.0))

    optPriceExact = P0T(T)*BS_Call_Put_Option_Price_Shifted(CP,L0,K,sigma,T,0.0,shift)
    plt.plot(K,optPriceMCV)
    plt.plot(K,optPriceExact,'--r')
    plt.grid()
    plt.xlabel('strike,K')
    plt.ylabel('option price')
    plt.legend(['Monte Carlo','Exact'])

    # Shift Effect on Option prices
    plt.figure(4)
    legend = []
    for shiftTemp in [0.2,0.3,0.4,0.5]:
        K = np.linspace(-shiftTemp,np.abs(L0)*6.0,25)
        optPriceExact = P0T(T)*BS_Call_Put_Option_Price_Shifted(CP,L0,K,sigma,T,0.0,shiftTemp)
        plt.plot(K,optPriceExact)
        legend.append('shift={0}'.format(shiftTemp))
    plt.grid()
    plt.xlabel('strike,K')
    plt.ylabel('option price')
    plt.legend(legend)

mainCalculation()

max_stars_repo_path: Lecture 07-Pricing of Swaptions and Negative Interest Rates/Materials/ShiftedLognormal.py

reason: import numpy,import scipy