text:
import numpy as np
import torch
import os
from collections import OrderedDict
from torch.autograd import Variable
import itertools
import util.util as util
from util.image_pool import ImagePool
from .base_model import BaseModel
from . import networks
import sys

from yolo.utils.datasets import pad_to_square, resize
import torchvision.transforms as transforms
from torchvision.utils import save_image
import torch.nn.functional as F
from data.crop_bboxes import bbox_rescale, resize_img, crop_object_by_bbox, save_traffic_cone, get_label, img_path2label_path

def mse_loss(input, target):
    return torch.sum((input - target)**2) / input.data.nelement()

def yolo_data_transform(img, img_size=416, normalized_labels=True):
        # ---------
        #  Image
        # ---------

        # img_path = self.img_files[index % len(self.img_files)].rstrip()

        # # Extract image as PyTorch tensor
        # img = transforms.ToTensor()(Image.open(img_path).convert('RGB'))

        # Handle images with less than three channels
        # if len(img.shape) != 3:
        #     img = img.unsqueeze(0)
        #     img = img.expand((3, img.shape[1:]))
        #img = transforms.ToTensor()(img)
        #print('img.shape = ',img.shape)
        _, _, h, w = img.shape
        h_factor, w_factor = (h, w) if normalized_labels else (1, 1)
        # Pad to square resolution
        img, pad = pad_to_square(img, 0)

        return resize(img, img_size)

class CycleGANModel(BaseModel):
    def name(self):
        return 'CycleGANModel'

    def initialize(self, opt):
        BaseModel.initialize(self, opt)

        nb = opt.batchSize
        size = opt.fineSize
        self.input_A = self.Tensor(nb, opt.input_nc, size, size)
        self.input_B = self.Tensor(nb, opt.output_nc, size, size)

        # load/define networks
        # The naming conversion is different from those used in the paper
        # Code (paper): G_A (G), G_B (F), D_A (D_Y), D_B (D_X)

        self.netG_A = networks.define_G(opt.input_nc, opt.output_nc,
                                        opt.ngf, opt.which_model_netG, opt.norm, not opt.no_dropout, self.gpu_ids)
        self.netG_B = networks.define_G(opt.output_nc, opt.input_nc,
                                        opt.ngf, opt.which_model_netG, opt.norm, not opt.no_dropout, self.gpu_ids)

        if self.isTrain:
            use_sigmoid = opt.no_lsgan
            self.netD_A = networks.define_D(opt.output_nc, opt.ndf,
                                            opt.which_model_netD,
                                            opt.n_layers_D, opt.norm, use_sigmoid, self.gpu_ids)
            self.netD_B = networks.define_D(opt.input_nc, opt.ndf,
                                            opt.which_model_netD,
                                            opt.n_layers_D, opt.norm, use_sigmoid, self.gpu_ids)
            self.netFeat = networks.define_feature_network(opt.which_model_feat, self.gpu_ids)

            # load pretrained YOLO3
            self.YOLO3 = networks.define_pretrained_yolo3(opt.yolo_weight, opt.model_def, self.gpu_ids)

        if not self.isTrain or opt.continue_train:
            which_epoch = opt.which_epoch
            self.load_network(self.netG_A, 'G_A', which_epoch)
            self.load_network(self.netG_B, 'G_B', which_epoch)
            if self.isTrain:
                self.load_network(self.netD_A, 'D_A', which_epoch)
                self.load_network(self.netD_B, 'D_B', which_epoch)

        if self.isTrain:
            self.old_lr = opt.lr
            self.fake_A_pool = ImagePool(opt.pool_size)
            self.fake_B_pool = ImagePool(opt.pool_size)
            # define loss functions
            self.criterionGAN = networks.GANLoss(use_lsgan=not opt.no_lsgan, tensor=self.Tensor)
            self.criterionCycle = torch.nn.L1Loss()
            self.criterionIdt = torch.nn.L1Loss()
            self.criterionFeat = mse_loss
            # yolo data transform
            self.yolo_data_transform = yolo_data_transform

            # utils
            #self.crop_img_by_bbox = crop_img_by_bbox

            # initialize optimizers
            self.optimizer_G = torch.optim.Adam(itertools.chain(self.netG_A.parameters(), self.netG_B.parameters()),
                                                lr=opt.lr, betas=(opt.beta1, 0.999))
            self.optimizer_D_A = torch.optim.Adam(self.netD_A.parameters(), lr=opt.lr, betas=(opt.beta1, 0.999))
            self.optimizer_D_B = torch.optim.Adam(self.netD_B.parameters(), lr=opt.lr, betas=(opt.beta1, 0.999))

        print('---------- Networks initialized -------------')
        networks.print_network(self.netG_A)
        networks.print_network(self.netG_B)
        if self.isTrain:
            networks.print_network(self.netD_A)
            networks.print_network(self.netD_B)
        print('-----------------------------------------------')

    def set_input(self, input):
        '''
        see unaligned_dataset.py

        return {'A': A_img, 'B': B_img,
                'A_paths': A_path, 'B_paths': B_path,
                'A_targets': A_targets}
        '''
        AtoB = self.opt.which_direction == 'AtoB'
        input_A = input['A' if AtoB else 'B'] # see unaligned_dataset.py
        input_B = input['B' if AtoB else 'A']
        self.input_A.resize_(input_A.size()).copy_(input_A)
        self.input_B.resize_(input_B.size()).copy_(input_B)
        self.image_paths = input['A_paths' if AtoB else 'B_paths']
        #print('\nA path: ', self.image_paths)
        self.A_targets = input['A_targets'] # prepared yolo taragets of A (synthesis traffic cones)
        #self.B_targets = input['B_targets'] # get bboxes of B_img to crop traffic cones
        #self.targets = input['target_batch'] # see unaligned_dataset.py & dataloader YOLOv3

        # For bbox crop
        self.A_path = input['A_paths']
        self.B_path = input['B_paths']

    def forward(self):
        self.real_A = Variable(self.input_A)
        self.real_B = Variable(self.input_B)

    def test(self):
        self.real_A = Variable(self.input_A, volatile=True)
        self.fake_B = self.netG_A.forward(self.real_A) # use it
        self.rec_A = self.netG_B.forward(self.fake_B)

        self.real_B = Variable(self.input_B, volatile=True)
        self.fake_A = self.netG_B.forward(self.real_B)
        self.rec_B = self.netG_A.forward(self.fake_A)

    # get image paths
    def get_image_paths(self):
        return self.image_paths

    def backward_D_basic(self, netD, real, fake):
        # Real
        pred_real = netD.forward(real)
        loss_D_real = self.criterionGAN(pred_real, True)
        # Fake
        pred_fake = netD.forward(fake.detach())
        loss_D_fake = self.criterionGAN(pred_fake, False)
        # Combined loss
        loss_D = (loss_D_real + loss_D_fake) * 0.5
        # backward
        loss_D.backward()
        return loss_D

    def backward_D_A(self):
        fake_B = self.fake_B_pool.query(self.fake_B)
        self.loss_D_A = self.backward_D_basic(self.netD_A, self.real_B, fake_B)

    def backward_D_B(self):
        fake_A = self.fake_A_pool.query(self.fake_A)
        self.loss_D_B = self.backward_D_basic(self.netD_B, self.real_A, fake_A)

    def backward_G(self):
        lambda_idt = self.opt.identity
        lambda_A = self.opt.lambda_A
        lambda_B = self.opt.lambda_B

        lambda_feat_AfB = self.opt.lambda_feat_AfB
        lambda_feat_BfA = self.opt.lambda_feat_BfA

        lambda_feat_fArecB = self.opt.lambda_feat_fArecB
        lambda_feat_fBrecA = self.opt.lambda_feat_fBrecA

        lambda_feat_ArecA = self.opt.lambda_feat_ArecA
        lambda_feat_BrecB = self.opt.lambda_feat_BrecB

        # Identity loss
        if lambda_idt > 0:
            # G_A should be identity if real_B is fed.
            self.idt_A = self.netG_A.forward(self.real_B)
            self.loss_idt_A = self.criterionIdt(self.idt_A, self.real_B) * lambda_B * lambda_idt
            # G_B should be identity if real_A is fed.
            self.idt_B = self.netG_B.forward(self.real_A)
            self.loss_idt_B = self.criterionIdt(self.idt_B, self.real_A) * lambda_A * lambda_idt
        else:
            self.loss_idt_A = 0
            self.loss_idt_B = 0

        # GAN loss
        # D_A(G_A(A))
        self.fake_B = self.netG_A.forward(self.real_A) # fake_B
        pred_fake = self.netD_A.forward(self.fake_B)
        self.loss_G_A = self.criterionGAN(pred_fake, True)
        # D_B(G_B(B))
        self.fake_A = self.netG_B.forward(self.real_B)
        pred_fake = self.netD_B.forward(self.fake_A)
        self.loss_G_B = self.criterionGAN(pred_fake, True)
        # Forward cycle loss
        self.rec_A = self.netG_B.forward(self.fake_B)

        # gamma = 1.
        # l_rec_A =  .2126 * self.rec_A[:,0]**gamma + .7152 * self.rec_A[:,1]**gamma + .0722 * self.rec_A[:,2]**gamma
        # l_real_A =  .2126 * self.real_A[:,0]**gamma + .7152 * self.real_A[:,1]**gamma + .0722 * self.real_A[:,2]**gamma

        self.loss_cycle_A = self.criterionCycle(self.rec_A, self.real_A) * lambda_A

        # Backward cycle loss
        self.rec_B = self.netG_A.forward(self.fake_A)

        # gamma = 1.
        # l_rec_B =  .2126 * self.rec_B[:,0]**gamma + .7152 * self.rec_B[:,1]**gamma + .0722 * self.rec_B[:,2]**gamma
        # l_real_B =  .2126 * self.real_B[:,0]**gamma + .7152 * self.real_B[:,1]**gamma + .0722 * self.real_B[:,2]**gamma

        self.loss_cycle_B = self.criterionCycle(self.rec_B, self.real_B) * lambda_B

        # print ('self.netFeat(self.real_A).parameters()', self.netFeat(self.real_A).parameters())
        # print ('self.netFeat(self.fake_B).parameters()', self.netFeat(self.fake_B).parameters())
        # print ('self.criterionFeat(self.netFeat(self.real_A), self.netFeat(self.fake_B)).parameters()', self.criterionFeat(self.netFeat(self.real_A), self.netFeat(self.fake_B)).parameters())

        # Perceptual Loss on Image Level
        #print('\nself.fake_B info: ', self.fake_B.shape)

        self.feat_loss_AfB = self.criterionFeat(self.netFeat(self.real_A), self.netFeat(self.fake_B)) * lambda_feat_AfB
        self.feat_loss_BfA = self.criterionFeat(self.netFeat(self.real_B), self.netFeat(self.fake_A)) * lambda_feat_BfA

        self.feat_loss_fArecB = self.criterionFeat(self.netFeat(self.fake_A), self.netFeat(self.rec_B)) * lambda_feat_fArecB
        self.feat_loss_fBrecA = self.criterionFeat(self.netFeat(self.fake_B), self.netFeat(self.rec_A)) * lambda_feat_fBrecA

        self.feat_loss_ArecA = self.criterionFeat(self.netFeat(self.real_A), self.netFeat(self.rec_A)) * lambda_feat_ArecA
        self.feat_loss_BrecB = self.criterionFeat(self.netFeat(self.real_B), self.netFeat(self.rec_B)) * lambda_feat_BrecB

        self.feat_loss = self.feat_loss_AfB + self.feat_loss_BfA + self.feat_loss_fArecB + self.feat_loss_fBrecA + self.feat_loss_ArecA + self.feat_loss_BrecB

        ###### yolov3 loss (regression + classification) #####
        # print('\nself.fake_B = ', self.fake_B.shape)
        self.fake_B_transformed = self.yolo_data_transform(self.fake_B)

        Tensor = torch.cuda.FloatTensor if torch.cuda.is_available() else torch.FloatTensor
        device = torch.device("cuda" if torch.cuda.is_available() else "cpu")

        self.fake_B_transformed = Variable(self.fake_B_transformed.type(Tensor), requires_grad=False)
        self.A_targets = Variable(self.A_targets.to(device), requires_grad=False)

        # print('\nfake_B_transformed = ', self.fake_B_transformed.shape)
        # print('\nself.A_targets shape = ', self.A_targets[0].shape)

        # print('\nself.A_targets = ', self.A_targets[0])

        ###### Perceptual Loss on Object Level #####

        # Crop traffic cones from images using bboxes
        #print('self.fake_B_transformed: ', self.fake_B_transformed.shape)
        # self.list_cones_syn = self.crop_img_by_bbox(self.fake_B_transformed, self.A_targets)
        # self.list_cones_real = self.crop_img_by_bbox(self.real_B, self.B_targets)
        #print('A_path: ', self.A_path[0].rstrip())
        #with torch.no_grad():
        self.A_img_path = self.A_path[0].rstrip()
        self.A_label_path = img_path2label_path(self.A_img_path, synthetic=True)
        #print('\nself.A_label_path: ', self.A_label_path)
        self.A_boxes = get_label(self.A_label_path)
        self.A_boxes_rescaled = bbox_rescale(self.A_boxes, 416)
        self.fake_B_resized = resize(self.fake_B, 416)
        #print('\nself.fake_B_resized: ', self.fake_B_resized.shape)
        self.fake_B_cones = crop_object_by_bbox(self.fake_B_resized, self.A_boxes_rescaled)

        self.B_img_path = self.B_path[0].rstrip()
        self.B_label_path = img_path2label_path(self.B_img_path, synthetic=False)
        #print('\nself.B_label_path: ', self.B_label_path)
        self.B_boxes = get_label(self.B_label_path)
        self.B_boxes_rescaled = bbox_rescale(self.B_boxes, 416)
        self.B_resized = resize(self.input_B, 416)
        #print('\nself.B_resized: ', self.B_resized.shape)
        self.real_B_cones = crop_object_by_bbox(self.B_resized, self.B_boxes_rescaled)

        # Extract traffic cones features

        # self.feat_cone_syns = torch.stack([self.netFeat(cone_syn) for cone_syn in self.list_cones_syn], dim=0)
        # self.feat_cone_syn_1d = self.feat_cone_syns.mean(dim=0)

        # self.feat_cone_reals = torch.stack([self.netFeat(cone_real) for cone_real in self.list_cones_real], dim=0)
        # self.feat_cone_real_1d = self.feat_cone_reals.mean(dim=0)

        fake_B_cones_feat_ = []
        for cone_syn in self.fake_B_cones:
            #print('\ncone fake info: ', cone_syn.shape)
            cone_syn_resized = resize(cone_syn, 360)
            #print('\ncone_syn_reseized fake B: ', cone_syn_resized.shape)
            fake_B_cones_feat_.append(self.netFeat(cone_syn_resized))
        self.fake_B_cones_features = torch.stack(fake_B_cones_feat_, dim=0)
        self.fake_B_cones_feature_mean = self.fake_B_cones_features.mean(dim=0)
        #print('\nself.fake_B_cones_feature_mean: ', self.fake_B_cones_feature_mean.shape)

        real_B_cones_feat_ = []
        # cropped_cone_path = '/home/basic/cyclegan/data/tmp_cones/'
        # i = 0
        for cone_real in self.real_B_cones:
            #i += 1
            #print('\ncone_real B: ', cone_real.shape)
            cone_real_resized = resize(cone_real, 360)
            # save_traffic_cone(cone_real_resized[0], cropped_cone_path, i)
            # print('\ncone_real_resized: ', cone_syn_resized[0].shape)
            real_B_cones_feat_.append(self.netFeat(cone_real_resized))
            #self.real_B_cones_features = torch.stack([self.netFeat(cone_real_resized)], dim=0)
        self.real_B_cones_features = torch.stack(real_B_cones_feat_, dim=0)
        self.real_B_cones_features_mean = self.real_B_cones_features.mean(dim=0)

        # Perceptual Loss at Object Level
        # self.feat_loss_object_level = self.criterionFeat(self.feat_cone_syn_1d, self.feat_cone_real_1d)
        self.feat_loss_object_level = self.criterionFeat(self.fake_B_cones_feature_mean, self.real_B_cones_features_mean)

        #with torch.no_grad():
            #self.yolo_loss, _ = self.YOLO3(self.fake_B_transformed, self.A_targets[0]) # create self.targets, this only worked for batchSize = 1
        self.yolo_loss, _ = self.YOLO3(self.fake_B_transformed, self.A_targets) # create self.targets

        # combined loss
        self.loss_G = self.loss_G_A + self.loss_G_B + self.loss_cycle_A \
                    + self.loss_cycle_B + self.loss_idt_A + self.loss_idt_B \
                    + self.feat_loss \
                    + self.yolo_loss * 10 \
                    + self.feat_loss_object_level * 20

        self.loss_G.backward()

    def optimize_parameters(self):
        # forward
        self.forward()
        # G_A and G_B
        self.optimizer_G.zero_grad()
        self.backward_G()
        self.optimizer_G.step()
        # D_A
        self.optimizer_D_A.zero_grad()
        self.backward_D_A()
        self.optimizer_D_A.step()
        # D_B
        self.optimizer_D_B.zero_grad()
        self.backward_D_B()
        self.optimizer_D_B.step()

    def get_current_errors(self):
        D_A = self.loss_D_A.data#[0]
        G_A = self.loss_G_A.data#[0]
        Cyc_A = self.loss_cycle_A.data#[0]
        D_B = self.loss_D_B.data#[0]
        G_B = self.loss_G_B.data#[0]
        Cyc_B = self.loss_cycle_B.data#[0]
        yolo_loss = self.yolo_loss.data
        loss_G = self.loss_G.data
        if self.opt.identity > 0.0:
            idt_A = self.loss_idt_A.data#[0]
            idt_B = self.loss_idt_B.data#[0]
            return OrderedDict([('D_A', D_A), ('G_A', G_A), ('Cyc_A', Cyc_A), ('idt_A', idt_A),
                                ('D_B', D_B), ('G_B', G_B), ('Cyc_B', Cyc_B), ('idt_B', idt_B), ('yolo_loss', yolo_loss), ('loss_G', loss_G)])
        else:
            return OrderedDict([('D_A', D_A), ('G_A', G_A), ('Cyc_A', Cyc_A),
                                ('D_B', D_B), ('G_B', G_B), ('Cyc_B', Cyc_B), ('yolo_loss', yolo_loss), ('loss_G', loss_G)])

    def get_current_visuals(self):
        real_A = util.tensor2im(self.real_A.data)
        fake_B = util.tensor2im(self.fake_B.data)
        rec_A = util.tensor2im(self.rec_A.data)
        real_B = util.tensor2im(self.real_B.data)
        fake_A = util.tensor2im(self.fake_A.data)
        rec_B = util.tensor2im(self.rec_B.data)
        if self.opt.identity > 0.0:
            idt_A = util.tensor2im(self.idt_A.data)
            idt_B = util.tensor2im(self.idt_B.data)
            return OrderedDict([('real_A', real_A), ('fake_B', fake_B), ('rec_A', rec_A), ('idt_B', idt_B),
                                ('real_B', real_B), ('fake_A', fake_A), ('rec_B', rec_B), ('idt_A', idt_A)])
        else:
            return OrderedDict([('real_A', real_A), ('fake_B', fake_B), ('rec_A', rec_A),
                                ('real_B', real_B), ('fake_A', fake_A), ('rec_B', rec_B)])

    def save(self, label):
        self.save_network(self.netG_A, 'G_A', label, self.gpu_ids)
        self.save_network(self.netD_A, 'D_A', label, self.gpu_ids)
        self.save_network(self.netG_B, 'G_B', label, self.gpu_ids)
        self.save_network(self.netD_B, 'D_B', label, self.gpu_ids)

    def update_learning_rate(self):
        lrd = self.opt.lr / self.opt.niter_decay
        lr = self.old_lr - lrd
        for param_group in self.optimizer_D_A.param_groups:
            param_group['lr'] = lr
        for param_group in self.optimizer_D_B.param_groups:
            param_group['lr'] = lr
        for param_group in self.optimizer_G.param_groups:
            param_group['lr'] = lr

        print('update learning rate: %f -> %f' % (self.old_lr, lr))
        self.old_lr = lr

max_stars_repo_name: basicskywards/cyclegan-yolo

max_issues_repo_name: basicskywards/cyclegan-yolo

max_forks_repo_name: basicskywards/cyclegan-yolo


text:
import open3d as o3d
import numpy as np
import python.open3d_tutorial as o3dtut
import copy

# http://www.open3d.org/docs/release/tutorial/geometry/mesh.html

########################################################################################################################
# 1. Mesh
########################################################################################################################
'''
Open3D has a data structure for 3D triangle meshes called TriangleMesh.
The code below shows how to read a triangle mesh from a ply file and print its vertices and triangles.

The TriangleMesh class has a few data fields such as vertices and triangles.
Open3D provides direct memory access to these fields via numpy.
'''
print("Testing mesh in Open3D...")
mesh = o3dtut.get_knot_mesh()
print(mesh)
print('Vertices:')
print(np.asarray(mesh.vertices))
print('Triangles:')
print(np.asarray(mesh.triangles))

########################################################################################################################
# 2. Visualize a 3D mesh
########################################################################################################################
print("Try to render a mesh with normals (exist: " +
      str(mesh.has_vertex_normals()) + ") and colors (exist: " +
      str(mesh.has_vertex_colors()) + ")")
o3d.visualization.draw_geometries([mesh],
                                  width=800, height=600,
                                  window_name='Visualize a 3D mesh')
print("A mesh with no normals and no colors does not look good.")
'''
You can rotate and move the mesh but it is painted with uniform gray color and does not look “3d”.
The reason is that the current mesh does not have normals for vertices or faces.
So uniform color shading is used instead of a more sophisticated Phong shading.
'''

########################################################################################################################
# 3. Surface normal estimation
########################################################################################################################
print("Computing normal and rendering it.")
mesh.compute_vertex_normals()
print(np.asarray(mesh.triangle_normals))
o3d.visualization.draw_geometries([mesh],
                                  width=800, height=600,
                                  window_name='Surface normal estimation')

########################################################################################################################
# 4. Crop mesh
########################################################################################################################
'''
We remove half of the surface by directly operating on the triangle and triangle_normals data fields of the mesh.
This is done via numpy.
'''
print("We make a partial mesh of only the first half triangles.")
mesh1 = copy.deepcopy(mesh)
mesh1.triangles = o3d.utility.Vector3iVector(np.asarray(mesh1.triangles)[:len(mesh1.triangles) // 2, :])
mesh1.triangle_normals = o3d.utility.Vector3dVector(
    np.asarray(mesh1.triangle_normals)[:len(mesh1.triangle_normals) // 2, :])
print(mesh1.triangles)
o3d.visualization.draw_geometries([mesh1],
                                  width=800, height=600,
                                  window_name='Crop mesh')

########################################################################################################################
# 5. Paint mesh
########################################################################################################################
# paint_uniform_color paints the mesh with a uniform color. The color is in RGB space, [0, 1] range.
print("Painting the mesh")
mesh1.paint_uniform_color([1, 0.706, 0])
o3d.visualization.draw_geometries([mesh1],
                                  width=800, height=600,
                                  window_name='Paint mesh')

########################################################################################################################
# 6. Mesh properties
########################################################################################################################
'''
https://zh.wikipedia.org/wiki/%E6%B5%81%E5%BD%A2
流形（英語：Manifolds）是可以局部欧几里得空间化的一个拓扑空间，是欧几里得空间中的曲线、曲面等概念的推广。欧几里得空间就是最简单的流
形的实例。

地球表面这样的球面则是一个稍微复杂的例子。一般的流形可以通过把许多平直的片折弯并粘连而成。

流形在数学中用于描述几何形体，它们为研究形体的可微性提供了一个自然的平台。物理上，经典力学的相空间和构造广义相对论的时空模型的四维伪
黎曼流形都是流形的实例。位形空间中也可以定义流形。环面就是双摆的位形空间。

一般可以把几何形体的拓扑结构看作是完全“柔软”的，因为所有变形（同胚）会保持拓扑结构不变；而把解析几何结构看作是“硬”的，因为整体的结构
都是固定的。例如，当一个多项式在 {displaystyle (0,1)}(0,1) 区间的取值确定了，则其在整个实数范围的值都被固定，可见局部的变动会导致
全局的变化。光滑流形可以看作是介于两者之间的模型：其无穷小的结构是“硬”的，而整体结构则是“柔软”的。这也许是中文译名“流形”的原因（整体
的形态可以流动）。该译名由著名数学家和数学教育学家江泽涵引入。这样，流形的硬度使它能够容纳微分结构，而它的软度使得它可以作为很多需要
独立的局部扰动的数学和物理的模型。

A triangle mesh has several properties that can be tested with Open3D. One important property is the manifold property,
where we can test the triangle mesh if it is edge manifold is_edge_manifold and if it is is_vertex_manifold. A triangle
mesh is edge manifold, if each edge is bounding either one or two triangles. The function is_edge_manifold has the bool
parameter allow_boundary_edges that defines if boundary edges should be allowed. Further, a triangle mesh is vertex
manifold if the star of the vertex is edge‐manifold and edge‐connected, e.g., two or more faces connected only by a
vertex and not by an edge.

Another property is the test of self-intersection.
The function is_self_intersecting returns True if there exists a triangle in the mesh that is intersecting another mesh.
A watertight mesh can be defined as a mesh that is edge manifold, vertex manifold and not self intersecting. The
function is_watertight implements this check in Open3D.

We also can test the triangle mesh, if it is orientable, i.e. the triangles can be oriented in such a way that all
normals point towards the outside. The corresponding function in Open3D is called is_orientable.

The code below tests a number of triangle meshes against those properties and visualizes the results. Non-manifold
edges are shown in red, boundary edges in green, non-manifold vertices are visualized as green points, and
self-intersecting triangles are shown in pink.
'''

def check_properties(name, mesh):
    mesh.compute_vertex_normals()

    edge_manifold = mesh.is_edge_manifold(allow_boundary_edges=True)
    edge_manifold_boundary = mesh.is_edge_manifold(allow_boundary_edges=False)
    vertex_manifold = mesh.is_vertex_manifold()
    self_intersecting = mesh.is_self_intersecting()
    watertight = mesh.is_watertight()
    orientable = mesh.is_orientable()

    print(name)
    print(f"  edge_manifold:          {edge_manifold}")
    print(f"  edge_manifold_boundary: {edge_manifold_boundary}")
    print(f"  vertex_manifold:        {vertex_manifold}")
    print(f"  self_intersecting:      {self_intersecting}")
    print(f"  watertight:             {watertight}")
    print(f"  orientable:             {orientable}")

    geoms = [mesh]
    if not edge_manifold:
        edges = mesh.get_non_manifold_edges(allow_boundary_edges=True)
        geoms.append(o3dtut.edges_to_lineset(mesh, edges, (1, 0, 0)))
    if not edge_manifold_boundary:
        edges = mesh.get_non_manifold_edges(allow_boundary_edges=False)
        geoms.append(o3dtut.edges_to_lineset(mesh, edges, (0, 1, 0)))
    if not vertex_manifold:
        verts = np.asarray(mesh.get_non_manifold_vertices())
        pcl = o3d.geometry.PointCloud(
            points=o3d.utility.Vector3dVector(np.asarray(mesh.vertices)[verts]))
        pcl.paint_uniform_color((0, 0, 1))
        geoms.append(pcl)
    if self_intersecting:
        intersecting_triangles = np.asarray(mesh.get_self_intersecting_triangles())
        intersecting_triangles = intersecting_triangles[0:1]
        intersecting_triangles = np.unique(intersecting_triangles)
        print("  # visualize self-intersecting triangles")
        triangles = np.asarray(mesh.triangles)[intersecting_triangles]
        edges = [
            np.vstack((triangles[:, i], triangles[:, j]))
            for i, j in [(0, 1), (1, 2), (2, 0)]
        ]
        edges = np.hstack(edges).T
        edges = o3d.utility.Vector2iVector(edges)
        geoms.append(o3dtut.edges_to_lineset(mesh, edges, (1, 0, 1)))
    o3d.visualization.draw_geometries(geoms, mesh_show_back_face=True,
                                      width=800, height=600,
                                      window_name='Mesh properties: ' + name)

check_properties('Knot', o3dtut.get_knot_mesh())
check_properties('Moebius', o3d.geometry.TriangleMesh.create_moebius(twists=1))
check_properties("non-manifold edge", o3dtut.get_non_manifold_edge_mesh())
check_properties("non-manifold vertex", o3dtut.get_non_manifold_vertex_mesh())
check_properties("open box", o3dtut.get_open_box_mesh())
check_properties("intersecting_boxes", o3dtut.get_intersecting_boxes_mesh())
########################################################################################################################
# 7. Mesh filtering - Average filter
########################################################################################################################
'''
The simplest filter is the average filter. A given vertex vi is given by the average of the adjacent vertices N
This filter can be used to denoise meshes as demonstrated in the code below. The parameter number_of_iterations in the
function filter_smooth_simple defines the how often the filter is applied to the mesh.
'''
print('create noisy mesh')
mesh_in = o3dtut.get_knot_mesh()
vertices = np.asarray(mesh_in.vertices)
noise = 5
vertices += np.random.uniform(0, noise, size=vertices.shape)
mesh_in.vertices = o3d.utility.Vector3dVector(vertices)
mesh_in.compute_vertex_normals()
o3d.visualization.draw_geometries([mesh_in],
                                  width=800, height=600,
                                  window_name='Average filter with noise')

print('filter with average with 1 iteration')
mesh_out = mesh_in.filter_smooth_simple(number_of_iterations=1)
mesh_out.compute_vertex_normals()
o3d.visualization.draw_geometries([mesh_out],
                                  width=800, height=600,
                                  window_name='Average filter with noise 1 iteration')

print('filter with average with 5 iterations')
mesh_out = mesh_in.filter_smooth_simple(number_of_iterations=5)
mesh_out.compute_vertex_normals()
o3d.visualization.draw_geometries([mesh_out],
                                  width=800, height=600,
                                  window_name='Average filter with noise 5 iteration')
########################################################################################################################
# 8. Mesh filtering - Laplacian
########################################################################################################################
'''
Another important mesh filter is the Laplacian defined as vi=vi⋅λ∑wnvn−vi,

where λ is the strength of the filter and wn are normalized weights that relate to the distance of the neighboring
vertices. The filter is implemented in filter_smooth_laplacian and has the parameters number_of_iterations and lambda.
'''
print('filter with Laplacian with 10 iterations')
mesh_out = mesh_in.filter_smooth_laplacian(number_of_iterations=10)
mesh_out.compute_vertex_normals()
o3d.visualization.draw_geometries([mesh_out],
                                  width=800, height=600,
                                  window_name='Laplacian with noise 10 iteration')

print('filter with Laplacian with 50 iterations')
mesh_out = mesh_in.filter_smooth_laplacian(number_of_iterations=50)
mesh_out.compute_vertex_normals()
o3d.visualization.draw_geometries([mesh_out],
                                  width=800, height=600,
                                  window_name='Laplacian with noise 50 iteration')

########################################################################################################################
# 9. Mesh filtering - Taubin filter
########################################################################################################################
'''
The problem with the average and Laplacian filter is that they lead to a shrinkage of the triangle mesh. [Taubin1995]
showed that the application of two Laplacian filters with different λ parameters can prevent the mesh shrinkage. The
filter is implemented in filter_smooth_taubin.
'''
print('filter with Taubin with 10 iterations')
mesh_out = mesh_in.filter_smooth_taubin(number_of_iterations=10)
mesh_out.compute_vertex_normals()
o3d.visualization.draw_geometries([mesh_out],
                                  width=800, height=600,
                                  window_name='Taubin with noise 10 iteration')

print('filter with Taubin with 100 iterations')
mesh_out = mesh_in.filter_smooth_taubin(number_of_iterations=100)
mesh_out.compute_vertex_normals()
o3d.visualization.draw_geometries([mesh_out],
                                  width=800, height=600,
                                  window_name='Taubin with noise 100 iteration')

########################################################################################################################
# 10. Sampling
########################################################################################################################
'''
Open3D includes functions to sample point clouds from a triangle mesh. The simplest method is sample_points_uniformly
that uniformly samples points from the 3D surface based on the triangle area. The parameter number_of_points defines
how many points are sampled from the triangle surface.
'''
mesh = o3d.geometry.TriangleMesh.create_sphere()
mesh.compute_vertex_normals()
o3d.visualization.draw_geometries([mesh],
                                  width=800, height=600,
                                  window_name='Mesh with Normal')

pcd = mesh.sample_points_uniformly(number_of_points=500)
o3d.visualization.draw_geometries([pcd],
                                  width=800, height=600,
                                  window_name='Meshed Point Cloud')

mesh = o3dtut.get_bunny_mesh()
mesh.compute_vertex_normals()
o3d.visualization.draw_geometries([mesh],
                                  width=800, height=600,
                                  window_name='Mesh with Normal')

pcd = mesh.sample_points_uniformly(number_of_points=500)
o3d.visualization.draw_geometries([pcd],
                                  width=800, height=600,
                                  window_name='Meshed Point Cloud')
'''
Uniform sampling can yield clusters of points on the surface, while a method called Poisson disk sampling can evenly
distribute the points on the surface. The method sample_points_poisson_disk implements sample elimination. It starts
with a sampled point cloud and removes points to satisfy the sampling criterion. The method supports two options to
provide the initial point cloud:
1. Default via the parameter init_factor: The method first samples uniformly a point cloud from the mesh with
   init_factor x number_of_points and uses this for the elimination.
2. One can provide a point cloud and pass it to the sample_points_poisson_disk method. Then, this point cloud is used
for elimination.
'''
mesh = o3d.geometry.TriangleMesh.create_sphere()
pcd = mesh.sample_points_poisson_disk(number_of_points=500, init_factor=5)
o3d.visualization.draw_geometries([pcd],
                                  width=800, height=600,
                                  window_name='Mesh with Normal')

pcd = mesh.sample_points_uniformly(number_of_points=2500)
pcd = mesh.sample_points_poisson_disk(number_of_points=500, pcl=pcd)
o3d.visualization.draw_geometries([pcd],
                                  width=800, height=600,
                                  window_name='Meshed Point Cloud')

mesh = o3dtut.get_bunny_mesh()
pcd = mesh.sample_points_poisson_disk(number_of_points=500, init_factor=5)
o3d.visualization.draw_geometries([pcd],
                                  width=800, height=600,
                                  window_name='Mesh with Normal')

pcd = mesh.sample_points_uniformly(number_of_points=2500)
pcd = mesh.sample_points_poisson_disk(number_of_points=500, pcl=pcd)
o3d.visualization.draw_geometries([pcd],
                                  width=800, height=600,
                                  window_name='Meshed Point Cloud')
########################################################################################################################
# 11. Mesh subdivision
########################################################################################################################
'''
In mesh subdivision, we divide each triangle into a number of smaller triangles. In the simplest case, we compute the
midpoint of each side per triangle and divide the triangle into four smaller triangles. This is implemented in the
subdivide_midpoint function. The 3D surface and area stays the same, but the number of vertices and triangles increases.
The parameter number_of_iterations defines how many times this process should be repeated.
'''
mesh = o3d.geometry.TriangleMesh.create_box()
mesh.compute_vertex_normals()
print(f'The mesh has {len(mesh.vertices)} vertices and {len(mesh.triangles)} triangles')

o3d.visualization.draw_geometries([mesh], zoom=0.8, mesh_show_wireframe=True,
                                  width=800, height=600,
                                  window_name='Raw Mesh')

mesh = mesh.subdivide_midpoint(number_of_iterations=1)
print(f'After subdivision it has {len(mesh.vertices)} vertices and {len(mesh.triangles)} triangles')
o3d.visualization.draw_geometries([mesh], zoom=0.8, mesh_show_wireframe=True,
                                  width=800, height=600,
                                  window_name='Mesh Subdivision')
'''
Open3D implements an additional subdivision method based on [Loop1987]. The method is based on a quartic box spline,
which generates C2 continuous limit surfaces everywhere except at extraordinary vertices where they are C1 continuous.
This leads to smoother corners.
'''
mesh = o3d.geometry.TriangleMesh.create_sphere()
mesh.compute_vertex_normals()
print(f'The mesh has {len(mesh.vertices)} vertices and {len(mesh.triangles)} triangles')
o3d.visualization.draw_geometries([mesh], zoom=0.8, mesh_show_wireframe=True,
                                  width=800, height=600,
                                  window_name='Raw Mesh')

mesh = mesh.subdivide_loop(number_of_iterations=2)
print(f'After subdivision it has {len(mesh.vertices)} vertices and {len(mesh.triangles)} triangles')
o3d.visualization.draw_geometries([mesh], zoom=0.8, mesh_show_wireframe=True,
                                  width=800, height=600,
                                  window_name='Mesh Subdivision')

mesh = o3dtut.get_knot_mesh()
mesh.compute_vertex_normals()
print(f'The mesh has {len(mesh.vertices)} vertices and {len(mesh.triangles)} triangles')
o3d.visualization.draw_geometries([mesh], zoom=0.8, mesh_show_wireframe=True,
                                  width=800, height=600,
                                  window_name='Raw Mesh')

mesh = mesh.subdivide_loop(number_of_iterations=1)
print(f'After subdivision it has {len(mesh.vertices)} vertices and {len(mesh.triangles)} triangles')
o3d.visualization.draw_geometries([mesh], zoom=0.8, mesh_show_wireframe=True,
                                  width=800, height=600,
                                  window_name='Mesh Subdivision')

########################################################################################################################
# 12. Mesh simplification
########################################################################################################################
'''
Sometimes we want to represent a high-resolution mesh with fewer triangles and vertices, but the low-resolution mesh
should still be close to the high-resolution mesh. For this purpose Open3D implements a number of mesh simplification
methods.
'''
########################################################################################################################
# 13. Mesh simplification - Vertex clustering
########################################################################################################################
'''
The vertex clustering method pools all vertices that fall into a voxel of a given size to a single vertex. The method is
implemented in simplify_vertex_clustering and has as parameters voxel_size that defines the size of the voxel grid and
contraction that defines how the vertices are pooled. o3d.geometry.SimplificationContraction.Average computes a simple
average.
'''
mesh_in = o3dtut.get_bunny_mesh()
print(f'Input mesh has {len(mesh_in.vertices)} vertices and {len(mesh_in.triangles)} triangles')
o3d.visualization.draw_geometries([mesh_in],
                                  width=800, height=600,
                                  window_name='Vertex clustering')

voxel_size = max(mesh_in.get_max_bound() - mesh_in.get_min_bound()) / 32
print(f'voxel_size = {voxel_size:e}')
mesh_smp = mesh_in.simplify_vertex_clustering(voxel_size=voxel_size,
                                              contraction=o3d.geometry.SimplificationContraction.Average)
print(f'Simplified mesh has {len(mesh_smp.vertices)} vertices and {len(mesh_smp.triangles)} triangles')
o3d.visualization.draw_geometries([mesh_smp],
                                  width=800, height=600,
                                  window_name='Vertex clustering: ' + str(voxel_size))

voxel_size = max(mesh_in.get_max_bound() - mesh_in.get_min_bound()) / 16
print(f'voxel_size = {voxel_size:e}')
mesh_smp = mesh_in.simplify_vertex_clustering(voxel_size=voxel_size,
                                              contraction=o3d.geometry.SimplificationContraction.Average)
print(f'Simplified mesh has {len(mesh_smp.vertices)} vertices and {len(mesh_smp.triangles)} triangles')
o3d.visualization.draw_geometries([mesh_smp],
                                  width=800, height=600,
                                  window_name='Vertex clustering: ' + str(voxel_size))
########################################################################################################################
# 14. Mesh decimation
########################################################################################################################
'''
Another category of mesh simplification methods is mesh decimation that operates in incremental steps. We select a
single triangle that minimizes an error metric and removes it. This is repeated until a required number of triangles is
achieved. Open3D implements simplify_quadric_decimation that minimizes error quadrics (distances to neighboring planes).
The parameter target_number_of_triangles defines the stopping critera of the decimation algorithm.
'''
mesh_smp = mesh_in.simplify_quadric_decimation(target_number_of_triangles=6500)
print(f'Simplified mesh has {len(mesh_smp.vertices)} vertices and {len(mesh_smp.triangles)} triangles')
o3d.visualization.draw_geometries([mesh_smp],
                                  width=800, height=600,
                                  window_name='Mesh decimation: ' + str(6500))

mesh_smp = mesh_in.simplify_quadric_decimation(target_number_of_triangles=1700)
print(f'Simplified mesh has {len(mesh_smp.vertices)} vertices and {len(mesh_smp.triangles)} triangles')
o3d.visualization.draw_geometries([mesh_smp],
                                  width=800, height=600,
                                  window_name='Mesh decimation: ' + str(1700))

########################################################################################################################
# 15. Connected components
########################################################################################################################
'''
The result of various reconstruction methods.
Open3D implements a connected components algorithm cluster_connected_triangles that assigns each triangle to a cluster
of connected triangles. It returns for each triangle the index of the cluster in triangle_clusters, and per cluster the
number of triangles in cluster_n_triangles and the surface area of the cluster in cluster_area.

This is useful in for instance RGBD Integration, which is not always a single triangle mesh, but a number of meshes.
Some of the smaller parts are due to noise and we most likely want to remove them.

The code below shows the application of cluster_connected_triangles and how it can be used to remove spurious triangles.
'''
print("Generate data")
mesh = o3dtut.get_bunny_mesh().subdivide_midpoint(number_of_iterations=2)
vert = np.asarray(mesh.vertices)
min_vert, max_vert = vert.min(axis=0), vert.max(axis=0)
for _ in range(30):
    cube = o3d.geometry.TriangleMesh.create_box()
    cube.scale(0.005, center=cube.get_center())
    cube.translate(
        (
            np.random.uniform(min_vert[0], max_vert[0]),
            np.random.uniform(min_vert[1], max_vert[1]),
            np.random.uniform(min_vert[2], max_vert[2]),
        ),
        relative=False,
    )
    mesh += cube
mesh.compute_vertex_normals()
print("Show input mesh")
o3d.visualization.draw_geometries([mesh],
                                  width=800, height=600,
                                  window_name='Input Mesh')

print("Cluster connected triangles")
with o3d.utility.VerbosityContextManager(o3d.utility.VerbosityLevel.Debug) as cm:
    triangle_clusters, cluster_n_triangles, cluster_area = (mesh.cluster_connected_triangles())
triangle_clusters = np.asarray(triangle_clusters)
cluster_n_triangles = np.asarray(cluster_n_triangles)
cluster_area = np.asarray(cluster_area)

print("Show mesh with small clusters removed")
mesh_0 = copy.deepcopy(mesh)
triangles_to_remove = cluster_n_triangles[triangle_clusters] < 100
mesh_0.remove_triangles_by_mask(triangles_to_remove)
o3d.visualization.draw_geometries([mesh_0],
                                  width=800, height=600,
                                  window_name='small clusters')

print("Show largest cluster")
mesh_1 = copy.deepcopy(mesh)
largest_cluster_idx = cluster_n_triangles.argmax()
triangles_to_remove = triangle_clusters != largest_cluster_idx
mesh_1.remove_triangles_by_mask(triangles_to_remove)
o3d.visualization.draw_geometries([mesh_1],
                                  width=800, height=600,
                                  window_name='largest cluster')

ext: py

reason: import numpy


text:
import math
import numbers
import operator
import unittest
import numpy as np
from functools import reduce

import datamol as dm

def random_fn(*args, op="mul", **kwargs):
    """Perform random functions on a list"""
    all_values = [x for x in args if isinstance(x, numbers.Number)]
    all_values += [x for x in kwargs.values() if isinstance(x, numbers.Number)]
    op_fn = getattr(operator, op, None)
    if op_fn is None:
        op_fn = getattr(math, op)
        return op_fn(all_values[0])
    return reduce(op_fn, all_values)

class TestJobs(unittest.TestCase):
    def test_sequential(self):
        jobrunner = dm.JobRunner(n_jobs=None, progress=False)
        # practically do nothing (add a single value with nothing)
        o1 = jobrunner(random_fn, [9, 25, 1024], op="add")
        self.assertEqual(o1, [9, 25, 1024])

        # take the sqrt
        o2 = jobrunner(random_fn, [9, 25, 1024], op="sqrt")
        self.assertEqual(o2, [3, 5, 32])

        # multiply all inputs
        o3 = jobrunner(random_fn, [(1, 2, 3), (4, 5, 6), (3, 4, 0)], arg_type="args", op="mul")
        self.assertEqual(o3, [6, 4 * 5 * 6, 0])

        # do the same thing but with kwargs
        o4 = jobrunner(
            random_fn,
            iter([dict(a=1, b=2, c=3), dict(a=4, b=5, c=6), dict(a=3, b=4, c=0)]),
            arg_type="kwargs",
            op="mul",
        )
        self.assertEqual(o4, [6, 4 * 5 * 6, 0])

        o5 = jobrunner(random_fn, np.asarray([9, 25, 1024]), op="add")
        self.assertEqual(o5, [9, 25, 1024])

    def test_parallel(self):
        jobrunner1 = dm.JobRunner(n_jobs=4, progress=True)  # use loky backend
        o1 = jobrunner1(random_fn, [9, 25, 1024], op="add")
        self.assertEqual(o1, [9, 25, 1024])

        o5 = jobrunner1(random_fn, np.asarray([9, 25, 1024]), op="add")
        self.assertEqual(o5, [9, 25, 1024])

        o3 = jobrunner1(random_fn, [(1, 2, 3), (4, 5, 6), (3, 4, 0)], arg_type="args", op="mul")
        self.assertEqual(o3, [6, 4 * 5 * 6, 0])

        # use threads instead, no progress
        jobrunner2 = dm.JobRunner(n_jobs=2, progress=False, prefer="threads")
        o2 = jobrunner2(random_fn, [9, 25, 1024], op="sqrt")
        self.assertEqual(o2, [3, 5, 32])

        o4 = jobrunner2(
            random_fn,
            iter([dict(a=1, b=2, c=3), dict(a=4, b=5, c=6), dict(a=3, b=4, c=0)]),
            arg_type="kwargs",
            op="mul",
        )
        self.assertEqual(o4, [6, 4 * 5 * 6, 0])

    def test_seq_vs_parallel(self):

        # test parallel vs sequential
        jobrunner = dm.JobRunner(n_jobs=4, progress=False)  # use loky backend
        o_seq = jobrunner.sequential(
            random_fn, [(1, 2, 3), (4, 5, 6), (3, 4, 0)], arg_type="args", op="mul"
        )
        o_par = jobrunner.parallel(
            random_fn, [(1, 2, 3), (4, 5, 6), (3, 4, 0)], arg_type="args", op="mul"
        )
        self.assertEqual(o_seq, o_par)

    def test_parallelized(self):
        def fn(x):
            return x ** 2

        results = dm.parallelized(
            fn,
            [{"x": i} for i in range(10)],
            scheduler="processes",
            n_jobs=None,
            arg_type="kwargs",
            progress=True,
        )
        assert results == [0, 1, 4, 9, 16, 25, 36, 49, 64, 81]

        results = dm.parallelized(
            fn,
            [[i] for i in range(10)],
            scheduler="processes",
            n_jobs=None,
            arg_type="args",
            progress=True,
        )
        assert results == [0, 1, 4, 9, 16, 25, 36, 49, 64, 81]

        results = dm.parallelized(
            fn,
            range(10),
            scheduler="processes",
            n_jobs=None,
            progress=False,
        )
        assert results == [0, 1, 4, 9, 16, 25, 36, 49, 64, 81]

    def test_job_kwargs(self):
        def fn(x):
            return x ** 2

        results = dm.parallelized(
            fn,
            [{"x": i} for i in range(10)],
            scheduler="processes",
            n_jobs=None,
            arg_type="kwargs",
            progress=True,
            verbose=100,
        )
        assert results == [0, 1, 4, 9, 16, 25, 36, 49, 64, 81]

    def test_tqdm_kwargs(self):
        def fn(x):
            return x ** 2

        results = dm.parallelized(
            fn,
            [{"x": i} for i in range(10)],
            scheduler="processes",
            n_jobs=None,
            arg_type="kwargs",
            progress=True,
            tqdm_kwargs=dict(desc="My progress bar"),
        )
        assert results == [0, 1, 4, 9, 16, 25, 36, 49, 64, 81]

    def test_with_batch_size(self):
        def _fn(n):
            return n * 3

        def _fn_return_none(n):
            return None

        results = dm.utils.parallelized(
            _fn,
            range(997),
            n_jobs=-1,
            progress=True,
            batch_size=10,
        )
        assert len(results) == 997

        results = dm.utils.parallelized(
            _fn_return_none,
            range(997),
            n_jobs=-1,
            progress=True,
            batch_size=10,
        )
        assert len(results) == 997


text:
"""
Created on Jul 3, 2019

@author: Hazem
"""

import numpy as np
import time
from numba import njit, guvectorize, float64
import scipy.optimize as opt
from matplotlib import pyplot as plt
import pdb
import csv

# Set
t = np.arange(1, 101)
NT = len(t)

# Parameters

fosslim = (
    6000  # Maximum cumulative extraction fossil fuels (GtC); denoted by CCum
)
tstep = 5  # Years per Period
ifopt = 0  # Indicator where optimized is 1 and base is 0

# Preferences

elasmu = 1.45  #  Elasticity of marginal utility of consumption
prstp = 0.015  #   Initial rate of social time preference per year

# ** Population and technology
gama = 0.300  #   Capital elasticity in production function         /.300 /
pop0 = 7403  # Initial world population 2015 (millions)          /7403 /
popadj = 0.134  #  Growth rate to calibrate to 2050 pop projection  /0.134/
popasym = 11500  # Asymptotic population (millions)                 /11500/
dk = 0.100  #     Depreciation rate on capital (per year)           /.100 /
q0 = 105.5  #     Initial world gross output 2015 (trill 2010 USD) /105.5/
k0 = 223  #     Initial capital value 2015 (trill 2010 USD)        /223  /
a0 = 5.115  #     Initial level of total factor productivity       /5.115/
ga0 = 0.076  #    Initial growth rate for TFP per 5 years          /0.076/
dela = 0.005  #   Decline rate of TFP per 5 years                  /0.005/

# ** Emissions parameters
gsigma1 = -0.0152  # Initial growth of sigma (per year)            /-0.0152/
dsig = -0.001  #   Decline rate of decarbonization (per period)    /-0.001 /
eland0 = 2.6  #  Carbon emissions from land 2015 (GtCO2 per year)   / 2.6   /
deland = 0.115  # Decline rate of land emissions (per period)        / .115  /
e0 = 35.85  #    Industrial emissions 2015 (GtCO2 per year)       /35.85  /
miu0 = 0.03  #   Initial emissions control rate for base case 2015  /.03    /

# ** Carbon cycle
# * Initial Conditions
mat0 = 851  #  Initial Concentration in atmosphere 2015 (GtC)       /851  /
mu0 = 460  #  Initial Concentration in upper strata 2015 (GtC)     /460  /
ml0 = 1740  #  Initial Concentration in lower strata 2015 (GtC)    /1740 /
mateq = 588  # mateq Equilibrium concentration atmosphere  (GtC)    /588  /
mueq = 360  # mueq Equilibrium concentration in upper strata (GtC) /360  /
mleq = 1720  # mleq Equilibrium concentration in lower strata (GtC) /1720 /

# * Flow paramaters, denoted by Phi_ij in the model
b12 = 0.12  #    Carbon cycle transition matrix                     /.12  /
b23 = 0.007  #   Carbon cycle transition matrix                    /0.007/
# * These are for declaration and are defined later
b11 = None  # Carbon cycle transition matrix
b21 = None  # Carbon cycle transition matrix
b22 = None  # Carbon cycle transition matrix
b32 = None  # Carbon cycle transition matrix
b33 = None  # Carbon cycle transition matrix
sig0 = None  # Carbon intensity 2010 (kgCO2 per output 2005 USD 2010)

# ** Climate model parameters
t2xco2 = 3.1  # Equilibrium temp impact (oC per doubling CO2)    / 3.1 /
fex0 = 0.5  #   2015 forcings of non-CO2 GHG (Wm-2)              / 0.5 /
fex1 = 1.0  #   2100 forcings of non-CO2 GHG (Wm-2)              / 1.0 /
tocean0 = 0.0068  # Initial lower stratum temp change (C from 1900) /.0068/
tatm0 = 0.85  #  Initial atmospheric temp change (C from 1900)    /0.85/
c1 = 0.1005  #     Climate equation coefficient for upper level  /0.1005/
c3 = 0.088  #     Transfer coefficient upper to lower stratum    /0.088/
c4 = 0.025  #     Transfer coefficient for lower level           /0.025/
fco22x = 3.6813  # eta in the model; Eq.22 : Forcings of equilibrium CO2 doubling (Wm-2)   /3.6813 /

# ** Climate damage parameters
a10 = 0  #     Initial damage intercept                         /0   /
a20 = None  #     Initial damage quadratic term
a1 = 0  #      Damage intercept                                 /0   /
a2 = 0.00236  #      Damage quadratic term                     /0.00236/
a3 = 2.00  #      Damage exponent                              /2.00   /

# ** Abatement cost
expcost2 = 2.6  # Theta2 in the model, Eq. 10 Exponent of control cost function             / 2.6  /
pback = 550  #   Cost of backstop 2010$ per tCO2 2015          / 550  /
gback = 0.025  #   Initial cost decline backstop cost per period / .025/
limmiu = 1.2  #  Upper limit on control rate after 2150        / 1.2 /
tnopol = 45  #  Period before which no emissions controls base  / 45   /
cprice0 = 2  # Initial base carbon price (2010$ per tCO2)      / 2    /
gcprice = 0.02  # Growth rate of base carbon price per year     /.02  /

# ** Scaling and inessential parameters
# * Note that these are unnecessary for the calculations
# * They ensure that MU of first period's consumption =1 and PV cons = PV utilty
scale1 = 0.0302455265681763  #    Multiplicative scaling coefficient           /0.0302455265681763 /
scale2 = -10993.704  #    Additive scaling coefficient       /-10993.704/;

# * Parameters for long-run consistency of carbon cycle
# (Question)
b11 = 1 - b12
b21 = b12 * mateq / mueq
b22 = 1 - b21 - b23
b32 = b23 * mueq / mleq
b33 = 1 - b32

# * Further definitions of parameters
a20 = a2
sig0 = e0 / (q0 * (1 - miu0))  # From Eq. 14
lam = fco22x / t2xco2  # From Eq. 25

l = np.zeros(NT)
l[0] = pop0  # Labor force
al = np.zeros(NT)
al[0] = a0
gsig = np.zeros(NT)
gsig[0] = gsigma1
sigma = np.zeros(NT)
sigma[0] = sig0
ga = ga0 * np.exp(-dela * 5 * (t - 1))  # TFP growth rate dynamics, Eq. 7
pbacktime = pback * (1 - gback) ** (t - 1)  # Backstop price
etree = eland0 * (1 - deland) ** (t - 1)  # Emissions from deforestration
rr = 1 / ((1 + prstp) ** (tstep * (t - 1)))  # Eq. 3
# The following three equations define the exogenous radiative forcing; used in Eq. 23
forcoth = np.full(NT, fex0)
forcoth[0:18] = forcoth[0:18] + (1 / 17) * (fex1 - fex0) * (t[0:18] - 1)
forcoth[18:NT] = forcoth[18:NT] + (fex1 - fex0)
optlrsav = (
    (dk + 0.004) / (dk + 0.004 * elasmu + prstp) * gama
)  # Optimal long-run savings rate used for transversality (Question)
cost1 = np.zeros(NT)
cumetree = np.zeros(NT)
cumetree[0] = 100
cpricebase = cprice0 * (1 + gcprice) ** (5 * (t - 1))

@njit("(float64[:], int32)")
def InitializeLabor(il, iNT):
    for i in range(1, iNT):
        il[i] = il[i - 1] * (popasym / il[i - 1]) ** popadj

@njit("(float64[:], int32)")
def InitializeTFP(ial, iNT):
    for i in range(1, iNT):
        ial[i] = ial[i - 1] / (1 - ga[i - 1])

@njit("(float64[:], int32)")
def InitializeGrowthSigma(igsig, iNT):
    for i in range(1, iNT):
        igsig[i] = igsig[i - 1] * ((1 + dsig) ** tstep)

@njit("(float64[:], float64[:],float64[:],int32)")
def InitializeSigma(isigma, igsig, icost1, iNT):
    for i in range(1, iNT):
        isigma[i] = isigma[i - 1] * np.exp(igsig[i - 1] * tstep)
        icost1[i] = pbacktime[i] * isigma[i] / expcost2 / 1000

@njit("(float64[:], int32)")
def InitializeCarbonTree(icumetree, iNT):
    for i in range(1, iNT):
        icumetree[i] = icumetree[i - 1] + etree[i - 1] * (5 / 3.666)

"""
Functions of the model
"""

"""
First: Functions related to emissions of carbon and weather damages
"""

# Retuns the total carbon emissions; Eq. 18
@njit("float64(float64[:],int32)")
def fE(iEIND, index):
    return iEIND[index] + etree[index]

# Eq.14: Determines the emission of carbon by industry EIND
@njit("float64(float64[:],float64[:],float64[:],int32)")
def fEIND(iYGROSS, iMIU, isigma, index):
    return isigma[index] * iYGROSS[index] * (1 - iMIU[index])

# Cumulative industrial emission of carbon
@njit("float64(float64[:],float64[:],int32)")
def fCCA(iCCA, iEIND, index):
    return iCCA[index - 1] + iEIND[index - 1] * 5 / 3.666

# Cumulative total carbon emission
@njit("float64(float64[:],float64[:],int32)")
def fCCATOT(iCCA, icumetree, index):
    return iCCA[index] + icumetree[index]

# Eq. 22: the dynamics of the radiative forcing
@njit("float64(float64[:],int32)")
def fFORC(iMAT, index):
    return fco22x * np.log(iMAT[index] / 588.000) / np.log(2) + forcoth[index]

# Dynamics of Omega; Eq.9
@njit("float64(float64[:],int32)")
def fDAMFRAC(iTATM, index):
    return a1 * iTATM[index] + a2 * iTATM[index] ** a3

# Calculate damages as a function of Gross industrial production; Eq.8
@njit("float64(float64[:],float64[:],int32)")
def fDAMAGES(iYGROSS, iDAMFRAC, index):
    return iYGROSS[index] * iDAMFRAC[index]

# Dynamics of Lambda; Eq. 10 - cost of the reudction of carbon emission (Abatement cost)
@njit("float64(float64[:],float64[:],float64[:],int32)")
def fABATECOST(iYGROSS, iMIU, icost1, index):
    return iYGROSS[index] * icost1[index] * iMIU[index] ** expcost2

# Marginal Abatement cost
@njit("float64(float64[:],int32)")
def fMCABATE(iMIU, index):
    return pbacktime[index] * iMIU[index] ** (expcost2 - 1)

# Price of carbon reduction
@njit("float64(float64[:],int32)")
def fCPRICE(iMIU, index):
    return pbacktime[index] * (iMIU[index]) ** (expcost2 - 1)

# Eq. 19: Dynamics of the carbon concentration in the atmosphere
@njit("float64(float64[:],float64[:],float64[:],int32)")
def fMAT(iMAT, iMU, iE, index):
    if index == 0:
        return mat0
    else:
        return (
            iMAT[index - 1] * b11
            + iMU[index - 1] * b21
            + iE[index - 1] * 5 / 3.666
        )

# Eq. 21: Dynamics of the carbon concentration in the ocean LOW level
@njit("float64(float64[:],float64[:],int32)")
def fML(iML, iMU, index):
    if index == 0:
        return ml0
    else:
        return iML[index - 1] * b33 + iMU[index - 1] * b23

# Eq. 20: Dynamics of the carbon concentration in the ocean UP level
@njit("float64(float64[:],float64[:],float64[:],int32)")
def fMU(iMAT, iMU, iML, index):
    if index == 0:
        return mu0
    else:
        return (
            iMAT[index - 1] * b12 + iMU[index - 1] * b22 + iML[index - 1] * b32
        )

# Eq. 23: Dynamics of the atmospheric temperature
@njit("float64(float64[:],float64[:],float64[:],int32)")
def fTATM(iTATM, iFORC, iTOCEAN, index):
    if index == 0:
        return tatm0
    else:
        return iTATM[index - 1] + c1 * (
            iFORC[index]
            - (fco22x / t2xco2) * iTATM[index - 1]
            - c3 * (iTATM[index - 1] - iTOCEAN[index - 1])
        )

# Eq. 24: Dynamics of the ocean temperature
@njit("float64(float64[:],float64[:],int32)")
def fTOCEAN(iTATM, iTOCEAN, index):
    if index == 0:
        return tocean0
    else:
        return iTOCEAN[index - 1] + c4 * (
            iTATM[index - 1] - iTOCEAN[index - 1]
        )

"""
Second: Function related to economic variables
"""

# The total production without climate losses denoted previously by YGROSS
@njit("float64(float64[:],float64[:],float64[:],int32)")
def fYGROSS(ial, il, iK, index):
    return ial[index] * ((il[index] / 1000) ** (1 - gama)) * iK[index] ** gama

# The production under the climate damages cost
@njit("float64(float64[:],float64[:],int32)")
def fYNET(iYGROSS, iDAMFRAC, index):
    return iYGROSS[index] * (1 - iDAMFRAC[index])

# Production after abatement cost
@njit("float64(float64[:],float64[:],int32)")
def fY(iYNET, iABATECOST, index):
    return iYNET[index] - iABATECOST[index]

# Consumption Eq. 11
@njit("float64(float64[:],float64[:],int32)")
def fC(iY, iI, index):
    return iY[index] - iI[index]

# Per capita consumption, Eq. 12
@njit("float64(float64[:],float64[:],int32)")
def fCPC(iC, il, index):
    return 1000 * iC[index] / il[index]

# Saving policy: investment
@njit("float64(float64[:],float64[:],int32)")
def fI(iS, iY, index):
    return iS[index] * iY[index]

# Capital dynamics Eq. 13
@njit("float64(float64[:],float64[:],int32)")
def fK(iK, iI, index):
    if index == 0:
        return k0
    else:
        return (1 - dk) ** tstep * iK[index - 1] + tstep * iI[index - 1]

# Interest rate equation; Eq. 26 added in personal notes
@njit("float64(float64[:],int32)")
def fRI(iCPC, index):
    return (1 + prstp) * (iCPC[index + 1] / iCPC[index]) ** (
        elasmu / tstep
    ) - 1

# Periodic utility: A form of Eq. 2
@njit("float64(float64[:],float64[:],int32)")
def fCEMUTOTPER(iPERIODU, il, index):
    return iPERIODU[index] * il[index] * rr[index]

# The term between brackets in Eq. 2
@njit("float64(float64[:],float64[:],int32)")
def fPERIODU(iC, il, index):
    return ((iC[index] * 1000 / il[index]) ** (1 - elasmu) - 1) / (
        1 - elasmu
    ) - 1

# utility function
@guvectorize([(float64[:], float64[:])], "(n), (m)")
def fUTILITY(iCEMUTOTPER, resUtility):
    resUtility[0] = tstep * scale1 * np.sum(iCEMUTOTPER) + scale2

"""
In this part we implement the objective function
"""

# * Control rate limits
MIU_lo = np.full(NT, 0.01)
MIU_up = np.full(NT, limmiu)
MIU_up[0:29] = 1
MIU_lo[0] = miu0
MIU_up[0] = miu0
MIU_lo[MIU_lo == MIU_up] = 0.99999 * MIU_lo[MIU_lo == MIU_up]
bnds1 = []
for i in range(NT):
    bnds1.append((MIU_lo[i], MIU_up[i]))
# * Control variables
lag10 = t > NT - 10
S_lo = np.full(NT, 1e-1)
S_lo[lag10] = optlrsav
S_up = np.full(NT, 0.9)
S_up[lag10] = optlrsav
S_lo[S_lo == S_up] = 0.99999 * S_lo[S_lo == S_up]
bnds2 = []
for i in range(NT):
    bnds2.append((S_lo[i], S_up[i]))

# Arbitrary starting values for the control variables:
S_start = np.full(NT, 0.2)
S_start[S_start < S_lo] = S_lo[S_start < S_lo]
S_start[S_start > S_up] = S_lo[S_start > S_up]
MIU_start = 0.99 * MIU_up
MIU_start[MIU_start < MIU_lo] = MIU_lo[MIU_start < MIU_lo]
MIU_start[MIU_start > MIU_up] = MIU_up[MIU_start > MIU_up]

K = np.zeros(NT)
YGROSS = np.zeros(NT)
EIND = np.zeros(NT)
E = np.zeros(NT)
CCA = np.zeros(NT)
CCATOT = np.zeros(NT)
MAT = np.zeros(NT)
ML = np.zeros(NT)
MU = np.zeros(NT)
FORC = np.zeros(NT)
TATM = np.zeros(NT)
TOCEAN = np.zeros(NT)
DAMFRAC = np.zeros(NT)
DAMAGES = np.zeros(NT)
ABATECOST = np.zeros(NT)
MCABATE = np.zeros(NT)
CPRICE = np.zeros(NT)
YNET = np.zeros(NT)
Y = np.zeros(NT)
I = np.zeros(NT)
C = np.zeros(NT)
CPC = np.zeros(NT)
RI = np.zeros(NT)
PERIODU = np.zeros(NT)
CEMUTOTPER = np.zeros(NT)

# The objective function
# It returns the utility as scalar
def fOBJ(
    x,
    sign,
    iI,
    iK,
    ial,
    il,
    iYGROSS,
    isigma,
    iEIND,
    iE,
    iCCA,
    iCCATOT,
    icumetree,
    iMAT,
    iMU,
    iML,
    iFORC,
    iTATM,
    iTOCEAN,
    iDAMFRAC,
    iDAMAGES,
    iABATECOST,
    icost1,
    iMCABATE,
    iCPRICE,
    iYNET,
    iY,
    iC,
    iCPC,
    iPERIODU,
    iCEMUTOTPER,
    iRI,
    iNT,
):

    iMIU = x[0:NT]
    iS = x[NT : (2 * NT)]

    for i in range(iNT):
        iK[i] = fK(iK, iI, i)
        iYGROSS[i] = fYGROSS(ial, il, iK, i)
        iEIND[i] = fEIND(iYGROSS, iMIU, isigma, i)
        iE[i] = fE(iEIND, i)
        iCCA[i] = fCCA(iCCA, iEIND, i)
        iCCATOT[i] = fCCATOT(iCCA, icumetree, i)
        iMAT[i] = fMAT(iMAT, iMU, iE, i)
        iML[i] = fML(iML, iMU, i)
        iMU[i] = fMU(iMAT, iMU, iML, i)
        iFORC[i] = fFORC(iMAT, i)
        iTATM[i] = fTATM(iTATM, iFORC, iTOCEAN, i)
        iTOCEAN[i] = fTOCEAN(iTATM, iTOCEAN, i)
        iDAMFRAC[i] = fDAMFRAC(iTATM, i)
        iDAMAGES[i] = fDAMAGES(iYGROSS, iDAMFRAC, i)
        iABATECOST[i] = fABATECOST(iYGROSS, iMIU, icost1, i)
        iMCABATE[i] = fMCABATE(iMIU, i)
        iCPRICE[i] = fCPRICE(iMIU, i)
        iYNET[i] = fYNET(iYGROSS, iDAMFRAC, i)
        iY[i] = fY(iYNET, iABATECOST, i)
        iI[i] = fI(iS, iY, i)
        iC[i] = fC(iY, iI, i)
        iCPC[i] = fCPC(iC, il, i)
        iPERIODU[i] = fPERIODU(iC, il, i)
        iCEMUTOTPER[i] = fCEMUTOTPER(iPERIODU, il, i)
        iRI = fRI(iCPC, i)
    with open("parameter_check.csv", "w") as myfile:
        for parameter in [
            iK,
            iYGROSS,
            iEIND,
            iE,
            iCCA,
            iCCATOT,
            iMAT,
            iML,
            iMU,
            iFORC,
            iTATM,
            iTOCEAN,
            iDAMFRAC,
            iDAMAGES,
            iABATECOST,
            iMCABATE,
            iCPRICE,
            iYNET,
            iY,
            iI,
            iC,
            iCPC,
            iPERIODU,
            iCEMUTOTPER,
        ]:
            wr = csv.writer(myfile, delimiter="\t", quoting=csv.QUOTE_ALL)
            wr.writerow(parameter)
    resUtility = np.zeros(1)
    fUTILITY(iCEMUTOTPER, resUtility)

    return sign * resUtility[0]

# For the optimal allocation of x, calculates the whole system variables
def Optimality(
    x,
    iI,
    iK,
    ial,
    il,
    iYGROSS,
    isigma,
    iEIND,
    iE,
    iCCA,
    iCCATOT,
    icumetree,
    iMAT,
    iMU,
    iML,
    iFORC,
    iTATM,
    iTOCEAN,
    iDAMFRAC,
    iDAMAGES,
    iABATECOST,
    icost1,
    iMCABATE,
    iCPRICE,
    iYNET,
    iY,
    iC,
    iCPC,
    iPERIODU,
    iCEMUTOTPER,
    iRI,
    iNT,
):

    iMIU = x[0:NT]
    iS = x[NT : (2 * NT)]

    for i in range(iNT):
        iK[i] = fK(iK, iI, i)
        iYGROSS[i] = fYGROSS(ial, il, iK, i)
        iEIND[i] = fEIND(iYGROSS, iMIU, isigma, i)
        iE[i] = fE(iEIND, i)
        iCCA[i] = fCCA(iCCA, iEIND, i)
        iCCATOT[i] = fCCATOT(iCCA, icumetree, i)
        iMAT[i] = fMAT(iMAT, iMU, iE, i)
        iML[i] = fML(iML, iMU, i)
        iMU[i] = fMU(iMAT, iMU, iML, i)
        iFORC[i] = fFORC(iMAT, i)
        iTATM[i] = fTATM(iTATM, iFORC, iTOCEAN, i)
        iTOCEAN[i] = fTOCEAN(iTATM, iTOCEAN, i)
        iDAMFRAC[i] = fDAMFRAC(iTATM, i)
        iDAMAGES[i] = fDAMAGES(iYGROSS, iDAMFRAC, i)
        iABATECOST[i] = fABATECOST(iYGROSS, iMIU, icost1, i)
        iMCABATE[i] = fMCABATE(iMIU, i)
        iCPRICE[i] = fCPRICE(iMIU, i)
        iYNET[i] = fYNET(iYGROSS, iDAMFRAC, i)
        iY[i] = fY(iYNET, iABATECOST, i)
        iI[i] = fI(iS, iY, i)
        iC[i] = fC(iY, iI, i)
        iCPC[i] = fCPC(iC, il, i)
        iPERIODU[i] = fPERIODU(iC, il, i)
        iCEMUTOTPER[i] = fCEMUTOTPER(iPERIODU, il, i)
        iRI[i] = fRI(iCPC, i)

    resUtility = np.zeros(1)
    fUTILITY(iCEMUTOTPER, resUtility)

    return (
        resUtility[0],
        iI,
        iK,
        ial,
        il,
        iYGROSS,
        isigma,
        iEIND,
        iE,
        iCCA,
        iCCATOT,
        icumetree,
        iMAT,
        iMU,
        iML,
        iFORC,
        iTATM,
        iTOCEAN,
        iDAMFRAC,
        iDAMAGES,
        iABATECOST,
        icost1,
        iMCABATE,
        iCPRICE,
        iYNET,
        iY,
        iC,
        iCPC,
        iPERIODU,
        iCEMUTOTPER,
        iRI,
    )

def PlotFigures():

    figTATM = plt.figure()
    plt.plot(TT, TATM)
    figTATM.suptitle(
        "Increase temperature of the atmosphere (TATM)", fontsize=20
    )
    plt.xlabel("Years", fontsize=16)
    plt.ylabel("Degrees C from 1900", fontsize=16)

    figTOCEAN = plt.figure()
    plt.plot(TT, TOCEAN)
    figTOCEAN.suptitle(
        "Increase temperature of the ocean (TOCEAN)", fontsize=20
    )
    plt.xlabel("Years", fontsize=16)
    plt.ylabel("Degrees C from 1900", fontsize=16)

    figMU = plt.figure()
    plt.plot(TT, MU)
    figMU.suptitle(
        "Carbon concentration increase in shallow oceans (MU)", fontsize=20
    )
    plt.xlabel("Years", fontsize=16)
    plt.ylabel("GtC from 1750", fontsize=16)

    figML = plt.figure()
    plt.plot(TT, ML)
    figML.suptitle(
        "Carbon concentration increase in lower oceans (ML)", fontsize=20
    )
    plt.xlabel("Years", fontsize=16)
    plt.ylabel("GtC from 1750", fontsize=16)

    figDAM = plt.figure()
    plt.plot(TT, DAMAGES)
    figDAM.suptitle("Damages", fontsize=20)
    plt.xlabel("Years", fontsize=16)
    plt.ylabel("trillions 2010 USD per year", fontsize=16)

    figDAMFRAC = plt.figure()
    plt.plot(TT, DAMFRAC)
    figDAMFRAC.suptitle("Damages as fraction of gross output", fontsize=20)
    plt.xlabel("Years", fontsize=16)
    plt.ylabel("", fontsize=16)

    figCOSTRED = plt.figure()
    plt.plot(TT, ABATECOST)
    figCOSTRED.suptitle("Cost of emissions reductions", fontsize=20)
    plt.xlabel("Years", fontsize=16)
    plt.ylabel("trillions 2010 USD per year", fontsize=16)

    figMarg = plt.figure()
    plt.plot(TT, MCABATE)
    figMarg.suptitle("Marginal abatement cost", fontsize=20)
    plt.xlabel("Years", fontsize=16)
    plt.ylabel("2010 USD per ton CO2", fontsize=16)

    figMIU = plt.figure()
    plt.plot(TT, result.x[0:NT])
    figMIU.suptitle("Carbon emission control rate", fontsize=20)
    plt.xlabel("Years", fontsize=16)
    plt.ylabel("Rate", fontsize=16)

    figE = plt.figure()
    plt.plot(TT, E)
    figE.suptitle("Total CO2 emission", fontsize=20)
    plt.xlabel("Years", fontsize=16)
    plt.ylabel("GtCO2 per year", fontsize=16)

    figMAT = plt.figure()
    plt.plot(TT, MAT)
    figMAT.suptitle(
        "Carbon concentration increase in the atmosphere", fontsize=20
    )
    plt.xlabel("Years", fontsize=16)
    plt.ylabel("GtC from 1750", fontsize=16)

    figFORC = plt.figure()
    plt.plot(TT, FORC)
    figFORC.suptitle("Increase in radiative forcing", fontsize=20)
    plt.xlabel("Years", fontsize=16)
    plt.ylabel("watts per m2 from 1900", fontsize=16)

    figRI = plt.figure()
    plt.plot(TT, RI)
    figRI.suptitle("Real interest rate", fontsize=20)
    plt.xlabel("Years", fontsize=16)
    plt.ylabel("Rate per annum", fontsize=16)

    figC = plt.figure()
    plt.plot(TT, C)
    figC.suptitle("Consumption", fontsize=20)
    plt.xlabel("Years", fontsize=16)
    plt.ylabel("trillions 2010 USD per year", fontsize=16)

    figY = plt.figure()
    plt.plot(TT, Y)
    figY.suptitle("Gross product net of abatement and damages", fontsize=20)
    plt.xlabel("Years", fontsize=16)
    plt.ylabel("trillions 2010 USD per year", fontsize=16)

    figYGROSS = plt.figure()
    plt.plot(TT, YGROSS)
    figYGROSS.suptitle("World gross product", fontsize=20)
    plt.xlabel("Years", fontsize=16)
    plt.ylabel("trillions 2010 USD per year", fontsize=16)

    figYGROSSbyY = plt.figure()
    plt.plot(TT, YGROSS - Y)
    figYGROSSbyY.suptitle("Abatement and damages costs", fontsize=20)
    plt.xlabel("Years", fontsize=16)
    plt.ylabel("trillions 2010 USD per year", fontsize=16)

    figS = plt.figure()
    plt.plot(TT, result.x[NT : (2 * NT)])
    figS.suptitle("Saving rate", fontsize=20)
    plt.xlabel("Years", fontsize=16)
    plt.ylabel("rate", fontsize=16)

    figI = plt.figure()
    plt.plot(TT, I)
    figI.suptitle("Investment (I)", fontsize=20)
    plt.xlabel("Years", fontsize=16)
    plt.ylabel("trillions 2010 USD per year", fontsize=16)

    plt.show()

if __name__ == "__main__":

    TT = np.linspace(2000, 2500, 100, dtype=np.int32)

    InitializeLabor(l, NT)
    InitializeTFP(al, NT)
    InitializeGrowthSigma(gsig, NT)
    InitializeSigma(sigma, gsig, cost1, NT)
    InitializeCarbonTree(cumetree, NT)

    x_start = np.concatenate([MIU_start, S_start])
    bnds = bnds1 + bnds2
    x = np.ones(2 * NT) / 4
    print(
        fOBJ(
            x,
            1.0,
            I,
            K,
            al,
            l,
            YGROSS,
            sigma,
            EIND,
            E,
            CCA,
            CCATOT,
            cumetree,
            MAT,
            MU,
            ML,
            FORC,
            TATM,
            TOCEAN,
            DAMFRAC,
            DAMAGES,
            ABATECOST,
            cost1,
            MCABATE,
            CPRICE,
            YNET,
            Y,
            C,
            CPC,
            PERIODU,
            CEMUTOTPER,
            RI,
            NT,
        )
    )

reason: import numpy,import scipy,from numba


text:
import numpy as np
import nibabel as nib
import os
import sys
import pdb
import glob

def load_data(path):
    data_dim, data, data_valid = load_ni(path)

    data = data.reshape(data_dim)
    data = np.moveaxis(data,1,2)
    data = np.moveaxis(data,2,0)
    data_dim = data.shape
    data = data.reshape((data_dim[0]*data_dim[1]*data_dim[2]))

    return data, data_dim

def load_ni(path):
    filename, file_extension = os.path.splitext(path)
    if file_extension=='.mnc':
        data_dim, data, data_valid = load_minc(path)
    else:
        data_dim, data, data_valid = load_nii(path)
    return data_dim, data, data_valid

def load_minc(path):
    ni_img = nib.load(path)
    ni_data = ni_img.get_data()
    ni_data_dim = ni_data.shape
    data_dim = ni_data_dim
    data = ni_data.reshape((ni_data_dim[0]*ni_data_dim[1]*ni_data_dim[2]))
    data_valid = True

    return data_dim, data, data_valid

def load_nii(path):
    ni_img = nib.load(path)
    ni_data = ni_img.get_data()
    ni_data_dim = ni_data.shape
    data_dim = ni_data_dim[0:3]
    if (len(ni_data_dim)==4):
       data = ni_data[:,:,:,0].reshape((ni_data_dim[0]*ni_data_dim[1]*ni_data_dim[2]))
    else:
       data = ni_data.reshape((ni_data_dim[0]*ni_data_dim[1]*ni_data_dim[2]))
    data_valid = True

    return data_dim, data, data_valid

files = glob.glob('/mnt/hippocampus/mri_data/brainweb/*.bc.nrm.nii')

for fid in files:
    in_path = fid
    in_filename, in_file_extension = os.path.splitext(in_path)
    out_path = in_filename + '.bin'

    data, data_dim = load_data(in_path)

    print data_dim
    with open(out_path,'wb') as f:
        f.write(np.asarray(data_dim,dtype=np.uint32).tobytes())
        f.write(np.asarray(data,dtype=np.float32).tobytes())

fid = '/mnt/hippocampus/mri_data/brainweb/phantom_1.0mm_normal_crisp.clean_ref.nii'
in_path = fid
in_filename, in_file_extension = os.path.splitext(in_path)
out_path = in_filename + '.bin'
data, data_dim = load_data(in_path)
print data_dim
with open(out_path,'wb') as f:
    f.write(np.asarray(data_dim,dtype=np.uint32).tobytes())
    f.write(np.asarray(data,dtype=np.float32).tobytes())

max_stars_repo_name: summerblind/neurotrace

max_issues_repo_name: summerblind/neurotrace

max_forks_repo_name: summerblind/neurotrace


text:
import numpy as np

SLOPE = (1, 3)

def part1(input_str):
    input_split = input_str.split('\n')
    data = np.zeros((len(input_split), len(input_split[0])), dtype=str)
    for row, line in enumerate(input_split):
        data[row, :] = np.array([char for char in line])

    out = ''
    for row in range(data.shape[0]):
        for col in range(data.shape[1]):
            out += data[row, col]
        out += '\n'
    out = out[:-1]
    assert input_str == out

    travel_path = [(0, 0)]
    while travel_path[-1][0] < data.shape[0] - SLOPE[0]:
        travel_path.append(tuple(sum(x) for x in zip(travel_path[-1], SLOPE)))
        if travel_path[-1][1] >= data.shape[1]:
            travel_path[-1] = (travel_path[-1][0], travel_path[-1][1] - data.shape[1])
    output = data.copy()
    for idx in travel_path[2:]:
        output[idx] = 'X' if data[idx] == '#' else 'O'
    num_trees = len(output[output == 'X'])
    return num_trees


text:
from typing import Dict
import numpy as np

class CcmClickModel:
    """
    An implementation of the CCM click model, see the paper for the details and references.
    """

    def __init__(self, click_relevance: Dict[int, float], stop_relevance: Dict[int, float], name: str, depth: int):
        """
        :param click_relevance: A mapping from a relevance label to the probability of the doc being clicked (assuming
        it is examined)
        :param stop_relevance: A mapping from a relevance label to the probability of the user stopping after clicking
        on this document
        :param name: Name of the click model instance
        :param depth: How deep the user examines the result page
        """
        self.click_relevance = click_relevance
        self.stop_relevance = stop_relevance
        self.name = name
        self.depth = depth

    def __call__(self, relevances: np.ndarray, random_state: np.random.RandomState) -> np.ndarray:
        """
        Generates an indicator array of the click events for the ranked documents with relevance labels encoded in
        `relevances`.
        :param relevances: Relevance labels of the documents
        :param random_state: Random generator state
        :return: Indicator array of the clicks on the documents

        As an example, consider a model of a user who always clicks on a highly relevant result and immediately stops
        >>> model = CcmClickModel(click_relevance={0: 0.0, 1: 0.0, 2: 1.0},
        ...                       stop_relevance={0: 0.0, 1: 0.0, 2: 1.0}, name="Model", depth=10)
        >>> # With the result list with highly relevant docs on positions 2 and 4,
        >>> doc_relevances = np.array([1, 0, 2, 0, 2, 0])
        >>> # We expect the user to click on the 3rd document, as it the first highly relevant:
        >>> model(doc_relevances, np.random.RandomState(1)).tolist()
        [0.0, 0.0, 1.0, 0.0, 0.0, 0.0]
        """
        n_docs = relevances.shape[0]
        click_label = np.zeros(n_docs)
        for i in range(min(self.depth, n_docs)):
            r = relevances[i]

            p_click = self.click_relevance[r]
            p_stop = self.stop_relevance[r]

            click_simu = random_state.uniform()
            stop_simu = random_state.uniform()

            if click_simu < p_click:
                click_label[i] = 1
                if stop_simu < p_stop:
                    break
        return click_label

# add PBM click model
class PbmClickModel:
    """
    An implementation of the PBM click model.
    """

    def __init__(self, click_relevance: Dict[int, float], stop_relevance: Dict[int, float], name: str, depth: int):
        """
        :param click_relevance: A mapping from a relevance label to the probability of the doc being clicked (assuming
        it is examined)
        :param stop_relevance: A mapping from a relevance label to the probability of the user stopping after clicking
        on this document
        :param name: Name of the click model instance
        :param depth: How deep the user examines the result page
        """
        self.click_relevance = click_relevance
        self.stop_relevance = stop_relevance
        self.name = name
        self.depth = depth

    def __call__(self, relevances: np.ndarray, random_state: np.random.RandomState, eta = 1) -> np.ndarray:
        """
        Generates an indicator array of the click events for the ranked documents with relevance labels encoded in
        `relevances`.
        :param relevances: Relevance labels of the documents
        :param random_state: Random generator state
        :return: Indicator array of the clicks on the documents
        """
        n_docs = relevances.shape[0]
        result = np.zeros(n_docs)
        for i in range(min(self.depth, n_docs)):
            r = relevances[i]
            propensity = (1.0 / (i + 1)) ** eta # position bias
            p_click = self.click_relevance[r] * propensity

            if random_state.uniform() < p_click:
                result[i] = 1
        return result


text:
# -*- coding: utf-8 -*-
import numpy as np

from vispy.scene.visuals import Image
from vispy.testing import (requires_application, TestingCanvas,
                           assert_image_equal, run_tests_if_main)

@requires_application()
def test_image():
    """Test image visual"""
    size = (100, 50)
    with TestingCanvas(size=size, bgcolor='w') as c:
        for three_d in (True, False):
            shape = size[::-1] + ((3,) if three_d else ())
            data = np.random.rand(*shape)
            image = Image(data, cmap='grays', clim=[0, 1])
            c.draw_visual(image)
            if three_d:
                expected = data
            else:
                expected = np.tile(data[:, :, np.newaxis], (1, 1, 3))
            assert_image_equal("screenshot", expected)

run_tests_if_main()

max_stars_repo_name: izaid/vispy

max_issues_repo_name: izaid/vispy

max_forks_repo_name: izaid/vispy


text:
"""
:Authors: - Wilker Aziz
"""
import dgm4nlp.device as dev
config = dev.tf_config('0', allow_growth=True)

import os
import numpy as np
import logging
import tempfile
import shutil
from datetime import datetime
from dgm4nlp.annealing import AnnealingSchedule
from embedalign.trainer import test_embedalign
from embedalign.infnet import get_infnet

import tensorflow as tf

np.random.seed(42)

logging.basicConfig(level=logging.INFO, format='%(asctime)s %(message)s', datefmt='%Y-%m-%d %H:%M:%S')

# 1. Get a unique working directory and save script for reproducibility
base_dir = 'debug/hansards'
os.makedirs(base_dir, exist_ok=True)  # make sure base_dir exists
output_dir = tempfile.mkdtemp(prefix=datetime.now().strftime("%y-%m-%d.%Hh%Mm%Ss."), dir=base_dir)
logging.info('Workspace: %s', output_dir)
shutil.copy(os.path.abspath(__file__), output_dir)

test_embedalign(
    training_x='hansards/shuffled.e',
    training_y='hansards/shuffled.f',
    val_x='hansards/trial.en-fr.en',
    val_y='hansards/trial.en-fr.fr',
    val_naacl='hansards/trial.en-fr.naacl',
    test_x='hansards/test.en-fr.en',
    test_y='hansards/test.en-fr.fr',
    test_naacl='hansards/test.en-fr.naacl',
    output_dir=output_dir,
    # preprocessing
    nb_words=[30000, 30000],
    longest_sequence=[50, 50],
    # architecture
    dx=128,
    dh=128,
    q=get_infnet('rnn-z', dz=100),
    nb_softmax_samples=[1000, 1000],
    softmax_approximation=['botev-batch', 'botev-batch'],
    softmax_q_inv_temperature=[0., 0.],
    # optimiser
    nb_epochs=20,
    batch_size=100,
    optimizer=tf.train.AdamOptimizer(learning_rate=0.001),
    grad_clipping=[-1., 1.],
    annealing_z=AnnealingSchedule(initial=0., final=1., step=0.01, nb_updates=100),
    config=config
)

max_stars_repo_name: uva-slpl/embedalign

max_issues_repo_name: uva-slpl/embedalign

max_forks_repo_name: uva-slpl/embedalign


text:
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Simple feature that checks whether this is a tweet as part of a thread
@author: marcelklehr
"""

import numpy as np
import re
from code.feature_extraction.feature_extractor import FeatureExtractor

class Threads(FeatureExtractor):
    """class for extracting the character-based length as a feature"""

    # constructor
    def __init__(self, input_column):
        super().__init__([input_column], "{0}_threads".format(input_column))

    # don't need to fit, so don't overwrite _set_variables()

    # compute the word length based on the inputs
    def _get_values(self, inputs):
        result = np.array(inputs[0].str.contains(r'🧵|^\s?[0-9]/[0-9]?\s|\s[0-9]/[0-9]?\s?$', flags=re.IGNORECASE, regex=True))
        result = result.reshape(-1, 1)
        return result

max_stars_repo_name: louiskhub/TweetViralityClassifier

max_issues_repo_name: louiskhub/TweetViralityClassifier

max_forks_repo_name: louiskhub/TweetViralityClassifier


text:
# Copyright 2021 Huawei Technologies Co., Ltd
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
# ============================================================================
"""train posenet"""
import ast
import argparse
import os
import shutil
import numpy as np
from mindspore.common import set_seed
from mindspore import Tensor, load_checkpoint, load_param_into_net, export, context
from mindspore.context import ParallelMode
from mindspore.communication.management import init, get_rank
from mindspore.nn import Adagrad
from mindspore.train.model import Model
from mindspore.train.callback import ModelCheckpoint, CheckpointConfig, LossMonitor, TimeMonitor
import mindspore.common.dtype as ms
from src.posenet import PoseNet
from src.config import common_config, KingsCollege, StMarysChurch
from src.dataset import create_posenet_dataset
from src.loss import PosenetWithLoss

set_seed(1)

parser = argparse.ArgumentParser(description='Posenet train.')
parser.add_argument("--run_distribute", type=ast.literal_eval, default=False, help="Run distribute, default is false.")
parser.add_argument('--dataset', type=str, default='KingsCollege',
                    choices=['KingsCollege', 'StMarysChurch'], help='Name of dataset.')
parser.add_argument('--device_num', type=int, default=1, help='Number of device.')
# 模型输出目录
parser.add_argument('--train_url', type=str, default='', help='the path model saved')
# 数据集目录
parser.add_argument('--data_url', type=str, default='', help='the training data')
# 抽取出来的超参配置
parser.add_argument('--pre_trained', type=ast.literal_eval, default=False, help='Pretrained checkpoint path')
parser.add_argument('--device_id', type=int, default=None, help='device id of GPU or Ascend. (Default: None)')
parser.add_argument('--device_target', type=str, default='Ascend', choices=("Ascend", "GPU", "CPU"),
                    help="Device target, support Ascend, GPU and CPU.")
parser.add_argument('--max_steps', type=int, default=30000, help='max_steps')
parser.add_argument('--save_checkpoint_epochs', type=int, default=5, help='save_checkpoint_epochs')
parser.add_argument('--keep_checkpoint_max', type=int, default=10, help='keep_checkpoint_max')
parser.add_argument('--save_checkpoint', type=ast.literal_eval, default=True, help='save_checkpoint')
parser.add_argument("--file_name", type=str, default="posenet", help="output file name.")
parser.add_argument('--is_modelarts', type=ast.literal_eval, default=True, help='Train in Modelarts.')
args_opt = parser.parse_args()
CACHE_TRAINING_URL = "/cache/training/"
CACHE = "/cache/"
src = "/"
local_data_path = '/cache/data/'
if not os.path.isdir(CACHE_TRAINING_URL):
    os.makedirs(CACHE_TRAINING_URL)

if __name__ == '__main__':
    cfg = common_config
    if args_opt.dataset == "KingsCollege":
        dataset_cfg = KingsCollege
    elif args_opt.dataset == "StMarysChurch":
        dataset_cfg = StMarysChurch

    device_target = args_opt.device_target
    context.set_context(mode=context.GRAPH_MODE, device_target=device_target)
    if args_opt.run_distribute:
        if device_target == "Ascend":
            device_id = int(os.getenv('DEVICE_ID'))
            context.set_context(device_id=device_id, enable_auto_mixed_precision=True)
            context.set_auto_parallel_context(parallel_mode=ParallelMode.DATA_PARALLEL,
                                              gradients_mean=True,
                                              auto_parallel_search_mode="recursive_programming")
            init()
        elif device_target == "GPU":
            init()
            context.set_auto_parallel_context(device_num=args_opt.device_num,
                                              parallel_mode=ParallelMode.DATA_PARALLEL,
                                              gradients_mean=True,
                                              auto_parallel_search_mode="recursive_programming")
    else:
        if args_opt.device_id is not None:
            context.set_context(device_id=args_opt.device_id)
        else:
            context.set_context(device_id=cfg.device_id)

    train_dataset_path = dataset_cfg.dataset_path
    if args_opt.is_modelarts:
        import moxing as mox
        mox.file.copy_parallel(args_opt.data_url, CACHE)
    if args_opt.dataset == "KingsCollege":
        mindrecord_file_name = "KingsCollege_posenet_train.mindrecord"
    elif args_opt.dataset == "StMarysChurch":
        mindrecord_file_name = "StMarysChurch_posenet_train.mindrecord"
    mindrecord_file = os.path.join(CACHE, mindrecord_file_name)
    dataset = create_posenet_dataset(mindrecord_file, batch_size=dataset_cfg.batch_size,
                                     device_num=args_opt.device_num, is_training=True)
    step_per_epoch = dataset.get_dataset_size()

    net_with_loss = PosenetWithLoss(args_opt.pre_trained)
    opt = Adagrad(params=net_with_loss.trainable_params(),
                  learning_rate=dataset_cfg.lr_init,
                  weight_decay=dataset_cfg.weight_decay)
    model = Model(net_with_loss, optimizer=opt)

    time_cb = TimeMonitor(data_size=step_per_epoch)
    loss_cb = LossMonitor()
    cb = [time_cb, loss_cb]
    if args_opt.save_checkpoint:
        config_ck = CheckpointConfig(save_checkpoint_steps=args_opt.save_checkpoint_epochs * step_per_epoch,
                                     keep_checkpoint_max=args_opt.keep_checkpoint_max)
        if args_opt.is_modelarts:
            save_checkpoint_path = CACHE_TRAINING_URL
            if args_opt.device_num == 1:
                ckpt_cb = ModelCheckpoint(prefix='train_posenet_' + args_opt.dataset,
                                          directory=save_checkpoint_path,
                                          config=config_ck)
                cb += [ckpt_cb]
            if args_opt.device_num > 1 and get_rank() % 8 == 0:
                ckpt_cb = ModelCheckpoint(prefix='train_posenet_' + args_opt.dataset,
                                          directory=save_checkpoint_path,
                                          config=config_ck)
                cb += [ckpt_cb]
        else:
            save_checkpoint_path = cfg.checkpoint_dir
            if not os.path.isdir(save_checkpoint_path):
                os.makedirs(save_checkpoint_path)

            if args_opt.device_num == 1:
                ckpt_cb = ModelCheckpoint(prefix='train_posenet_' + args_opt.dataset,
                                          directory=save_checkpoint_path,
                                          config=config_ck)
                cb += [ckpt_cb]
            if args_opt.device_num > 1 and get_rank() % 8 == 0:
                ckpt_cb = ModelCheckpoint(prefix='train_posenet_' + args_opt.dataset,
                                          directory=save_checkpoint_path,
                                          config=config_ck)
                cb += [ckpt_cb]

    epoch_size = args_opt.max_steps // args_opt.device_num // step_per_epoch
    model.train(1, dataset, callbacks=cb)

    net = PoseNet()
    file_name1 = "train_posenet_KingsCollege-1_16.ckpt"
    assert cfg.checkpoint_dir is not None, "cfg.checkpoint_dir is None."
    param_dict = load_checkpoint(os.path.join(CACHE_TRAINING_URL, file_name1))
    load_param_into_net(net, param_dict)
    input_arr = Tensor(np.ones([1, 3, 224, 224]), ms.float32)
    export(net, input_arr, file_name=args_opt.file_name, file_format='AIR')
    shutil.copy('posenet.air', CACHE_TRAINING_URL)
    if args_opt.is_modelarts:
        mox.file.copy_parallel(src_url=CACHE_TRAINING_URL, dst_url=args_opt.train_url)

max_stars_repo_name: leelige/mindspore

max_issues_repo_name: leelige/mindspore

max_forks_repo_name: leelige/mindspore


text:
import numpy as np
from pyxsim.event_list import ConvolvedEventList
from soxs.instrument import \
    AuxiliaryResponseFile, \
    RedistributionMatrixFile
from soxs.utils import parse_prng
from yt.funcs import issue_deprecation_warning

class InstrumentSimulator(object):
    def __init__(self, name, arf_file, rmf_file):
        """
        Construct an instrument simulator.

        Parameters
        ----------
        inst_name : string
            The string corresponding to the name of the SOXS
            instrument specification.
        """
        self.name = name
        self.arf_file = arf_file
        self.rmf_file = rmf_file

    _arf = None
    @property
    def arf(self):
        if self._arf is None:
            self._arf = AuxiliaryResponseFile(self.arf_file)
        return self._arf

    _rmf = None
    @property
    def rmf(self):
        if self._rmf is None:
            self._rmf = RedistributionMatrixFile(self.rmf_file)
        return self._rmf

    def __call__(self, events, prng=None):
        """
        Calling method for :class:`~pyxsim.instruments.InstrumentSimulator`.

        Parameters
        ----------
        events : :class:`~pyxsim.events.EventList`
            An EventList instance of unconvolved events.
        prng : integer or :class:`~numpy.random.RandomState` object
            A pseudo-random number generator. Typically will only be specified
            if you have a reason to generate the same set of random numbers, such as for a
            test. Default is to use the :mod:`numpy.random` module.
        """
        issue_deprecation_warning("The pyXSIM built-in instrument simulators "
                                  "have been deprecated and will be removed "
                                  "in a future release!")
        if "pi" in events or "pha" in events:
            raise RuntimeError("These events have already been convolved with a response!!")
        prng = parse_prng(prng)
        flux = np.sum(events["eobs"]).to("erg") / \
               events.parameters["exp_time"]/events.parameters["area"]
        exp_time = events.parameters["exp_time"]
        emin = events["eobs"].min().value
        emax = events["eobs"].max().value
        new_events = {}
        new_events.update(events.events)
        new_events["energy"] = new_events.pop("eobs")
        new_events = self.arf.detect_events(new_events, exp_time, flux,
                                            [emin, emax], prng=prng)
        new_events = self.rmf.scatter_energies(new_events, prng=prng)
        new_events["eobs"] = new_events.pop("energy")
        chantype = self.rmf.header["CHANTYPE"].lower()
        new_events[chantype] = new_events.pop(self.rmf.header["CHANTYPE"])
        parameters = {}
        parameters.update(events.parameters)
        parameters["channel_type"] = chantype
        parameters["mission"] = self.rmf.header.get("MISSION", "")
        parameters["instrument"] = self.rmf.header["INSTRUME"]
        parameters["telescope"] = self.rmf.header["TELESCOP"]
        parameters["arf"] = self.arf.filename
        parameters["rmf"] = self.rmf.filename
        return ConvolvedEventList(new_events, parameters)

# Specific instrument approximations

ACIS_S = InstrumentSimulator("acis-s", "aciss_aimpt_cy19.arf",
                             "aciss_aimpt_cy19.rmf")
ACIS_I = InstrumentSimulator("acis-i", "acisi_aimpt_cy19.arf",
                             "acisi_aimpt_cy19.rmf")
Hitomi_SXS = InstrumentSimulator("hitomi_sxs", "hitomi_sxs_ptsrc.arf",
                                 "hitomi_sxs.rmf")
Athena_WFI = InstrumentSimulator("athena_wfi",
                                 "athena_wfi_1469_onaxis_w_filter_v20150326.arf",
                                 "athena_wfi_rmf_v20150326.rmf")
Athena_XIFU = InstrumentSimulator("athena_xifu",
                                  "athena_xifu_1469_onaxis_pitch249um_v20160401.arf",
                                  "athena_xifu_rmf_v20160401.rmf")
Lynx_Imager = InstrumentSimulator("lynx_hdxi", "xrs_hdxi_3x10.arf", "xrs_hdxi.rmf")
Lynx_Calorimeter = InstrumentSimulator("lynx_lxm", "xrs_mucal_3x10_3.0eV.arf", "xrs_mucal.rmf")

max_stars_repo_name: MaggieLieu/pyxsim

max_issues_repo_name: MaggieLieu/pyxsim

max_forks_repo_name: MaggieLieu/pyxsim