text:
###############################################################################
# MIT License
#
# Copyright (c) 2020 Daniel Vitor Ruiz
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.
###############################################################################

import os
import cv2
import numpy as np
from PIL import Image
import time
import random

# for fancy parameterization
from argparse import ArgumentParser

def parse_args():
    parser = ArgumentParser(
        description='Compute resulting image using IDA techinique for the DUTS dataset')

    parser.add_argument(
        '-obj_path', '--obj_path',
        type=str, default="/home/bakrinski/workspace/datasets/duts/training/images/",
        help='OBJ_FOLDER_IMG input images path'
    )

    parser.add_argument(
        '-obj_mask_path', '--obj_mask_path',
        type=str, default="/home/bakrinski/workspace/datasets/duts/training/colored_masks/",
        help='OBJ_FOLDER_MASK input masks path'
    )

    parser.add_argument(
        '-bg_path', '--bg_path',
        type=str, default="/home/bakrinski/workspace/datasets/duts/training/inpaintedv2_nogrid/",
        help='BG_FOLDER_IMG background images path'
    )

    parser.add_argument(
        '-index_obj_path', '--index_obj_path',
        type=str, default="dataset.txt",
        help='LIST_OF_N_OBJECTS filepath for the file containing per line a indice, e.g. "dataset.txt" resulting from computeKnn.py'
    )

    parser.add_argument(
        '-index_bg_path', '--index_bg_path',
        type=str, default="indices_cosine.txt",
        help='LIST_OF_INDICES filepath for the file containing per line a indice, e.g. "indices_cosine.txt" resulting from computeKnn.py'
    )

    parser.add_argument(
        '-out_path', '--out_path',
        type=str, default="output/",
        help='output path containing a folder named images and masks, e.g."output/" '
    )

    parser.add_argument(
        '-seed', '--seed',
        type=int, default=22,
        help='seed number for the pseudo-random computation'
    )

    parser.add_argument(
        '-size', '--size',
        type=int, default=10553,
        help='number of images in the dataset'
    )

    parser.add_argument(
        '-n_bgs', '--n_bgs',
        type=int, default=1,
        help='N_OF_BACKGROUNDS'
    )

    parser.add_argument(
        '-n_ops', '--n_ops',
        type=int, default=1,
        help='N_OF_OPS'
    )

    parser.add_argument(
        '-dataset_file', '--dataset_file',
        type=str, default="/home/bakrinski/workspace/datasets/duts/training/paths.txt",
        help='DATASET_FILE'
    )

    return parser.parse_args()

# SETTINGS

# CALL PARSER
args = parse_args()
#

OBJ_FOLDER_IMG = args.obj_path
OBJ_FOLDER_MASK = args.obj_mask_path
BG_FOLDER_IMG = args.bg_path
OUTPUT_FOLDER_IMG = args.out_path + "images/"
OUTPUT_FOLDER_MASK = args.out_path + "masks/"
LIST_OF_N_OBJECTS = args.index_obj_path
N_OBJECT = args.size
N_OF_BACKGROUNDS = args.n_bgs
N_OF_OPS = args.n_ops
LIST_OF_INDICES = args.index_bg_path
DATASET_FILE = args.dataset_file

kernelErode = np.ones((3, 3), np.uint8)

maxH = 512
maxW = 512

random.seed(args.seed)
np.random.seed(args.seed)
# noise_scale = np.random.uniform(low=0.975, high=1.025, size=N_OBJECT)
# noise_x = np.random.uniform(low=0.0, high=1.0, size=N_OBJECT)
# noise_y = np.random.uniform(low=0.0, high=1.0, size=N_OBJECT)
#

def computeMaximalDistanceForTranslation(obj, obj_mask, bg, newYmax, newYmin, newXmax, newXmin, newOrigin, border, M):
    import torch
    from scipy.spatial import distance
    import skimage.feature as ft

    DEVICE = "cpu"
    NBINS = 64

    def getHistogramsTorch(img):

        img = Image.fromarray((img).astype(np.uint8))

        h, s, v = img.convert('HSV').split()
        h = np.array(h)
        s = np.array(s)
        v = np.array(v)

        histH = torch.histc(torch.from_numpy(h).float().to(
            DEVICE), bins=NBINS, min=0.0, max=255.0)
        histS = torch.histc(torch.from_numpy(s).float().to(
            DEVICE), bins=NBINS, min=0.0, max=255.0)
        histV = torch.histc(torch.from_numpy(v).float().to(
            DEVICE), bins=NBINS, min=0.0, max=255.0)

        imgGray = img.convert('L')

        # # settings for LBP
        # # radius = 3
        # # n_points = 8 * radius
        # # METHOD = 'uniform'
        # # lbp = ft.local_binary_pattern(imgGray, n_points, radius, METHOD)
        lbp = ft.local_binary_pattern(imgGray, 24, 3, 'uniform')

        histLBP = torch.histc(torch.from_numpy(lbp).float().to(
            DEVICE), bins=NBINS, min=0.0, max=255.0)

        hist = torch.stack((histH, histS, histV, histLBP))
        hist = hist.view(-1)

        # img.save("tmp/bgTmp"+str((time.time()))+".png")

        return hist.numpy()

    def getHistogramsWithMaskTorch(img, mask):

        img = Image.fromarray((img).astype(np.uint8))
        mask = Image.fromarray((mask).astype(np.uint8))

        if img.mode == "L":
            img = img.convert('RGB')

        if mask.mode != "L":
            mask = mask.convert('L')

        h, s, v = img.convert('HSV').split()
        h = np.array(h) * mask
        s = np.array(s) * mask
        v = np.array(v) * mask

        histH = torch.histc(torch.from_numpy(h).float().to(
            DEVICE), bins=NBINS, min=0.0, max=255.0)
        histS = torch.histc(torch.from_numpy(s).float().to(
            DEVICE), bins=NBINS, min=0.0, max=255.0)
        histV = torch.histc(torch.from_numpy(v).float().to(
            DEVICE), bins=NBINS, min=0.0, max=255.0)

        imgGray = img.convert('L')

        # # settings for LBP
        # # radius = 3
        # # n_points = 8 * radius
        # # METHOD = 'uniform'
        # # lbp = ft.local_binary_pattern(imgGray, n_points, radius, METHOD)
        lbp = ft.local_binary_pattern(imgGray, 24, 3, 'uniform')

        histLBP = torch.histc(torch.from_numpy(lbp).float().to(
            DEVICE), bins=NBINS, min=0.0, max=255.0)

        hist = torch.stack((histH, histS, histV, histLBP))
        hist = hist.view(-1)

        # imgGray.save("tmp/graytmp"+str(int(time.time()))+".png")
        # mask.save("tmp/masktmp"+str(int(time.time()))+".png")

        return hist.numpy()

    # print("-------------")
    fd = getHistogramsWithMaskTorch(obj, obj_mask)

    # imgtmp = Image.fromarray((bg).astype(np.uint8))
    # imgtmp.save("tmp/bgFullTmp"+str((time.time()))+".png")

    bg_h, bg_w, _ = bg.shape
    obj_h, obj_w, _ = obj.shape
    # print("bg.shape=", bg.shape)
    # print("obj.shape=", obj.shape)

    h = int(np.ceil(bg_h / obj_h))
    w = int(np.ceil(bg_w / obj_w))
    # print("h=", h)
    # print("w=", w)

    cell = np.zeros((h, w, NBINS * 4))
    # print(cell.shape)

    disMax = 0.0
    startXmax = 0.0
    endXmax = 0.0
    startYmax = 0.0
    endYmax = 0.0

    for i in range(0, h):
        for j in range(0, w):

            starty = min((i) * obj_h, bg_h - obj_h)
            endy = min((i + 1) * obj_h, bg_h)

            startx = min((j) * obj_w, bg_w - obj_w)
            endx = min((j + 1) * obj_w, bg_w)

            if((endx - startx) < obj_w):
                startx = bg_w - obj_w
                endx = bg_w
            if((endy - starty) < obj_h):
                starty = bg_h - obj_h
                endy = bg_h

            # print("starty=",starty,"endy=",endy,"startx=",startx,"endx=",endx)

            # print(getHistogramsTorch(bg[starty:endy,startx:endx]))
            cell[i][j] = getHistogramsTorch(bg[starty:endy, startx:endx])
            dis = distance.cosine(cell[i][j], fd)
            if(dis > disMax):
                disMax = dis
                startXmax = startx
                endXmax = endx
                startYmax = starty
                endYmax = endy
            # print(dis)

    # print("Max=", disMax, "imax=", iMax, "jmax=", jMax)
    # print("-------------")

    # Begin translate LOGIC
    # newOrigin[0]#y
    # newOrigin[1]#x

    posY = startYmax + newOrigin[0]
    posX = startXmax + newOrigin[1]

    # y
    M[1][2] += posY - newYmin
    # x
    M[0][2] += posX - newXmin

    # End translate LOGIC

    return M

def randomTranslateInside(newYmax, newYmin, newXmax, newXmin, newOrigin, border, M):
    noise_x = np.random.uniform(low=0.0, high=1.0)
    noise_y = np.random.uniform(low=0.0, high=1.0)
    # check if bbox can move in y
    if((newYmax - newYmin) < border[0]):
        # check the direction of free space
        if((newYmax) < newOrigin[0] + border[0]):
            if((newYmin) > newOrigin[0]):
                freeSpacePos = (newOrigin[0] + border[0]) - newYmax
                freeSpaceNeg = newYmin - newOrigin[0]

                luck = np.random.randint(low=0, high=2)
                if(luck == 0):
                    M[1][2] += np.floor(noise_y * freeSpacePos)
                else:
                    M[1][2] -= np.floor(noise_y * freeSpaceNeg)

            else:
                freeSpace = (newOrigin[0] + border[0]) - newYmax
                M[1][2] += np.floor(noise_y * freeSpace)
        else:
            if((newYmin) > newOrigin[0]):
                freeSpace = newYmin - newOrigin[0]
                M[1][2] -= np.floor(noise_y * freeSpace)

    if((newXmax - newXmin) < border[1]):
        # check the direction of free space
        if((newXmax) < newOrigin[1] + border[1]):
            if((newXmin) > newOrigin[1]):
                freeSpacePos = (newOrigin[1] + border[1]) - newXmax
                freeSpaceNeg = newXmin - newOrigin[1]

                luck = np.random.randint(low=0, high=2)
                if(luck == 0):
                    M[0][2] += np.floor(noise_x * freeSpacePos)
                else:
                    M[0][2] -= np.floor(noise_x * freeSpaceNeg)

            else:
                freeSpace = (newOrigin[1] + border[1]) - newXmax
                M[0][2] += np.floor(noise_x * freeSpace)
        else:
            if((newXmin) > newOrigin[1]):
                freeSpace = newXmin - newOrigin[1]
                M[0][2] -= np.floor(noise_x * freeSpace)
    return M

def geometricOp2(resizedImg, resizedMask, bgOriginalshape, op, globalIndex, bgOriginalImg):
    #######################################################
    diffH = int((resizedImg.shape[0] - bgOriginalshape[0]) / 2)
    diffW = int((resizedImg.shape[1] - bgOriginalshape[1]) / 2)
    ####
    # compute bounding box
    ymin, ymax, xmin, xmax = bbox(resizedMask)

    # compute obj size
    propX = (xmax - xmin)
    propY = (ymax - ymin)

    # compute area
    areaOBJ = propX * propY
    areaIMG = bgOriginalshape[0] * bgOriginalshape[1]

    # compute obj size in relation to bg
    prop = areaOBJ / areaIMG
    ###

    # scale definition
    op = globalIndex % 3

    if(op == 0):
        # beta = 0.05 * noise_scale[globalIndex]
        beta = np.random.uniform(low=0.075, high=0.10)
    if(op == 1):
        # beta = 0.15 * noise_scale[globalIndex]
        beta = np.random.uniform(low=0.10, high=0.20)
    if(op == 2):
        # beta = 0.25 * noise_scale[globalIndex]
        beta = np.random.uniform(low=0.20, high=0.30)
    # if(op == 3):
    #     beta = 0.35 * noise_scale[globalIndex]

    scale = np.sqrt((beta * areaIMG) / areaOBJ)
    ###

    diffx = ((xmax - xmin) / 2)
    diffy = ((ymax - ymin) / 2)
    centerx = xmin + diffx
    centery = ymin + diffy

    pts1 = np.float32([[xmin, ymin], [xmax, ymin], [xmin, ymax]])

    newXmin = centerx - diffx * scale
    newXmax = centerx + diffx * scale

    newYmin = centery - diffy * scale
    newYmax = centery + diffy * scale

    # LOGIC HERE
    newOrigin = [diffH, diffW]
    border = [bgOriginalshape[0], bgOriginalshape[1]]

    # check if the aspect of the object is the same as the bg
    obj_orientation = -1
    bg_orientation = -1

    if(diffx >= diffy):
        obj_orientation = 0
    else:
        obj_orientation = 1

    if(bgOriginalshape[1] >= bgOriginalshape[0]):
        bg_orientation = 0
    else:
        bg_orientation = 1

    # check if can fit
    if((newYmax - newYmin <= border[0])and(newXmax - newXmin <= border[1])):
        # ok then it can fit
        # but does it need translation?

        pts2 = np.float32(
            [[newXmin, newYmin], [newXmax, newYmin], [newXmin, newYmax]])

        M = cv2.getAffineTransform(pts1, pts2)

        # origin of object must be >= newOrigin
        if(newYmin <= newOrigin[0]):
            local_diff_y = newOrigin[0] - newYmin
            M[1][2] += (local_diff_y)

        if(newXmin <= newOrigin[1]):
            local_diff_x = newOrigin[1] - newXmin
            M[0][2] += (local_diff_x)

        # maxdim must be <= border with the correct origin
        if(newYmax >= (border[0] + newOrigin[0])):
            local_diff_y = newYmax - (border[0] + newOrigin[0])
            M[1][2] -= (local_diff_y)

        if(newXmax >= (border[1] + newOrigin[1])):
            local_diff_x = newXmax - (border[1] + newOrigin[1])
            M[0][2] -= (local_diff_x)

        newXmin = xmin * M[0][0] + ymin * M[0][1] + M[0][2]
        newXmax = xmax * M[0][0] + ymax * M[0][1] + M[0][2]

        newYmin = xmin * M[1][0] + ymin * M[1][1] + M[1][2]
        newYmax = xmax * M[1][0] + ymax * M[1][1] + M[1][2]

        newXminTmp = min(newXmin, newXmax)
        newXmaxTmp = max(newXmin, newXmax)

        newYminTmp = min(newYmin, newYmax)
        newYmaxTmp = max(newYmin, newYmax)

        newXmin = newXminTmp
        newXmax = newXmaxTmp

        newYmin = newYminTmp
        newYmax = newYmaxTmp

        # NEW
        # resizedImgTmp = cv2.warpAffine(
        #     resizedImg, M, (bgOriginalshape[1], bgOriginalshape[0]), flags=cv2.INTER_LINEAR)
        # resizedMaskTmp = cv2.warpAffine(
        #     resizedMask, M, (bgOriginalshape[1], bgOriginalshape[0]), flags=cv2.INTER_NEAREST)
        resizedImgTmp = cv2.warpAffine(
            resizedImg, M, (maxW, maxH), flags=cv2.INTER_LINEAR)
        resizedMaskTmp = cv2.warpAffine(
            resizedMask, M, (maxW, maxH), flags=cv2.INTER_NEAREST)
        # print([int(newYmin),int(newYmax),int(newXmin),int(newXmax)])
        M = computeMaximalDistanceForTranslation(resizedImgTmp[int(newYmin):int(newYmax), int(newXmin):int(newXmax)], resizedMaskTmp[int(
            newYmin):int(newYmax), int(newXmin):int(newXmax)], bgOriginalImg, newYmax, newYmin, newXmax, newXmin, newOrigin, border, M)
        #######

        # M = randomTranslateInside(
        #     newYmax, newYmin, newXmax, newXmin, newOrigin, border, M)

        resizedImg = cv2.warpAffine(
            resizedImg, M, (maxW, maxH), flags=cv2.INTER_LINEAR)
        resizedMask = cv2.warpAffine(
            resizedMask, M, (maxW, maxH), flags=cv2.INTER_NEAREST)
        # resizedImg = cv2.warpAffine(
        #     resizedImg, M, (bgOriginalshape[1], bgOriginalshape[0]), flags=cv2.INTER_LINEAR)
        # resizedMask = cv2.warpAffine(
        #     resizedMask, M, (bgOriginalshape[1], bgOriginalshape[0]), flags=cv2.INTER_NEAREST)
    else:
        # it cannot fit
        # resize
        if(obj_orientation == bg_orientation):
            print("same orientatation")
            print("scale=",np.sqrt((beta * areaIMG) / areaOBJ))

            # scale must consider translation
            scale = min((border[0]) / (ymax - ymin),
                        (border[1]) / (xmax - xmin))
            scale = scale*0.5
            #
            newXmin = centerx - diffx * scale
            newXmax = centerx + diffx * scale

            newYmin = centery - diffy * scale
            newYmax = centery + diffy * scale

            pts2 = np.float32(
                [[newXmin, newYmin], [newXmax, newYmin], [newXmin, newYmax]])

            M = cv2.getAffineTransform(pts1, pts2)

            # origin of object must be >= newOrigin
            if(newYmin <= newOrigin[0]):
                local_diff_y = newOrigin[0] - newYmin
                M[1][2] += (local_diff_y)

            if(newXmin <= newOrigin[1]):
                local_diff_x = newOrigin[1] - newXmin
                M[0][2] += (local_diff_x)

            # maxdim must be <= border with the correct origin
            if(newYmax >= (border[0] + newOrigin[0])):
                local_diff_y = newYmax - (border[0] + newOrigin[0])
                M[1][2] -= (local_diff_y)

            if(newXmax >= (border[1] + newOrigin[1])):
                local_diff_x = newXmax - (border[1] + newOrigin[1])
                M[0][2] -= (local_diff_x)

            newXmin = xmin * M[0][0] + ymin * M[0][1] + M[0][2]
            newXmax = xmax * M[0][0] + ymax * M[0][1] + M[0][2]

            newYmin = xmin * M[1][0] + ymin * M[1][1] + M[1][2]
            newYmax = xmax * M[1][0] + ymax * M[1][1] + M[1][2]

            newXminTmp = min(newXmin, newXmax)
            newXmaxTmp = max(newXmin, newXmax)

            newYminTmp = min(newYmin, newYmax)
            newYmaxTmp = max(newYmin, newYmax)

            newXmin = newXminTmp
            newXmax = newXmaxTmp

            newYmin = newYminTmp
            newYmax = newYmaxTmp

            # NEW
            # resizedImgTmp = cv2.warpAffine(
            #     resizedImg, M, (bgOriginalshape[1], bgOriginalshape[0]), flags=cv2.INTER_LINEAR)
            # resizedMaskTmp = cv2.warpAffine(
            #     resizedMask, M, (bgOriginalshape[1], bgOriginalshape[0]), flags=cv2.INTER_NEAREST)
            resizedImgTmp = cv2.warpAffine(
                resizedImg, M, (maxW, maxH), flags=cv2.INTER_LINEAR)
            resizedMaskTmp = cv2.warpAffine(
                resizedMask, M, (maxW, maxH), flags=cv2.INTER_NEAREST)
            # print([int(newYmin),int(newYmax),int(newXmin),int(newXmax)])
            M = computeMaximalDistanceForTranslation(resizedImgTmp[int(newYmin):int(newYmax), int(newXmin):int(newXmax)], resizedMaskTmp[int(
                newYmin):int(newYmax), int(newXmin):int(newXmax)], bgOriginalImg, newYmax, newYmin, newXmax, newXmin, newOrigin, border, M)
            #######

            # M = randomTranslateInside(
            #     newYmax, newYmin, newXmax, newXmin, newOrigin, border, M)

            resizedImg = cv2.warpAffine(
                resizedImg, M, (maxW, maxH), flags=cv2.INTER_LINEAR)
            resizedMask = cv2.warpAffine(
                resizedMask, M, (maxW, maxH), flags=cv2.INTER_NEAREST)

        else:
            print("different orientatation")
            # print("scale=",np.sqrt((beta * areaIMG) / areaOBJ))
            # check if a rotated obj fits

            idxmod = np.random.randint(low=0, high=2)
            if(idxmod == 0):
                degrees = -90
            if(idxmod == 1):
                degrees = 90

            M = cv2.getRotationMatrix2D(((maxW / 2), (maxH / 2)), degrees, 1)

            newXmin = xmin * M[0][0] + ymin * M[0][1] + M[0][2]
            newXmax = xmax * M[0][0] + ymax * M[0][1] + M[0][2]

            newYmin = xmin * M[1][0] + ymin * M[1][1] + M[1][2]
            newYmax = xmax * M[1][0] + ymax * M[1][1] + M[1][2]

            newXminTmp = min(newXmin, newXmax)
            newXmaxTmp = max(newXmin, newXmax)

            newYminTmp = min(newYmin, newYmax)
            newYmaxTmp = max(newYmin, newYmax)

            newXmin = newXminTmp
            newXmax = newXmaxTmp

            newYmin = newYminTmp
            newYmax = newYmaxTmp

            # scale must consider translation
            scale = min((border[0]) / (newYmax - newYmin),
                        (border[1]) / (newXmax - newXmin))
            scale=scale*0.5
            #

            M[0][0] *= scale
            M[0][1] *= scale
            M[1][0] *= scale
            M[1][1] *= scale

            newXmin = xmin * M[0][0] + ymin * M[0][1] + M[0][2]
            newXmax = xmax * M[0][0] + ymax * M[0][1] + M[0][2]

            newYmin = xmin * M[1][0] + ymin * M[1][1] + M[1][2]
            newYmax = xmax * M[1][0] + ymax * M[1][1] + M[1][2]

            newXminTmp = min(newXmin, newXmax)
            newXmaxTmp = max(newXmin, newXmax)

            newYminTmp = min(newYmin, newYmax)
            newYmaxTmp = max(newYmin, newYmax)

            newXmin = newXminTmp
            newXmax = newXmaxTmp

            newYmin = newYminTmp
            newYmax = newYmaxTmp

            # origin of object must be >= newOrigin
            if(newYmin <= newOrigin[0]):
                local_diff_y = newOrigin[0] - newYmin
                M[1][2] += (local_diff_y)

            if(newXmin <= newOrigin[1]):
                local_diff_x = newOrigin[1] - newXmin
                M[0][2] += (local_diff_x)

            # maxdim must be <= border with the correct origin
            if(newYmax >= (border[0] + newOrigin[0])):
                local_diff_y = newYmax - (border[0] + newOrigin[0])
                M[1][2] -= (local_diff_y)

            if(newXmax >= (border[1] + newOrigin[1])):
                local_diff_x = newXmax - (border[1] + newOrigin[1])
                M[0][2] -= (local_diff_x)

            newXmin = xmin * M[0][0] + ymin * M[0][1] + M[0][2]
            newXmax = xmax * M[0][0] + ymax * M[0][1] + M[0][2]

            newYmin = xmin * M[1][0] + ymin * M[1][1] + M[1][2]
            newYmax = xmax * M[1][0] + ymax * M[1][1] + M[1][2]

            newXminTmp = min(newXmin, newXmax)
            newXmaxTmp = max(newXmin, newXmax)

            newYminTmp = min(newYmin, newYmax)
            newYmaxTmp = max(newYmin, newYmax)

            newXmin = newXminTmp
            newXmax = newXmaxTmp

            newYmin = newYminTmp
            newYmax = newYmaxTmp

            # NEW
            # resizedImgTmp = cv2.warpAffine(
            #     resizedImg, M, (bgOriginalshape[1], bgOriginalshape[0]), flags=cv2.INTER_LINEAR)
            # resizedMaskTmp = cv2.warpAffine(
            #     resizedMask, M, (bgOriginalshape[1], bgOriginalshape[0]), flags=cv2.INTER_NEAREST)
            resizedImgTmp = cv2.warpAffine(
                resizedImg, M, (maxW, maxH), flags=cv2.INTER_LINEAR)
            resizedMaskTmp = cv2.warpAffine(
                resizedMask, M, (maxW, maxH), flags=cv2.INTER_NEAREST)
            # print([int(newYmin),int(newYmax),int(newXmin),int(newXmax)])
            M = computeMaximalDistanceForTranslation(resizedImgTmp[int(newYmin):int(newYmax), int(newXmin):int(newXmax)], resizedMaskTmp[int(
                newYmin):int(newYmax), int(newXmin):int(newXmax)], bgOriginalImg, newYmax, newYmin, newXmax, newXmin, newOrigin, border, M)
            #######

            # M = randomTranslateInside(
            #     newYmax, newYmin, newXmax, newXmin, newOrigin, border, M)

            resizedImg = cv2.warpAffine(
                resizedImg, M, (maxW, maxH), flags=cv2.INTER_LINEAR)
            resizedMask = cv2.warpAffine(
                resizedMask, M, (maxW, maxH), flags=cv2.INTER_NEAREST)

        # print("ABORT! it cannot fit")
        # exit(1)
    ####
    # cv2.rectangle(resizedMask, (int(newXmin), int(newYmin)),
    #               (int(newXmax), int(newYmax)), (255, 255, 255), 1)
    #######################################################
    return resizedImg, resizedMask

def bbox(img):
    rows = np.any(img, axis=1)
    cols = np.any(img, axis=0)
    rmin, rmax = np.where(rows)[0][[0, -1]]
    cmin, cmax = np.where(cols)[0][[0, -1]]

    return rmin, rmax, cmin, cmax

def resize_with_pad(image, height, width):

    def get_padding_size(image, height, width):
        # h, w, _ = image.shape
        h = image.shape[0]
        w = image.shape[1]
        # print("h=",h,"w=",w)
        top, bottom, left, right = (0, 0, 0, 0)
        if h < height:
            dh = height - h
            top = dh // 2
            bottom = dh - top
        if w < width:
            dw = width - w
            left = dw // 2
            right = dw - left
        else:
            pass
        return top, bottom, left, right

    top, bottom, left, right = get_padding_size(image, height, width)
    BLACK = [0, 0, 0]
    constant = cv2.copyMakeBorder(
        image, top, bottom, left, right, cv2.BORDER_CONSTANT, value=BLACK)

    return constant

def resizeToOrgImg(bgOriginalshape, new):
    if(bgOriginalshape[0] < new.shape[0]):
        diffH = int((new.shape[0] - bgOriginalshape[0]) / 2)
        new = new[diffH:bgOriginalshape[0] + diffH, :, :]

    if(bgOriginalshape[1] < new.shape[1]):
        diffW = int((new.shape[1] - bgOriginalshape[1]) / 2)
        new = new[:, diffW:bgOriginalshape[1] + diffW, :]

    return new

def resizeToOrg(bgOriginalshape, new, newMask):
    if(bgOriginalshape[0] < new.shape[0]):
        diffH = int((new.shape[0] - bgOriginalshape[0]) / 2)
        new = new[diffH:bgOriginalshape[0] + diffH, :, :]
        newMask = newMask[diffH:bgOriginalshape[0] + diffH, :, :]

    if(bgOriginalshape[1] < new.shape[1]):
        diffW = int((new.shape[1] - bgOriginalshape[1]) / 2)
        new = new[:, diffW:bgOriginalshape[1] + diffW, :]
        newMask = newMask[:, diffW:bgOriginalshape[1] + diffW, :]

    return new, newMask

def loadResizedBG(index, datasetList):
    # bgName = "MSRA10K_image_{:06d}.png".format(index)
    bgName = datasetList[index][0]
    #special case
    bgName= bgName.replace(".jpg","_inpaint.jpg")

    bgFile = Image.open(BG_FOLDER_IMG + bgName)

    if bgFile.mode != "RGB":
        bgFile = bgFile.convert('RGB')

    bg = np.array(bgFile)
    bgOriginalshape = bg.shape
    resizedBg = resize_with_pad(bg, height=maxH, width=maxW)
    # resizedBg = bg
    bgFile.close()
    return resizedBg, bgOriginalshape

def load_list(file):
    mylist = []
    with open(file) as f:
        line = f.readline()
        while True:
            if not line:
                break
            line = line.split(" ")
            line[1] = line[1].replace("\n", "")

            tmp = line[0].split("/")
            line[0] = tmp[len(tmp) - 1]

            tmp = line[1].split("/")
            line[1] = tmp[len(tmp) - 1]

            mylist.append([line[0], line[1]])
            line = f.readline()
    return mylist

def main(op, multipleBgs, outPath):

    # read LIST_OF_N_OBJECTS
    arrOBJ = np.zeros(N_OBJECT, np.int)
    f = open(LIST_OF_N_OBJECTS, "r")
    for i in range(0, N_OBJECT):
        line = f.readline()
        args = line.split(" ")
        arrOBJ[i] = int(args[0])
    f.close()
    ###

    # read LIST_OF_N_OBJECTS
    arrBG = np.zeros((N_OBJECT, N_OF_BACKGROUNDS), np.int)
    f = open(LIST_OF_INDICES, "r")
    for i in range(0, N_OBJECT):
        line = f.readline()
        if line == '\n':
            arrOBJ[i] = -1
        else:
            args = line.split(" ")
            for bgindex in range(0, N_OF_BACKGROUNDS):
                arrBG[i][bgindex] = int(args[bgindex])
    f.close()
    ###

    # read DATASET_FILE
    datasetFile = DATASET_FILE
    datasetList = load_list(datasetFile)
    ###

    realI = 0

    for i in range(0, N_OBJECT, 1):
        if(arrOBJ[i] != -1):
            # imgName = "MSRA10K_image_{:06d}.jpg".format(arrOBJ[i])
            imgName = datasetList[i][0]
            imFile = Image.open(OBJ_FOLDER_IMG + imgName)

            if imFile.mode != "RGB":
                imFile = imFile.convert('RGB')

            img = np.array(imFile)

            maskName = datasetList[i][1]
            # maskName = imgName.replace(".jpg", ".png")
            # maskName = maskName.replace("image", "mask")

            maskFile = Image.open(OBJ_FOLDER_MASK + maskName)

            if maskFile.mode != "L":
                maskFile = maskFile.convert('L')

            mask = np.array(maskFile)

            mask = np.tile(mask[:, :, None], [1, 1, 3])

            resizedImg = resize_with_pad(img, height=maxH, width=maxW)
            resizedMask = resize_with_pad(mask, height=maxH, width=maxW)

            # resizedImg = img
            # resizedMask = mask

            imFile.close()
            maskFile.close()
            # print(stamp)

            resizedImgArr = [[None] * (N_OF_OPS)] * N_OF_BACKGROUNDS
            resizedMaskArr = [[None] * (N_OF_OPS)] * N_OF_BACKGROUNDS
            # print(resizedImgArr)

            resizedBg = [None] * (N_OF_BACKGROUNDS)
            bgOriginalshape = [None] * (N_OF_BACKGROUNDS)

            blur = [[None] * (N_OF_OPS)] * N_OF_BACKGROUNDS
            inv_blur = [[None] * (N_OF_OPS)] * N_OF_BACKGROUNDS
            new = [[None] * (N_OF_OPS)] * N_OF_BACKGROUNDS
            result = [[None] * (N_OF_OPS)] * N_OF_BACKGROUNDS
            resizedMaskFinal = [[None] * (N_OF_OPS)] * N_OF_BACKGROUNDS

            for bgindex in range(0, N_OF_BACKGROUNDS):
                resizedBg[bgindex], bgOriginalshape[bgindex] = loadResizedBG(
                    arrBG[i][bgindex], datasetList)

                # resizedImg = resize_with_pad(img, height=bgOriginalshape[bgindex][0], width=bgOriginalshape[bgindex][1])
                # resizedMask = resize_with_pad(mask, height=bgOriginalshape[bgindex][0], width=bgOriginalshape[bgindex][1])
                # print("bgs=",bgOriginalshape[bgindex][0],bgOriginalshape[bgindex][1])
                # print("resizedImg.shape=",resizedImg.shape)
                # print("resizedMask.shape=",resizedMask.shape)

                # calcule ops per bgs
                for opindex in range(0, N_OF_OPS):
                    globalIndex = (
                        ((realI * N_OF_BACKGROUNDS) + bgindex) * N_OF_OPS) + opindex
                    # print(globalIndex)
                    resizedImgArr[bgindex][opindex], resizedMaskArr[bgindex][opindex] = geometricOp2(
                        resizedImg, resizedMask, bgOriginalshape[bgindex], opindex, globalIndex, resizeToOrgImg(bgOriginalshape[bgindex], resizedBg[bgindex]))

                    # internalLoop
                    # BEGIN Smooth border copy

                    # disable border blur
                    # resizedMaskTmp = cv2.dilate(
                    #     resizedMaskArr[bgindex][opindex], kernelErode, iterations=1)
                    resizedMaskTmp = resizedMaskArr[bgindex][opindex]
                    blur[bgindex][opindex] = cv2.blur(resizedMaskTmp, (3, 3))

                    # blur[bgindex][opindex] = resizedMaskArr[bgindex][opindex]
                    blur[bgindex][opindex] = (
                        blur[bgindex][opindex] / 255) * 0.95
                    inv_blur[bgindex][opindex] = 1 - blur[bgindex][opindex]

                    new[bgindex][opindex] = blur[bgindex][opindex] * resizedImgArr[bgindex][opindex] + \
                        inv_blur[bgindex][opindex] * resizedBg[bgindex]
                    # END Smooth border copy

                    new[bgindex][opindex], resizedMaskArr[bgindex][opindex] = resizeToOrg(
                        bgOriginalshape[bgindex], new[bgindex][opindex], resizedMaskArr[bgindex][opindex])

                    #########################################################

                    result[bgindex][opindex] = Image.fromarray(
                        (new[bgindex][opindex]).astype(np.uint8))

                    resizedMaskFinal[bgindex][opindex] = Image.fromarray(
                        (resizedMaskArr[bgindex][opindex]).astype(np.uint8))

                    stamp = "{:06d}_{:06d}_{:03d}".format(
                        arrOBJ[i], arrBG[i][bgindex], opindex)

                    result[bgindex][opindex].save(OUTPUT_FOLDER_IMG +
                                                  "duts_image_" + stamp+ ".jpg")
                    resizedMaskFinal[bgindex][opindex].save(OUTPUT_FOLDER_MASK
                                                            + "duts_image_" + stamp+ ".png")

                    print(stamp)
                    #########################################################

            realI += 1

if __name__ == '__main__':

    if(args.n_bgs > 1):
        main(0, True, args.out_path)
    else:
        main(0, False, args.out_path)

max_stars_repo_name: VRI-UFPR/IDA

max_issues_repo_name: VRI-UFPR/IDA

max_forks_repo_name: VRI-UFPR/IDA


text:
import os
import random
import scipy
import math
import HC as hc
import SA as sa
import MCMCMC as mc3
import GA as ga

def DistanceMatrix(coords):
    '''Take a set of coordinates and calculate a matrix of
    Euclidean distances between the points.'''

    # we can assume that xs and ys are the same length
    stops = coords.shape[1]
    distMat = scipy.zeros((stops, stops))

    # this will be symmetric, so we only need to calculate
    # the upper triangular
    for i in range(stops):
        for j in range(i + 1, stops):
            xdist = coords[0, i] - coords[0, j]
            ydist = coords[1, i] - coords[1, j]
            distMat[i, j] = math.sqrt(xdist**2 + ydist**2)
    # add the transpose to make it symmetric
    distMat = distMat + distMat.transpose()
    return distMat

def GenerateMap(stops, fname=None, seed=None):
    '''Generate a map with "stops" stops for the salesman
    to traverse. Write coordinates to file if "fname" is
    specified. Return the distance matrix for all coordinates.'''
    random.seed(seed)
    # randomly place stop coordinates in the unit square
    xs = [random.uniform(0, 1) for x in range(stops)]
    ys = [random.uniform(0, 1) for x in range(stops)]
    coords = scipy.array([xs, ys])
    # calculate matrix of distances
    distMat = DistanceMatrix(coords)
    if fname is not None:
        scipy.savetxt(fname, coords)

    return distMat

def TSP(stops, Alg, steps, param, seed=None,
        coordfile='xycoords.txt'):
    '''A wrapper function that attempts to optimize the traveling
    salesperson problem using a specified algorithm. If coordfile
    exists, a preexisting set of coordinates will be used. Otherwise,
    a new set of "stops" coordinates will be generated for the person to
    traverse, and will be written to the specified file.'''

    # Create the distance matrix, which will be used to calculate
    # the fitness of a given path
    if os.path.isfile(coordfile):
        coords = scipy.genfromtxt(coordfile)
        distMat = DistanceMatrix(coords)
    else:
        distMat = GenerateMap(stops, fname=coordfile, seed=seed)

    if Alg == 'HC':
        # param is the number of solutions to try per step
        bestSol, fitHistory = hc.HillClimber(steps, param, distMat, seed)
    elif Alg == 'SA':
        # param is a placeholder
        bestSol, fitHistory = sa.SimulatedAnnealing(
            steps, param, distMat, seed)
    elif Alg == 'MC3':
        # param is the number of chains
        bestSol, fitHistory = mc3.MCMCMC(steps, param, distMat, seed)
    elif Alg == 'GA':
        # param is the population size
        bestSol, fitHistory = ga.GeneticAlgorithm(steps, param, distMat, seed)
    else:
        raise ValueError('Algorithm must be "HC", "SA", "MC3", or "GA".')

    outfname = coordfile + '-' + Alg + '-' + \
        str(steps) + '-' + str(param) + '.txt'
    scipy.savetxt(outfname, scipy.array(bestSol), fmt='%i')
    return bestSol, fitHistory

ext: py

max_stars_repo_name: mendiang/GoodEnoughAlgs

max_issues_repo_name: mendiang/GoodEnoughAlgs

max_forks_repo_name: mendiang/GoodEnoughAlgs

reason: import scipy


text:
import numpy as np
import os
from scipy.interpolate import interp1d
import configparser as cp
from colossus.halo.mass_defs import changeMassDefinition
from colossus.cosmology.cosmology import setCosmology
from colossus.halo.concentration import concentration
from scipy.special import erfc
from astropy.cosmology import FlatLambdaCDM

def mass_conversion(m200m, redshift, cosmology, mass_is_log=True):
    """Convert m200m to m500c

    Parameters
    ----------
    m200m : array_like
        Halo mass(es) calculated in a radius such that the mean density is 200 times
        the mean density of the universe.
    redshift : float
        Redshift of the halos used to calculate the concentration and perform
        the mass conversion
    cosmology : dict
        Cosmology parameters being used
    mass_is_log : bool
        Flag to tell the script if it should intake and output ln(M) or M

    Returns
    -------
    output : array_like
        Halo mass(es) calculated in a radius such that the mean density is 500 times
        the critical density of the universe.

    Notes
    -----
    We assume that every halo is at the same redshift.
    """
    setCosmology('myCosmo', cosmology)

    if mass_is_log:
        m200m = np.exp(m200m)

    m500c = changeMassDefinition(m200m, concentration(m200m, '200m', redshift),
                                 redshift, '200m', '500c')[0]

    if mass_is_log:
        m500c = np.log(m500c)

    return m500c

def load_simple_projection_model(addpath=''):
    """Read in the simple projection model

    The model is Gaussian and this function reads in the variance in the observed
    richness given the true richness.

    Parameters
    ----------
    addpath : str
        The function looks in a specific directory for the file. This parameter
        prepends path onto the default location.

    Returns
    -------
    lamlist : array_like
        A grid in true richness
    varlist : array_like
        The variance in observed richness given true richness on the grid given
        by lamlist
    """
    lamlist, varlist = np.loadtxt('{}inputs/p_lam_obs_z0.2.txt'.format(addpath),
                                  unpack=True, skiprows=1)
    return lamlist, varlist

def load_full_projection_model(addpath=''):
    """Read in the full projection model

    Parameters
    ----------
    addpath : str
        The function looks in a specific directory for the file. This parameter
        prepends path onto the default location.

    Returns
    -------
    projection_model : dict
        Dictionary containing lists to the 5 projection model parameters and true richness
        on a grid
    """
    input_model = np.loadtxt('{}projection_model/extended_projection_model.dat'.format(addpath))
    lam, fprj, fmsk, mu, sigma, tau = input_model
    return {'tau': tau, 'mu': mu, 'sigma': sigma, 'fprj': fprj, 'fmsk': fmsk, 'lam': lam}

def load_projection_model_interpolations(addpath=''):
    """Interpolate the full projection model

    Parameters
    ----------
    addpath : str
        The function looks in a specific directory for the file. This parameter
        prepends path onto the default location.

    Returns
    -------
    projection_model : dict
        Dictionary containing cubic splines for each of the 5 projection model parameters
    """
    lss_model = load_full_projection_model(addpath)

    fprj = interp1d(lss_model['lam'], lss_model['fprj'], kind='cubic')
    fmsk = interp1d(lss_model['lam'], lss_model['fmsk'], kind='cubic')
    mu = interp1d(lss_model['lam'], lss_model['mu'], kind='cubic')
    sigma = interp1d(lss_model['lam'], lss_model['sigma'], kind='cubic')
    tau = interp1d(lss_model['lam'], lss_model['tau'], kind='cubic')

    projection_model = {'tau': tau, 'mu': mu, 'sigma': sigma, 'fprj': fprj, 'fmsk': fmsk}
    return projection_model

def evaluate_full_projection_params(model, true_richness):
    """Given a true richness, determine the values of the projection model parameters

    Parameters
    ----------
    model : dict
        Dictionary of cubic splines for each of the projection model parameters
    true_richness : float
        The true richess that you wish to evaluate the projection model at

    Returns
    -------
    output : dict
        Dictionary of the model parameters evaluated at the requested true richness.
    """
    fprj = model['fprj'](true_richness)
    fmsk = model['fmsk'](true_richness)
    tau = model['tau'](true_richness)
    sigma = model['sigma'](true_richness)
    mu = model['mu'](true_richness)

    return {'fprj': fprj, 'fmsk': fmsk, 'mu': mu, 'sigma': sigma,
            'tau': tau, 'lamtrue': true_richness}

def evaluate_full_projection_model(lamobs, model):
    """Calculate P(observed_richness|true_richness)

    Parameters
    ----------
    lamobs : float
        Observed richness
    model : dict
        Model parameters and true richness at the true richness

    Returns
    -------
    output : float
        P(observed_richness|true_richness)
    """
    # Model params
    fprj, fmsk, lamtrue = model['fprj'], model['fmsk'], model['lamtrue']
    mu, sigma, tau = model['mu'], model['sigma'], model['tau']

    # Some stuff to make the calculation easier
    sig2 = sigma*sigma
    aa = np.exp(0.5*tau*(2.*mu + tau*sig2 - 2.*lamobs))
    root2siginv = 1./(np.sqrt(2.)*sigma)

    # The 4 terms in the model
    t1 = (1.-fmsk)*(1.-fprj)*np.exp(-(lamobs - mu)*(lamobs - mu)/(2*sig2))/np.sqrt(2.*np.pi*sig2)
    t2 = 0.5*((1.-fmsk)*fprj*tau + fmsk*fprj/lamtrue)*aa*erfc((mu + tau*sig2 - lamobs)*root2siginv)
    t3 = (fmsk*0.5/lamtrue)*(erfc((mu - lamobs - lamtrue)*root2siginv)
                             - erfc((mu - lamobs)*root2siginv))
    t4 = (fmsk*fprj*0.5/lamtrue)*np.exp(-1.*tau*lamtrue)*aa*erfc((mu + tau*sig2 - lamtrue
                                                                 - lamobs)*root2siginv)
    return t1+t2+t3-t4

class SurveyInfo(object):
    """Hold all of the survey details

    Attributes
    ----------
    area : float
        Survey area in sq deg
    zmin : float
        Minimum survey redshift
    zmax : float
        Maximum survey redshift
    cosmo : dict
        Cosmological model
    c : float
        The speed of light
    survey_volume : float
        3d volume of survey in Mpc^3
    solid_angle : float
        Survey area in steradians
    """
    def __init__(self, cfgin, cosmology):
        self.area = cfgin['Survey Info'].getfloat('survey_area')
        self.zmin = cfgin['Survey Info'].getfloat('zlo')
        self.zmax = cfgin['Survey Info'].getfloat('zhi')
        self.cosmo = cosmology
        self.c = 300000.
        self.solid_angle = self.area*np.pi**2/(180.**2)
        self._calc_comoving_volume()

    def _calc_comoving_volume(self):
        cosmo = FlatLambdaCDM(H0=self.cosmo['H0'], Om0=self.cosmo['Om0'])
        fullsky_volume = cosmo.comoving_volume(self.zmax)-cosmo.comoving_volume(self.zmin)
        survey_volume = fullsky_volume.value*(self.solid_angle/(4.*np.pi))
        self.survey_volume = survey_volume

    def reset_survey_area(self, area):
        self.area = area
        self.solid_angle = self.area*np.pi**2/(180.**2)
        self._calc_comoving_volume()

class Cluster(object):
    """Hold all of the details for a single galaxy cluster

    Attributes
    ----------
    lam : float
        Richness
    lamE : float
        Uncertainty in richness
    mgas : float
        Gas mass
    mgasE : float
        Uncertainty in gas mass
    lnm200m : float
        ln(M200m)
    lnm500c : float
        ln(M500c)
    """
    def __init__(self, richness, richness_error, mgas, mgas_error, rescales):
        self.lam = richness/rescales['rich']
        self.lamE = richness_error/rescales['rich']
        self.mgas = mgas/rescales['mgas']
        self.mgasE = mgas_error/rescales['mgas']

    def include_mass(self, lnm200m, lnm500c):
        self.lnm200m = lnm200m
        self.lnm500c = lnm500c

    def __lt__(self, other):
        return self.lam < other.lam

    def __repr__(self):
        return "%f %f %f %f\n" % (self.lam, self.lamE, self.mgas, self.mgasE)

    def printout(self):
        return [1.0, self.lam, self.lamE, self.mgas, self.mgasE]

    def printout_extra(self):
        return [1.0, self.lam, self.lamE, self.mgas, self.mgasE, self.lnm200m, self.lnm500c]

class Realization:
    """Hold a single survey realization.

    The intended use here is to allow the analysis to be run on several survey
    realizations at once.

    Attributes
    ----------
    cllist : array_like
        List of Cluster objects
    lam_n : float
        The richness of the N richest object for our selection function
        analysis
    lam_n_obs : float
        Observational uncertainty in lam_n
    """
    def __init__(self, cllist, lam_n, lam_n_obs):
        self.lam_n = lam_n
        self.lam_n_obs = lam_n_obs
        self.cllist = cllist

def read_hmf_files(addpath, cfgin, gencut=-1.):
    """Load the halo mass function

    My halo mass function comes from a cosmosis modules and is precalculated.
    This could probably be rewritten to calculate it on the fly using pyccl or
    some other tool.

    Parameters
    ----------
    addpath : str
        If not storing this in the default location, prepend onto the path
    cfgin : ConfigParser
        Config details for the run
    gencut : float, optional
        Use a different mass cut if I am generating a mock
    """
    log10m_cut = cfgin['General'].getfloat('log10Mcut')
    zhmf = cfgin['Survey Info'].getfloat('zhmf')
    cosmology = cfgin['General']['cosmology_name']
    path = "{}inputs/hmf/hmf_{}_cosmology/mass_function/".format(addpath, cosmology)

    # Find the index of the HMF we want to work with
    zlist = np.loadtxt(path+'z.txt', skiprows=1)
    zind = np.where(np.abs(zlist-zhmf) < 0.00003)[0][0]
    print("Loading HMF at redshift %f" % zlist[zind])

    # Extract the appropriate HMF
    with open(path+'dndmdz.txt', 'r') as fin:
        lines = fin.readlines()[1:]
        hmf = [float(row.split(' ')[zind]) for row in lines]

    # Load masses and make mass cut
    mlist = np.loadtxt(path+'mass.txt', skiprows=1)
    if gencut > 0:
        mass_cut = np.where((mlist >= 10.**gencut) & (mlist <= 10.**16))
    else:
        mass_cut = np.where((mlist >= 10.**log10m_cut) & (mlist <= 10.**16))

    # Recast into dn/(dz dlnm)
    lnmlist = np.asarray([np.log(m) for m in mlist])
    dndlnmdz = np.asarray([m*dndmdz for m, dndmdz in zip(mlist, hmf)])

    return lnmlist[mass_cut], dndlnmdz[mass_cut]

def load_clusters(cfgin, note, rescales, altpath=''):
    """Load the list of clusters being used for analysis

    Parameters
    ----------
    cfgin : ConfigParser
        Configuration details for the run
    note : int
        Additional details for the input file names. Used to tell realizations apart
    rescales : array_like
        Quantities to rescale the cluster measurements to ~unity
    altpath : str
        If clusters are stored somewhere besides the default path

    Returns
    -------
    real_list : array_like
        Realization objects loaded in
    lnpivot : float
        The median cluster mass (or estimate of such) to use as a default mass pivot
    """
    simtype = cfgin['General']['simtype']
    catalog = cfgin['General']['synth_catalog']
    n_clusters = cfgin['SurveyInfo'].getint('Nclusters')

    if simtype == 'synth10':
        real_list = load_synth10(catalog, note, rescales)
        return real_list, 0.0

    # Build the filename based on the type of simulation
    if simtype == 'real':
        cluster_fname = 'inputs/richestClustersInRedmapper_MantzXray_MatteoErrors.dat'
    elif len(altpath) > 5:
        cluster_fname = altpath
    elif simtype == 'synth':
        cluster_fname = './catalogs/mock/{}/richest{}_{}.dat'.format(catalog, n_clusters, note)

    # Actually load the clusters
    cllist, lnpivot, lam_n, lam_n_obs = load_clusters_from_file(cluster_fname, rescales)
    real_list = [Realization(cllist, lam_n, lam_n_obs)]

    return real_list, lnpivot

def load_clusters_from_file(fname, rescales):
    """Load the clusters from a file into Cluster objects

    Parameters
    ----------
    fname : str
        Path to the cluster catalog
    rescales : array_like
        Quantities to rescale the cluster measurements to ~unity

    Returns
    -------
    cllist : array_like
        Like of Cluster objects
    lnpivot : float
        Default pivot, the median cluster mass
    lam_n : float
        Richness of the least rich cluster
    lam_n_obs : float
        Uncertainty in lam_n
    """
    # Load the clusters
    inputs = np.loadtxt(fname, skiprows=1, usecols=[1, 2, 3, 4])
    lamobs, lamobs_uncert, mgas, mgas_uncert = inputs
    cllist = [Cluster(lamobs[i], lamobs_uncert[i], mgas[i], mgas_uncert[i], rescales)
              for i in range(len(lamobs))]

    # Find the median pivot
    lnpivot = np.log(np.median(mgas)/0.3)

    # Find the least rich cluster
    lam_n, idx = min((lamN, idx) for (idx, lamN) in enumerate(lamobs))
    lam_n_obs = lamobs_uncert[idx]

    print('%i Clusters loaded.' % (len(cllist)))
    return cllist, lnpivot, lam_n, lam_n_obs

def load_synth10(catalog, note, rescales):
    """Load 10 realizations"""
    real_list = []
    for kk in range(10):
        cluster_fname = './catalogs/mock/%s/richest30_%i.dat' % (catalog, note*10+kk)
        print("Read cluster", cluster_fname)
        cllist, lnpivot, lam_n, lam_n_obs = load_clusters_from_file(cluster_fname, rescales)
        real_list.append(Realization(cllist, lam_n, lam_n_obs))
    return real_list

class PriorContainer(object):
    """Hold all of the prior information"""
    def __init__(self, cfgin):
        path = '/home/matthewkirby/cosmosis/modules/xray-likelihood/priors/'
        fname = cfgin['MCMC']['priorfile']
        priorsin = cp.ConfigParser()
        priorsin.read(os.path.join(path + fname))

        self.mg0 = priorsin['MgasPriors'].getfloat('mg0')
        self.mg0_uncert = priorsin['MgasPriors'].getfloat('mg0e')
        self.alphamg = priorsin['MgasPriors'].getfloat('alphamg')
        self.alphamg_uncert = priorsin['MgasPriors'].getfloat('alphamge')

        self.amp_mr_rel = priorsin['RichnessPriors'].getfloat('as')
        self.amp_uncert_mr_rel = priorsin['RichnessPriors'].getfloat('ase')
        self.slope_mr_rel = priorsin['RichnessPriors'].getfloat('gamma')
        self.slope_uncert_mr_rel = priorsin['RichnessPriors'].getfloat('gammae')

        self.quant = priorsin['RichnessPriors'].getfloat('quant')
        self.quant_uncert = priorsin['RichnessPriors'].getfloat('quante')
        self.alphalam = priorsin['RichnessPriors'].getfloat('alphalam')
        self.alphalam_uncert = priorsin['RichnessPriors'].getfloat('alphalame')

    def __repr__(self):
        return "mg0 %f +- %f\nalphamg %f +- %f\nAs %f +- %f\ngamma %f +- %f" % \
            (self.mg0, self.mg0_uncert, self.alphamg, self.alphamg_uncert, self.amp_mr_rel,
             self.amp_uncert_mr_rel, self.slope_mr_rel, self.slope_uncert_mr_rel)

def build_output_files(cfgin, note):
    """Build the paths to the output files

    Parameters
    ----------
    cfgin : ConfigParser
        Configuration details for the run
    note : int
        Realization specific path details

    Returns
    -------
    out_chain : str
        Output file name for the chain
    out_like : str
        Output file name for the likelihood
    """
    simtype = cfgin['General']['simtype']
    prior = cfgin['MCMC']['prior']
    if simtype == 'synth' or simtype == 'synth10':
        catalog = cfgin['General']['synth_catalog']
        out_chain = './outputs/{}/chains/mcmc{}_{}prior_richest{}.out'.format(
            catalog, note, prior, cfgin['SurveyInfo'].getint('nclusters'))
        out_like = './outputs/{}/chains/lnL{}_{}priors.out'.format(catalog, note, prior)

    elif simtype == 'real':
        out_chain = './outputs/real_runs/chains/mcmc_{}priors.out'.format(prior)
        out_like = './outputs/real_runs/chains/lnL_{}priors.out'.format(prior)

    else:
        raise TypeError("Simulation type {} is an invalid type.".format(simtype))

    return out_chain, out_like

def select_pivots(cfgin, lnpivot):
    """Select the mass pivots that will be used in the analysis

    Parameters
    ----------
    cfgin : ConfigParser
        Configuration details for the run
    lnpivot : float
        Default pivot to use

    Returns
    -------
    pivots : array_like
        The two pivots to use in the analysis
    """
    lam_pivot = cfgin['General'].getfloat('rich_pivot')
    mgas_pivot = cfgin['General'].getfloat('mgas_pivot')
    pivots = {'rich': 0.0, 'mgas': 0.0}

    if abs(lam_pivot+1.) <= 0.003:
        pivots['rich'] = lnpivot
    else:
        pivots['rich'] = np.log(lam_pivot)
    if abs(mgas_pivot+1.) <= 0.003:
        pivots['mgas'] = lnpivot
    else:
        pivots['mgas'] = np.log(mgas_pivot)

    return pivots

def load_truth(fname):
    """Load the true model"""
    model = cp.ConfigParser()
    model.read(fname)
    params = [model['Model'].getfloat('r'),
              model['Model'].getfloat('mgas0'),
              model['Model'].getfloat('alphamg'),
              model['Model'].getfloat('s0mg'),
              model['Model'].getfloat('lam0'),
              model['Model'].getfloat('alphalam'),
              model['Model'].getfloat('s0lam'),
              model['Model'].getfloat('dhmf')]

    return params

def find_mock_lamobs_max(simtype, addpath=''):
    """Find the maximum relevant observed richness as a function of true richness

    This is used as the upper bound on the integrals.
    """
    lamtrue_grid = np.loadtxt('{}projection_model/cdfs/lamtrue.dat'.format(addpath))
    lamobs_grid = np.loadtxt('{}projection_model/cdfs/lamobs.dat'.format(addpath))
    cdf_grid = np.loadtxt('{}projection_model/cdfs/cdflist.dat'.format(addpath))
    maxlamobs = []

    if simtype == 'real':
        return [lamtrue_grid, np.ones(len(lamtrue_grid))*1999.]

    for i in range(len(lamtrue_grid)):
        cdf = cdf_grid[i]
        idxs = np.array(np.where(cdf < 1.0)[0])
        lamobs_grid_cut = lamobs_grid[idxs]
        maxlamobs.append(lamobs_grid_cut[-1])

    return [lamtrue_grid, maxlamobs]

def load_cosmology(filename):
    """Load cosmology from a file into a dictionary

    Parameters
    ----------
    filename : str
        Full path and file name for the cosmology file

    Returns
    -------
    output : dict
        Cosmology parameters
    """
    cosmoin = cp.ConfigParser()
    cosmoin.read(filename)
    cosmology = {'name': cosmoin['Cosmology']['name'],
                 'flat': bool(cosmoin['Cosmology']['flat']),
                 'H0': cosmoin['Cosmology'].getfloat('h0'),
                 'Om0': cosmoin['Cosmology'].getfloat('omm'),
                 'Ob0': cosmoin['Cosmology'].getfloat('omb'),
                 'sigma8': cosmoin['Cosmology'].getfloat('sigma8'),
                 'ns': cosmoin['Cosmology'].getfloat('ns')}

    return cosmology

reason: import numpy,from scipy,from astropy


text:
import numpy

from chainer import cuda
from chainer import function
from chainer.utils import type_check

class Where(function.Function):

    """Choose elements depending on condition."""

    def check_type_forward(self, in_types):
        type_check.expect(in_types.size() == 3)
        c_type, x_type, y_type = in_types

        type_check.expect(
            c_type.dtype == numpy.bool_,
            x_type.dtype == y_type.dtype,
            x_type.shape == c_type.shape,
            y_type.shape == c_type.shape,
        )

    def forward(self, inputs):
        xp = cuda.get_array_module(*inputs)
        condition, x, y = inputs
        return xp.where(condition, x, y),

    def backward(self, inputs, grads):
        xp = cuda.get_array_module(*inputs)
        condition, x, y = inputs
        gx = xp.where(condition, grads[0], 0)
        gy = xp.where(condition, 0, grads[0])
        return None, gx, gy

def where(condition, x, y):
    """Choose elements depending on condition.

    This function choose values depending on a given ``condition``.
    All ``condition``, ``x``, and ``y`` must have the same shape.

    Args:
        condition (~chainer.Variable): Variable containing the condition.
            Only boolean array is permitted.
        x (~chainer.Variable): Variable chosen when ``condition`` is ``True``.
        y (~chainer.Variable): Variable chosen when ``condition`` is ``False``.

    Returns:
        ~chainer.Variable: Variable containing chosen values.
    """

    return Where()(condition, x, y)

reason: import numpy


text:
# coding=utf-8
# Copyright 2019 The Google Research Authors.
#
# 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.

"""Configurations for MAML training (Reinforcement Learning).

See maml_rl.py for usage examples.
An easy task to get started with is: RL_MINITAUR_POINT_CONFIG_CIRCLE.
"""

# b/128310658.
# pytype: skip-file

from __future__ import absolute_import
from __future__ import division
from __future__ import print_function

import functools
import random
import numpy as np
import tensorflow as tf

from norml import networks
from norml import policies
from norml.envs import cartpole_sensor_bias_env
from norml.envs import halfcheetah_motor_env
from norml.envs import move_point_env

def _early_termination_avg(rewards, num_steps, avg_reward):
  """Early termination based on average reward."""
  flat_reward = np.array(rewards).ravel()
  len_ok = len(flat_reward) >= num_steps
  val_ok = np.mean(flat_reward[-num_steps:]) >= avg_reward
  return len_ok and val_ok

MOVE_POINT_ROTATE_MAML = dict(
    random_seed=random.randint(0, 1000000),
    num_outer_iterations=1000,
    task_generator=functools.partial(
        move_point_env.MovePointEnv,
        start_pos=(0, 0),
        end_pos=(1, 0),
        goal_reached_distance=-1,
        trial_length=10),
    task_env_modifiers=[{
        '_action_rotation': i
    } for i in np.linspace(-np.pi, np.pi, 5000)],
    network_generator=networks.FullyConnectedNetworkGenerator(
        dim_input=2,
        dim_output=2,
        layer_sizes=(
            50,
            50,
        ),
        activation_fn=tf.nn.tanh),
    input_dims=2,
    pol_log_std_init=-3.,
    output_dims=2,
    reward_disc=0.9,
    learn_offset=False,
    policy=policies.GaussianPolicy,
    tasks_batch_size=10,
    num_inner_rollouts=25,
    outer_optimizer_algo=tf.train.AdamOptimizer,
    advantage_function='returns-values',
    whiten_values=False,
    always_full_rollouts=False,
    inner_lr_init=0.02,
    outer_lr_init=7e-3,
    outer_lr_decay=True,
    first_order=False,
    learn_inner_lr=True,
    learn_inner_lr_tensor=True,
    fixed_tasks=False,
    ppo=True,
    ppo_clip_value=0.2,
    max_num_batch_env=1000,
    max_rollout_len=10,
    log_every=10,
)

MOVE_POINT_ROTATE_MAML_OFFSET = MOVE_POINT_ROTATE_MAML.copy()
MOVE_POINT_ROTATE_MAML_OFFSET.update(
    learn_offset=True,
    inner_lr_init=0.1,
    outer_lr_init=3e-3,
    pol_log_std_init=-3.)

MOVE_POINT_ROTATE_MAML_LAF = MOVE_POINT_ROTATE_MAML.copy()
MOVE_POINT_ROTATE_MAML_LAF.update(
    learn_inner_lr=False,
    learn_inner_lr_tensor=False,
    learn_advantage_function_inner=True,
    advantage_generator=networks.FullyConnectedNetworkGenerator(
        dim_input=2 * 2 + 2,
        dim_output=1,
        layer_sizes=(
            50,
            50,
        ),
        activation_fn=tf.nn.tanh),
    inner_lr_init=0.7,
    outer_lr_init=6e-4,
    pol_log_std_init=-3.25)

MOVE_POINT_ROTATE_NORML = MOVE_POINT_ROTATE_MAML_LAF.copy()
MOVE_POINT_ROTATE_NORML.update(
    learn_offset=True,
    inner_lr_init=10.,
    outer_lr_init=6e-3,
    pol_log_std_init=-0.75)

MOVE_POINT_ROTATE_SPARSE_MAML = MOVE_POINT_ROTATE_MAML.copy()
MOVE_POINT_ROTATE_SPARSE_MAML.update(
    max_rollout_len=100,
    task_generator=functools.partial(
        move_point_env.MovePointEnv,
        start_pos=(0, 0),
        end_pos=(1, 0),
        goal_reached_distance=0.1,
        trial_length=100,
        sparse_reward=True),
    inner_lr_init=1e-4,
    outer_lr_init=2e-3,
    pol_log_std_init=-1.25)

MOVE_POINT_ROTATE_SPARSE_MAML_OFFSET = MOVE_POINT_ROTATE_MAML_OFFSET.copy()
MOVE_POINT_ROTATE_SPARSE_MAML_OFFSET.update(
    max_rollout_len=100,
    task_generator=functools.partial(
        move_point_env.MovePointEnv,
        start_pos=(0, 0),
        end_pos=(1, 0),
        goal_reached_distance=0.1,
        trial_length=100,
        sparse_reward=True),
    inner_lr_init=7.,
    outer_lr_init=2e-3,
    pol_log_std_init=-0.5)

MOVE_POINT_ROTATE_SPARSE_MAML_LAF = MOVE_POINT_ROTATE_MAML_LAF.copy()
MOVE_POINT_ROTATE_SPARSE_MAML_LAF.update(
    max_rollout_len=100,
    task_generator=functools.partial(
        move_point_env.MovePointEnv,
        start_pos=(0, 0),
        end_pos=(1, 0),
        goal_reached_distance=0.1,
        trial_length=100,
        sparse_reward=True),
    inner_lr_init=2e-5,
    outer_lr_init=1e-3,
    pol_log_std_init=-1.)

MOVE_POINT_ROTATE_SPARSE_NORML = MOVE_POINT_ROTATE_NORML.copy()
MOVE_POINT_ROTATE_SPARSE_NORML.update(
    max_rollout_len=100,
    task_generator=functools.partial(
        move_point_env.MovePointEnv,
        start_pos=(0, 0),
        end_pos=(1, 0),
        goal_reached_distance=0.1,
        trial_length=100,
        sparse_reward=True),
    inner_lr_init=9.,
    outer_lr_init=2.6e-3,
    pol_log_std_init=-0.6)

CARTPOLE_SENSOR_DR = dict(
    random_seed=random.randint(0, 1000000),
    num_outer_iterations=1000,
    task_generator=functools.partial(
        cartpole_sensor_bias_env.CartpoleSensorBiasEnv),
    task_env_modifiers=[{
        '_angle_observation_bias': theta
    } for theta in np.linspace(-np.pi / 18, np.pi / 18, 5000)],
    network_generator=networks.FullyConnectedNetworkGenerator(
        dim_input=4,
        dim_output=1,
        layer_sizes=(
            50,
            50,
        ),
        activation_fn=tf.nn.tanh),
    input_dims=4,
    pol_log_std_init=-4.0,
    output_dims=1,
    reward_disc=0.97,
    learn_offset=False,
    policy=policies.GaussianPolicy,
    tasks_batch_size=10,
    num_inner_rollouts=25,
    outer_optimizer_algo=tf.train.AdamOptimizer,
    advantage_function='returns-values',
    whiten_values=False,
    always_full_rollouts=False,
    inner_lr_init=0.,
    outer_lr_init=2e-4,
    outer_lr_decay=True,
    first_order=False,
    learn_inner_lr=False,
    learn_inner_lr_tensor=False,
    fixed_tasks=False,
    ppo=True,
    ppo_clip_value=0.2,
    max_num_batch_env=1000,
    max_rollout_len=500,
    log_every=10,
)

CARTPOLE_SENSOR_MAML = CARTPOLE_SENSOR_DR.copy()
CARTPOLE_SENSOR_MAML.update(
    learn_inner_lr=True,
    learn_inner_lr_tensor=True,
    inner_lr_init=1e-2,
    outer_lr_init=1e-2,
    pol_log_std_init=-0.5)

CARTPOLE_SENSOR_MAML_OFFSET = CARTPOLE_SENSOR_MAML.copy()
CARTPOLE_SENSOR_MAML_OFFSET.update(
    learn_offset=True,
    inner_lr_init=1e-1,
    outer_lr_init=4e-4,
    pol_log_std_init=-3.5)

CARTPOLE_SENSOR_MAML_LAF = CARTPOLE_SENSOR_MAML.copy()
CARTPOLE_SENSOR_MAML_LAF.update(
    learn_advantage_function_inner=True,
    advantage_generator=networks.FullyConnectedNetworkGenerator(
        dim_input=2 * 4 + 1,
        dim_output=1,
        layer_sizes=(
            50,
            50,
        ),
        activation_fn=tf.nn.tanh),
    inner_lr_init=1e-5,
    outer_lr_init=3e-3,
    pol_log_std_init=-0.5)

CARTPOLE_SENSOR_NORML = CARTPOLE_SENSOR_MAML_LAF.copy()
CARTPOLE_SENSOR_NORML.update(
    learn_offset=True,
    inner_lr_init=7e-4,
    outer_lr_init=3e-4,
    pol_log_std_init=-3.5)

HALFCHEETAH_MOTOR_DR = dict(
    random_seed=random.randint(0, 1000000),
    num_outer_iterations=1000,
    task_generator=functools.partial(halfcheetah_motor_env.HalfcheetahMotorEnv),
    task_env_modifiers=[{
        '_swap_action': True
    }, {
        '_swap_action': False
    }] * 2,
    network_generator=networks.FullyConnectedNetworkGenerator(
        dim_input=14,
        dim_output=6,
        layer_sizes=(100,),
        activation_fn=tf.identity),
    pol_log_std_init=-1.61,
    input_dims=14,
    output_dims=6,
    reward_disc=0.99,
    learn_offset=False,
    policy=policies.GaussianPolicy,
    tasks_batch_size=4,
    num_inner_rollouts=50,
    outer_optimizer_algo=tf.train.AdamOptimizer,
    advantage_function='returns-values',
    whiten_values=True,
    always_full_rollouts=False,  # also learn from failures?
    inner_lr_init=0.,
    outer_lr_init=0.0012,
    outer_lr_decay=True,
    first_order=False,
    learn_inner_lr=False,
    learn_inner_lr_tensor=False,
    fixed_tasks=False,
    log_every=10,
    ppo=True,
    ppo_clip_value=0.2,
    max_rollout_len=1000,
    max_num_batch_env=300,
)

HALFCHEETAH_MOTOR_MAML = HALFCHEETAH_MOTOR_DR.copy()
HALFCHEETAH_MOTOR_MAML.update(
    learn_inner_lr=True,
    learn_inner_lr_tensor=True,
    learn_offset=False,
    inner_lr_init=4.5e-4,
    outer_lr_init=8.7e-5,
    pol_log_std_init=-1.17)

HALFCHEETAH_MOTOR_MAML_OFFSET = HALFCHEETAH_MOTOR_MAML.copy()
HALFCHEETAH_MOTOR_MAML_OFFSET.update(
    learn_offset=True,
    inner_lr_init=5e-5,
    outer_lr_init=1.5e-4,
    pol_log_std_init=-0.8)

HALFCHEETAH_MOTOR_MAML_LAF = HALFCHEETAH_MOTOR_MAML.copy()
HALFCHEETAH_MOTOR_MAML_LAF.update(
    learn_advantage_function_inner=True,
    advantage_generator=networks.FullyConnectedNetworkGenerator(
        dim_input=14 * 2 + 6,
        dim_output=1,
        layer_sizes=(50,),
        activation_fn=tf.nn.relu),
    learn_inner_lr=False,
    learn_inner_lr_tensor=False,
    inner_lr_init=3e-5,
    outer_lr_init=5e-4,
    pol_log_std_init=-1.5)

HALFCHEETAH_MOTOR_NORML = HALFCHEETAH_MOTOR_MAML_LAF.copy()
HALFCHEETAH_MOTOR_NORML.update(
    learn_offset=True,
    inner_lr_init=3e-5,
    outer_lr_init=5e-4,
    pol_log_std_init=-1.9)

max_issues_repo_name: RubensZimbres/google-research

max_forks_repo_name: RubensZimbres/google-research


text:
from .. import arr
import numpy as np
from math import factorial

cfg = {
    'name': 'EEG - Frequency Analysis',
    'group': 'Medical algorithms',
    'desc': 'Frequency Analysis of EEG.',
    'reference': '',
    'overlap': 58,
    'interval': 60,
    'inputs': [{'name': 'eeg', 'type': 'wav'}],
    'outputs': [
        {'name': 'TOTPOW', 'type': 'num', 'unit': 'dB', 'min': 0, 'max': 100},
        {'name': 'SEF', 'type': 'num', 'unit': 'Hz', 'min': 0, 'max': 30},
        {'name': 'MF', 'type': 'num', 'unit': 'Hz', 'min': 0, 'max': 30},

        {'name': 'DELTA', 'type': 'num', 'unit': '%', 'min': 0, 'max': 100},
        {'name': 'THETA', 'type': 'num', 'unit': '%', 'min': 0, 'max': 100},
        {'name': 'ALPHA', 'type': 'num', 'unit': '%', 'min': 0, 'max': 100},
        {'name': 'BETA', 'type': 'num', 'unit': '%', 'min': 0, 'max': 100},
        {'name': 'GAMMA', 'type': 'num', 'unit': '%', 'min': 0, 'max': 100}
        ]
    }

# http://scipy.github.io/old-wiki/pages/Cookbook/SavitzkyGolay
def savitzky_golay(y, window_size, order, deriv=0, rate=1):
     try:
         window_size = np.abs(np.int(window_size))
         order = np.abs(np.int(order))
     except ValueError:
         raise ValueError("window_size and order have to be of type int")

     if window_size % 2 != 1 or window_size < 1:
         raise TypeError("window_size size must be a positive odd number")

     if window_size < order + 2:
         raise TypeError("window_size is too small for the polynomials order")

     order_range = range(order+1)
     half_window = (window_size -1) // 2
     # precompute coefficients
     b = np.mat([[k**i for i in order_range] for k in range(-half_window, half_window+1)])
     m = np.linalg.pinv(b).A[deriv] * rate**deriv * factorial(deriv)
     # pad the signal at the extremes with values taken from the signal itself
     firstvals = y[0] - np.abs(y[1:half_window+1][::-1] - y[0])
     lastvals = y[-1] + np.abs(y[-half_window-1:-1][::-1] - y[-1])
     y = np.concatenate((firstvals, y, lastvals))
     return np.convolve(m[::-1], y, mode='valid')

def smooth(y):
    #return butter_bandpass(y, 0.5, 50, 128)
    #return lowess(y)
    return savitzky_golay(y, window_size=91, order=3)

def fromhz(f, fres):
    # if type(f) is np.array:
    #     return (f / fres).astype(int)
    return int(f / fres)

def tohz(i, fres):
    return fres * i

def repcols(v, nreps):
    return np.tile(v[:, None], nreps)

def run(inp, opt, cfg):
    data = arr.interp_undefined(inp['eeg']['vals'])
    data -= smooth(np.array(data))
    srate = int(inp['eeg']['srate'])
    nfft = srate * 2  # srate * epoch size
    fres = srate / nfft  # frequency resolution (hz)

    # frequency domain analysis
    EPOCH_SIZE = int(srate * 2)
    STRIDE_SIZE = int(srate * 0.5)
    ps = []
    for epoch_start in range(0, len(data) - EPOCH_SIZE + 1, STRIDE_SIZE):  # 0.5초 마다 겹침
        epoch_w = data[epoch_start:epoch_start + EPOCH_SIZE]  # 2초 epoch
        epoch_w = (epoch_w - np.mean(epoch_w)) * np.blackman(EPOCH_SIZE)  # detrend and windowing
        dft = np.fft.fft(epoch_w)[:srate]  # 실수를 fft 했으므로 절반만 필요하다
        dft[0] = 0  # dc 성분은 지움
        ps.append(2 * np.abs(dft) ** 2) # 파워의 절대값인데 절반 날렸으므로
    ps = np.mean(np.array(ps), axis=0)
    pssum = np.cumsum(ps)  # cummulative sum
    pssum = pssum[1:]
    totpow = pssum[fromhz(30, fres)]
    sef = tohz(np.argmax(pssum > 0.95 * totpow), fres)
    mf = tohz(np.argmax(pssum > 0.5 * totpow), fres)

    delta = pssum[fromhz(4, fres) - 1] / pssum[-1] * 100
    theta = (pssum[fromhz(8, fres) - 1] - pssum[fromhz(4, fres)]) / pssum[-1] * 100
    alpha = (pssum[fromhz(12, fres) - 1] - pssum[fromhz(8, fres)]) / pssum[-1] * 100
    beta = (pssum[fromhz(30, fres) - 1] - pssum[fromhz(12, fres)]) / pssum[-1] * 100
    gamma = (pssum[-1] - pssum[fromhz(30, fres)]) / pssum[-1] * 100

    # pttmax_list.append()
    # pttdmax_list.append({'dt': dmax_dt, 'val': (dmax_dt - rpeak_dt) * 1000})
    # pttmin_list.append({'dt': min_dt, 'val': (min_dt - rpeak_dt) * 1000})
    #
    return [
        [{'dt': cfg['interval'], 'val': 10 * np.log10(totpow)}],
        [{'dt': cfg['interval'], 'val': sef}],
        [{'dt': cfg['interval'], 'val': mf}],
        [{'dt': cfg['interval'], 'val': delta}],
        [{'dt': cfg['interval'], 'val': theta}],
        [{'dt': cfg['interval'], 'val': alpha}],
        [{'dt': cfg['interval'], 'val': beta}],
        [{'dt': cfg['interval'], 'val': gamma}]
    ]
    # pttmin_list,
    # pttdmax_list,
    # arr.get_samples(ecg_data, ecg_srate, ecg_rlist),
    # pttmax_list]

max_stars_repo_name: vitaldb/pyvital

max_issues_repo_name: vitaldb/pyvital

max_forks_repo_name: vitaldb/pyvital


text:
"""
The main sampling script.

Takes pixel positions (x, y) and the number of components (npeaks)
from the command line (via sys.argv[]), sets the priors, the likelihood
function, preforms some sanity checks, and fires up MultiNest through
pymultinest.
"""

import os
import sys
import warnings
import mpi4py
import numpy as np
from astropy import log
from astropy.io import fits
# Those import warnings are annoying
with warnings.catch_warnings():
    warnings.simplefilter('ignore')
    import pymultinest
    from pyspeckit.spectrum.models.ammonia import cold_ammonia_model
    from pyspecnest.ammonia import get_nh3_model
    from pyspecnest.chaincrunch import pars_xy, lnZ_xy, get_zero_evidence

# All the I/O functions now reside here
import opencube
# Configuration for line modelling setup (gets passed to pyspeckit)
from config import line_names, npars
# Path settings
from config import name_id, proj_dir, file_Zs
from config import default_yx, default_npeaks
# Kwargs for digesting and saving spectral cube (making it on the
# fly every time we need a spectrum is too slow!)
from config import cube_save_kwargs
# Finally, MultiNest settings and priors
from config import n_live_points, sampling_efficiency, get_priors_xoff_wrapped

# Compatibility with Python 2.7
try:
    FileNotFoundError
except NameError:
    FileNotFoundError = IOError

def mpi_rank():
    """
    Returns the rank of the calling process.
    """
    comm = mpi4py.MPI.COMM_WORLD
    rank = comm.Get_rank()

    return rank

# pop culture references always deserve their own function
def i_am_root():
    """
    Checks if the running subprocess is of rank 0
    """
    try:
        return True if mpi_rank() == 0 else False
    except AttributeError:
        # not running MPI
        return True

try:
    npeaks = int(sys.argv[1])
except IndexError:
    npeaks = default_npeaks
    log.info("npeaks not specified, setting to {}".format(npeaks))

try:
    yx = int(sys.argv[2]), int(sys.argv[3])
except IndexError:
    yx = default_yx
    log.info("xy-pixel not specified, setting to {}".format(yx[::-1]))

try:
    plotting = bool(int(sys.argv[4]))
except IndexError:
    plotting = 0

if not plotting:
    plot_fit, plot_corner, show_fit, show_corner = False, False, False, False
else: # defaults a for non-batch run
    plot_fit = True
    plot_corner = True
    show_fit = True
    show_corner = True
    from chainconsumer import ChainConsumer  # Optional if no plotting is done
    import matplotlib.pyplot as plt
    plt.rc('text', usetex=True)

y, x = yx
sp = opencube.get_spectrum(x, y, **cube_save_kwargs)

fittype_fmt = 'cold_ammonia_x{}'
fitmodel = cold_ammonia_model
sp.specfit.Registry.add_fitter(fittype_fmt.format(npeaks), npars=npars,
                               function=fitmodel(line_names=line_names))
# is this still needed?
opencube.update_model(sp, fittype_fmt.format(npeaks))

# needed because of this:
# https://github.com/pyspeckit/pyspeckit/issues/179
sp.specfit.fitter.npeaks = npeaks
# npeaks > 1 seems to break because of fitter.parnames is not mirrored
if len(sp.specfit.fitter.parnames) == npars and npeaks > 1:
    sp.specfit.fitter.parnames *= npeaks

# TODO: make a wrapper function for pyspecnest instead!
priors = get_priors_xoff_wrapped(npeaks)

nh3_model = get_nh3_model(sp, line_names, sp.error,
                          priors=priors, npeaks=npeaks)

# Safeguard - check some common causes of failure before scheduling
# a job that would just throw tens of thousands of errors at us
no_valid_chans = not np.any(np.isfinite(nh3_model.ydata))
sanity_check = np.isfinite(
    nh3_model.log_likelihood([15, 5, 15, 0.2, 7, 0.5] * npeaks,
                             nh3_model.npars, nh3_model.dof))
if no_valid_chans or not sanity_check:
    # This should fail if, e.g., the errors are not finite
    log.error("no valid pixels at x={}; y={}. Aborting.".format(*yx[::-1]))
    sys.exit()

# The first process gets to make the directory structure!
output_dir = os.path.join(proj_dir, 'nested-sampling/')
fig_dir = os.path.join(output_dir, 'figs/')
suffix = 'x{1}y{0}'.format(*yx)
chains_dir = '{}chains/{}_{}'.format(output_dir, name_id, suffix)
if not os.path.exists(chains_dir):
    try: # hacks around a race condition
        os.makedirs(chains_dir)
    except OSError as e:
        if e.errno != 17:
            raise

chains_dir = '{}/{}-'.format(chains_dir, npeaks)

# Run MultiNest on the model+priors specified
pymultinest.run(nh3_model.xoff_symmetric_log_likelihood,
                nh3_model.prior_uniform, nh3_model.npars,
                outputfiles_basename=chains_dir,
                verbose=True, n_live_points=n_live_points,
                sampling_efficiency=sampling_efficiency)

# The remainder of the script is not essential for sampling, and can be safely
# moved out into a script of its own.
if i_am_root() and plot_fit:
    # parse the results as sensible output
    from pyspecnest.chaincrunch import analyzer_xy
    a = analyzer_xy(x, y, npeaks, output_dir=output_dir,
                    name_id=name_id, npars=npars)

    a_lnZ = a.get_stats()['global evidence']
    log.info('ln(Z) for model with {} line(s) = {:.1f}'.format(npeaks, a_lnZ))

    try:
        lnZ0 = fits.getdata(file_Zs)[0]
    except (FileNotFoundError, OSError) as e:
        cubes = opencube.make_cube_shh()
        lnZ0 = get_zero_evidence(data=cubes.cube, rms=cubes.errorcube,
                                 normalize=False)

    Zs = lnZ_xy(list(np.arange(npeaks+1)), x=x, y=y, output_dir=output_dir,
                name_id=name_id, silent=True, lnZ0=(lnZ0[y, x], 0))

    log.info('ln(Z{}/Z{}) = {:.2f}'.format(npeaks, npeaks-1,
                                           Zs[npeaks] - Zs[npeaks-1]))
    if npeaks > 1:
        log.info('ln(Z{}/Z{}) = {:.2f}'.format(npeaks, 0, Zs[npeaks] - Zs[0]))

if plot_fit and i_am_root():
    sp.plotter(errstyle='fill')

    mle_pars = pars_xy(x=x, y=y, npars=npars, npeaks=npeaks,
                       output_dir=output_dir, name_id=name_id)

    mle_parinfo = sp.specfit.fitter._make_parinfo(mle_pars, npeaks=npeaks)[0]

    try:
        sp.specfit.plot_fit(xarr=sp.xarr, pars=mle_parinfo,
                            show_components=True)
    except TypeError:
        # eh? does it want pars or parinfo?
        sp.specfit.plot_fit(xarr=sp.xarr, pars=mle_pars, show_components=True)

    # annotate the Bayes factors
    plt.annotate('ln(Z{}/Z{}) = {:.2f}'.format(npeaks, npeaks-1,
                 Zs[npeaks] - Zs[npeaks-1]), xy=(0.05, 0.90),
                 xycoords='axes fraction')
    if npeaks > 1:
        plt.annotate('ln(Z{}/Z{}) = {:.2f}'.format(npeaks, 0,
                     Zs[npeaks] - Zs[0]), xy=(0.05, 0.85),
                     xycoords='axes fraction')
    if show_fit:
        plt.show()

    fig_name = "{}-fit-{}-x{}".format(name_id, suffix, npeaks)
    plt.savefig(os.path.join(fig_dir, fig_name + ".pdf"))

if plot_corner and i_am_root():
    mle_multinest = pars_xy(x=x, y=y, npars=npars, npeaks=npeaks,
                            output_dir=output_dir, name_id=name_id)
    unfrozen_slice = nh3_model.get_nonfixed_slice(a.data.shape, axis=1)
    c = ChainConsumer()
    parameters = nh3_model.get_names(latex=True, no_fixed=True)
    c.add_chain(a.data[:, 2:][unfrozen_slice], parameters=parameters)
    c.configure(statistics="max", summary=True)
    fig = c.plotter.plot(figsize="column")
    fig.get_size_inches()
    fig.set_size_inches(9, 7)

    fig_name = "{}-corner-{}-x{}".format(name_id, suffix, npeaks)
    plt.savefig(fig_dir + fig_name + ".pdf")

    if show_corner:
        plt.show()

reason: import numpy,from astropy


text:
# RUN: SUPPORT_LIB=%mlir_runner_utils_dir/libmlir_c_runner_utils%shlibext \
# RUN:   %PYTHON %s | FileCheck %s

import ctypes
import numpy as np
import os

import mlir.all_passes_registration

from mlir import ir
from mlir import runtime as rt
from mlir import execution_engine
from mlir import passmanager

from mlir.dialects import sparse_tensor as st
from mlir.dialects import builtin
from mlir.dialects.linalg.opdsl import lang as dsl

@dsl.linalg_structured_op
def matmul_dsl(
    A=dsl.TensorDef(dsl.T, dsl.S.M, dsl.S.K),
    B=dsl.TensorDef(dsl.T, dsl.S.K, dsl.S.N),
    C=dsl.TensorDef(dsl.T, dsl.S.M, dsl.S.N, output=True)):
  C[dsl.D.m, dsl.D.n] += A[dsl.D.m, dsl.D.k] * B[dsl.D.k, dsl.D.n]

def build_SpMM(attr: st.EncodingAttr):
  """Build SpMM kernel.

  This method generates a linalg op with for matrix multiplication using
  just the Python API. Effectively, a generic linalg op is constructed
  that computes C(i,j) += A(i,k) * B(k,j) for annotated matrix A.
  """
  module = ir.Module.create()
  f64 = ir.F64Type.get()
  a = ir.RankedTensorType.get([3, 4], f64, attr)
  b = ir.RankedTensorType.get([4, 2], f64)
  c = ir.RankedTensorType.get([3, 2], f64)
  arguments = [a, b, c]
  with ir.InsertionPoint(module.body):

    @builtin.FuncOp.from_py_func(*arguments)
    def spMxM(*args):
      return matmul_dsl(args[0], args[1], outs=[args[2]])

  return module

def boilerplate(attr: st.EncodingAttr):
  """Returns boilerplate main method.

  This method sets up a boilerplate main method that takes three tensors
  (a, b, c), converts the first tensor a into s sparse tensor, and then
  calls the sparse kernel for matrix multiplication. For convenience,
  this part is purely done as string input.
  """
  return f"""
func @main(%ad: tensor<3x4xf64>, %b: tensor<4x2xf64>, %c: tensor<3x2xf64>) -> tensor<3x2xf64>
  attributes {{ llvm.emit_c_interface }} {{
  %a = sparse_tensor.convert %ad : tensor<3x4xf64> to tensor<3x4xf64, {attr}>
  %0 = call @spMxM(%a, %b, %c) : (tensor<3x4xf64, {attr}>,
                                  tensor<4x2xf64>,
                                  tensor<3x2xf64>) -> tensor<3x2xf64>
  return %0 : tensor<3x2xf64>
}}
"""

def build_compile_and_run_SpMM(attr: st.EncodingAttr, support_lib: str,
                               compiler):
  # Build.
  module = build_SpMM(attr)
  func = str(module.operation.regions[0].blocks[0].operations[0].operation)
  module = ir.Module.parse(func + boilerplate(attr))

  # Compile.
  compiler(module)
  engine = execution_engine.ExecutionEngine(
      module, opt_level=0, shared_libs=[support_lib])

  # Set up numpy input and buffer for output.
  a = np.array(
      [[1.1, 0.0, 0.0, 1.4], [0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 3.3, 0.0]],
      np.float64)
  b = np.array([[1.0, 2.0], [4.0, 3.0], [5.0, 6.0], [8.0, 7.0]], np.float64)
  c = np.zeros((3, 2), np.float64)

  mem_a = ctypes.pointer(ctypes.pointer(rt.get_ranked_memref_descriptor(a)))
  mem_b = ctypes.pointer(ctypes.pointer(rt.get_ranked_memref_descriptor(b)))
  mem_c = ctypes.pointer(ctypes.pointer(rt.get_ranked_memref_descriptor(c)))
  # Allocate a MemRefDescriptor to receive the output tensor.
  # The buffer itself is allocated inside the MLIR code generation.
  ref_out = rt.make_nd_memref_descriptor(2, ctypes.c_double)()
  mem_out = ctypes.pointer(ctypes.pointer(ref_out))

  # Invoke the kernel and get numpy output.
  # Built-in bufferization uses in-out buffers.
  # TODO: replace with inplace comprehensive bufferization.
  engine.invoke('main', mem_out, mem_a, mem_b, mem_c)

  # Sanity check on computed result.
  expected = np.matmul(a, b);
  c = rt.ranked_memref_to_numpy(mem_out[0])
  if np.allclose(c, expected):
    pass
  else:
    quit(f'FAILURE')

class SparseCompiler:
  """Sparse compiler passes."""

  def __init__(self, options: str):
    pipeline = (
        f'builtin.func(linalg-generalize-named-ops,linalg-fuse-elementwise-ops),'
        f'sparsification{{{options}}},'
        f'sparse-tensor-conversion,'
        f'builtin.func(linalg-bufferize,convert-linalg-to-loops,convert-vector-to-scf),'
        f'convert-scf-to-std,'
        f'func-bufferize,'
        f'tensor-constant-bufferize,'
        f'builtin.func(tensor-bufferize,std-bufferize,finalizing-bufferize),'
        f'convert-vector-to-llvm{{reassociate-fp-reductions=1 enable-index-optimizations=1}},'
        f'lower-affine,'
        f'convert-memref-to-llvm,'
        f'convert-std-to-llvm,'
        f'reconcile-unrealized-casts')
    self.pipeline = pipeline

  def __call__(self, module: ir.Module):
    passmanager.PassManager.parse(self.pipeline).run(module)

def main():
  support_lib = os.getenv('SUPPORT_LIB')
  assert support_lib is not None, 'SUPPORT_LIB is undefined'
  if not os.path.exists(support_lib):
    raise FileNotFoundError(errno.ENOENT, os.strerror(errno.ENOENT), support_lib)

  # CHECK-LABEL: TEST: testSpMM
  print('\nTEST: testSpMM')
  with ir.Context() as ctx, ir.Location.unknown():
    count = 0
    # Loop over various ways to compile and annotate the SpMM kernel with
    # a *single* sparse tensor. Note that we deliberate do not exhaustively
    # search the full state space to reduce runtime of the test. It is
    # straightforward to adapt the code below to explore more combinations.
    par = 0
    vec = 0
    vl = 1
    e = False
    opt = (f'parallelization-strategy={par} '
           f'vectorization-strategy={vec} '
           f'vl={vl} enable-simd-index32={e}')
    levels = [[st.DimLevelType.dense, st.DimLevelType.dense],
              [st.DimLevelType.dense, st.DimLevelType.compressed],
              [st.DimLevelType.compressed, st.DimLevelType.dense],
              [st.DimLevelType.compressed, st.DimLevelType.compressed]]
    orderings = [
        ir.AffineMap.get_permutation([0, 1]),
        ir.AffineMap.get_permutation([1, 0])
    ]
    bitwidths = [0]
    for level in levels:
      for ordering in orderings:
        for pwidth in bitwidths:
          for iwidth in bitwidths:
            attr = st.EncodingAttr.get(level, ordering, pwidth, iwidth)
            compiler = SparseCompiler(options=opt)
            build_compile_and_run_SpMM(attr, support_lib, compiler)
            count = count + 1
    # CHECK: Passed 8 tests
    print('Passed ', count, 'tests')

if __name__ == '__main__':
  main()

max_stars_repo_name: yuriykoch/llvm

max_issues_repo_name: yuriykoch/llvm

max_forks_repo_name: yuriykoch/llvm


text:
import numpy as np
from icp import ConformalPredictor
from sklearn.metrics import classification_report, matthews_corrcoef, accuracy_score, precision_recall_fscore_support, f1_score
from sklearn.model_selection import KFold
from transformers import AutoTokenizer, AutoModelForSequenceClassification, TrainingArguments, Trainer
import plotly.express as px
import plotly.graph_objects as go
import matplotlib.pyplot as plt
import seaborn as sns
import pandas as pd

class CrossConformalPredictor():

	def __init__(self, X_train, y_train, X_test, y_test, validation_set, model_name, tokenizer, labels, n_folds=3, mondrian=True):
		self.X_train = X_train
		self.y_train = y_train
		self.X_test = X_test
		self.y_test = y_test
		self.validation_set = validation_set
		self.model_name = model_name
		self.tokenizer = tokenizer
		self.labels = labels
		self.num_labels = len(self.labels)
		self.p_values = np.zeros((len(y_test), len(labels)))
		self.n_folds = n_folds
		self.mondrian = mondrian
		self.proper_trains = DatasetDict()
		self.calis = DatasetDict()
		self.preds = {}
		kf = KFold(n_splits=self.n_folds, shuffle=True, random_state=1986)
		for f, (pt, cal) in enumerate(kf.split(data['train'])):
			    self.proper_trains[f] = X_train.select(pt)
			    self.calis[f] = X_train.select(cal)

	def train(self, train_args):
		for i in range(self.n_folds):
		    preds[i] = {}

		    model = AutoModelForSequenceClassification.from_pretrained(self.model_name, num_labels=self.num_labels)
		    model.to(device)

		    trainer = Trainer(
		        model,
		        args,
		        train_dataset=proper_trains[i],
		        eval_dataset=self.validation_set,
		        tokenizer=self.tokenizer,
		        compute_metrics=compute_metrics
		    )
		    trainer.train()

		    # conformal prediction
		    y_cal = calis[i]['label'].numpy()
		    y_test = data['test']['label'].numpy()
		    cal_predictions = trainer.predict(calis[i])
		    y_cal_proba = cal_predictions.predictions
		    y_test_proba = trainer.predict(data['test']).predictions
		    labels = [i for i in range(self.num_labels)]

		    preds[i]['y_cal'] = y_cal
		    preds[i]['y_cal_proba'] = y_cal_proba
		    preds[i]['y_test'] = y_test
		    preds[i]['y_test_proba'] = y_test_proba

	def compute_metrics(eval_pred):
	    predictions, labels = eval_pred
	    preds = np.argmax(predictions, axis=1)
	    precision, recall, f1_macro, _ = precision_recall_fscore_support(labels, preds, average='macro')
	    acc = accuracy_score(labels, preds)
	    f1_weighted = f1_score(labels, preds, average='weighted')
	    return {
	        'accuracy': acc,
	        'f1_macro': f1_macro,
	        'f1_weighted': f1_weighted,
	        'precision': precision,
	        'recall': recall
	    }

	def fit(self, ncf_measure):
		for i in range(n_folds):
		    preds[i]['icp'] = ConformalPredictor(preds[i]['y_cal_proba'], preds[i]['y_cal'],
		                                      preds[i]['y_test_proba'], preds[i]['y_test'],
		                                      self.labels, mondrian=self.mondrian)
		    preds[i]['icp'].fit(ncf_measure)
		    preds[i]['ranks'] = preds[i]['icp'].ranks - 1  # necessary for cross-conformal

		    rz = np.sum(np.array([preds[i]['ranks'] for i in range(self.n_folds)]), axis=0)  # sum ranks over all folds

		    if is_mondrian:
			    for label in labels:
			        self.p_values[:, label] = (rz[:, label] + 1) / (int(sum(X_train['label'] == label)) + 1)
			else:
			    self.p_values = (rz + 1) / (len(self.X_train) + 1)

	def point_prediction_performance(self, digits=2):
		y_hat_test = np.argmax(self.p_values, axis=1)
		print(classification_report(self.y_test, y_hat_test, digits=digits))

	def average_false_pvalue(self):
		y_test_false = np.logical_not(self.y_test).astype(int)
		false_p_values = np.take_along_axis(self.p_values, y_test_false[:, None], axis=1)
		return np.average(false_p_values)

	def top_high_confidence(self, n=10):
		second_largest = np.sort(self.p_values, axis=1)[:, -2]
		return np.argsort(second_largest)[:n]  # it's in increasing order

	def top_low_confidence(self, n=10):
		second_largest = np.sort(self.p_values, axis=1)[:, -2]
		return np.argsort(second_largest)[-n:]  # it's in increasing order

	def top_high_credibility(self, n=10):
		return np.argsort(np.max(self.p_values, axis=1))[-n:]

	def top_low_credibility(self, n=10):
		return np.argsort(np.max(self.p_values, axis=1))[:n]

	def plot_validity(self, labelwise=False):  # TODO: granularity;
		significance = []
		error_rates = []
		pd_labels = []

		if labelwise:
			labels = self.labels
		else:
			labels = ['all']

		for label in labels:
			for i in range(0, 101, 2):
				eps = i/100
				significance.append(eps)
				if labelwise:
					regions = self.p_values[self.y_test==label, :] > eps
					true_label_not_present = np.invert(regions[np.arange(len(regions)), self.y_test[self.y_test==label]])
					error_rate = np.sum(true_label_not_present) / len(self.y_test[self.y_test==label])
				else:
					regions = self.p_values > eps
					true_label_not_present = np.invert(regions[np.arange(len(regions)), self.y_test])
					error_rate = np.sum(true_label_not_present) / len(self.y_test)
				error_rates.append(error_rate)
				pd_labels.append(str(label))

		df = pd.DataFrame({'significance':significance, 'error rate':error_rates, 'label':pd_labels})

		fig = px.scatter(df, x='significance', y='error rate', color='label', width=500, height=400,  # 700, 500
			template='none',
			color_discrete_sequence=px.colors.qualitative.T10
		)
		fig.add_trace(
			go.Scatter(x=significance,
				y=significance,
				opacity=0.2,
				line=dict(
					color='gray',
					dash='dash'
					),
				showlegend=False
				)
		)
		fig.update_layout(
			legend=dict(
			    yanchor="top",
			    y=0.94,
			    xanchor="left",
			    x=0.1),
			# font_family='Officina Sans ITC Book'
		)
		fig.update_xaxes(nticks=10)

		# fig = sns.scatterplot(data=df, x='significance', y='error rate', hue='label')

		return fig


text:
import numpy as np
import torch
from torch import nn
from pytorch_lightning.utilities.cli import instantiate_class

from src.common_module import DomainClassifier
from src.system.source_only import DABase

class DANN(DABase):
    def __init__(self, *args, hidden_dim: int = 1024, gamma: int = 10, **kwargs):
        super(DANN, self).__init__(*args, **kwargs)
        self.gamma = gamma
        self.dc = DomainClassifier(kwargs['embed_dim'], hidden_dim)
        self.criterion_dc = nn.CrossEntropyLoss()

    def training_step(self, batch, batch_idx, optimizer_idx=None, child_compute_already=None):
        (x_s, y_s), (x_t, y_t) = batch
        if child_compute_already:
            (embed_s, y_hat_s), (embed_t, y_hat_t) = child_compute_already
        else:
            embed_s, y_hat_s = self.get_feature(x_s)
            embed_t, y_hat_t = self.get_feature(x_t)

        loss_dc = self.compute_dc_loss(embed_s, embed_t, y_hat_s, y_hat_t)
        loss_cls = self.criterion(y_hat_s, y_s)
        loss = loss_cls+loss_dc

        metric = self.train_metric(y_hat_s, y_s)
        self.log_dict({f'train/loss': loss})
        self.log_dict(metric)
        return loss

    def compute_dc_loss(self, embed_s, embed_t, y_hat_s, y_hat_t):
        y_hat_dc = self.dc(torch.cat([embed_s, embed_t]), self.get_alpha())
        y_dc = torch.cat([torch.zeros_like(y_hat_s[:, 0]), torch.ones_like(y_hat_t[:, 0])]).long()
        loss_dc = self.criterion_dc(y_hat_dc, y_dc)
        return loss_dc

    def get_feature(self, x, domain=None):
        feature = self.backbone(x)
        embed = self.bottleneck(feature)
        y_hat = self.fc(embed)
        return embed, y_hat

    def get_alpha(self):
        return 2. / (1. + np.exp(-self.gamma * self.global_step / (self.num_step * self.max_epochs))) - 1

    def configure_optimizers(self):
        optimizer = instantiate_class([
            {'params': self.backbone.parameters(), 'lr': self.optimizer_init_config['init_args']['lr'] * 0.1},
            {'params': self.bottleneck.parameters()},
            {'params': self.fc.parameters()},
            {'params': self.dc.parameters()},
        ], self.optimizer_init_config)

        lr_scheduler = {'scheduler': instantiate_class(optimizer, self.update_and_get_lr_scheduler_config()),
                        'interval': 'step'}
        return {'optimizer': optimizer, 'lr_scheduler': lr_scheduler}