text:
import numpy as np
from tqdm import trange
from lunch import items_list
from mess_env import serve, est_waste, input_vec, generate_menu, sigmoid
from lunch import dals, sabzis, sweets

import keras.layers as layers
from keras.models import Model, load_model
from keras.optimizers import Adam
import keras.backend as K
from keras.initializers import glorot_normal

import matplotlib.pyplot as plt

def make_policy_model(hidden_layer_neurons, lr):
    num_inputs = len(items_list) # All food items
    num_outputs = len(items_list) + 1 # All food items plus wastage
    inp = layers.Input(shape=np.array(items_list).shape, name='Inputs')
    adv = layers.Input(shape=[1], name='Advantages')
    x1 = layers.Dense(hidden_layer_neurons,
                    activation='relu',
                    use_bias = False,
                    kernel_initializer=glorot_normal(seed=0),
                    name="Hidden1")(inp)
    x2 = layers.Dense(hidden_layer_neurons,
                    activation='relu',
                    use_bias = False,
                    kernel_initializer=glorot_normal(seed=0),
                    name="Hidden2")(x1)
    out = layers.Dense(num_outputs,
                    activation='relu',
                    use_bias = False,
                    kernel_initializer=glorot_normal(seed=0),
                    name="Out")(x2)
    model_train = Model(inputs=[inp,adv], outputs=[out])
    model_train.compile(loss='mean_squared_error', optimizer=Adam(lr))
    model_predict = Model(inputs=[inp], outputs=out)
    return model_train, model_predict

def discount_rewards(r, gamma=0.99):
    prior = 0
    out = []
    for val in r:
        new_val = val + prior * gamma
        out.append(new_val)
        prior = new_val
    return np.array(out[::-1])

def score_model(model, num_tests):
    scores = []
    for test in range(num_tests):
        reward_sum = 0
        for _ in range(7): # Run for a week
            state = np.reshape(input_vec(serve()),[1,num_inputs])
            predicted = model.predict(state)[0]
            action = np.zeros(num_outputs)[:-1]
            action[np.argmax(predicted[:len(dals)])] = 1
            action[len(dals)+np.argmax(predicted[len(dals):len(dals)+len(sabzis)])] = 1
            action[len(dals)+len(sabzis)+np.argmax(predicted[len(sabzis):len(sabzis)+len(sweets)])] = 1
            suggested_menu = generate_menu(action)
            waste = est_waste(suggested_menu)
            reward_sum += 1/waste
        scores.append(reward_sum)
    return np.mean(scores)

def run_model(num_weeks=10):
    training_model, prediction_model = make_policy_model(30, 1e-2)
    reward = reward_sum = 0
    losses = []
    scores = []

    # Placeholders for our observations, outputs and rewards
    states = np.empty(0).reshape(0,num_inputs)
    actions = np.empty(0).reshape(0,num_outputs)
    rewards = np.empty(0).reshape(0,1)
    discounted_rewards = np.empty(0).reshape(0,1)

    for day in trange(7*num_weeks):
        state = np.reshape(input_vec(serve()),[1,num_inputs])
        predicted = prediction_model.predict(state)[0]

        # choosing the proposed move
        action = np.zeros(num_outputs)
        action[np.random.choice(range(len(predicted[:len(dals)])),p=predicted[:len(dals)]/np.sum(predicted[:len(dals)]))] = 1
        action[len(dals) + np.random.choice(range(len(predicted[len(dals):len(dals) + len(sabzis)])),p=(predicted[len(dals):len(dals) + len(sabzis)])/np.sum(predicted[len(dals):len(dals) + len(sabzis)]))] = 1
        action[len(dals) + len(sabzis) + np.random.choice(range(len(predicted[len(dals) + len(sabzis) :len(dals) + len(sabzis) + len(sweets)])),p=(predicted[len(dals) + len(sabzis) :len(dals) + len(sabzis) + len(sweets)])/np.sum(predicted[len(dals) + len(sabzis) :len(dals) + len(sabzis) + len(sweets)]))] = 1
        action[-1] =  predicted[-1]

        states = np.vstack([states, state])
        actions = np.vstack([actions, action])
        suggested_menu = generate_menu(action[:-1])
        waste = est_waste(suggested_menu)
    #     reward = 1/(sigmoid(waste)+0.01)
        reward = 1/waste
        reward_sum += reward
        rewards = np.vstack([rewards, reward])

        if (day + 1) % 7 == 0: # End of week

            discounted_rewards_episode = discount_rewards(rewards, gamma)
            discounted_rewards = np.vstack([discounted_rewards, discounted_rewards_episode])
            rewards = np.empty(0).reshape(0,1)

            discounted_rewards -= discounted_rewards.mean()
            discounted_rewards /= discounted_rewards.std()

            loss = training_model.train_on_batch([states,discounted_rewards],actions)
            losses.append(loss)

            states = np.empty(0).reshape(0,num_inputs)
            actions = np.empty(0).reshape(0,num_outputs)
            discounted_rewards = np.empty(0).reshape(0,1)

        if (day + 1) % print_every == 0:
            score = score_model(prediction_model,50)
            scores.append(score)

    return losses, scores

if __name__=='__main__':
    # Constants defining our neural network
    hidden_layer_neurons = 60
    num_inputs = len(items_list) # All food items
    num_outputs = len(items_list) + 1 # All food items plus wastage
    gamma = .99
    print_every = 14
    batch_size = 7

    num_weeks = 10
    lr = 1e-2
    goal = 0.15

    losses, scores = run_model(500)
    plt.clf()
    plt.plot(losses)
    plt.savefig('losses.png')

    meal = serve()
    state = np.reshape(input_vec(meal),[1,num_inputs])
    prediction_model = load_model('predict.h5')
    prediction_model.compile(loss='mean_squared_error', optimizer=Adam(lr))
    predicted = prediction_model.predict(state)[0]

    food = generate_menu(predicted[:-1])
    print(meal)
    print( food)
    print(predicted[-1])
    print(est_waste(meal))

max_stars_repo_name: deadbeatfour/chowhound

max_issues_repo_name: deadbeatfour/chowhound

max_forks_repo_name: deadbeatfour/chowhound


text:
import os
import cv2
import numpy as np
from tqdm import tqdm
import ersa_utils
import processBlock

class GammaAdjust(processBlock.BasicProcess):
    """
    Adjust images' brightness with O = I^(1/G)
    """
    def __init__(self, gamma, path=None, ds_name=None, name='gamma_adjust'):
        """
        :param path: path to save the adjusted images
        :param ds_name: where to save data if path is None
        :param gamma: gamma to adjust images' brightness
        """
        self.gamma = gamma
        pe_name = '{}_gamma{}'.format(name, self.gamma)
        if path is None:
            assert ds_name is not None
            path = ersa_utils.get_block_dir('data', ['preprocess', ds_name, name])
        super().__init__(pe_name, path, func=self.process)

    def process(self, **kwargs):
        """
        Extract the patches
        :param kwargs:
            file_list: list of files
        :return:
        """
        pbar = tqdm(kwargs['file_list'])
        for files in pbar:
            pbar.set_description('Adjusting {} with gamma={}'.format(os.path.basename(files), self.gamma))
            img = ersa_utils.load_file(files)
            img = self.adjust_gamma(img, self.gamma)
            old_name = os.path.basename(files)
            tile_name = '{}_gamma{}.{}'.format(old_name.split('.')[0], ersa_utils.float2str(self.gamma),
                                               old_name.split('.')[1])
            tile_name = os.path.join(self.path, tile_name)
            ersa_utils.save_file(tile_name, img.astype(np.uint8))

    @staticmethod
    def adjust_gamma(img, gamma):
        """
        Adjust the gamma for given image
        :param img:
        :param gamma:
        :return:
        """
        inv_gamma = 1.0 / gamma
        table = np.array([((i / 255.0) ** inv_gamma) * 255 for i in np.arange(0, 256)]).astype('uint8')
        img_adjust = cv2.LUT(img, table)
        return img_adjust

ext: py

max_stars_repo_name: bohaohuang/ersa

max_issues_repo_name: bohaohuang/ersa

max_forks_repo_name: bohaohuang/ersa

reason: import numpy


text:
"""
match rdi with cnf
"""
import numpy as np
import os, re
from pathlib import Path
import pandas as pd
import dash
import dash_table
import dash_core_components as dcc
import dash_html_components as html
import plotly.express as px
import plotly.graph_objects as go
from mongoengine import connect
#connect('cnf')
from models import (
    CNFFoodName, CNFConversionFactor, CNFNutrientAmount,
        CNFYieldAmount, CNFRefuseAmount, CNFNutrientName
)

# todo: use this for cumulation
rdi_nutrients = {
    'water': "", "fat": "", 'fiber': "", 'linoleicAcid': "", 'alphaLinolenicAcid': "",
    'carbohydrate': "", 'protein': "",
    'calcium': "", 'chromium': "", 'copper': "", 'fluoride': "", 'iodine': "", 'iron': "",
    'magnesium': "",
    'manganese': "", 'molybdenum': "", 'phosphorus': "", 'selenium': "", 'zinc': "",
    'potassium': "",
    'sodium': "", 'chloride': "",
    'vitaminA': "", 'vitaminC': "", 'vitaminD': "", 'vitaminE': "", 'vitaminK': "", 'thiamin': "",
    'riboflavin': "", 'niacin': "", 'vitaminB6': "", 'folate': "", 'vitaminB12': "",
    'pantothenicAcid': "", 'biotin': "", 'choline': ""
}

def make_food_to_id_dict():
    """
    call this right away
    """
    food_name_query_set = CNFFoodName.objects(description__exists=True)
    food_names_arr = []
    food_ids_arr = []

    for food in food_name_query_set:
        '''
        # dict with description as key, values: id, group, source
        food_dict[food['description']] = [food['id'], food['food_group'],
                                          food['food_source'],
                                          food['scientific_name']]
        '''
        food_names_arr.append(food['description'])
        food_ids_arr.append(food['id'])

    #assert len(food_names_arr) == len(food_ids)

    food_to_id_dict = dict(zip(food_names_arr, food_ids_arr))

    return food_to_id_dict, food_names_arr, food_ids_arr

#global helper fn
def get_unit_names(food_id):
    """
    class_collection: one of the CNF collections
    get the unit names and return after cleaning and keeping only unique
    """
    # get units
    conversions = CNFConversionFactor.objects.filter(food=str(food_id))
    unit_names = []
    for c in conversions:
        unit_names.append(c.measure.name)
    #print(f'unit_names {unit_names}')
    # take out the numbers and whitespace
    # todo: pattern has to include all non-alpha
    pattern = '[0-9]'
    temp_arr = [re.sub(pattern, '', x) for x in unit_names]
    unit_names = [x.strip() for x in temp_arr]
    # only keep the unique measure names
    measure_set = set(unit_names)
    unit_names = list(measure_set)

    return unit_names

def make_foodgroup_df(food_id):
    """
    ("Group"), html.Th("Source"), html.Th("Scientific Name")])
    """

    food = CNFFoodName.objects.get(id=str(food_id))
    food_grp = food.food_group.name
    food_src = food.food_source.description
    food_sci_name = "n/a"
    if food.scientific_name:
        food_sci_name = food.scientific_name
    foodgroup_df = pd.DataFrame(
        {'Group': [food_grp],
         "Source": [food_src],
         "Scientific Name": [food_sci_name]
         }
    )
    return foodgroup_df

def make_conversions_df(food_id):
    """
    param
    """
    food = CNFFoodName.objects.get(id=str(food_id))
    conversions = CNFConversionFactor.objects.filter(food=food)
    measure_names = []
    measure_vals = []
    for c in conversions:
        measure_names.append(c.measure.name)
        measure_vals.append(c.value)

    conversions_df = pd.DataFrame(
        {'Name': measure_names,
         'Multiplier': measure_vals
         })

    return conversions_df

#global helper fn
def make_nutrients_df(food_id):
    food = CNFFoodName.objects.get(id=str(food_id))
    nutrients = CNFNutrientAmount.objects.filter(food=food, nutrient_value__gt=0)

    nutrient_names = []
    nutrient_vals = []
    nutrient_units = []

    for n in nutrients:
        nutrient_names.append(n.nutrient_name.name)
        nutrient_vals.append(n.nutrient_value)
        nutrient_units.append(n.nutrient_name.unit)

    nutrients_df = pd.DataFrame(
        {"Name": nutrient_names,
         "Value": nutrient_vals,
         "Units": nutrient_units
         })

    return nutrients_df

def get_conversions_multiplier(conversions_df, units):
    curr_multiplier = ''
    measure_num = -1
    measure_name = ''
    for index, row in conversions_df.iterrows():
        if units in str(row['Name']):
            # get number
            measure_name = str(row['Name'])
            measure_num = float(re.findall("\d+", measure_name)[0])
            curr_multiplier = float(row['Multiplier'])
            break
    return curr_multiplier, measure_num

def mult_nutrients_df(nutrients_df,curr_multiplier, measure_num, amt):
    new_nutrients_df = nutrients_df.copy()
    new_nutrients_df.set_index('Name', inplace=True, drop=False)
    multiplier_factor = float((float(amt)) / measure_num)
    new_multiplier = multiplier_factor * float(curr_multiplier)

    # multiply all nutrients by new_multiplier
    for index, row in nutrients_df.iterrows():
        row_val = float(row['Value']) * new_multiplier
        #nutrients_df.row['Value'] = str(new_row_val)
        nut_name = row['Name']
        new_nutrients_df.loc[[nut_name], ['Value']] = str(row_val)

    return new_nutrients_df


text:
import baseline2_model
import torch
import torchvision
import torchvision.transforms as transforms
import torch.optim as optim
import numpy as np

import matplotlib.pyplot as plt
import cv2
import time

# PyTorch TensorBoard support
from torch.utils.tensorboard import SummaryWriter
from datetime import datetime
VIDEO_PATH = '../../Anticipating-Accidents/dataset/videos/testing/positive/'
def detection_collate(batch):
    """Custom collate fn for dealing with batches of images that have a different
    number of associated object annotations (bounding boxes).
    Arguments:
        batch: (tuple) A tuple of tensor images and lists of annotations
    Return:
        A tuple containing:
            1) (tensor) batch of images stacked on their 0 dim
            2) (list of tensors) annotations for a given image are stacked on
                                 0 dim
    """
    data = []
    labels = []
    file_name = []
    bbox = []
    for sample in batch:
        a,b,c,d = sample
        data.append(a)
        labels.append(b)
        file_name.append(c)
        bbox.append(d)
    data = np.array(data)
    labels = np.array(labels)
    return torch.from_numpy(data).to('cuda'), torch.from_numpy(labels).to('cuda'),file_name,bbox

def training(batch_size,n_epoch,learning_rate):
    net = baseline2_model.Baseline_SA(False,features_size=8192).to('cuda')
    net.train()
    criterion = baseline2_model.custom_loss(100)
    optimizer = optim.Adam(net.parameters(), lr=learning_rate)
    trainset = baseline2_model.SADataset('/mnt/sdb/Dataset/SA_Maskformer',True)
    trainloader = torch.utils.data.DataLoader(trainset, batch_size=batch_size,
                                          shuffle=True, num_workers=0,collate_fn=detection_collate)
    # testset = baseline2_model.SADataset('/mnt/sdb/Dataset/SA',False)
    # testloader = torch.utils.data.DataLoader(testset, batch_size=batch_size,
                                        #   shuffle=True, num_workers=0,collate_fn=detection_collate)
    total_batch = len(trainloader)
    best_vloss = 100000.0
    for epoch in range(n_epoch):  # loop over the dataset multiple times
        print("Epoch %d/%d"%(epoch+1,n_epoch))
        running_loss = 0.0
        for i, data in enumerate(trainloader):
            # get the inputs; data is a list of [inputs, labels]
            inputs,labels,_,_ = data
            # zero the parameter gradients
            optimizer.zero_grad()
            # forward + backward + optimize
            pred, prob, pred_loss= net(inputs)
            loss = criterion(pred_loss, labels)
            loss.backward()
            # print statistics
            running_loss += loss.item()
            # print(loss.item(),running_loss)
            # running_loss = 0.0
            print("Batch : %d/%d" % (i+1,total_batch),end='\r')
            optimizer.step()
        # optimizer.zero_grad()
        print('loss: %f'%(running_loss/total_batch))
        if running_loss < best_vloss:
            best_vloss = running_loss
            model_path = 'baseline2_maskformer/model_{}'.format(epoch+1)
            print("Saving model..")
            torch.save(net.state_dict(), model_path)
    return
def demo(batch_size,model_path):
    net = baseline2_model.Baseline_SA(False)
    net.load_state_dict(torch.load(model_path))
    net.to('cuda')
    net.eval()
    testset = baseline2_model.SADataset('/mnt/sdb/Dataset/SA',False)
    testloader = torch.utils.data.DataLoader(testset, batch_size=batch_size,
                                          shuffle=True, num_workers=0,collate_fn=detection_collate)
    with torch.no_grad():
        for data in testloader:
            inputs,labels,file_names,det = data
            preds, probs, _= net(inputs)
            for i,label in enumerate(labels):
                if label:
                    pred = preds[i].cpu().numpy()
                    weight = probs[i].cpu().numpy()
                    plt.figure(figsize=(14,5))
                    plt.plot(pred[:90,1],linewidth=3.0)
                    plt.ylim(0, 1)
                    plt.ylabel('Probability')
                    plt.xlabel('Frame')
                    plt.show()
                    file_name = file_names[i]
                    bboxes = det[i]
                    new_weight = weight*255
                    counter = 0
                    print(file_name)
                    cap = cv2.VideoCapture(VIDEO_PATH+str(file_name))
                    ret, frame = cap.read()
                    while(ret):
                        attention_frame = np.zeros((frame.shape[0],frame.shape[1]),dtype = np.uint8)
                        now_weight = new_weight[counter,:]
                        new_bboxes = bboxes[counter,:]
                        # index = np.argsort(now_weight)
                        # print(now_weight)
                        for num_box in range(20):
                            if now_weight[num_box]/255.0>0.4:
                                cv2.rectangle(frame,(int(new_bboxes[num_box,0]),int(new_bboxes[num_box,1])),(int(new_bboxes[num_box,2]),int(new_bboxes[num_box,3])),(0,255,0),3)
                            else:
                                cv2.rectangle(frame,(int(new_bboxes[num_box,0]),int(new_bboxes[num_box,1])),(int(new_bboxes[num_box,2]),int(new_bboxes[num_box,3])),(255,0,0),2)
                            font = cv2.FONT_HERSHEY_SIMPLEX
                            cv2.putText(frame,str(round(now_weight[num_box]/255.0*10000)/10000),(int(new_bboxes[num_box,0]),int(new_bboxes[num_box,1])), font, 0.5,(0,0,255),1,cv2.LINE_AA)
                            attention_frame[int(new_bboxes[num_box,1]):int(new_bboxes[num_box,3]),int(new_bboxes[num_box,0]):int(new_bboxes[num_box,2])] = now_weight[num_box]
                            # cv2.putText(frame,str(num_box),(int(new_bboxes[num_box,3]),int(new_bboxes[num_box,2])), font, 1,(0,0,255),1,cv2.LINE_AA)
                        attention_frame = cv2.applyColorMap(attention_frame, cv2.COLORMAP_HOT)
                        dst = cv2.addWeighted(frame,0.6,attention_frame,0.4,0)
                        cv2.putText(dst,str(counter+1),(10,30), font, 1,(255,255,255),3)
                        cv2.imshow('result',dst)
                        time.sleep(0.05)
                        c = cv2.waitKey(50)
                        ret, frame = cap.read()
                        if c == ord('q') and c == 27 and ret:
                            break
                        counter += 1
                    cap.release()

                cv2.destroyAllWindows()

if __name__ == '__main__':
    learning_rate = 0.001
    batch_size = 10
    n_epoch = 40
    model_path = 'baseline2/model_40'
    training(batch_size,n_epoch,learning_rate)
    # demo(batch_size,model_path)

max_stars_repo_name: HCIS-Lab/Interaction-benchmark

max_issues_repo_name: HCIS-Lab/Interaction-benchmark

max_forks_repo_name: HCIS-Lab/Interaction-benchmark


text:
from control4.core.mdp import MDP
import theano, theano.tensor as TT

class SymbolicMDP(MDP):

    """
    MDP for which call() is defined by symbolic operation (in theano)
    """

    def __init__(self):
        x = TT.matrix(name="x",dtype=self.input_dtype("x"))
        u = TT.matrix(name="u",dtype=self.input_dtype("u"))
        outputs = self.symbolic_call(x,u)
        self._call = theano.function([x,u],outputs)

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

    # New

    def symbolic_call(self,x,u):
        raise NotImplementedError

    def output_names(self):
        raise NotImplementedError

    # Inherited

    def initialize_mdp_arrays(self):
        raise NotImplementedError

    def input_info(self):
        raise NotImplementedError

    def output_info(self):
        raise NotImplementedError

    def plot(self, mdp_arrs, pol_arrs):
        raise NotImplementedError

    def cost_names(self):
        raise NotImplementedError

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

    def call(self,input_arrs):
        outputs = self._call(input_arrs["x"],input_arrs["u"])
        return dict(zip(self.output_names(),outputs))

max_stars_repo_name: sfpd/rlreloaded

max_issues_repo_name: sfpd/rlreloaded

max_forks_repo_name: sfpd/rlreloaded

reason: import theano


text:
#!/usr/bin/env python2
# -*- coding: utf-8 -*-
"""
Created on Tue Sep 24 10:09:58 2019

@author: daniel
"""
from pylab import *
import pickle
import os
from pathlib import Path
import numpy as np
import networkx as nx
import scipy.stats as sts
import scipy.spatial as sps

from statsmodels.stats.multitest import multipletests
import networkx as nx

import matplotlib.pyplot as plt
plt.style.use('ggplot')

from itertools import chain, combinations
from itertools import chain
from itertools import groupby

def replication_w_error(p_in, p_del, env, model_reactions,phenotype, r_dict, pool):
    process = np.random.choice(['insert','delete', 'remain' ], p=np.array([p_in, p_del, 1-(p_in+p_del)]))
    if process =='insert':
        reaction = np.random.choice(pool)
        mod = model_reactions.copy()
        mod.add(reaction)
        reacs = [r_dict[i] for i in mod]
        m=Model(reacs)
        pc = m.is_growing(env, phenotype)
    elif process=='delete':
        reaction = np.random.choice(list(model_reactions))
        mod = model_reactions.copy()
        mod.discard(reaction)
        reacs = [r_dict[i] for i in mod]
        m=Model(reacs)
        pc = m.is_growing(env, phenotype)
    else:
        mod = model_reactions.copy()
        reacs = [r_dict[i] for i in mod]
        m=Model(reacs)
        pc = m.is_growing(env, phenotype)

    return mod,pc

def moran_process(population, panEFM, env, iter_n, p_in,p_del, phenotype, r_dict):
    panEFM = list(panEFM.keys())
    pop = {}
    for i in range(population):
        pop[i] = set(np.random.choice(panEFM))
    pop_i = pop.copy()
    for i in range(iter_n):

        pool = list(chain.from_iterable(list(pop.values())))
        selected = np.random.choice(list(pop.keys()), size=2)
        m,pc = replication_w_error(p_in, p_del, env, pop[selected[1]],phenotype, r_dict, pool)
        if pc:
            pop[selected[0]] = m

    return pop_i, pop

class Reaction:
    '''
    reaction object. Reactions have a name, reactants, and products.
    Can be reversible. Stoichiometry is one-to-one in the toy model.
    '''

    def __init__(self, name, reactants, products, reversibility=False):
        self.name = name
        self.reactants = set(reactants)
        self.products = set(products)
        self.reversibility = reversibility

    def is_active(self, metabolite_pool):
        '''
        given a metabolite pool, will the reaction be used.

        Parameters
        ----------
        metabolite_pool : TYPE set
            set of metabolites. If the reactants (in any direction)
            are in the set, the reaction is active and produces products.

        Returns
        -------
        active : TYPE bool

        '''
        active=True
        mp = set(metabolite_pool)

        if len(mp.intersection(self.reactants))<len(self.reactants):
            active=False
        if self.reversibility and active==False:
            if len(mp.intersection(self.products))==len(self.products):
                active=True

        return active
    def build_string(self):
        '''
        build a string representation for the reaction
        '''
        if self.reversibility:
            symbol=' <--> '
        else:
            symbol=' --> '

        return ' + '.join(self.reactants) + symbol + ' + '.join(self.products)

class Model:
    '''
    Model object. Model has intracellular and extracellular reactions and metabolites;
    Extracellular metabolites are identified by "EX" in their names.
    '''
    def __init__(self, reactions,name='model'):
        self.name = name
        self.reactions = set(reactions)
        self.ic_products, self.ic_reactants, self.ex_products, self.ex_reactants, self.reactants, self.products, self.m_to_reactions= self.parse_reactions(self.reactions)

    def parse_reactions(self, reactions):
        ic_prod, ic_react, ex_prod, ex_react, prod, react = set(),set(),set(),set(),set(),set()
        m_to_reactions={}

        for i in reactions:
            for r in i.reactants:
                if r in m_to_reactions:
                    m_to_reactions[r].append(i)
                else:
                    m_to_reactions[r]=[i]
                react.add(r)
                if 'EX' in r:
                    ex_react.add(r)
                else:
                    ic_react.add(r)
            for p in i.products:
                if p in m_to_reactions:
                    m_to_reactions[p].append(i)
                else:
                    m_to_reactions[p]=[i]
                prod.add(p)
                if 'EX' in p:
                    ex_prod.add(p)
                else:
                    ic_prod.add(p)
        return ic_prod, ic_react, ex_prod, ex_react, prod,react, m_to_reactions

    def is_producing_phenotype(self, metabolite_pool, phenotype):

        mp = set(metabolite_pool)

        allmets = mp.copy()

        allmets=allmets.union(self.products)

        #print allmets
        active=[]

        for i in self.reactions:
            if i.is_active(allmets):
                active.append(i)

        t=[]
        while set(t)!=set(active):

            t=active[:]
            active=[]

            allmets = mp.copy()
            for i in t:
                allmets=allmets.union(i.products)

            for i in t:
                if i.is_active(allmets):
                    active.append(i)

        allmets=mp.copy()
        for i in active:
            allmets=allmets.union(i.products)

        #print allmets
        #print [i.name for i in active]
        if phenotype==allmets.intersection(phenotype):
            return True
        else:
            return False

    def is_growing(self, metabolite_pool, phenotype):

        mp = set(metabolite_pool)

        allmets = mp.copy()

        allmets=allmets.union(self.products)

        #print allmets
        active=[]

        for i in self.reactions:
            if i.is_active(allmets):
                active.append(i)

        t=[]
        while set(t)!=set(active):

            t=active[:]
            active=[]

            allmets = mp.copy()
            for i in t:
                allmets=allmets.union(i.products)

            for i in t:
                if i.is_active(allmets):
                    active.append(i)

        allmets=mp.copy()
        for i in active:
            allmets=allmets.union(i.products)

        #print allmets
        #print [i.name for i in active]
        if phenotype==allmets.intersection(phenotype):
            #final check!
            prods=set()
            reacts=set()
            for i in active:
                for p in i.products:
                    prods.add(p)
                for r in i.reactants:
                    reacts.add(r)
            prods=prods-phenotype
            prods = prods
            b=prods.difference(reacts.intersection(prods))

            if b==self.ex_products.intersection(b):
                return True
            else:
                remove=[]

                for i in b:
                    remove+=self.m_to_reactions[i]

                r =self.reactions-set(remove)

                m=Model(r)
                if m.is_producing_phenotype(metabolite_pool, phenotype):
                    return True
                else:
                    return False
        else:
            return False

def autoIncrement():
 global rec
 pStart = 1 #adjust start value, if req'd
 pInterval = 1 #adjust interval value, if req'd
 if (rec == 0):
  rec = pStart
 else:
  rec = rec + pInterval
 return str(rec).zfill(4)

def get_panEFM(reactions, env, phenotype, iterations):

    rd = {i.name:i for i in reactions}

    reaction_list = list(rd.keys())

    rl = reaction_list[:]

    panEFM={}
    envs={}
    for i in range(iterations):
        np.random.shuffle(rl)
        np.random.shuffle(rl)

        rlc=rl[:]

        rlt=rlc[:]

        for z in range(len(rl)):
            target = rlc[z]
            rlt.remove(target)

            r = [rd[w] for w in rlt]
            m = Model(r)
            #print rlt
            if not m.is_growing(env, phenotype):
                rlt.append(target)

        rlt.sort()
        r = [rd[w] for w in rlt]
        m = Model(r)

        if frozenset(rlt) in panEFM:
            panEFM[frozenset(rlt)]+=1
            if m.ex_reactants not in envs[frozenset(rlt)]:
                envs[frozenset(rlt)].append(m.ex_reactants)
        else:
            panEFM[frozenset(rlt)]=1
            envs[frozenset(rlt)] = [m.ex_reactants]

    return panEFM,envs

def get_freq(reaction_list, panEFM):
    f = {i:0 for i in reaction_list}

    for i_set in panEFM:

        for reaction in i_set:
            f[reaction]+=1./len(panEFM)
    return np.array([f[i] for i in reaction_list])

def get_average_env(env_vec, env_d):
    k = np.zeros(len(env_vec))
    denom = len(env_d)

    for i, t in enumerate(env_d):
        for j,m in enumerate(env_vec):
            if m in env_d[t][0]:
                k[j] += 1/denom
    return k

def cosine(a,b):
    return np.inner(a, b)/(np.linalg.norm(a)*np.linalg.norm(b))

def assign_to_class(vector, sig_level=0.05, mt=False):
    '''
    for environment_frequency - average
    -1: reaction is significantly less frequent in this environment
    0: reaction is neutral in this environment
    1: reaction is significantly more frequent in this environment

    '''
    if sum(abs(vector))==0:
        return np.zeros(len(vector))
    SE = np.std(vector)/np.sqrt(len(vector))
    z_norm_deviate =(vector-np.mean(vector))/SE

    sign = z_norm_deviate/np.abs(z_norm_deviate)

    pvs =sts.norm.sf(np.abs(z_norm_deviate))*2
    if mt:
        pvs = multipletests(pvs, alpha=sig_level, method='bonferroni')[1]

    pvs[pvs>=sig_level]=2
    pvs[(pvs<sig_level) & (sign==-1)]=3
    pvs[(pvs<sig_level) & (sign==1)]=4

    pvs[pvs==2]=0
    pvs[pvs==3]=-1
    pvs[pvs==4]=1

    return pvs

def assign_to_rank(vector, p_cutoff, n_cutoff):
    '''

    '''
    # SE = np.std(vector)/np.sqrt(len(vector))
    # z_norm_deviate =(vector-np.mean(vector))/SE
    ranks = np.zeros(len(vector))
    ranks[(vector>0) & (vector>=p_cutoff)] =1
    ranks[(vector<0) & (vector<=n_cutoff)]=-1
    return ranks

def build_association_network(association_d, reacs, mets):
    G=nx.MultiGraph()

    for i in reacs:
        G.add_node(i, type='reaction')
    for i in mets:
        G.add_node(i, type='metabolite')

    for i,r in enumerate(reacs):
        for z,m in enumerate(mets):
            if association_d[r][z]==-1:
                G.add_edge(r,m,type='negative')
            elif association_d[r][z]==1:
                G.add_edge(r,m,type='positive')
    return G

rec=0
mets  =['m'+autoIncrement() for i in range(23)]

ex_mets=['EX_' + mets[i] for i in range(10)]

rec=0
react_ids = ['r' +autoIncrement() for i in range(15)]

r1 = Reaction(react_ids[0], [ex_mets[0]], [mets[0], mets[1]])
r2 = Reaction(react_ids[1], [ex_mets[1]], [mets[0], mets[2]])
r3 = Reaction(react_ids[2], [ex_mets[2]], [mets[3]])
r4 = Reaction(react_ids[3], [ex_mets[3]], [mets[4]])
r5 = Reaction(react_ids[4], [ex_mets[4]], [mets[5]])
r6 = Reaction(react_ids[5], [ex_mets[5],ex_mets[6]], [mets[6]])
r7 = Reaction(react_ids[6], [ex_mets[7]], [mets[7]])
r8 = Reaction(react_ids[7], [ex_mets[8]], [mets[7]])
r9 = Reaction(react_ids[8], [mets[0]], [mets[8]])
r10 = Reaction(react_ids[9], [mets[1]], [ex_mets[9]])
r11 = Reaction(react_ids[10], [mets[6]], [mets[7]])
r12 = Reaction(react_ids[11], [mets[5],mets[7]], [mets[9]])
r13 = Reaction(react_ids[12], [mets[4]], [mets[10]])
r14 = Reaction(react_ids[13], [mets[2],mets[8]], [mets[11]])
r15 = Reaction(react_ids[14], [mets[3]], [mets[4], mets[2]])

reactions ={r1,r2,r3,r4,r5,r6,r7,r8,r9,r10,r11,r12,r13,r14,r15}

react_dict = {i.name:i for i in reactions}

k=list(react_dict.keys())
k.sort()

phenotype={mets[9], mets[10], mets[11]}

m=Model(reactions)

m.is_growing(ex_mets, phenotype)

r,e=get_panEFM(reactions, ex_mets, phenotype, 20000)

rd={}

emet = np.array(ex_mets)

rf = np.zeros((20000, len(k)))

envs=np.zeros((20000, len(ex_mets)))

used=[]
env_that_sup_growth = []

for i in range(20000):
    j = np.round(np.random.uniform(size=10))
    if list(j) not in used:
        used.append(list(j))

        env = emet[j.astype(np.bool)]
        if m.is_growing(env,phenotype):
            env_that_sup_growth.append(j)

            panEFM,environment = get_panEFM(reactions, env, phenotype, 1000)
            envs[i]=get_average_env(emet, environment)
            rf[i] = get_freq(k, panEFM)

#eliminate the environments with no growth

s_rf = np.sum(rf, axis=1)

quantile_freq=.90
quantile_met=.80
quantile_ass=.95

used_environment = envs[s_rf!=0].copy()

full_freq_m = rf[s_rf!=0].copy()
reactome = np.array(k)
transporter = emet.copy()

av_freq_m = np.mean(full_freq_m, axis=0)
diff_freq_m = full_freq_m-av_freq_m

#filter out noise and find reactions that are driven by the environment
m_diff_freq_m = np.max(np.abs(diff_freq_m), axis=0)
env_driven_reactome = reactome[m_diff_freq_m>.005]
diff_freq_m_envd = diff_freq_m.T[m_diff_freq_m>.005].T

diff_freq_pool = diff_freq_m_envd.flatten()
clss_freq_m = np.zeros(diff_freq_m_envd.shape)
for i,v in enumerate(diff_freq_m_envd):
    clss_freq_m[i] = v

#for the environment
av_used_env = np.mean(used_environment, axis=0)
diff_used_env = used_environment-av_used_env
#filter out noise and find metabolites that are driven by the environment
m_diff_used_env = np.max(np.abs(diff_used_env), axis=0)
driving_mets = transporter[m_diff_used_env>0.005]
diff_used_env_envd = diff_used_env.T[m_diff_used_env>0.005].T

diff_used_env_pool = diff_used_env_envd.flatten()
clss_used_env = np.zeros(diff_used_env_envd.shape)
for i,v in enumerate(diff_used_env_envd):
    clss_used_env[i] = v

s_clss_fm = np.sum(np.abs(clss_freq_m), axis=0)
s_clss_ue = np.sum(np.abs(clss_used_env), axis=0)

#env_driven_reactome
envd_reactions =env_driven_reactome[s_clss_fm!=0]

#driving_metabolites
dm = driving_mets.copy()
dm = dm[s_clss_ue!=0]

#profiles
envd_prof = clss_freq_m.T[s_clss_fm!=0].T
dm_prof = clss_used_env.T[s_clss_ue!=0].T

cosine_dict={}

for i, reac in enumerate(envd_prof.T):
    cosine_dict[envd_reactions[i]] = np.array([cosine(reac.flatten(), metab.flatten()) for metab in dm_prof.T])

cosine_pool = np.array(list(cosine_dict.values())).flatten()
pc = np.quantile(cosine_pool[cosine_pool>0], quantile_ass)
nc = np.quantile(cosine_pool[cosine_pool<0], 1-quantile_ass)
association_d={}

for i, reac in enumerate(envd_prof.T):
    v = cosine_dict[envd_reactions[i]]
    ncp=sqrt(len(v))*(mean(v**2)-0)/std(v**2)
    alpha=ncp/sqrt(len(v))
    pa = alpha/2
    na = -pa
    association_d[envd_reactions[i]] = assign_to_rank(v, pa,na)

g= build_association_network(association_d, envd_reactions, dm)

#to save the network

#nx.write_graphml(g, FileToSave)

reacs = list(react_dict.keys())
predicted = []
average = []
predicted_environments=[]
true_used_env=[]
panEFM_avg_env = []
evolved_freq = []

# evolve conditioned to an environment
for i in range(len(env_that_sup_growth)):
    print(i)
    j=env_that_sup_growth[i]
    env = emet[j.astype(np.bool)]
    panEFM,environment = get_panEFM(reactions, env, phenotype, 1000)
    envs=get_average_env(emet, environment)
    panEFM_avg_env.append(envs)
    pop_i, pop = moran_process(1000, r, env, 1000000,.05,.5,phenotype, react_dict)
    g1, g2 = [],[]

    for pp in pop_i:
        h=list(pop_i[pp].copy())
        h.sort()
        if h not in g1:
            g1.append(h)

    for pp in pop:
        h=list(pop[pp].copy())
        h.sort()
        if h not in g2:
            g2.append(h)

    #get the reaction_frequencies of evolved population

    f1,f2,f3=[],[],[]
    p1 = list(chain.from_iterable(list(g1)))
    p2 = list(chain.from_iterable(list(g2)))
    p3 = list(chain.from_iterable(list(r.keys())))
    for rr in reactome:
        f1.append(p1.count(rr)/len(g1))
    for rr in reactome:
        f2.append(p2.count(rr)/len(g2))
    for rr in reactome:
        f3.append(p3.count(rr)/9)
    f1=np.array(f1)
    f2=np.array(f2)
    f3=np.array(f3)

    evolved_freq.append(f2)

    #get the average metabolite usage from the evolved population
    used_mets=[]
    for mm in g2:
        reacs=[react_dict[z] for z in mm]
        m=Model(reacs)
        used_mets.append(m.ex_reactants)
    used_mets = list(chain.from_iterable(used_mets))
    mf=[]
    for mm in dm:
        mf.append(used_mets.count(mm)/len(g2))
    true_used_env.append(mf)

    from sklearn.linear_model import MultiTaskElasticNetCV as EN
    enet  = EN(cv=50, max_iter=100000)
    x = full_freq_m.T[m_diff_freq_m>.005].T
    y = used_environment.T[m_diff_used_env>0.005].T
    mod=enet.fit(x, y)
    p = mod.predict(f2[m_diff_freq_m>.005].reshape(1,-1))

    p=p.flatten()
    p = p+abs(min(p))
    p=p/max(p)

    c = [sts.pearsonr(mf,used_environment[ee][m_diff_used_env>0.005])[0] for ee in range(len(used_environment))]

    predicted.append(sts.pearsonr(p, mf)[0])

    average.append(mean(c))

    predicted_environments.append(p)

store = {'full_freq_m': full_freq_m, 'diff_freq_m':diff_freq_m, 'envd_prof':envd_prof, 'envd_reactions': envd_reactions,'used_environment':used_environment, 'ex_mets':ex_mets, 'dm':dm, 'diff_used_env':diff_used_env,         'dm_prof':dm_prof, 'x':x, 'y':y, 'predicted':predicted, 'average':average, 'predicted_environment':predicted_environments,'true_used_env':true_used_env, 'panEFM_avg_env': panEFM_avg_env, 'evolved_freq':evolved_freq, 'cosine_dict':cosine_dict, 'met_filter1':m_diff_used_env>0.005, 'met_filter2':s_clss_ue!=0, 'reac_filter1':m_diff_freq_m>.005, 'reac_filter2':s_clss_fm!=0, 'evolved_freq':evolved_freq}

data_folder = os.path.join(Path(os.getcwd()).parents[1], 'data')
pickle.dump(store, open(data_folder + '/pickles/toy_model.pkl', 'wb'))

max_stars_repo_name: danielriosgarza/NutritionOrNature

max_issues_repo_name: danielriosgarza/NutritionOrNature

max_forks_repo_name: danielriosgarza/NutritionOrNature

reason: import numpy,import scipy,from statsmodels,import networkx


text:
import numpy

data1 = []
data2 = []

i = 9
j = 19

for val in range(100):
    data1.append(val*i)
    data2.append(val*j)

print len(data1), data1
print len(data2), data2

nd_arr1 = numpy.array(zip(data1,data2))

max_stars_repo_name: mgadagin/PythonClass

max_issues_repo_name: mgadagin/PythonClass

max_forks_repo_name: mgadagin/PythonClass


text:
#!/usr/bin/env python3
# -*- coding: utf-8 -*-

"""
評価用のテストパターン作成ツール集

"""

import os
import cv2
import matplotlib.pyplot as plt
import numpy as np
from colour.colorimetry import CMFS, ILLUMINANTS
from colour.models import XYZ_to_xy, xy_to_XYZ, XYZ_to_RGB, RGB_to_XYZ
from colour.models import xy_to_xyY, xyY_to_XYZ, Lab_to_XYZ
from colour.models import BT709_COLOURSPACE
from colour.utilities import normalise_maximum
from colour import models
from colour import RGB_COLOURSPACES, COLOURCHECKERS
from scipy.spatial import Delaunay
from scipy.ndimage.filters import convolve
import math

import transfer_functions as tf

CMFS_NAME = 'CIE 1931 2 Degree Standard Observer'
D65_WHITE = ILLUMINANTS[CMFS_NAME]['D65']
YCBCR_CHECK_MARKER = [0, 0, 0]

UNIVERSAL_COLOR_LIST = ["#F6AA00", "#FFF100", "#03AF7A",
                        "#005AFF", "#4DC4FF", "#804000"]

def preview_image(img, order='rgb', over_disp=False):
    if order == 'rgb':
        cv2.imshow('preview', img[:, :, ::-1])
    elif order == 'bgr':
        cv2.imshow('preview', img)
    elif order == 'mono':
        cv2.imshow('preview', img)
    else:
        raise ValueError("order parameter is invalid")

    if over_disp:
        cv2.resizeWindow('preview', )
    cv2.waitKey(0)
    cv2.destroyAllWindows()

def equal_devision(length, div_num):
    """
    # 概要
    length を div_num で分割する。
    端数が出た場合は誤差拡散法を使って上手い具合に分散させる。
    """
    base = length / div_num
    ret_array = [base for x in range(div_num)]

    # 誤差拡散法を使った辻褄合わせを適用
    # -------------------------------------------
    diff = 0
    for idx in range(div_num):
        diff += math.modf(ret_array[idx])[0]
        if diff >= 1.0:
            diff -= 1.0
            ret_array[idx] = int(math.floor(ret_array[idx]) + 1)
        else:
            ret_array[idx] = int(math.floor(ret_array[idx]))

    # 計算誤差により最終点が +1 されない場合への対処
    # -------------------------------------------
    diff = length - sum(ret_array)
    if diff != 0:
        ret_array[-1] += diff

    # 最終確認
    # -------------------------------------------
    if length != sum(ret_array):
        raise ValueError("the output of equal_division() is abnormal.")

    return ret_array

def do_matrix(img, mtx):
    """
    img に対して mtx を適用する。
    """
    base_shape = img.shape

    r, g, b = img[..., 0], img[..., 1], img[..., 2]
    ro = r * mtx[0][0] + g * mtx[0][1] + b * mtx[0][2]
    go = r * mtx[1][0] + g * mtx[1][1] + b * mtx[1][2]
    bo = r * mtx[2][0] + g * mtx[2][1] + b * mtx[2][2]

    out_img = np.dstack((ro, go, bo)).reshape(base_shape)

    return out_img

def _get_cmfs_xy():
    """
    xy色度図のプロットのための馬蹄形の外枠のxy値を求める。

    Returns
    -------
    array_like
        xy coordinate for chromaticity diagram

    """
    # 基本パラメータ設定
    # ------------------
    cmf = CMFS.get(CMFS_NAME)
    d65_white = D65_WHITE

    # 馬蹄形のxy値を算出
    # --------------------------
    cmf_xy = XYZ_to_xy(cmf.values, d65_white)

    return cmf_xy

def get_primaries(name='ITU-R BT.2020'):
    """
    prmary color の座標を求める

    Parameters
    ----------
    name : str
        a name of the color space.

    Returns
    -------
    array_like
        prmaries. [[rx, ry], [gx, gy], [bx, by], [rx, ry]]

    """
    primaries = RGB_COLOURSPACES[name].primaries
    primaries = np.append(primaries, [primaries[0, :]], axis=0)

    rgb = np.array([[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]])

    return primaries, rgb

def xy_to_rgb(xy, name='ITU-R BT.2020', normalize='maximum', specific=None):
    """
    xy値からRGB値を算出する。
    いい感じに正規化もしておく。

    Parameters
    ----------
    xy : array_like
        xy value.
    name : string
        color space name.
    normalize : string
        normalize method. You can select 'maximum', 'specific' or None.

    Returns
    -------
    array_like
        rgb value. the value is normalized.
    """
    illuminant_XYZ = D65_WHITE
    illuminant_RGB = D65_WHITE
    chromatic_adaptation_transform = 'CAT02'
    large_xyz_to_rgb_matrix = get_xyz_to_rgb_matrix(name)
    if normalize == 'specific':
        xyY = xy_to_xyY(xy)
        xyY[..., 2] = specific
        large_xyz = xyY_to_XYZ(xyY)
    else:
        large_xyz = xy_to_XYZ(xy)

    rgb = XYZ_to_RGB(large_xyz, illuminant_XYZ, illuminant_RGB,
                     large_xyz_to_rgb_matrix,
                     chromatic_adaptation_transform)

    """
    そのままだとビデオレベルが低かったりするので、
    各ドット毎にRGB値を正規化＆最大化する。必要であれば。
    """
    if normalize == 'maximum':
        rgb = normalise_maximum(rgb, axis=-1)
    else:
        if(np.sum(rgb > 1.0) > 0):
            print("warning: over flow has occured at xy_to_rgb")
        if(np.sum(rgb < 0.0) > 0):
            print("warning: under flow has occured at xy_to_rgb")
        rgb[rgb < 0] = 0
        rgb[rgb > 1.0] = 1.0

    return rgb

def get_white_point(name):
    """
    white point を求める。CIE1931ベース。
    """
    if name != "DCI-P3":
        illuminant = RGB_COLOURSPACES[name].illuminant
        white_point = ILLUMINANTS[CMFS_NAME][illuminant]
    else:
        white_point = ILLUMINANTS[CMFS_NAME]["D65"]

    return white_point

def get_secondaries(name='ITU-R BT.2020'):
    """
    secondary color の座標を求める

    Parameters
    ----------
    name : str
        a name of the color space.

    Returns
    -------
    array_like
        secondaries. the order is magenta, yellow, cyan.

    """
    secondary_rgb = np.array([[1.0, 0.0, 1.0],
                              [1.0, 1.0, 0.0],
                              [0.0, 1.0, 1.0]])
    illuminant_XYZ = D65_WHITE
    illuminant_RGB = D65_WHITE
    chromatic_adaptation_transform = 'CAT02'
    rgb_to_xyz_matrix = get_rgb_to_xyz_matrix(name)
    large_xyz = RGB_to_XYZ(secondary_rgb, illuminant_RGB,
                           illuminant_XYZ, rgb_to_xyz_matrix,
                           chromatic_adaptation_transform)

    xy = XYZ_to_xy(large_xyz, illuminant_XYZ)

    return xy, secondary_rgb.reshape((3, 3))

# def plot_chromaticity_diagram(
#         rate=480/755.0*2, xmin=0.0, xmax=0.8, ymin=0.0, ymax=0.9, **kwargs):
#     # キーワード引数の初期値設定
#     # ------------------------------------
#     monitor_primaries = kwargs.get('monitor_primaries', None)
#     secondaries = kwargs.get('secondaries', None)
#     test_scatter = kwargs.get('test_scatter', None)
#     intersection = kwargs.get('intersection', None)

#     # プロット用データ準備
#     # ---------------------------------
#     xy_image = get_chromaticity_image(
#         xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
#     cmf_xy = _get_cmfs_xy()

#     bt709_gamut, _ = get_primaries(name=cs.BT709)
#     bt2020_gamut, _ = get_primaries(name=cs.BT2020)
#     dci_p3_gamut, _ = get_primaries(name=cs.P3_D65)
#     ap0_gamut, _ = get_primaries(name=cs.ACES_AP0)
#     ap1_gamut, _ = get_primaries(name=cs.ACES_AP1)
#     xlim = (min(0, xmin), max(0.8, xmax))
#     ylim = (min(0, ymin), max(0.9, ymax))

#     ax1 = pu.plot_1_graph(fontsize=20 * rate,
#                           figsize=((xmax - xmin) * 10 * rate,
#                                    (ymax - ymin) * 10 * rate),
#                           graph_title="CIE1931 Chromaticity Diagram",
#                           graph_title_size=None,
#                           xlabel=None, ylabel=None,
#                           axis_label_size=None,
#                           legend_size=18 * rate,
#                           xlim=xlim, ylim=ylim,
#                           xtick=[x * 0.1 + xmin for x in
#                                  range(int((xlim[1] - xlim[0])/0.1) + 1)],
#                           ytick=[x * 0.1 + ymin for x in
#                                  range(int((ylim[1] - ylim[0])/0.1) + 1)],
#                           xtick_size=17 * rate,
#                           ytick_size=17 * rate,
#                           linewidth=4 * rate,
#                           minor_xtick_num=2,
#                           minor_ytick_num=2)
#     ax1.plot(cmf_xy[..., 0], cmf_xy[..., 1], '-k', lw=3.5*rate, label=None)
#     ax1.plot((cmf_xy[-1, 0], cmf_xy[0, 0]), (cmf_xy[-1, 1], cmf_xy[0, 1]),
#              '-k', lw=3.5*rate, label=None)
#     ax1.plot(bt709_gamut[:, 0], bt709_gamut[:, 1],
#              c=UNIVERSAL_COLOR_LIST[0], label="BT.709", lw=2.75*rate)
#     ax1.plot(bt2020_gamut[:, 0], bt2020_gamut[:, 1],
#              c=UNIVERSAL_COLOR_LIST[1], label="BT.2020", lw=2.75*rate)
#     ax1.plot(dci_p3_gamut[:, 0], dci_p3_gamut[:, 1],
#              c=UNIVERSAL_COLOR_LIST[2], label="DCI-P3", lw=2.75*rate)
#     ax1.plot(ap1_gamut[:, 0], ap1_gamut[:, 1],
#              c=UNIVERSAL_COLOR_LIST[3], label="ACES AP1", lw=2.75*rate)
#     ax1.plot(ap0_gamut[:, 0], ap0_gamut[:, 1],
#              c=UNIVERSAL_COLOR_LIST[4], label="ACES AP0", lw=2.75*rate)
#     if monitor_primaries is not None:
#         ax1.plot(monitor_primaries[:, 0], monitor_primaries[:, 1],
#                  c="#202020", label="???", lw=3*rate)
#     if secondaries is not None:
#         xy, rgb = secondaries
#         ax1.scatter(xy[..., 0], xy[..., 1], s=700*rate, marker='s', c=rgb,
#                     edgecolors='#404000', linewidth=2*rate)
#     if test_scatter is not None:
#         xy, rgb = test_scatter
#         ax1.scatter(xy[..., 0], xy[..., 1], s=300*rate, marker='s', c=rgb,
#                     edgecolors='#404040', linewidth=2*rate)
#     if intersection is not None:
#         ax1.scatter(intersection[..., 0], intersection[..., 1],
#                     s=300*rate, marker='s', c='#CCCCCC',
#                     edgecolors='#404040', linewidth=2*rate)

#     ax1.imshow(xy_image, extent=(xmin, xmax, ymin, ymax))
#     plt.legend(loc='upper right')
#     plt.savefig('temp_fig.png', bbox_inches='tight')
#     plt.show()

def get_chromaticity_image(samples=1024, antialiasing=True, bg_color=0.9,
                           xmin=0.0, xmax=1.0, ymin=0.0, ymax=1.0):
    """
    xy色度図の馬蹄形の画像を生成する

    Returns
    -------
    ndarray
        rgb image.
    """

    """
    色域設定。sRGBだと狭くて少し変だったのでBT.2020に設定。
    若干色が薄くなるのが難点。暇があれば改良したい。
    """
    # color_space = models.BT2020_COLOURSPACE
    # color_space = models.S_GAMUT3_COLOURSPACE
    color_space = models.ACES_CG_COLOURSPACE

    # 馬蹄形のxy値を算出
    # --------------------------
    cmf_xy = _get_cmfs_xy()

    """
    馬蹄の内外の判別をするために三角形で領域分割する(ドロネー図を作成)。
    ドロネー図を作れば後は外積計算で領域の内外を判別できる（たぶん）。

    なお、作成したドロネー図は以下のコードでプロット可能。
    1点補足しておくと、```plt.triplot``` の第三引数は、
    第一、第二引数から三角形を作成するための **インデックス** のリスト
    になっている。[[0, 1, 2], [2, 4, 3], ...]的な。

    ```python
    plt.figure()
    plt.triplot(xy[:, 0], xy[:, 1], triangulation.simplices.copy(), '-o')
    plt.title('triplot of Delaunay triangulation')
    plt.show()
    ```
    """
    triangulation = Delaunay(cmf_xy)

    """
    ```triangulation.find_simplex()``` で xy がどのインデックスの領域か
    調べることができる。戻り値が ```-1``` の場合は領域に含まれないため、
    0以下のリストで領域判定の mask を作ることができる。
    """
    xx, yy\
        = np.meshgrid(np.linspace(xmin, xmax, samples),
                      np.linspace(ymax, ymin, samples))
    xy = np.dstack((xx, yy))
    mask = (triangulation.find_simplex(xy) < 0).astype(np.float)

    # アンチエイリアシングしてアルファチャンネルを滑らかに
    # ------------------------------------------------
    if antialiasing:
        kernel = np.array([
            [0, 1, 0],
            [1, 2, 1],
            [0, 1, 0],
        ]).astype(np.float)
        kernel /= np.sum(kernel)
        mask = convolve(mask, kernel)

    # ネガポジ反転
    # --------------------------------
    mask = 1 - mask[:, :, np.newaxis]

    # xy のメッシュから色を復元
    # ------------------------
    illuminant_XYZ = D65_WHITE
    illuminant_RGB = color_space.whitepoint
    chromatic_adaptation_transform = 'XYZ Scaling'
    large_xyz_to_rgb_matrix = color_space.XYZ_to_RGB_matrix
    xy[xy == 0.0] = 1.0  # ゼロ割対策
    large_xyz = xy_to_XYZ(xy)
    rgb = XYZ_to_RGB(large_xyz, illuminant_XYZ, illuminant_RGB,
                     large_xyz_to_rgb_matrix,
                     chromatic_adaptation_transform)

    """
    そのままだとビデオレベルが低かったりするので、
    各ドット毎にRGB値を正規化＆最大化する。
    """
    rgb[rgb == 0] = 1.0  # ゼロ割対策
    rgb = normalise_maximum(rgb, axis=-1)

    # mask 適用
    # -------------------------------------
    mask_rgb = np.dstack((mask, mask, mask))
    rgb *= mask_rgb

    # 背景色をグレーに変更
    # -------------------------------------
    bg_rgb = np.ones_like(rgb)
    bg_rgb *= (1 - mask_rgb) * bg_color

    rgb += bg_rgb

    rgb = rgb ** (1/2.2)

    return rgb

def get_csf_color_image(width=640, height=480,
                        lv1=np.uint16(np.array([1.0, 1.0, 1.0]) * 1023 * 0x40),
                        lv2=np.uint16(np.array([1.0, 1.0, 1.0]) * 512 * 0x40),
                        stripe_num=18):
    """
    長方形を複数個ズラして重ねることでCSFパターンっぽいのを作る。
    入力信号レベルは16bitに限定する。

    Parameters
    ----------
    width : numeric.
        width of the pattern image.
    height : numeric.
        height of the pattern image.
    lv1 : numeric
        video level 1. this value must be 10bit.
    lv2 : numeric
        video level 2. this value must be 10bit.
    stripe_num : numeric
        number of the stripe.

    Returns
    -------
    array_like
        a cms pattern image.
    """
    width_list = equal_devision(width, stripe_num)
    height_list = equal_devision(height, stripe_num)
    h_pos_list = equal_devision(width // 2, stripe_num)
    v_pos_list = equal_devision(height // 2, stripe_num)
    lv1_16bit = lv1
    lv2_16bit = lv2
    img = np.zeros((height, width, 3), dtype=np.uint16)

    width_temp = width
    height_temp = height
    h_pos_temp = 0
    v_pos_temp = 0
    for idx in range(stripe_num):
        lv = lv1_16bit if (idx % 2) == 0 else lv2_16bit
        temp_img = np.ones((height_temp, width_temp, 3), dtype=np.uint16)
        # temp_img *= lv
        temp_img[:, :] = lv
        ed_pos_h = h_pos_temp + width_temp
        ed_pos_v = v_pos_temp + height_temp
        img[v_pos_temp:ed_pos_v, h_pos_temp:ed_pos_h] = temp_img
        width_temp -= width_list[stripe_num - 1 - idx]
        height_temp -= height_list[stripe_num - 1 - idx]
        h_pos_temp += h_pos_list[idx]
        v_pos_temp += v_pos_list[idx]

    return img

def plot_xyY_color_space(name='ITU-R BT.2020', samples=1024,
                         antialiasing=True):
    """
    SONY の HDR説明資料にあるような xyY の図を作る。

    Parameters
    ----------
    name : str
        name of the target color space.

    Returns
    -------
    None

    """

    # 馬蹄の領域判別用データ作成
    # --------------------------
    primary_xy, _ = get_primaries(name=name)
    triangulation = Delaunay(primary_xy)

    xx, yy\
        = np.meshgrid(np.linspace(0, 1, samples), np.linspace(1, 0, samples))
    xy = np.dstack((xx, yy))
    mask = (triangulation.find_simplex(xy) < 0).astype(np.float)

    # アンチエイリアシングしてアルファチャンネルを滑らかに
    # ------------------------------------------------
    if antialiasing:
        kernel = np.array([
            [0, 1, 0],
            [1, 2, 1],
            [0, 1, 0],
        ]).astype(np.float)
        kernel /= np.sum(kernel)
        mask = convolve(mask, kernel)

    # ネガポジ反転
    # --------------------------------
    mask = 1 - mask[:, :, np.newaxis]

    # xy のメッシュから色を復元
    # ------------------------
    illuminant_XYZ = D65_WHITE
    illuminant_RGB = RGB_COLOURSPACES[name].whitepoint
    chromatic_adaptation_transform = 'CAT02'
    large_xyz_to_rgb_matrix = get_xyz_to_rgb_matrix(name)
    rgb_to_large_xyz_matrix = get_rgb_to_xyz_matrix(name)
    large_xyz = xy_to_XYZ(xy)
    rgb = XYZ_to_RGB(large_xyz, illuminant_XYZ, illuminant_RGB,
                     large_xyz_to_rgb_matrix,
                     chromatic_adaptation_transform)

    """
    そのままだとビデオレベルが低かったりするので、
    各ドット毎にRGB値を正規化＆最大化する。
    """
    rgb_org = normalise_maximum(rgb, axis=-1)

    # mask 適用
    # -------------------------------------
    mask_rgb = np.dstack((mask, mask, mask))
    rgb = rgb_org * mask_rgb
    rgba = np.dstack((rgb, mask))

    # こっからもういちど XYZ に変換。Yを求めるために。
    # ---------------------------------------------
    large_xyz2 = RGB_to_XYZ(rgb, illuminant_RGB, illuminant_XYZ,
                            rgb_to_large_xyz_matrix,
                            chromatic_adaptation_transform)

    # ログスケールに変換する準備
    # --------------------------
    large_y = large_xyz2[..., 1] * 1000
    large_y[large_y < 1] = 1.0

    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d')
    # ax.plot_wireframe(xy[..., 0], xy[..., 1], np.log10(large_y),
    #                   rcount=100, ccount=100)
    ax.plot_surface(xy[..., 0], xy[..., 1], np.log10(large_y),
                    rcount=64, ccount=64, facecolors=rgb_org)
    ax.set_xlabel("x")
    ax.set_ylabel("y")
    ax.set_zlabel("Y")
    ax.set_zticks([0, 1, 2, 3])
    ax.set_zticklabels([1, 10, 100, 1000])

    # chromatcity_image の取得。z=0 の位置に貼り付ける
    # ----------------------------------------------
    cie1931_rgb = get_chromaticity_image(samples=samples, bg_color=0.0)

    alpha = np.zeros_like(cie1931_rgb[..., 0])
    rgb_sum = np.sum(cie1931_rgb, axis=-1)
    alpha[rgb_sum > 0.00001] = 1
    cie1931_rgb = np.dstack((cie1931_rgb[..., 0], cie1931_rgb[..., 1],
                             cie1931_rgb[..., 2], alpha))
    zz = np.zeros_like(xy[..., 0])
    ax.plot_surface(xy[..., 0], xy[..., 1], zz,
                    facecolors=cie1931_rgb)

    plt.show()

def log_tick_formatter(val, pos=None):
    return "{:.0e}".format(10**val)

def get_3d_grid_cube_format(grid_num=4):
    """
    # 概要
    (0, 0, 0), (1, 0, 0), (0, 1, 0), (1, 1, 0), (0, 0, 1), ...
    みたいな配列を返す。
    CUBE形式の3DLUTを作成する時に便利。
    """

    base = np.linspace(0, 1, grid_num)
    ones_x = np.ones((grid_num, grid_num, 1))
    ones_y = np.ones((grid_num, 1, grid_num))
    ones_z = np.ones((1, grid_num, grid_num))
    r_3d = base[np.newaxis, np.newaxis, :] * ones_x
    g_3d = base[np.newaxis, :, np.newaxis] * ones_y
    b_3d = base[:, np.newaxis, np.newaxis] * ones_z
    r_3d = r_3d.flatten()
    g_3d = g_3d.flatten()
    b_3d = b_3d.flatten()

    return np.dstack((r_3d, g_3d, b_3d))

def quadratic_bezier_curve(t, p0, p1, p2, samples=1024):
    # x = ((1 - t) ** 2) * p0[0] + 2 * (1 - t) * t * p1[0]\
    #     + (t ** 2) * p2[0]
    # y = ((1 - t) ** 2) * p0[1] + 2 * (1 - t) * t * p1[1]\
    #     + (t ** 2) * p2[1]

    x = ((1 - t) ** 2) * p0[0] + 2 * (1 - t) * t * p1[0]\
        + (t ** 2) * p2[0]
    y = ((1 - t) ** 2) * p0[1] + 2 * (1 - t) * t * p1[1]\
        + (t ** 2) * p2[1]

    # ax1 = pu.plot_1_graph(fontsize=20,
    #                       figsize=(10, 8),
    #                       graph_title="Title",
    #                       graph_title_size=None,
    #                       xlabel="X Axis Label", ylabel="Y Axis Label",
    #                       axis_label_size=None,
    #                       legend_size=17,
    #                       xlim=None,
    #                       ylim=None,
    #                       xtick=None,
    #                       ytick=None,
    #                       xtick_size=None, ytick_size=None,
    #                       linewidth=3,
    #                       minor_xtick_num=None,
    #                       minor_ytick_num=None)
    # ax1.plot(x, y, label='aaa')
    # plt.legend(loc='upper left')
    # plt.show()

def gen_step_gradation(width=1024, height=128, step_num=17,
                       bit_depth=10, color=(1.0, 1.0, 1.0),
                       direction='h', debug=False):
    """
    # 概要
    階段状に変化するグラデーションパターンを作る。
    なお、引数の調整により正確に1階調ずつ変化するパターンも作成可能。

    # 注意事項
    正確に1階調ずつ変化するグラデーションを作る場合は
    ```step_num = (2 ** bit_depth) + 1```
    となるようにパラメータを指定すること。具体例は以下のExample参照。

    # Example
    ```
    grad_8 = gen_step_gradation(width=grad_width, height=grad_height,
                                step_num=257, bit_depth=8,
                                color=(1.0, 1.0, 1.0), direction='h')

    grad_10 = gen_step_gradation(width=grad_width, height=grad_height,
                                 step_num=1025, bit_depth=10,
                                 color=(1.0, 1.0, 1.0), direction='h')
    ```
    """
    max = 2 ** bit_depth

    # グラデーション方向設定
    # ----------------------
    if direction == 'h':
        pass
    else:
        temp = height
        height = width
        width = temp

    if (max + 1 != step_num):
        """
        1階調ずつの増加では無いパターン。
        末尾のデータが 256 や 1024 になるため -1 する。
        """
        val_list = np.linspace(0, max, step_num)
        val_list[-1] -= 1
    else:
        """
        正確に1階調ずつ変化するパターン。
        末尾のデータが 256 や 1024 になるため除外する。
        """
        val_list = np.linspace(0, max, step_num)[0:-1]
        step_num -= 1  # step_num は 引数で余計に +1 されてるので引く

        # 念のため1階調ずつの変化か確認
        # ---------------------------
        diff = val_list[1:] - val_list[0:-1]
        if (diff == 1).all():
            pass
        else:
            raise ValueError("calculated value is invalid.")

    # まずは水平1LINEのグラデーションを作る
    # -----------------------------------
    step_length_list = equal_devision(width, step_num)
    step_bar_list = []
    for step_idx, length in enumerate(step_length_list):
        step = [np.ones((length)) * color[c_idx] * val_list[step_idx]
                for c_idx in range(3)]
        if direction == 'h':
            step = np.dstack(step)
            step_bar_list.append(step)
            step_bar = np.hstack(step_bar_list)
        else:
            step = np.dstack(step).reshape((length, 1, 3))
            step_bar_list.append(step)
            step_bar = np.vstack(step_bar_list)

    # ブロードキャストを利用して2次元に拡張する
    # ------------------------------------------
    if direction == 'h':
        img = step_bar * np.ones((height, 1, 3))
    else:
        img = step_bar * np.ones((1, height, 3))

    # np.uint16 にコンバート
    # ------------------------------
    # img = np.uint16(np.round(img * (2 ** (16 - bit_depth))))

    if debug:
        preview_image(img, 'rgb')

    return img

def merge(img_a, img_b, pos=(0, 0)):
    """
    img_a に img_b をマージする。
    img_a にデータを上書きする。

    pos = (horizontal_st, vertical_st)
    """
    b_width = img_b.shape[1]
    b_height = img_b.shape[0]

    img_a[pos[1]:b_height+pos[1], pos[0]:b_width+pos[0]] = img_b

def merge_with_alpha(bg_img, fg_img, tf_str=tf.SRGB, pos=(0, 0)):
    """
    合成する。

    Parameters
    ----------
    bg_img : array_like(float, 3-channel)
        image data.
    fg_img : array_like(float, 4-channel)
        image data
    tf : strings
        transfer function
    pos : list(int)
        (pos_h, pos_v)
    """
    f_width = fg_img.shape[1]
    f_height = fg_img.shape[0]

    bg_merge_area = bg_img[pos[1]:f_height+pos[1], pos[0]:f_width+pos[0]]
    bg_linear = tf.eotf_to_luminance(bg_merge_area, tf_str)
    fg_linear = tf.eotf_to_luminance(fg_img, tf_str)
    alpha = fg_linear[:, :, 3:] / tf.PEAK_LUMINANCE[tf_str]

    out_linear = (1 - alpha) * bg_linear + fg_linear[:, :, :-1]
    out_merge_area = tf.oetf_from_luminance(out_linear, tf_str)
    bg_img[pos[1]:f_height+pos[1], pos[0]:f_width+pos[0]] = out_merge_area

    return bg_img

def dot_pattern(dot_size=4, repeat=4, color=np.array([1.0, 1.0, 1.0])):
    """
    dot pattern 作る。

    Parameters
    ----------
    dot_size : integer
        dot size.
    repeat : integer
        The number of high-low pairs.
    color : array_like
        color value.

    Returns
    -------
    array_like
        dot pattern image.

    """
    # 水平・垂直のピクセル数
    pixel_num = dot_size * 2 * repeat

    # High-Log の 論理配列を生成
    even_logic = [(np.arange(pixel_num) % (dot_size * 2)) - dot_size < 0]
    even_logic = np.dstack((even_logic, even_logic, even_logic))
    odd_logic = np.logical_not(even_logic)

    # 着色
    color = color.reshape((1, 1, 3))
    even_line = (np.ones((1, pixel_num, 3)) * even_logic) * color
    odd_line = (np.ones((1, pixel_num, 3)) * odd_logic) * color

    # V方向にコピー＆Even-Oddの結合
    even_block = np.repeat(even_line, dot_size, axis=0)
    odd_block = np.repeat(odd_line, dot_size, axis=0)
    pair_block = np.vstack((even_block, odd_block))

    img = np.vstack([pair_block for x in range(repeat)])

    return img

def complex_dot_pattern(kind_num=3, whole_repeat=2,
                        fg_color=np.array([1.0, 1.0, 1.0]),
                        bg_color=np.array([0.15, 0.15, 0.15])):
    """
    dot pattern 作る。

    Parameters
    ----------
    kind_num : integer
        作成するドットサイズの種類。
        例えば、kind_num=3 ならば、1dot, 2dot, 4dot のパターンを作成。
    whole_repeat : integer
        異なる複数種類のドットパターンの組数。
        例えば、kind_num=3, whole_repeat=2 ならば、
        1dot, 2dot, 4dot のパターンを水平・垂直に2組作る。
    fg_color : array_like
        foreground color value.
    bg_color : array_like
        background color value.
    reduce : bool
        HDRテストパターンの3840x2160専用。縦横を半分にする。

    Returns
    -------
    array_like
        dot pattern image.

    """
    max_dot_width = 2 ** kind_num
    img_list = []
    for size_idx in range(kind_num)[::-1]:
        dot_size = 2 ** size_idx
        repeat = max_dot_width // dot_size
        dot_img = dot_pattern(dot_size, repeat, fg_color)
        img_list.append(dot_img)
        img_list.append(np.ones_like(dot_img) * bg_color)
        # preview_image(dot_img)

    line_upper_img = np.hstack(img_list)
    line_upper_img = np.hstack([line_upper_img for x in range(whole_repeat)])
    line_lower_img = line_upper_img.copy()[:, ::-1, :]
    h_unit_img = np.vstack((line_upper_img, line_lower_img))

    img = np.vstack([h_unit_img for x in range(kind_num * whole_repeat)])
    # preview_image(img)
    # cv2.imwrite("hoge.tiff", np.uint8(img * 0xFF)[..., ::-1])

    return img

def make_csf_color_image(width=640, height=640,
                         lv1=np.array([940, 940, 940], dtype=np.uint16),
                         lv2=np.array([1023, 1023, 1023], dtype=np.uint16),
                         stripe_num=6):
    """
    長方形を複数個ズラして重ねることでCSFパターンっぽいのを作る。
    入力信号レベルは10bitに限定する。

    Parameters
    ----------
    width : numeric.
        width of the pattern image.
    height : numeric.
        height of the pattern image.
    lv1 : array_like
        video level 1. this value must be 10bit.
    lv2 : array_like
        video level 2. this value must be 10bit.
    stripe_num : numeric
        number of the stripe.

    Returns
    -------
    array_like
        a cms pattern image.
    """
    width_list = equal_devision(width, stripe_num)
    height_list = equal_devision(height, stripe_num)
    h_pos_list = equal_devision(width // 2, stripe_num)
    v_pos_list = equal_devision(height // 2, stripe_num)
    img = np.zeros((height, width, 3), dtype=np.uint16)

    width_temp = width
    height_temp = height
    h_pos_temp = 0
    v_pos_temp = 0
    for idx in range(stripe_num):
        lv = lv1 if (idx % 2) == 0 else lv2
        temp_img = np.ones((height_temp, width_temp, 3), dtype=np.uint16)
        temp_img = temp_img * lv.reshape((1, 1, 3))
        ed_pos_h = h_pos_temp + width_temp
        ed_pos_v = v_pos_temp + height_temp
        img[v_pos_temp:ed_pos_v, h_pos_temp:ed_pos_h] = temp_img
        width_temp -= width_list[stripe_num - 1 - idx]
        height_temp -= height_list[stripe_num - 1 - idx]
        h_pos_temp += h_pos_list[idx]
        v_pos_temp += v_pos_list[idx]

    # preview_image(img / 1023)

    return img

def make_tile_pattern(width=480, height=960, h_tile_num=4,
                      v_tile_num=4, low_level=(940, 940, 940),
                      high_level=(1023, 1023, 1023)):
    """
    タイル状の縞々パターンを作る
    """
    width_array = equal_devision(width, h_tile_num)
    height_array = equal_devision(height, v_tile_num)
    high_level = np.array(high_level, dtype=np.uint16)
    low_level = np.array(low_level, dtype=np.uint16)

    v_buf = []

    for v_idx, height in enumerate(height_array):
        h_buf = []
        for h_idx, width in enumerate(width_array):
            tile_judge = (h_idx + v_idx) % 2 == 0
            h_temp = np.zeros((height, width, 3), dtype=np.uint16)
            h_temp[:, :] = high_level if tile_judge else low_level
            h_buf.append(h_temp)

        v_buf.append(np.hstack(h_buf))
    img = np.vstack(v_buf)
    # preview_image(img/1024.0)
    return img

def get_marker_idx(img, marker_value):
    return np.all(img == marker_value, axis=-1)

def make_ycbcr_checker(height=480, v_tile_num=4):
    """
    YCbCr係数誤りを確認するテストパターンを作る。
    正直かなり汚い組み方です。雑に作ったパターンを悪魔合体させています。

    Parameters
    ----------
    height : numeric.
        height of the pattern image.
    v_tile_num : numeric
        number of the tile in the vertical direction.

    Note
    ----
    横長のパターンになる。以下の式が成立する。

    ```
    h_tile_num = v_tile_num * 2
    width = height * 2
    ```

    Returns
    -------
    array_like
        ycbcr checker image
    """

    cyan_img = make_tile_pattern(width=height, height=height,
                                 h_tile_num=v_tile_num,
                                 v_tile_num=v_tile_num,
                                 low_level=[0, 990, 990],
                                 high_level=[0, 1023, 1023])
    magenta_img = make_tile_pattern(width=height, height=height,
                                    h_tile_num=v_tile_num,
                                    v_tile_num=v_tile_num,
                                    low_level=[990, 0, 312],
                                    high_level=[1023, 0, 312])

    out_img = np.hstack([cyan_img, magenta_img])

    # preview_image(out_img/1023.0)

    return out_img

def plot_color_checker_image(rgb, rgb2=None, size=(1920, 1080),
                             block_size=1/4.5, padding=0.01):
    """
    ColorCheckerをプロットする

    Parameters
    ----------
    rgb : array_like
        RGB value of the ColorChecker.
        RGB's shape must be (24, 3).
    rgb2 : array_like
        It's a optional parameter.
        If You want to draw two different ColorCheckers,
        set the RGB value to this variable.
    size : tuple
        canvas size.
    block_size : float
        A each block's size.
        This value is ratio to height of the canvas.
    padding : float
        A padding to the block.

    Returns
    -------
    array_like
        A ColorChecker image.

    """
    IMG_HEIGHT = size[1]
    IMG_WIDTH = size[0]
    COLOR_CHECKER_SIZE = block_size
    COLOR_CHECKER_H_NUM = 6
    COLOR_CHECKER_V_NUM = 4
    COLOR_CHECKER_PADDING = 0.01
    # 基本パラメータ算出
    # --------------------------------------
    COLOR_CHECKER_H_NUM = 6
    COLOR_CHECKER_V_NUM = 4
    img_height = IMG_HEIGHT
    img_width = IMG_WIDTH
    patch_st_h = int(IMG_WIDTH / 2.0
                     - (IMG_HEIGHT * COLOR_CHECKER_SIZE
                        * COLOR_CHECKER_H_NUM / 2.0
                        + (IMG_HEIGHT * COLOR_CHECKER_PADDING
                           * (COLOR_CHECKER_H_NUM / 2.0 - 0.5)) / 2.0))
    patch_st_v = int(IMG_HEIGHT / 2.0
                     - (IMG_HEIGHT * COLOR_CHECKER_SIZE
                        * COLOR_CHECKER_V_NUM / 2.0
                        + (IMG_HEIGHT * COLOR_CHECKER_PADDING
                           * (COLOR_CHECKER_V_NUM / 2.0 - 0.5)) / 2.0))
    patch_width = int(img_height * COLOR_CHECKER_SIZE)
    patch_height = patch_width
    patch_space = int(img_height * COLOR_CHECKER_PADDING)

    # 24ループで1枚の画像に24パッチを描画
    # -------------------------------------------------
    img_all_patch = np.zeros((img_height, img_width, 3), dtype=np.uint8)
    for idx in range(COLOR_CHECKER_H_NUM * COLOR_CHECKER_V_NUM):
        v_idx = idx // COLOR_CHECKER_H_NUM
        h_idx = (idx % COLOR_CHECKER_H_NUM)
        patch = np.ones((patch_height, patch_width, 3))
        patch[:, :] = rgb[idx]
        st_h = patch_st_h + (patch_width + patch_space) * h_idx
        st_v = patch_st_v + (patch_height + patch_space) * v_idx
        img_all_patch[st_v:st_v+patch_height, st_h:st_h+patch_width] = patch

        # pt1 = (st_h, st_v)  # upper left
        pt2 = (st_h + patch_width, st_v)  # upper right
        pt3 = (st_h, st_v + patch_height)  # lower left
        pt4 = (st_h + patch_width, st_v + patch_height)  # lower right
        pts = np.array((pt2, pt3, pt4))
        sub_color = rgb[idx].tolist() if rgb2 is None else rgb2[idx].tolist()
        cv2.fillPoly(img_all_patch, [pts], sub_color)

    preview_image(img_all_patch)

    return img_all_patch

def get_log10_x_scale(
        sample_num=8, ref_val=1.0, min_exposure=-1, max_exposure=6):
    """
    Log10スケールのx軸データを作る。

    Examples
    --------
    >>> get_log2_x_scale(
    ...     sample_num=8, ref_val=1.0, min_exposure=-1, max_exposure=6)
    array([  1.0000e-01   1.0000e+00   1.0000e+01   1.0000e+02
             1.0000e+03   1.0000e+04   1.0000e+05   1.0000e+06])
    """
    x_min = np.log10(ref_val * (10 ** min_exposure))
    x_max = np.log10(ref_val * (10 ** max_exposure))
    x = np.linspace(x_min, x_max, sample_num)

    return 10.0 ** x

def get_log2_x_scale(
        sample_num=32, ref_val=1.0, min_exposure=-6.5, max_exposure=6.5):
    """
    Log2スケールのx軸データを作る。

    Examples
    --------
    >>> get_log2_x_scale(sample_num=10, min_exposure=-4.0, max_exposure=4.0)
    array([[  0.0625       0.11573434   0.214311     0.39685026   0.73486725
              1.36079      2.5198421    4.66611616   8.64047791  16.        ]])
    """
    x_min = np.log2(ref_val * (2 ** min_exposure))
    x_max = np.log2(ref_val * (2 ** max_exposure))
    x = np.linspace(x_min, x_max, sample_num)

    return 2.0 ** x

def shaper_func_linear_to_log2(
        x, mid_gray=0.18, min_exposure=-6.5, max_exposure=6.5):
    """
    ACESutil.Lin_to_Log2_param.ctl を参考に作成。
    https://github.com/ampas/aces-dev/blob/master/transforms/ctl/utilities/ACESutil.Lin_to_Log2_param.ctl

    Parameters
    ----------
    x : array_like
        linear data.
    mid_gray : float
        18% gray value on linear scale.
    min_exposure : float
        minimum value on log scale.
    max_exposure : float
        maximum value on log scale.

    Returns
    -------
    array_like
        log2 value that is transformed from linear x value.

    Examples
    --------
    >>> shaper_func_linear_to_log2(
    ...     x=0.18, mid_gray=0.18, min_exposure=-6.5, max_exposure=6.5)
    0.5
    >>> shaper_func_linear_to_log2(
    ...     x=np.array([0.00198873782209, 16.2917402385])
    ...     mid_gray=0.18, min_exposure=-6.5, max_exposure=6.5)
    array([  1.58232402e-13   1.00000000e+00])
    """
    # log2空間への変換。mid_gray が 0.0 となるように補正
    y = np.log2(x / mid_gray)

    # min, max の範囲で正規化。
    y_normalized = (y - min_exposure) / (max_exposure - min_exposure)

    y_normalized[y_normalized < 0] = 0

    return y_normalized

def shaper_func_log2_to_linear(
        x, mid_gray=0.18, min_exposure=-6.5, max_exposure=6.5):
    """
    ACESutil.Log2_to_Lin_param.ctl を参考に作成。
    https://github.com/ampas/aces-dev/blob/master/transforms/ctl/utilities/ACESutil.Log2_to_Lin_param.ctl

    Log2空間の補足は shaper_func_linear_to_log2() の説明を参照

    Examples
    --------
    >>> x = np.array([0.0, 1.0])
    >>> shaper_func_log2_to_linear(
    ...     x, mid_gray=0.18, min_exposure=-6.5, max_exposure=6.5)
    array([0.00198873782209, 16.2917402385])
    """
    x_re_scale = x * (max_exposure - min_exposure) + min_exposure
    y = (2.0 ** x_re_scale) * mid_gray
    # plt.plot(x, y)
    # plt.show()

    return y

def draw_straight_line(img, pt1, pt2, color, thickness):
    """
    直線を引く。OpenCV だと 8bit しか対応してないっぽいので自作。

    Parameters
    ----------
    img : array_like
        image data.
    pt1 : list(pos_h, pos_v)
        start point.
    pt2 : list(pos_h, pos_v)
        end point.
    color : array_like
        color
    thickness : int
        thickness.

    Returns
    -------
    array_like
        image data with line.

    Notes
    -----
    thickness のパラメータは pt1 の点から右下方向に効きます。
    pt1 を中心として太さではない事に注意。

    Examples
    --------
    >>> pt1 = (0, 0)
    >>> pt2 = (1920, 0)
    >>> color = (940, 940, 940)
    >>> thickness = 4
    >>> draw_straight_line(img, pt1, pt2, color, thickness)
    """
    # parameter check
    if (pt1[0] != pt2[0]) and (pt1[1] != pt2[1]):
        raise ValueError("invalid pt1, pt2 parameters")

    # check direction
    if pt1[0] == pt2[0]:
        thickness_direction = 'h'
    else:
        thickness_direction = 'v'

    if thickness_direction == 'h':
        for h_idx in range(thickness):
            img[pt1[1]:pt2[1], pt1[0] + h_idx, :] = color

    elif thickness_direction == 'v':
        for v_idx in range(thickness):
            img[pt1[1] + v_idx, pt1[0]:pt2[0], :] = color

def draw_outline(img, fg_color, outline_width):
    """
    img に対して外枠線を引く

    Parameters
    ----------
    img : array_like
        image data.
    fg_color : array_like
        color
    outline_width : int
        thickness.

    Returns
    -------
    array_like
        image data with line.

    Examples
    --------
    >>> img = np.zeros((1080, 1920, 3))
    >>> color = (940, 940, 940)
    >>> thickness = 2
    >>> draw_outline(img, color, thickness)
    """
    width = img.shape[1]
    height = img.shape[0]
    # upper left
    pt1 = (0, 0)
    pt2 = (width, 0)
    draw_straight_line(
        img, pt1, pt2, fg_color, outline_width)
    pt1 = (0, 0)
    pt2 = (0, height)
    draw_straight_line(
        img, pt1, pt2, fg_color, outline_width)
    # lower right
    pt1 = (width - outline_width, 0)
    pt2 = (width - outline_width, height)
    draw_straight_line(
        img, pt1, pt2, fg_color, outline_width)
    pt1 = (0, height - outline_width)
    pt2 = (width, height - outline_width)
    draw_straight_line(
        img, pt1, pt2, fg_color, outline_width)

def convert_luminance_to_color_value(luminance, transfer_function):
    """
    輝度[cd/m2] から code value の RGB値に変換する。
    luminance の単位は [cd/m2]。無彩色である。

    Examples
    --------
    >>> convert_luminance_to_color_value(100, tf.GAMMA24)
    >>> [ 1.0  1.0  1.0 ]
    >>> convert_luminance_to_color_value(100, tf.ST2084)
    >>> [ 0.50807842  0.50807842  0.50807842 ]
    """
    code_value = convert_luminance_to_code_value(
        luminance, transfer_function)
    return np.array([code_value, code_value, code_value])

def convert_luminance_to_code_value(luminance, transfer_function):
    """
    輝度[cd/m2] から code value に変換する。
    luminance の単位は [cd/m2]
    """
    return tf.oetf_from_luminance(luminance, transfer_function)

def calc_rad_patch_idx2(outmost_num=5, current_num=3):
    """
    以下のような、中心がGray、周りは CIELAB 空間の a*b*平面のカラーパッチの
    RGB値のリストを得る。
    https://user-images.githubusercontent.com/3609012/75444470-d3bc5600-59a6-11ea-962b-c315648782a9.png

    得られたデータは並べ替えが済んでいないため、calc_rad_patch_idx2() で
    得られる変換テーブルを使った変換が必要。
    本関数はまさにその変換を行う。
    """
    base = np.arange(outmost_num ** 2).reshape((outmost_num, outmost_num))
    # print(base)
    t_idx = (outmost_num - current_num) // 2
    trimmed = base[t_idx:t_idx+current_num, t_idx:t_idx+current_num]
    # print(trimmed)
    # print(np.arange(current_num**2).reshape((current_num, current_num)))

    half_num = current_num // 2
    conv_idx = []
    for idx in range(half_num):
        val = (current_num ** 2) // 2 + half_num - current_num * idx
        conv_idx.append(val)
    for idx in range(current_num)[::-1]:
        conv_idx.append(idx)
    for idx in range(1, current_num - 1):
        conv_idx.append(idx * current_num)
    for idx in range(current_num):
        val = (current_num ** 2) - current_num + idx
        conv_idx.append(val)
    for idx in range(1, half_num):
        val = (current_num ** 2) - 1 - idx * current_num
        conv_idx.append(val)

    conv_idx = trimmed.flatten()[conv_idx]

    return conv_idx

def _calc_rgb_from_same_lstar_radial_data(
        lstar, temp_chroma, current_num, color_space):
    """
    以下のような、中心がGray、周りは CIELAB 空間の a*b*平面のカラーパッチの
    RGB値のリストを得る。
    https://user-images.githubusercontent.com/3609012/75444470-d3bc5600-59a6-11ea-962b-c315648782a9.png

    得られたデータは並べ替えが済んでいないため、calc_rad_patch_idx2() で
    得られる変換テーブルを使った変換が必要。
    """
    current_patch_num = (current_num - 1) * 4 if current_num > 1 else 1
    rad = np.linspace(0, 2 * np.pi, current_patch_num, endpoint=False)
    ll = np.ones((current_patch_num)) * lstar
    aa = np.cos(rad) * temp_chroma
    bb = np.sin(rad) * temp_chroma
    lab = np.dstack((ll, aa, bb))
    large_xyz = Lab_to_XYZ(lab)
    rgb = XYZ_to_RGB(large_xyz, D65_WHITE, D65_WHITE,
                     color_space.XYZ_to_RGB_matrix)

    return np.clip(rgb, 0.0, 1.0)

def calc_same_lstar_radial_color_patch_data(
        lstar=58, chroma=32.5, outmost_num=9,
        color_space=BT709_COLOURSPACE,
        transfer_function=tf.GAMMA24):
    """
    以下のような、中心がGray、周りは CIELAB 空間の a*b*平面のカラーパッチの
    RGB値のリストを得る。
    https://user-images.githubusercontent.com/3609012/75444470-d3bc5600-59a6-11ea-962b-c315648782a9.png

    得られた RGB値のリストは最初のデータが画像左上の緑データ、
    最後のデータが画像右下の紫データとなるよう既に**並べ替え**が行われている。

    よってパッチをプロットする場合はRGB値リストの先頭から順にデータを取り出し、
    右下に向かって並べていけば良い。
    """
    patch_num = outmost_num ** 2
    transfer_function = tf.GAMMA24
    rgb_list = np.ones((patch_num, 3))

    current_num_list = range(1, outmost_num + 1, 2)
    chroma_list = np.linspace(0, chroma, len(current_num_list))
    for temp_chroma, current_num in zip(chroma_list, current_num_list):
        current_patch_num = (current_num - 1) * 4 if current_num > 1 else 1
        rgb = _calc_rgb_from_same_lstar_radial_data(
            lstar, temp_chroma, current_num, color_space)
        rgb = np.reshape(rgb, (current_patch_num, 3))
        rgb = tf.oetf(rgb, transfer_function)
        conv_idx = calc_rad_patch_idx2(
            outmost_num=outmost_num, current_num=current_num)
        for idx in range(current_patch_num):
            rgb_list[conv_idx[idx]] = rgb[idx]

    return rgb_list

def _plot_same_lstar_radial_color_patch_data(
        lstar=58, chroma=32.5, outmost_num=9,
        color_space=BT709_COLOURSPACE,
        transfer_function=tf.GAMMA24):
    patch_size = 1080 // outmost_num
    img = np.ones((1080, 1080, 3)) * 0.0
    rgb = calc_same_lstar_radial_color_patch_data(
        lstar=lstar, chroma=chroma, outmost_num=outmost_num,
        color_space=color_space, transfer_function=transfer_function)

    for idx in range(outmost_num ** 2):
        h_idx = idx % outmost_num
        v_idx = idx // outmost_num
        st_pos = (h_idx * patch_size, v_idx * patch_size)
        temp_img = np.ones((patch_size, patch_size, 3))\
            * rgb[idx][np.newaxis, np.newaxis, :]
        merge(img, temp_img, st_pos)

    cv2.imwrite("hoge2.tiff", np.uint16(np.round(img[:, :, ::-1] * 0xFFFF)))

def get_accelerated_x_1x(sample_num=64):
    """
    単調増加ではなく、加速度が 0→1→0 となるような x を作る

    Parameters
    ----------
    sample_num : int
        the number of the sample.

    Returns
    -------
    array_like
        accelerated value list

    Examples
    --------
    >>> x0 = np.linspace(0, 1, 8)
    >>> x1 = get_accelerated_x_1x(8)
    >>> print(x0)
    >>> [ 0.  0.142  0.285  0.428  0.571  0.714  0.857  1. ]
    >>> print(x1)
    >>> [ 0.  0.049  0.188  0.388  0.611  0.811  0.950  1. ]
    """
    rad = np.linspace(-0.5 * np.pi, 0.5 * np.pi, sample_num)
    x = (np.sin(rad) + 1) / 2

    return x

def get_accelerated_x_2x(sample_num=64):
    """
    単調増加ではなく、加速度が 0→1→0 となるような x を作る。
    加速度が `get_accelerated_x_1x` の2倍！！

    Parameters
    ----------
    sample_num : int
        the number of the sample.

    Returns
    -------
    array_like
        accelerated value list

    Examples
    --------
    >>> x0 = np.linspace(0, 1, 8)
    >>> x2 = get_accelerated_x_2x(8)
    >>> print(x0)
    >>> [ 0.  0.142  0.285  0.428  0.571  0.714  0.857  1. ]
    >>> print(x2)
    >>> [ 0.  0.006  0.084  0.328  0.671  0.915  0.993  1. ]
    """
    rad = np.linspace(-0.5 * np.pi, 0.5 * np.pi, sample_num)
    rad = np.sin(rad) * 0.5 * np.pi
    x = (np.sin(rad) + 1) / 2

    return x

def get_accelerated_x_4x(sample_num=64):
    """
    単調増加ではなく、加速度が 0→1→0 となるような x を作る。
    加速度が `get_accelerated_x_1x` の4倍！！

    Parameters
    ----------
    sample_num : int
        the number of the sample.

    Returns
    -------
    array_like
        accelerated value list
    """
    rad = np.linspace(-0.5 * np.pi, 0.5 * np.pi, sample_num)
    rad = np.sin(rad) * 0.5 * np.pi
    rad = np.sin(rad) * 0.5 * np.pi
    x = (np.sin(rad) + 1) / 2

    return x

def get_accelerated_x_8x(sample_num=64):
    """
    単調増加ではなく、加速度が 0→1→0 となるような x を作る。
    加速度が `get_accelerated_x_1x` の4倍！！

    Parameters
    ----------
    sample_num : int
        the number of the sample.

    Returns
    -------
    array_like
        accelerated value list
    """
    rad = np.linspace(-0.5 * np.pi, 0.5 * np.pi, sample_num)
    rad = np.sin(rad) * 0.5 * np.pi
    rad = np.sin(rad) * 0.5 * np.pi
    rad = np.sin(rad) * 0.5 * np.pi
    x = (np.sin(rad) + 1) / 2

    return x

def generate_color_checker_rgb_value(
        color_space=BT709_COLOURSPACE, target_white=D65_WHITE):
    """
    Generate the 24 RGB values of the color checker.

    Parameters
    ----------
    color_space : color space
        color space object in `colour` module.

    target_white : array_like
        the xy values of the white point of target color space.

    Returns
    -------
    array_like
        24 RGB values. This is linear. OETF is not applied.

    Examples
    --------
    >>> generate_color_checker_rgb_value(
    ...     color_space=colour.models.BT709_COLOURSPACE,
    ...     target_white=[0.3127, 0.3290])
    >>> [[ 0.17289286  0.08205728  0.05714562]
    >>>  [ 0.5680292   0.29250401  0.21951748]
    >>>  [ 0.10435534  0.19656108  0.32958666]
    >>>  [ 0.1008804   0.14839018  0.05327639]
    >>>  [ 0.22303549  0.2169701   0.43166537]
    >>>  [ 0.10715338  0.513512    0.41415978]
    >>>  [ 0.74639182  0.20020473  0.03081343]
    >>>  [ 0.05947812  0.10659045  0.39897686]
    >>>  [ 0.5673215   0.08485376  0.11945382]
    >>>  [ 0.11177253  0.04285397  0.14166202]
    >>>  [ 0.34250836  0.5062777   0.0557734 ]
    >>>  [ 0.79262553  0.35803886  0.025485  ]
    >>>  [ 0.01864598  0.05139665  0.28886469]
    >>>  [ 0.054392    0.29876719  0.07187681]
    >>>  [ 0.45628547  0.03075684  0.04092033]
    >>>  [ 0.85379178  0.56503558  0.01475575]
    >>>  [ 0.53533883  0.09006355  0.3047824 ]
    >>>  [-0.03662977  0.24753781  0.39824679]
    >>>  [ 0.91177068  0.91497623  0.89427332]
    >>>  [ 0.57973934  0.59203191  0.59370647]
    >>>  [ 0.35495537  0.36538027  0.36772001]
    >>>  [ 0.19009594  0.19180133  0.19316719]
    >>>  [ 0.08524707  0.08890587  0.09255774]
    >>>  [ 0.03038879  0.03118623  0.03279615]]
    """
    colour_checker_param = COLOURCHECKERS.get('ColorChecker 2005')
    # 今回の処理では必要ないデータもあるので xyY と whitepoint だけ抽出
    # -------------------------------------------------------------
    _name, data, whitepoint = colour_checker_param
    temp_xyY = []
    for key in data.keys():
        temp_xyY.append(data[key])
    temp_xyY = np.array(temp_xyY)
    large_xyz = xyY_to_XYZ(temp_xyY)
    rgb_white_point = D65_WHITE
    illuminant_XYZ = whitepoint   # ColorCheckerのオリジナルデータの白色点
    illuminant_RGB = rgb_white_point  # XYZ to RGB 変換後の白色点を設定
    chromatic_adaptation_transform = 'CAT02'
    large_xyz_to_rgb_matrix = color_space.XYZ_to_RGB_matrix

    rgb = XYZ_to_RGB(
        large_xyz, illuminant_XYZ, illuminant_RGB,
        large_xyz_to_rgb_matrix, chromatic_adaptation_transform)

    return rgb

def make_color_checker_image(rgb, width=1920, padding_rate=0.01):
    """
    6x4 の カラーチェッカーの画像を作る。
    Height は Width から自動計算される。padding_rate で少し値が変わる。
    """
    h_patch_num = 6
    v_patch_num = 4

    # 各種パラメータ計算
    each_padding = int(width * padding_rate + 0.5)
    h_padding_total = each_padding * (h_patch_num + 1)
    h_patch_width_total = width - h_padding_total
    patch_height = h_patch_width_total // h_patch_num
    height = patch_height * v_patch_num + each_padding * (v_patch_num + 1)
    patch_width_list = equal_devision(h_patch_width_total, h_patch_num)

    # パッチを並べる
    img = np.zeros((height, width, 3))
    for v_idx in range(v_patch_num):
        h_pos_st = each_padding
        v_pos_st = each_padding + v_idx * (patch_height + each_padding)
        for h_idx in range(h_patch_num):
            rgb_idx = v_idx * h_patch_num + h_idx
            pos = (h_pos_st, v_pos_st)
            patch_img = np.ones((patch_height, patch_width_list[h_idx], 3))\
                * rgb[rgb_idx]
            merge(img, patch_img, pos)
            h_pos_st += (patch_width_list[h_idx] + each_padding)

    return img

def calc_st_pos_for_centering(bg_size, fg_size):
    """
    Calculate start postion for centering.

    Parameters
    ----------
    bg_size : touple(int)
        (width, height) of the background image.

    fg_size : touple(int)
        (width, height) of the foreground image.

    Returns
    -------
    touple (int)
        (st_pos_h, st_pos_v)

    Examples
    --------
    >>> calc_st_pos_for_centering(bg_size=(1920, 1080), fg_size=(640, 480))
    >>> (640, 300)
    """
    bg_width = bg_size[0]
    bg_height = bg_size[1]

    fg_width = fg_size[0]
    fg_height = fg_size[1]

    st_pos_h = bg_width // 2 - fg_width // 2
    st_pos_v = bg_height // 2 - fg_height // 2

    return (st_pos_h, st_pos_v)

def get_size_from_image(img):
    """
    `calc_st_pos_for_centering()` の引数計算が面倒だったので関数化。
    """
    return (img.shape[1], img.shape[0])

if __name__ == '__main__':
    os.chdir(os.path.dirname(os.path.abspath(__file__)))
    # print(calc_rad_patch_idx(outmost_num=9, current_num=1))
    # _plot_same_lstar_radial_color_patch_data(
    #     lstar=58, chroma=32.5, outmost_num=7,
    #     color_space=BT709_COLOURSPACE,
    #     transfer_function=tf.GAMMA24)
    # calc_rad_patch_idx2(outmost_num=9, current_num=7)
    # print(convert_luminance_to_color_value(100, tf.ST2084))
    # print(generate_color_checker_rgb_value(target_white=[0.3127, 0.3290]))
    print(calc_st_pos_for_centering(bg_size=(1920, 1080), fg_size=(640, 480)))

reason: import numpy,from scipy