text:
############################ 改变模型数据写入需要改变路径#########################
#训练并测试
from __future__ import print_function

import os
from skimage.transform import resize
from skimage.io import imsave
import numpy as np
from keras.models import Model
from keras.layers import Input, concatenate, Conv2D, MaxPooling2D, Conv2DTranspose
from keras.optimizers import Adam
from keras.callbacks import ModelCheckpoint
from keras import backend as K# keras后端   https://keras.io/zh/backend/

K.set_image_data_format('channels_last')  # TF dimension ordering in this code

img_rows = 512# 图像长
img_cols = 512# 图像宽

smooth = 1.

def load_train_data():# 读矩阵
    imgs_train = np.load('imgs_train.npy')
    imgs_mask_train = np.load('imgs_mask_train.npy')
    return imgs_train, imgs_mask_train

def load_test_data():# 读矩阵
    imgs_test = np.load('imgs_test.npy')
    imgs_id = np.load('imgs_id_test.npy')
    return imgs_test, imgs_id

# 计算dice指标，判断分割好坏
def dice_coef(y_true, y_pred):# y_true为真实准确值，y_pred为预测值
    y_true_f = K.flatten(y_true)# 捋直
    y_pred_f = K.flatten(y_pred)# 捋直
    # K.sum不加axi（指定方向求和，返回对应方向向量）,则为全元素求和，返回一个数字
    intersection = K.sum(y_true_f * y_pred_f)# 求预测准确的结果（真实准确值和预测值的交集）
    # 原始公式：（2*预测准确值）/（真实准确值+预测值），越大效果越好
    return (2. * intersection + smooth) / (K.sum(y_true_f) + K.sum(y_pred_f) + smooth)

# 损失函数
def dice_coef_loss(y_true, y_pred):
    return -dice_coef(y_true, y_pred)

# ################################################ 原始u-net 网络模型  #############################################
def get_unet(pretrained_weights=None):
    print('使用的是 原始get_unet')
    weight=32
    nb_filter = [weight, weight*2, weight*4, weight*8, weight*16]

    inputs = Input((img_rows, img_cols, 1))
    conv1 = Conv2D(nb_filter[0], (3, 3), activation='relu', padding='same')(inputs)
    conv1 = Conv2D(nb_filter[0], (3, 3), activation='relu', padding='same')(conv1)
    pool1 = MaxPooling2D(pool_size=(2, 2))(conv1)

    conv2 = Conv2D(nb_filter[1], (3, 3), activation='relu', padding='same')(pool1)
    conv2 = Conv2D(nb_filter[1], (3, 3), activation='relu', padding='same')(conv2)
    pool2 = MaxPooling2D(pool_size=(2, 2))(conv2)

    conv3 = Conv2D(nb_filter[2], (3, 3), activation='relu', padding='same')(pool2)
    conv3 = Conv2D(nb_filter[2], (3, 3), activation='relu', padding='same')(conv3)
    pool3 = MaxPooling2D(pool_size=(2, 2))(conv3)

    conv4 = Conv2D(nb_filter[3], (3, 3), activation='relu', padding='same')(pool3)
    conv4 = Conv2D(nb_filter[3], (3, 3), activation='relu', padding='same')(conv4)
    pool4 = MaxPooling2D(pool_size=(2, 2))(conv4)

    conv5 = Conv2D(nb_filter[4], (3, 3), activation='relu', padding='same')(pool4)
    conv5 = Conv2D(nb_filter[4], (3, 3), activation='relu', padding='same')(conv5)

    up6 = concatenate([Conv2DTranspose(256, (2, 2), strides=(2, 2), padding='same')(conv5), conv4], axis=3)
    conv6 = Conv2D(nb_filter[3], (3, 3), activation='relu', padding='same')(up6)
    conv6 = Conv2D(nb_filter[3], (3, 3), activation='relu', padding='same')(conv6)

    up7 = concatenate([Conv2DTranspose(128, (2, 2), strides=(2, 2), padding='same')(conv6), conv3], axis=3)
    conv7 = Conv2D(nb_filter[2], (3, 3), activation='relu', padding='same')(up7)
    conv7 = Conv2D(nb_filter[2], (3, 3), activation='relu', padding='same')(conv7)

    up8 = concatenate([Conv2DTranspose(64, (2, 2), strides=(2, 2), padding='same')(conv7), conv2], axis=3)
    conv8 = Conv2D(nb_filter[1], (3, 3), activation='relu', padding='same')(up8)
    conv8 = Conv2D(nb_filter[1], (3, 3), activation='relu', padding='same')(conv8)

    up9 = concatenate([Conv2DTranspose(32, (2, 2), strides=(2, 2), padding='same')(conv8), conv1], axis=3)
    conv9 = Conv2D(nb_filter[0], (3, 3), activation='relu', padding='same')(up9)
    conv9 = Conv2D(nb_filter[0], (3, 3), activation='relu', padding='same')(conv9)

    conv10 = Conv2D(1, (1, 1), activation='sigmoid')(conv9)

    model = Model(inputs=[inputs], outputs=[conv10])

    # optimizer：优化器，如Adam
    # loss:计算损失，dice_coef_loss损失函数
    # metrics: 列表，包含评估模型在训练和测试时的性能的指标，dice_coef
    model.compile(optimizer=Adam(lr=1e-5), loss=dice_coef_loss, metrics=[dice_coef])# dice_coef_loss损失函数

    if (pretrained_weights):
        model.load_weights(pretrained_weights)

    return model

# ################################################ 原始u-net wide 网络模型  #############################################
def get_unetw(pretrained_weights=None):
    print('使用的是 原始get_unet wide')

    weight=38
    nb_filter = [weight,weight*2,weight*4,weight*8,weight*16]

    inputs = Input((img_rows, img_cols, 1))
    conv1 = Conv2D(nb_filter[0], (3, 3), activation='relu', padding='same')(inputs)
    conv1 = Conv2D(nb_filter[0], (3, 3), activation='relu', padding='same')(conv1)
    pool1 = MaxPooling2D(pool_size=(2, 2))(conv1)

    conv2 = Conv2D(nb_filter[1], (3, 3), activation='relu', padding='same')(pool1)
    conv2 = Conv2D(nb_filter[1], (3, 3), activation='relu', padding='same')(conv2)
    pool2 = MaxPooling2D(pool_size=(2, 2))(conv2)

    conv3 = Conv2D(nb_filter[2], (3, 3), activation='relu', padding='same')(pool2)
    conv3 = Conv2D(nb_filter[2], (3, 3), activation='relu', padding='same')(conv3)
    pool3 = MaxPooling2D(pool_size=(2, 2))(conv3)

    conv4 = Conv2D(nb_filter[3], (3, 3), activation='relu', padding='same')(pool3)
    conv4 = Conv2D(nb_filter[3], (3, 3), activation='relu', padding='same')(conv4)
    pool4 = MaxPooling2D(pool_size=(2, 2))(conv4)

    conv5 = Conv2D(nb_filter[4], (3, 3), activation='relu', padding='same')(pool4)
    conv5 = Conv2D(nb_filter[4], (3, 3), activation='relu', padding='same')(conv5)

    up6 = concatenate([Conv2DTranspose(256, (2, 2), strides=(2, 2), padding='same')(conv5), conv4], axis=3)
    conv6 = Conv2D(nb_filter[3], (3, 3), activation='relu', padding='same')(up6)
    conv6 = Conv2D(nb_filter[3], (3, 3), activation='relu', padding='same')(conv6)

    up7 = concatenate([Conv2DTranspose(128, (2, 2), strides=(2, 2), padding='same')(conv6), conv3], axis=3)
    conv7 = Conv2D(nb_filter[2], (3, 3), activation='relu', padding='same')(up7)
    conv7 = Conv2D(nb_filter[2], (3, 3), activation='relu', padding='same')(conv7)

    up8 = concatenate([Conv2DTranspose(64, (2, 2), strides=(2, 2), padding='same')(conv7), conv2], axis=3)
    conv8 = Conv2D(nb_filter[1], (3, 3), activation='relu', padding='same')(up8)
    conv8 = Conv2D(nb_filter[1], (3, 3), activation='relu', padding='same')(conv8)

    up9 = concatenate([Conv2DTranspose(32, (2, 2), strides=(2, 2), padding='same')(conv8), conv1], axis=3)
    conv9 = Conv2D(nb_filter[0], (3, 3), activation='relu', padding='same')(up9)
    conv9 = Conv2D(nb_filter[0], (3, 3), activation='relu', padding='same')(conv9)

    conv10 = Conv2D(1, (1, 1), activation='sigmoid')(conv9)

    model = Model(inputs=[inputs], outputs=[conv10])

    # optimizer：优化器，如Adam
    # loss:计算损失，dice_coef_loss损失函数
    # metrics: 列表，包含评估模型在训练和测试时的性能的指标，dice_coef
    model.compile(optimizer=Adam(lr=1e-5), loss=dice_coef_loss, metrics=[dice_coef])# dice_coef_loss损失函数

    if (pretrained_weights):
        model.load_weights(pretrained_weights)

    return model

# ################################################ 原始u-net++网络模型  #############################################

def get_unetpp(num_class=1, deep_supervision=False):
    print('使用的是原始unet++')
    nb_filter = [32,64,128,256,512]
    img_rows=256
    img_cols=256
    color_type=1
    bn_axis = 3

    img_input = Input(shape=(img_rows, img_cols, color_type), name='main_input')

    conv1_1 = Conv2D(nb_filter[0], (3, 3), activation='relu', padding='same')(img_input)
    conv1_1 = Conv2D(nb_filter[0], (3, 3), activation='relu', padding='same')(conv1_1)
    pool1 = MaxPooling2D((2, 2), strides=(2, 2), name='pool1')(conv1_1)

    conv2_1 = Conv2D(nb_filter[1], (3, 3), activation='relu', padding='same')(pool1)
    conv2_1 = Conv2D(nb_filter[1], (3, 3), activation='relu', padding='same')(conv2_1)
    pool2 = MaxPooling2D((2, 2), strides=(2, 2), name='pool2')(conv2_1)

    up1_2 = Conv2DTranspose(nb_filter[0], (2, 2), strides=(2, 2), name='up12', padding='same')(conv2_1)
    conv1_2 = concatenate([up1_2, conv1_1], name='merge12', axis=bn_axis)
    conv1_2 = Conv2D(nb_filter[0], (3, 3), activation='relu', padding='same')(conv1_2)
    conv1_2 = Conv2D(nb_filter[0], (3, 3), activation='relu', padding='same')(conv1_2)

    conv3_1 = Conv2D(nb_filter[2], (3, 3), activation='relu', padding='same')(pool2)
    conv3_1 = Conv2D(nb_filter[2], (3, 3), activation='relu', padding='same')(conv3_1)
    pool3 = MaxPooling2D((2, 2), strides=(2, 2), name='pool3')(conv3_1)

    up2_2 = Conv2DTranspose(nb_filter[1], (2, 2), strides=(2, 2), name='up22', padding='same')(conv3_1)
    conv2_2 = concatenate([up2_2, conv2_1], name='merge22', axis=bn_axis)
    conv2_2 = Conv2D(nb_filter[1], (3, 3), activation='relu', padding='same')(conv2_2)
    conv2_2 = Conv2D(nb_filter[1], (3, 3), activation='relu', padding='same')(conv2_2)

    up1_3 = Conv2DTranspose(nb_filter[0], (2, 2), strides=(2, 2), name='up13', padding='same')(conv2_2)
    conv1_3 = concatenate([up1_3, conv1_1, conv1_2], name='merge13', axis=bn_axis)
    conv1_3 = Conv2D(nb_filter[0], (3, 3), activation='relu', padding='same')(conv1_3)
    conv1_3 = Conv2D(nb_filter[0], (3, 3), activation='relu', padding='same')(conv1_3)

    conv4_1 = Conv2D(nb_filter[3], (3, 3), activation='relu', padding='same')(pool3)
    conv4_1 = Conv2D(nb_filter[3], (3, 3), activation='relu', padding='same')(conv4_1)
    pool4 = MaxPooling2D((2, 2), strides=(2, 2), name='pool4')(conv4_1)

    up3_2 = Conv2DTranspose(nb_filter[2], (2, 2), strides=(2, 2), name='up32', padding='same')(conv4_1)
    conv3_2 = concatenate([up3_2, conv3_1], name='merge32', axis=bn_axis)
    conv3_2 = Conv2D(nb_filter[2], (3, 3), activation='relu', padding='same')(conv3_2)
    conv3_2 = Conv2D(nb_filter[2], (3, 3), activation='relu', padding='same')(conv3_2)

    up2_3 = Conv2DTranspose(nb_filter[1], (2, 2), strides=(2, 2), name='up23', padding='same')(conv3_2)
    conv2_3 = concatenate([up2_3, conv2_1, conv2_2], name='merge23', axis=bn_axis)
    conv2_3 = Conv2D(nb_filter[1], (3, 3), activation='relu', padding='same')(conv2_3)
    conv2_3 = Conv2D(nb_filter[1], (3, 3), activation='relu', padding='same')(conv2_3)

    up1_4 = Conv2DTranspose(nb_filter[0], (2, 2), strides=(2, 2), name='up14', padding='same')(conv2_3)
    conv1_4 = concatenate([up1_4, conv1_1, conv1_2, conv1_3], name='merge14', axis=bn_axis)
    conv1_4 = Conv2D(nb_filter[0], (3, 3), activation='relu', padding='same')(conv1_4)
    conv1_4 = Conv2D(nb_filter[0], (3, 3), activation='relu', padding='same')(conv1_4)

    conv5_1 = standard_unit(pool4, stage='51', nb_filter=nb_filter[4])
    conv5_1 = Conv2D(nb_filter[4], (3, 3), activation='relu', padding='same')(pool4)
    conv5_1 = Conv2D(nb_filter[4], (3, 3), activation='relu', padding='same')(conv5_1)

    up4_2 = Conv2DTranspose(nb_filter[3], (2, 2), strides=(2, 2), name='up42', padding='same')(conv5_1)
    conv4_2 = concatenate([up4_2, conv4_1], name='merge42', axis=bn_axis)
    conv4_2 = Conv2D(nb_filter[3], (3, 3), activation='relu', padding='same')(conv4_2)
    conv4_2 = Conv2D(nb_filter[3], (3, 3), activation='relu', padding='same')(conv4_2)

    up3_3 = Conv2DTranspose(nb_filter[2], (2, 2), strides=(2, 2), name='up33', padding='same')(conv4_2)
    conv3_3 = concatenate([up3_3, conv3_1, conv3_2], name='merge33', axis=bn_axis)
    conv3_3 = Conv2D(nb_filter[2], (3, 3), activation='relu', padding='same')(conv3_3)
    conv3_3 = Conv2D(nb_filter[2], (3, 3), activation='relu', padding='same')(conv3_3)

    up2_4 = Conv2DTranspose(nb_filter[1], (2, 2), strides=(2, 2), name='up24', padding='same')(conv3_3)
    conv2_4 = concatenate([up2_4, conv2_1, conv2_2, conv2_3], name='merge24', axis=bn_axis)
    conv2_4 = Conv2D(nb_filter[1], (3, 3), activation='relu', padding='same')(conv2_4)
    conv2_4 = Conv2D(nb_filter[1], (3, 3), activation='relu', padding='same')(conv2_4)

    up1_5 = Conv2DTranspose(nb_filter[0], (2, 2), strides=(2, 2), name='up15', padding='same')(conv2_4)
    conv1_5 = concatenate([up1_5, conv1_1, conv1_2, conv1_3, conv1_4], name='merge15', axis=bn_axis)
    conv1_5 = Conv2D(nb_filter[0], (3, 3), activation='relu', padding='same')(conv1_5)
    conv1_5 = Conv2D(nb_filter[0], (3, 3), activation='relu', padding='same')(conv1_5)

    nestnet_output_1 = Conv2D(num_class, (1, 1), activation='sigmoid', name='output_1', kernel_initializer = 'he_normal', padding='same', kernel_regularizer=l2(1e-4))(conv1_2)
    nestnet_output_2 = Conv2D(num_class, (1, 1), activation='sigmoid', name='output_2', kernel_initializer = 'he_normal', padding='same', kernel_regularizer=l2(1e-4))(conv1_3)
    nestnet_output_3 = Conv2D(num_class, (1, 1), activation='sigmoid', name='output_3', kernel_initializer = 'he_normal', padding='same', kernel_regularizer=l2(1e-4))(conv1_4)
    nestnet_output_4 = Conv2D(num_class, (1, 1), activation='sigmoid', name='output_4', kernel_initializer = 'he_normal', padding='same', kernel_regularizer=l2(1e-4))(conv1_5)

    model = Model(input=img_input, output=[nestnet_output_4])

    # optimizer：优化器，如Adam
    # loss:计算损失，dice_coef_loss损失函数
    # metrics: 列表，包含评估模型在训练和测试时的性能的指标，dice_coef
    model.compile(optimizer=Adam(lr=1e-5), loss=dice_coef_loss, metrics=[dice_coef])# dice_coef_loss损失函数

    return model

def preprocess(imgs):# 矩阵变换
    # np shape函数为对应矩阵维度的**，类似于长宽高，返回元组，这里[0]为对应图片数量
    imgs_p = np.ndarray((imgs.shape[0], img_rows, img_cols), dtype=np.uint8)# 创建矩阵（对应图片数量*256*256）
    for i in range(imgs.shape[0]):
        imgs_p[i] = resize(imgs[i], (img_cols, img_rows), preserve_range=True)# 图片缩放：512 -> 256

    imgs_p = imgs_p[..., np.newaxis]# 增加一个维度,3维 -> 4维
    return imgs_p

def train_and_predict():
    print('Loading and preprocessing train data...')
    imgs_train, imgs_mask_train = load_train_data()# 读矩阵（对应图片数量*图片长（512）*图片宽（512））

    imgs_train = preprocess(imgs_train)# 矩阵变换
    imgs_mask_train = preprocess(imgs_mask_train)
    imgs_train = imgs_train.astype('float32')
    mean = np.mean(imgs_train)  # 平均值mean for data centering
    std = np.std(imgs_train)  # 矩阵全局标准差std for data normalization

    imgs_train -= mean
    imgs_train /= std

    imgs_mask_train = imgs_mask_train.astype('float32')
    imgs_mask_train /= 255.  # 矩阵归一化scale masks to [0, 1]

    # 图片信息矩阵输入并开始训练
    #model = get_unet('/opt/bin/unet_weights.h5') #重复训练则读取本文件继续加载训练,加载好train1和train2

    #model = get_unet()  # 原始unet
    model = get_unetw()  # 原始unetw
    #model = get_unetpp() # 原始unet++

    # 图片保存路径
    #pred_dir = 'preds_unet'
    pred_dir = 'preds_unetw_38'
    #pred_dir = 'preds_unetpp'

    #save_res = '/opt/bin/unet_weights.h5'############################ 改变模型数据写入需要改变路径#########################
    save_res = '/opt/bin/unet_weightsw_38.h5'
    #save_res = '/opt/bin/unet_weightspp.h5'

    #'''
	#在训练生成 .h5 文件后可以直接注释该段选用其他测试集进行测试
    model_checkpoint = ModelCheckpoint(save_res, monitor='val_loss', save_best_only=True)
    print('Fitting model...')
    model.fit(imgs_train,
              imgs_mask_train,
              batch_size=16,
              nb_epoch=60, # 训练次数
              verbose=1,
              shuffle=True,
              validation_split=0.2,
              callbacks=[model_checkpoint])
    #'''

    # 获取测试集信息矩阵
    imgs_test, imgs_id_test = load_test_data()
    imgs_test = preprocess(imgs_test)

    imgs_test = imgs_test.astype('float32')
    imgs_test -= mean
    imgs_test /= std

    model.load_weights(save_res)
    # 开始测试
    print('Predicting masks on test data...')

    imgs_mask_test = model.predict(imgs_test, verbose=1)
    np.save('imgs_mask_test.npy', imgs_mask_test)

    if not os.path.exists(pred_dir):
        os.mkdir(pred_dir)
    for image, image_id in zip(imgs_mask_test, imgs_id_test):
        image = (image[:, :, 0] * 255.).astype(np.uint8)
        imsave(os.path.join(pred_dir, str(image_id) + '_mask.png'), image)

if __name__ == '__main__':
    train_and_predict()


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

import networkx as nx

class MapAnalyzer(object):

    def __init__(self, webpage_currency, output_encoding):

        self.webpage_currency = webpage_currency
        self.output_encoding  = output_encoding

    def list_flights(self, flights_map, price_limit, flights_limit, departure_point):

        print("Flights from {0}:".format(flights_map.points_names[departure_point].encode(self.output_encoding)))

        costs, paths = nx.single_source_dijkstra(flights_map.map, departure_point)
        self._find_cycle_from_dep_point(flights_map.map, departure_point, costs, paths, price_limit)

        price_w = len(str(max(costs)))

        for point in sorted(costs, key=lambda y: costs[y]):
            print("{0:{1}} {2} - {3}".format(costs[point], price_w, self.webpage_currency, flights_map.points_names[point].encode(self.output_encoding))),

            p = paths[point][1:-1]

            if len(p)>0 :
                print('\t( via'),
                for c in p:
                    print(flights_map.points_names[c].encode(self.output_encoding)),
                print(")"),

            print("")

    @staticmethod
    def  _find_cycle_from_dep_point(map, point, costs, paths, price_limit):

        costs[point] = float("inf")

        if map.edge[point].has_key(point) and map.edge[point][point]["weight"] <= price_limit:
            costs[point] = map.edge[point][point]["weight"]

        map.add_node("point_copy")

        for p in map.edge[point].copy():
            map.add_edge(p, "point_copy", weight=map.edge[point][p]["weight"])

        copy_cost, copy_path = nx.single_source_dijkstra(map, point, "point_copy")

        if copy_cost["point_copy"] < costs[point] and copy_cost["point_copy"] <= price_limit:
            costs[point] = copy_cost["point_copy"]
            paths[point] = copy_path["point_copy"]

        map.remove_node("point_copy")

        if costs[point] == float("inf"):
            costs.pop(point, None)

ext: py

reason: import networkx


text:
"""Wrap the details of an inversion.

Hide most of the implementation details.
Take prior and parameters, return posterior.

The interface is still experimental at this point.  I am likely to
remove functions and parameters as I find better ways to organize
things.  The code still works fine as an example.
"""
from distutils.version import LooseVersion
import subprocess
import datetime
import sys
import os
try:
    from shlex import quote
except ImportError:
    from pipes import quote
from pwd import getpwuid

import numpy as np
import dateutil.tz

import pandas as pd
import xarray

import atmos_flux_inversion.optimal_interpolation
import atmos_flux_inversion.correlations as inv_corr
import atmos_flux_inversion.covariances
import atmos_flux_inversion.util

HOURS_PER_DAY = 24
OBSERVATION_TEMPORAL_CORRELATION_FUNCTION = (  # noqa: C0103
    inv_corr.ExponentialCorrelation)
UTC = dateutil.tz.tzutc()
UDUNITS_DATE = "%Y-%m-%d %H:%M:%S%z"
ACDD_DATE = "%Y-%m-%dT%H:%M:%S%z"
CALENDAR = "standard"
RUN_DATE = datetime.datetime.now(tz=UTC)
_PACKAGE_INFO = None

def invert_uniform(prior_fluxes, observations,
                   observation_operator,
                   prior_correlation_length, prior_correlation_structure,
                   prior_correlation_time_days, prior_correlation_time_hours,
                   prior_flux_stds,
                   observation_correlation_time,
                   method=(
                       atmos_flux_inversion.optimal_interpolation.fold_common
                   ),
                   output_uncertainty_frequency="MS"):
    """Perform an inversion.

    Assumes error correlations are invariant in time and space and
    that the error variances are constant in time.  Also assumes the
    observation correlations between towers is zero.

    Parameters
    ----------
    prior_fluxes: xarray.Dataarray[flux_time, y, x]
        The prior or background estimate of the fluxes
    observations: xarray.Dataarray[obs_time, site]
        The observed mole fractions
    observation_operator: xarray.Dataarray[obs_time, site, flux_time, y, x]
        The linearized operator mapping a flux distribution to
        observations.
    prior_correlation_length: float
        The lengthscale for the spatial correlations in the prior.
    prior_correlation_structure: inv_corr.DistanceCorrelationFunction
        The structure of the spatial correlations in the prior.
    prior_correlation_time_days: float
        The correlation timescale for fluxes across different days.
    prior_correlation_time_hours: float
        The correlation timescale for fluxes within the same day.
    prior_flux_stds: xarray.Dataarray[y, x]
        The standard deviations of the prior flux distribution at each
        point in space.
    observation_correlation_time: float
        The correlation timescale for consecutive observations at a
        single site.
    method: function
        The method to use.  Must match the signature of
        :func:`atmos_flux_inversion.optimal_interpolation.simple`
    output_uncertainty_frequency: str
        One of the frequencies from
        `<http://pandas.pydata.org/pandas-docs/stable/user_guide/timeseries.html#offset-aliases>`_
        or "season".  Often "M", output_uncertainty_frequency, "Y",
        "YS", or "season" make the most sense.  Used to provide
        lower-temporal-resolution posterior uncertainties.

    Returns
    -------
    xarray.Dataset
        The results of the inversion.
        Contents:

        prior[flux_time, y, x]
            The prior mean for the inversion.
        increment[flux_time, y, x]
            The change from the prior to the posterior mean estimates.
        posterior[flux_time, y, x]
            The posterior mean for the inversion.
        post_cov[red_flux_time_adj, y_adj, x_adj, red_flux_time, y, x]
            The analytic uncertainty for the posterior, expressed as a
            covariance matrix on a reduced-resolution domain.  If the
            prior and observation likelihood are Gaussian, this will
            be exact.
    """
    y_index_name = [name for name in prior_fluxes.coords
                    if "y" in name.lower()][0]
    x_index_name = y_index_name.replace("y", "x")
    flux_time_index_name = [name for name in prior_fluxes.coords
                            if "time" in name.lower()][0]
    site_index_names = [name for name in observations.coords
                        if "site" in name.lower()]
    obs_time_index_names = [name for name in observations.coords
                            if "time" in name.lower()]

    obs_is_multi = False
    for index_name, index in observations.indexes.items():
        if isinstance(index, pd.MultiIndex):
            obs_is_multi = True
            obs_time_index_names.extend([name for name in index.names
                                         if "time" in name])
            site_index_names.extend([name for name in index.names
                                     if "site" in name])
            observations = observations.rename(
                {index_name: "observations"})
            observation_operator = observation_operator.rename(
                {index_name: "observations"})

    if len(site_index_names) == 1:
        site_index_name = site_index_names[0]
    elif "site" in site_index_names:
        site_index_name = "site"

    obs_time_index_name = obs_time_index_names[0]

    x_index = prior_fluxes[x_index_name]
    y_index = prior_fluxes[y_index_name]
    obs_time_index = observations.coords[obs_time_index_name]
    flux_time_index = prior_fluxes.indexes[flux_time_index_name]

    site_names = list(set(observations.coords[site_index_name].values))

    if not obs_is_multi:
        observations_for_inversion = observations.stack(dict(
            observations=[obs_time_index_name, site_index_name]
        ))
    else:
        observations_for_inversion = observations
    long_obs_times = observations_for_inversion.coords[obs_time_index_name]
    long_site = observations_for_inversion.coords[site_index_name]

    flux_time_adjoint_index = flux_time_index.copy()
    flux_time_adjoint_index.name += "_adjoint"

    dx = abs(x_index[1] - x_index[0])
    obs_interval = abs(obs_time_index[1] - obs_time_index[0])
    flux_interval = abs(flux_time_index[1] - flux_time_index[0])

    spatial_correlations = (
        atmos_flux_inversion.correlations.
        HomogeneousIsotropicCorrelation.from_function(
            prior_correlation_structure(prior_correlation_length / dx),
            prior_fluxes.shape[1:]))
    sqrt_spatial_variances = prior_flux_stds.transpose(
        y_index_name, x_index_name)
    spatial_covariances = (
        atmos_flux_inversion.covariances.CorrelationStandardDeviation(
            spatial_correlations, sqrt_spatial_variances.data
        )
    )

    hour_correlations = (
        atmos_flux_inversion.correlations.
        HomogeneousIsotropicCorrelation.from_function(
            atmos_flux_inversion.correlations.ExponentialCorrelation(
                pd.Timedelta(hours=prior_correlation_time_hours) /
                flux_interval),
            pd.Timedelta(hours=HOURS_PER_DAY) // flux_interval))
    day_correlations = (
        atmos_flux_inversion.correlations.make_matrix(
            atmos_flux_inversion.correlations.ExponentialCorrelation(
                prior_correlation_time_days),
            prior_fluxes.shape[0]))
    temporal_covariances = atmos_flux_inversion.util.kron(
        day_correlations.dot(np.eye(*day_correlations.shape)),
        hour_correlations.dot(np.eye(*hour_correlations.shape))
    )

    temporal_covariance_dataset = xarray.DataArray(
        dims=(flux_time_index_name, flux_time_adjoint_index.name),
        data=temporal_covariances,
        attrs=dict(long_name="temporal_covariances"),
        coords={flux_time_index_name: flux_time_index,
                flux_time_adjoint_index.name: flux_time_adjoint_index}
    )

    reduced_temp_cov_ds = temporal_covariance_dataset.resample(**{
        flux_time_adjoint_index.name: output_uncertainty_frequency
    }).sum(flux_time_adjoint_index.name).resample(**{
        flux_time_index_name: output_uncertainty_frequency
    }).sum(flux_time_index_name)

    reduced_obs_op = observation_operator.resample(**{
        flux_time_index_name: output_uncertainty_frequency
    }).sum(flux_time_index_name)

    prior_covariances = (
        atmos_flux_inversion.util.kronecker_product(temporal_covariances,
                                                    spatial_covariances)
    )
    reduced_prior_covariances = (
        atmos_flux_inversion.util.kron(
            reduced_temp_cov_ds.data,
            spatial_covariances.dot(np.eye(*spatial_covariances.shape))
        )
    )

    observation_temporal_correlations = (
        OBSERVATION_TEMPORAL_CORRELATION_FUNCTION(
            observation_correlation_time)
    )(abs(
        long_obs_times.values[:, np.newaxis] -
        long_obs_times.values[np.newaxis, :]
    ).astype("m8[h]").astype(int))
    observation_site_correlations = (
        long_site.values[:, np.newaxis] == long_site.values[np.newaxis, :])
    observation_correlations = (
        observation_temporal_correlations * observation_site_correlations)
    observation_covariance = observation_correlations

    if not obs_is_multi:
        obs_stack_dict = dict(
            observations=[obs_time_index_name, site_index_name],
            fluxes=[flux_time_index_name, y_index_name, x_index_name]
        )
    else:
        obs_stack_dict = dict(
            fluxes=[flux_time_index_name, y_index_name, x_index_name]
        )

    posterior_fluxes, posterior_covariance = method(
        prior_fluxes.data.reshape(prior_covariances.shape[0]),
        prior_covariances,
        observations.data.reshape(observation_covariance.shape[0]),
        observation_covariance,
        # Question: should this transpose obs_op to ensure the proper
        # dimension order?
        observation_operator.stack(obs_stack_dict).transpose(
            "observations", "fluxes"),
        reduced_prior_covariances,
        reduced_obs_op.stack(obs_stack_dict).transpose(
            "observations", "fluxes"),
    )
    posterior_fluxes = posterior_fluxes.reshape(prior_fluxes.shape)

    posterior_var_atts = prior_fluxes.attrs.copy()
    posterior_var_atts.update(dict(
        long_name="posterior_fluxes",
        description=(
            "posterior mean fluxes obtained using atmos_flux_inversion package"
        ),
        origin="atmos_flux_inversion package",
        units=prior_fluxes.attrs["units"],
    ))
    increment_var_atts = prior_fluxes.attrs.copy()
    increment_var_atts.update(dict(
        long_name="flux_increment",
        units=prior_fluxes.attrs["units"],
        description="Change from prior to posterior using dask",
        origin="OI and dask",
    ))
    posterior_covariance_attributes = prior_fluxes.attrs.copy()
    posterior_covariance_attributes.update(dict(
        long_name="inversion_posterior_covariance",
        units="({units:s})^2".format(units=prior_fluxes.attrs["units"]),
        description=("Posterior flux covariance"
                     "obtained using atmos_flux_inversion package"),
        origin="atmos_flux_inversion package",
    ))
    posterior_global_atts = global_attributes_dict()
    posterior_global_atts.update(dict(
        title="Bayesian flux inversion results",
        summary="Results of a Bayesian flux inversion",
        product_version="v0.0.0.dev0",
        cdm_data_type="grid",
        source="Inversion package",
        # Descriptions of the inversion that may be hard to get later
        prior_flux_long_name=prior_fluxes.attrs.get("long_name", ""),
        correlation_length=prior_correlation_length,
        correlation_time_day=prior_correlation_time_days,
        correlation_time_hour=prior_correlation_time_hours,
        spatial_correlation_structure=prior_correlation_structure.__name__,
        observation_interval=obs_interval,
        observation_sites_used=site_names,
    ))

    posterior_ds = xarray.Dataset(
        dict(
            posterior=(
                prior_fluxes.dims,
                posterior_fluxes.reshape(prior_fluxes.shape),
                posterior_var_atts
            ),
            prior=prior_fluxes,
            increment=(
                prior_fluxes.dims,
                posterior_fluxes - prior_fluxes,
                increment_var_atts
            ),
            posterior_covariance=(
                ["{dim:s}_adjoint".format(dim=dim)
                 if "time" not in dim else
                 "reduced_{dim:s}_adjoint".format(dim=dim)
                 for dim in prior_fluxes.dims] +
                [dim if "time" not in dim else
                 "reduced_{dim:s}".format(dim=dim)
                 for dim in prior_fluxes.dims],
                posterior_covariance.reshape(
                    (reduced_temp_cov_ds.shape[0],
                     len(y_index),
                     len(x_index)) * 2
                ),
                posterior_covariance_attributes
            ),
        ),
        prior_fluxes.coords,
        posterior_global_atts)

    coord_list = set(posterior_ds.coords)
    for coord in coord_list:
        if "time" in coord:
            continue
        posterior_ds.coords["{coord:s}_adjoint".format(coord=coord)] = (
            posterior_ds.coords[coord]
        )

    for coord in reduced_temp_cov_ds.coords:
        posterior_ds.coords["reduced_{coord:s}".format(coord=coord)] = (
            reduced_temp_cov_ds.coords[coord]
        )

    return posterior_ds

def get_installed_modules():
    """Get the list of installed modules.

    Returns
    -------
    list of list of str
        List of currently installed packages, using conda if that is
        available or pip if it is not.  List of two-element strings
        giving name-version pairs.
    """
    global _PACKAGE_INFO
    if _PACKAGE_INFO is not None:
        return _PACKAGE_INFO

    try:  # pragma: no branch
        with open(os.devnull, "wb") as stderr:
            output = subprocess.check_output(
                ["conda", "list", "--export", "--no-show-channel-urls"],
                universal_newlines=True, stderr=stderr)
        package_info = [line.split("=")[:2] for line in output.split("\n")
                        if line and not line.startswith("#")]
        python_versions = [pair[1] for pair in package_info
                           if pair[0] == "python"]

        if ((LooseVersion(python_versions[0]) ==
             LooseVersion(sys.version.split()[0]))):
            _PACKAGE_INFO = package_info
            return package_info
    except (subprocess.CalledProcessError, OSError):  # pragma: no cover
        # Testing with python3 from conda and python2 from the system
        # will not enter this branch.
        pass  # pragma: no cover

    # Conda not present or not providing current executable
    pip_args = []
    if sys.executable is not None:  # pragma: no branch
        pip_args = [sys.executable, "-m"]
    pip_args.extend(["pip", "freeze"])
    try:  # pragma: no branch
        output = subprocess.check_output(
            pip_args,
            universal_newlines=True
        )
        package_info = [line.split("==") for line in output.split("\n")
                        if line and not line.startswith("-")]

        _PACKAGE_INFO = package_info
        return package_info
    except OSError:  # pragma: no cover
        # Avoid crashes on rare platforms lacking pip
        _PACKAGE_INFO = []    # pragma: no cover
        return _PACKAGE_INFO  # pragma: no cover

def global_attributes_dict():
    """Set global attributes required by conventions.

    Currently CF-1.6 and ACDD-1.3.

    Still needs title, summary, source, creator_institution,
    product_version, references, cdm_data_type, institution,
    geospatial_vertical_{min,max,positive,units}, ...

    Users may want to overwrite the history and creator_name
    attributes, as these are fancy guesses.

    Returns
    -------
    global_atts: dict
        attribute_name: attribute_value mapping.
    """
    return dict(
        Conventions="CF-1.6 ACDD-1.3",
        standard_name_vocabulary="CF Standard Name Table v32",
        history=("{now:{date_fmt:s}}: Created by {progname:s}"
                 "with command line: {cmd_line:s}").format(
            now=RUN_DATE, date_fmt=UDUNITS_DATE, progname=sys.argv[0],
            cmd_line=" ".join(quote(arg) for arg in sys.argv)),
        date_created="{now:{date_fmt:s}}".format(
            now=RUN_DATE, date_fmt=ACDD_DATE),
        date_modified="{now:{date_fmt:s}}".format(
            now=RUN_DATE, date_fmt=ACDD_DATE),
        date_metadata_modified="{now:{date_fmt:s}}".format(
            now=RUN_DATE, date_fmt=ACDD_DATE),
        creator_name=getpwuid(os.getuid())[0],
        installed_modules=[
            "{name:s}={version:s}".format(name=name, version=version)
            for name, version in get_installed_modules()],
    )

reason: import numpy


text:
# encoding: utf-8
# module pandas._libs.window
# from C:\Python27\lib\site-packages\pandas\_libs\window.pyd
# by generator 1.147
# no doc

# imports
import numpy as np # C:\Python27\lib\site-packages\numpy\__init__.pyc
import __builtin__ as __builtins__ # <module '__builtin__' (built-in)>

# functions

def ewma(*args, **kwargs): # real signature unknown
    """
    Compute exponentially-weighted moving average using center-of-mass.

        Parameters
        ----------
        vals : ndarray (float64 type)
        com : float64
        adjust: int
        ignore_na: int
        minp: int

        Returns
        -------
        y : ndarray
    """
    pass

def ewmcov(*args, **kwargs): # real signature unknown
    """
    Compute exponentially-weighted moving variance using center-of-mass.

        Parameters
        ----------
        input_x : ndarray (float64 type)
        input_y : ndarray (float64 type)
        com : float64
        adjust: int
        ignore_na: int
        minp: int
        bias: int

        Returns
        -------
        y : ndarray
    """
    pass

def get_window_indexer(*args, **kwargs): # real signature unknown
    """
    return the correct window indexer for the computation

        Parameters
        ----------
        values: 1d ndarray
        win: integer, window size
        minp: integer, minimum periods
        index: 1d ndarray, optional
            index to the values array
        closed: string, default None
            {'right', 'left', 'both', 'neither'}
            window endpoint closedness. Defaults to 'right' in
            VariableWindowIndexer and to 'both' in FixedWindowIndexer
        floor: optional
            unit for flooring the unit
        use_mock: boolean, default True
            if we are a fixed indexer, return a mock indexer
            instead of the FixedWindow Indexer. This is a type
            compat Indexer that allows us to use a standard
            code path with all of the indexers.

        Returns
        -------
        tuple of 1d int64 ndarrays of the offsets & data about the window
    """
    pass

def roll_count(*args, **kwargs): # real signature unknown
    pass

def roll_generic(*args, **kwargs): # real signature unknown
    pass

def roll_kurt(*args, **kwargs): # real signature unknown
    pass

def roll_max(*args, **kwargs): # real signature unknown
    """
    Moving max of 1d array of any numeric type along axis=0 ignoring NaNs.

        Parameters
        ----------
        values: numpy array
        window: int, size of rolling window
        minp: if number of observations in window
              is below this, output a NaN
        index: ndarray, optional
           index for window computation
        closed: 'right', 'left', 'both', 'neither'
                make the interval closed on the right, left,
                both or neither endpoints
    """
    pass

def roll_mean(*args, **kwargs): # real signature unknown
    pass

def roll_median_c(*args, **kwargs): # real signature unknown
    pass

def roll_min(*args, **kwargs): # real signature unknown
    """
    Moving max of 1d array of any numeric type along axis=0 ignoring NaNs.

        Parameters
        ----------
        values: numpy array
        window: int, size of rolling window
        minp: if number of observations in window
              is below this, output a NaN
        index: ndarray, optional
           index for window computation
    """
    pass

def roll_quantile(*args, **kwargs): # real signature unknown
    """ O(N log(window)) implementation using skip list """
    pass

def roll_skew(*args, **kwargs): # real signature unknown
    pass

def roll_sum(*args, **kwargs): # real signature unknown
    pass

def roll_var(*args, **kwargs): # real signature unknown
    """ Numerically stable implementation using Welford's method. """
    pass

def roll_window(*args, **kwargs): # real signature unknown
    """ Assume len(weights) << len(values) """
    pass

def _check_minp(*args, **kwargs): # real signature unknown
    """
    Parameters
        ----------
        win: int
        minp: int or None
        N: len of window
        floor: int, optional
            default 1

        Returns
        -------
        minimum period
    """
    pass

def __pyx_unpickle_Enum(*args, **kwargs): # real signature unknown
    pass

def __pyx_unpickle_FixedWindowIndexer(*args, **kwargs): # real signature unknown
    pass

def __pyx_unpickle_MockFixedWindowIndexer(*args, **kwargs): # real signature unknown
    pass

def __pyx_unpickle_VariableWindowIndexer(*args, **kwargs): # real signature unknown
    pass

def __pyx_unpickle_WindowIndexer(*args, **kwargs): # real signature unknown
    pass

# classes

class WindowIndexer(object):
    # no doc
    def get_data(self, *args, **kwargs): # real signature unknown
        pass

    def __init__(self, *args, **kwargs): # real signature unknown
        pass

    @staticmethod # known case of __new__
    def __new__(S, *more): # real signature unknown; restored from __doc__
        """ T.__new__(S, ...) -> a new object with type S, a subtype of T """
        pass

    def __reduce__(self, *args, **kwargs): # real signature unknown
        pass

    def __setstate__(self, *args, **kwargs): # real signature unknown
        pass

class FixedWindowIndexer(WindowIndexer):
    """
    create a fixed length window indexer object
        that has start & end, that point to offsets in
        the index object; these are defined based on the win
        arguments

        Parameters
        ----------
        values: ndarray
            values data array
        win: int64_t
            window size
        minp: int64_t
            min number of obs in a window to consider non-NaN
        index: object
            index of the values
        floor: optional
            unit for flooring the unit
        left_closed: bint
            left endpoint closedness
        right_closed: bint
            right endpoint closedness
    """
    def __init__(self, *args, **kwargs): # real signature unknown
        pass

    @staticmethod # known case of __new__
    def __new__(S, *more): # real signature unknown; restored from __doc__
        """ T.__new__(S, ...) -> a new object with type S, a subtype of T """
        pass

    def __reduce__(self, *args, **kwargs): # real signature unknown
        pass

    def __setstate__(self, *args, **kwargs): # real signature unknown
        pass

class MockFixedWindowIndexer(WindowIndexer):
    """
    We are just checking parameters of the indexer,
        and returning a consistent API with fixed/variable
        indexers.

        Parameters
        ----------
        values: ndarray
            values data array
        win: int64_t
            window size
        minp: int64_t
            min number of obs in a window to consider non-NaN
        index: object
            index of the values
        floor: optional
            unit for flooring
        left_closed: bint
            left endpoint closedness
        right_closed: bint
            right endpoint closedness
    """
    def __init__(self, *args, **kwargs): # real signature unknown
        pass

    @staticmethod # known case of __new__
    def __new__(S, *more): # real signature unknown; restored from __doc__
        """ T.__new__(S, ...) -> a new object with type S, a subtype of T """
        pass

    def __reduce__(self, *args, **kwargs): # real signature unknown
        pass

    def __setstate__(self, *args, **kwargs): # real signature unknown
        pass

class VariableWindowIndexer(WindowIndexer):
    """
    create a variable length window indexer object
        that has start & end, that point to offsets in
        the index object; these are defined based on the win
        arguments

        Parameters
        ----------
        values: ndarray
            values data array
        win: int64_t
            window size
        minp: int64_t
            min number of obs in a window to consider non-NaN
        index: ndarray
            index of the values
        left_closed: bint
            left endpoint closedness
            True if the left endpoint is closed, False if open
        right_closed: bint
            right endpoint closedness
            True if the right endpoint is closed, False if open
        floor: optional
            unit for flooring the unit
    """
    def build(self, *args, **kwargs): # real signature unknown
        pass

    def __init__(self, *args, **kwargs): # real signature unknown
        pass

    @staticmethod # known case of __new__
    def __new__(S, *more): # real signature unknown; restored from __doc__
        """ T.__new__(S, ...) -> a new object with type S, a subtype of T """
        pass

    def __reduce__(self, *args, **kwargs): # real signature unknown
        pass

    def __setstate__(self, *args, **kwargs): # real signature unknown
        pass

# variables with complex values

interpolation_types = {
    'higher': 2,
    'linear': 0,
    'lower': 1,
    'midpoint': 4,
    'nearest': 3,
}

__test__ = {}


text:
import numpy as np
import torch
import torch.nn.functional as F
import kornia
from typing import Union, Tuple, Sequence
from scipy.spatial.transform import Rotation
from typing import Optional

class Simulator:
    def __init__(self):
        fname = "/home/kejie/repository/fast_sdf/dist-model.txt"
        data = np.loadtxt(fname, comments='%', skiprows=5)
        dist = np.empty([80, 80, 5])

        for y in range(0, 80):
            for x in range(0, 80):
                idx = (y * 80 + x) * 23 + 3
                if (data[idx:idx + 5] < 8000).all():
                    dist[y, x, :] = 0
                else:
                    dist[y, x, :] = data[idx + 15: idx + 20]

        self.model = dist

    def undistort(self, x, y, z):

        i2 = int((z + 1) / 2)
        i1 = i2 - 1
        a = (z - (i1 * 2 + 1)) / 2
        x = int(x / 8)
        y = int(y / 6)
        f = (1 - a) * self.model[y, x, min(max(i1, 0), 4)
                                 ] + a * self.model[y, x, min(i2, 4)]

        if f == 0:
            return 0
        else:
            return z / f

    def simulate(self, clean_depth):

        # a = np.array(Image.open(inputpng)).astype(np.float32) / 1000.0
        img_h, img_w = clean_depth.shape
        b = np.copy(clean_depth)
        it = np.nditer(clean_depth, flags=[
                       'multi_index'], op_flags=['writeonly'])

        while not it.finished:

            # pixel shuffle
            x = min(
                max(round(it.multi_index[1] + np.random.normal(0, 0.25)), 0), img_w-1)
            y = min(
                max(round(it.multi_index[0] + np.random.normal(0, 0.25)), 0), img_h-1)

            # downsample
            d = b[y - y % 2, x - x % 2]

            # distortion
            d = self.undistort(x, y, d)

            # quantization and high freq noise
            if d == 0:
                it[0] = 0
            else:
                it[0] = 35.130 * 8 / \
                    round((35.130 / d + np.random.normal(0, 0.027778)) * 8)

            it.iternext()
        return clean_depth

def get_homogeneous(
    pts: Union['np.ndarray', 'torch.tensor']
    ) -> Union['np.ndarray', 'torch.tensor']:
    """ convert [(b), N, 3] pts to homogeneous coordinate

    Args:
        pts ([(b), N, 3] Union['np.ndarray', 'torch.tensor']): input point cloud

    Returns:
        homo_pts ([(b), N, 4] Union['np.ndarray', 'torch.tensor']): output point
            cloud

    Raises:
        ValueError: if the input tensor/array is not with the shape of [b, N, 3]
            or [N, 3]
        TypeError: if input is not either tensor or array
    """

    batch = False
    if len(pts.shape) == 3:
        batch = True
    elif len(pts.shape) == 2:
        pts = pts
    else:
        raise ValueError("only accept [b, n_pts, 3] or [n_pts, 3]")

    if isinstance(pts, torch.Tensor):
        ones = torch.ones_like(pts[..., :1])
        homo_pts = torch.cat([pts, ones], dim=-1)
    elif isinstance(pts, np.ndarray):
        ones = np.ones_like(pts[..., :1])
        homo_pts = np.concatenate([pts, ones], axis=-1)
    else:
        raise TypeError("wrong data type")
    return homo_pts

def get_aabb(pc: 'np.ndarray') -> 'np.ndarray':
    """ get aabb of a point cloud

    Args:
        pc ([N, 3] np.ndarray): input point cloud

    Returns:
        aabb ([2, 3] np.ndarray): a 3D bbox represent by
            [[x_min, y_min, z_min], [x_max, y_max, z_max]]
    """

    x_min, y_min, z_min = np.min(pc, axis=0)
    x_max, y_max, z_max = np.max(pc, axis=0)
    aabb = np.array([[x_min, y_min, z_min], [x_max, y_max, z_max]])
    return aabb

# def get_aabb(pc: 'np.ndarray', img_w: int, img_h: int) -> 'np.ndarray':
#     """ get aabb of a point cloud

#     Args:
#         pc ([N, 2] np.ndarray): input point cloud

#     Returns:
#         aabb ([2, 2] np.ndarray): a 2D bbox represent by
#             [[x_min, y_min], [x_max, y_max]]
#     """

#     x_min, y_min = np.min(pc, axis=0)
#     x_max, y_max = np.max(pc, axis=0)
#     x_min = max(0, x_min)
#     y_min = max(0, y_min)
#     x_max = min(img_w, x_max)
#     y_max = min(img_h, y_max)
#     aabb = np.array([[x_min, y_min], [x_max, y_max]])
#     return aabb

def depth2xyz(depth, intr_mat):
    """ convert depth map to xyz map

    Args:
        depth ([H, W] np.ndarray): depth map

    Returns:
        xyz ([H, W, 3] np.ndarray): xyz map
    """

    height, width = depth.shape
    fx, fy, cx, cy = intr_mat[0, 0], intr_mat[1, 1], intr_mat[0, 2], intr_mat[1, 2]

    urange = (
        np.arange(width, dtype=np.float32).reshape(1, -1).repeat(height, 0) - cx
    ) / fx
    vrange = (
        np.arange(height, dtype=np.float32).reshape(-1, 1).repeat(width, 1) - cy
    ) / fy
    xyz = np.stack([urange, vrange, np.ones(urange.shape)], axis=-1)
    xyz = xyz * depth.reshape(height, width, 1)
    return xyz

def angle2class(angles, num_classes=30):
    """ convert angles between [0, 180] to class index for classification

    Args:
        angles (np.ndarray): angle in radian

    Returns:
        out_class (np.ndarray): angle is converted to class, the number of which
            is defined in num_classes
    """
    y = torch.sin(angles)
    x = torch.cos(angles)
    angles = torch.atan2(y, x) / np.pi * 180.
    angles = torch.where(angles<0, angles + 180, angles)
    out_class = angles // (180 / num_classes)
    assert (out_class >= 0).all()
    assert (out_class <= num_classes).all()
    out_class = np.clip(out_class, a_min=0, a_max=num_classes-1)
    return out_class

def iou_2d(bboxA: 'np.ndarray', bboxB: 'np.ndarray') -> float:
    """ calculate IoU between two 2D bboxes

    Args:
        bboxA ([2, 2] np.ndarray): input bbox A in AABB format
        bboxB ([2, 2] np.ndarray): input bbox B in AABB format

    Returns:
        IoU (float): output IoU
    """

    x_min = max(bboxA[0, 0], bboxB[0, 0])
    y_min = max(bboxA[0, 1], bboxB[0, 1])
    x_max = min(bboxA[1, 0], bboxB[1, 0])
    y_max = min(bboxA[1, 1], bboxB[1, 1])

    inter_area = max(0, (x_max - x_min)) * max(0, (y_max - y_min))
    area_A = np.prod(bboxA[1] - bboxA[0])
    area_B = np.prod(bboxB[1] - bboxB[0])
    IoU = inter_area / (area_A + area_B - inter_area)
    assert IoU <= 1 and IoU >= 0, "invalid IoU value"
    return IoU

def iou_3d(bboxA: 'np.ndarray', bboxB: 'np.ndarray') -> float:
    """ calculate 3D IoU between two 3D bboxes

    Args:
        bboxA ([2, 3] np.ndarray): input bbox A in AABB format
        bboxB ([2, 3] np.ndarray): input bbox B in AABB format

    Returns:
        IoU (float): 3D IoU
    """

    x_min = max(bboxA[0, 0], bboxB[0, 0])
    y_min = max(bboxA[0, 1], bboxB[0, 1])
    z_min = max(bboxA[0, 2], bboxB[0, 2])
    x_max = min(bboxA[1, 0], bboxB[1, 0])
    y_max = min(bboxA[1, 1], bboxB[1, 1])
    z_max = min(bboxA[1, 2], bboxB[1, 2])

    inter_volume = max(0, (x_max - x_min)) * max(0, (y_max - y_min)) * max(0, (z_max - z_min))
    volume_A = np.prod(bboxA[1] - bboxA[0])
    volume_B = np.prod(bboxB[1] - bboxB[0])
    IoU = inter_volume / (volume_A + volume_B - inter_volume)
    assert IoU <= 1 and IoU >= 0, "invalid IoU value"
    return IoU

def giou_3d(bboxA: 'np.ndarray', bboxB: 'np.ndarray') -> float:
    """ calculate generalized 3D IoU between two 3D bboxes

    Args:
        bboxA ([2, 3] np.ndarray): input bbox A in AABB format
        bboxB ([2, 3] np.ndarray): input bbox B in AABB format

    Returns:
        IoU (float): 3D Generalized IoU
    """

    x_min = max(bboxA[0, 0], bboxB[0, 0])
    y_min = max(bboxA[0, 1], bboxB[0, 1])
    z_min = max(bboxA[0, 2], bboxB[0, 2])
    x_max = min(bboxA[1, 0], bboxB[1, 0])
    y_max = min(bboxA[1, 1], bboxB[1, 1])
    z_max = min(bboxA[1, 2], bboxB[1, 2])

    inter_volume = max(0, (x_max - x_min)) * max(0, (y_max - y_min)) * max(0, (z_max - z_min))
    volume_A = np.prod(bboxA[1] - bboxA[0])
    volume_B = np.prod(bboxB[1] - bboxB[0])
    volume_union = (volume_A + volume_B - inter_volume)

    iou = iou_3d(bboxA, bboxB)

    x_min = min(bboxA[0, 0], bboxB[0, 0])
    y_min = min(bboxA[0, 1], bboxB[0, 1])
    z_min = min(bboxA[0, 2], bboxB[0, 2])
    x_max = max(bboxA[1, 0], bboxB[1, 0])
    y_max = max(bboxA[1, 1], bboxB[1, 1])
    z_max = max(bboxA[1, 2], bboxB[1, 2])

    volume_complete = (x_max - x_min) * (y_max - y_min) * (z_max - z_min)
    giou = iou - (volume_complete - volume_union) / volume_complete
    return giou

def get_corner_by_dims(dimensions) -> np.ndarray:
    """get 8 corner points of 3D bbox defined by self.dimensions

    Returns:
        a np.ndarray with shape [8,3] to represent 8 corner points'
        position of the 3D bounding box.
    """

    w, h, l = dimensions[0], dimensions[1], dimensions[2]
    x_corners = [w/2, -w/2, -w/2, w/2, w/2, -w/2, -w/2, w/2]
    y_corners = [h/2, h/2, h/2, h/2, -h/2, -h/2, -h/2, -h/2]
    z_corners = [l/2, l/2, -l/2, -l/2, l/2, l/2, -l/2, -l/2]
    corner_pts = np.array([x_corners, y_corners, z_corners], dtype=np.float32).T
    return corner_pts

def scale_vertices_to_bbox(
    vertices: 'np.ndarray', bbox_dims: 'np.ndarray'
    ) -> 'np.ndarray':
    """scale the vertices such that they are tightly bounded by the 3D bbox

    Args:
        vertices ([N, 3] np.ndarray): input vertices
        bbox_dims ([3] np.ndarray): bbox dimension in x, y, z axis

    Returns:
        vertices: the scaled vertices
    """

    vertices[:, 0] *= (bbox_dims[0] / (np.max(vertices[:, 0]) - np.min(vertices[:, 0])))
    vertices[:, 1] *= (bbox_dims[1] / (np.max(vertices[:, 1]) - np.min(vertices[:, 1])))
    vertices[:, 2] *= (bbox_dims[2] / (np.max(vertices[:, 2]) - np.min(vertices[:, 2])))

    return vertices

def unproject(pixel, depth, intr_mat):
    """ unproject from pixels and depths to 3D

    Args:
        pixel: [n, 2]
        depth: [n]
    """
    fx = intr_mat[0, 0]
    fy = intr_mat[1, 1]
    cx = intr_mat[0, 2]
    cy = intr_mat[1, 2]
    pts = np.concatenate([pixel, np.ones_like(pixel)[:, :1]], axis=1)
    pts[:, 0] = (pts[:, 0] - cx) / fx
    pts[:, 1] = (pts[:, 1] - cy) / fy
    pts = pts * depth[:, None]
    return pts

def projection(pts, intr_mat, keep_z=False):
    """perspective projection

    Args:
        pts ([(b), N, 3] or [(b), N, 4] np.ndarray or torch.tensor): 3D points
        intr_mat ([(b), 3, 3] or [(b), 3, 4] np.ndarray or torch.tensor): intrinsic
            matrix

    Returns:
        pts ([(b), N, 3], np.ndarray or torch.tensor): projected points
    """

    batch = False
    if len(pts.shape) == 3:
        assert len(intr_mat.shape) == 3, "intr_mat shape needs to match pts"
        batch = True
    elif len(pts.shape) == 2:
        assert len(intr_mat.shape) == 2, "intr_mat shape needs to match pts"
    else:
        ValueError("only accept [b, n_pts, 3] or [n_pts, 3]")
    if batch:
        if isinstance(pts, torch.Tensor):
            intr_mat = intr_mat.transpose(1, 2)
        else:
            intr_mat = intr_mat.transpose(0, 2, 1)
    else:
        intr_mat = intr_mat.T
    pts = pts @ intr_mat
    if isinstance(pts, torch.Tensor):
        z = torch.ones_like(pts[..., -1])
    else:
        z = np.ones_Like(pts[..., -1])
    if batch:
        if keep_z:
            z = pts[:, :, -1]
        pts = pts / pts[:, :, -1:]
        pts[:, :, -1] *= z
    else:
        if keep_z:
            z = pts[:, -1]
        pts = pts / pts[:, -1:]
        pts[:, -1] *= z
    return pts

def pad_transform_matrix(mat: 'np.ndarray') -> 'np.ndarray':
    """ pad a [3, 4] transform matrix to a [4, 4] matrix

    Args:
        mat ([3, 4] np.ndarray): the input [3, 4] matrix
    Returns:
        mat ([4, 4] np.ndarray): the output [4, 4] matrix
    """

    if mat.shape[0] < 4:
        pad = np.zeros((1, 4), dtype=np.float32)
        pad[0,-1] = 1
        return np.concatenate([mat, pad], axis=0)
    else:
        return mat

def rgbd_to_colored_pc(
    rgb: 'np.ndarray',
    depth: 'np.ndarray',
    fx: float,
    fy: float,
    cx: float,
    cy: float,
    cap: float = 200) -> Tuple['np.ndarray', 'np.ndarray']:
    """ convert a pair of rgb and depth iamge to a 3D colored point cloud

    Args:
        rgb ([H, W, 3] np.ndarray): rgb image
        depth ([H, W] np.ndarray): depth image
        fx, fy, cx, cy (float, float, float, float): camera intrinsic matrix
        cap (float): depth capping value

    Returns:
        a tuple containing:
            points ([N, 3] np.ndarray): 3D point positions
            colors ([N, 3] np.ndarray): color for each point
    """

    rgb_height, rgb_width, _ = rgb.shape
    X, Y = np.meshgrid(np.arange(rgb_width), np.arange(rgb_height))
    xyz_rgb = np.concatenate(
        [X[:, :, None], Y[:, :, None], depth[:, :, None], rgb],
        axis=2
    )
    xyz_rgb[:, :, 0] = (xyz_rgb[:, :, 0] - cx) * xyz_rgb[:, :, 2] / fx
    xyz_rgb[:, :, 1] = (xyz_rgb[:, :, 1] - cy) * xyz_rgb[:, :, 2] / fy
    points = xyz_rgb[:, :, :3].reshape(-1, 3)
    colors = xyz_rgb[:, :, 3:].reshape(-1, 3) / 255.
    cap_ind = np.logical_and((points[:, 2] < cap), (points[:, 2] > 0))
    points = points[cap_ind]
    colors = colors[cap_ind]
    return points, colors

def geodesic_distance(R1: 'np.ndarray', R2: 'np.ndarray') -> float:
    '''Returns the geodesic distance between two rotation matrices.

    Args:
        R1 ([3, 3] np.ndarray): input rotation matrix
        R2 ([3, 3] np.ndarray): input rotation matrix

    Returns:
        delta_theta (float): geodesic distance between the input rotation
            matrices
    '''

    delta_R = np.dot(R1, R2.T)
    rotvec = Rotation.from_dcm(delta_R).as_rotvec()
    delta_theta = np.linalg.norm(rotvec)
    return delta_theta

def pts_in_box(pts: 'np.ndarray', img_shape: 'np.ndarray') -> 'np.ndarray':
    """ check projected points are within image frame

    Args:
        pts ([N, 2] np.ndarray): a set of 2D points on image plane
        img_shape (aabb): bbox_size [x_min, y_min, x_max, y_max]
    Return:
        a boolean array of shape [N] indicating whether a point is within
        image frame
    """

    img_shape = img_shape.reshape(2, 2)
    larger_x_min = pts[:, 0] > img_shape[0, 0]
    smaller_x_max = pts[:, 0] < img_shape[1, 0]
    larger_y_min = pts[:, 1] > img_shape[0, 1]
    smaller_y_max = pts[:, 1] < img_shape[1, 1]
    return (larger_x_min * smaller_x_max * \
        larger_y_min * smaller_y_max)

def get_normalize_factor(pts):
    norm = np.max(np.linalg.norm(pts, axis=1))
    return norm

def create_meshgrid(
    height: int,
    width: int,
    normalized_coordinates: bool = True,
    device: Optional[torch.device] = torch.device('cpu'),
    dtype: torch.dtype = torch.float32,
) -> torch.Tensor:
    """Generates a coordinate grid for an image.

    When the flag `normalized_coordinates` is set to True, the grid is
    normalized to be in the range [-1,1] to be consistent with the pytorch
    function grid_sample.
    http://pytorch.org/docs/master/nn.html#torch.nn.functional.grid_sample

    Args:
        height (int): the image height (rows).
        width (int): the image width (cols).
        normalized_coordinates (bool): whether to normalize
          coordinates in the range [-1, 1] in order to be consistent with the
          PyTorch function grid_sample.
        device (torch.device): the device on which the grid will be generated. Default: cpu.
        dtype (torch.dtype): the data type of the generated gird. Default: float32.

    Return:
        torch.Tensor: returns a grid tensor with shape :math:`(1, H, W, 2)`.
    """
    xs: torch.Tensor = torch.linspace(0, width - 1, width, device=device, dtype=dtype)
    ys: torch.Tensor = torch.linspace(0, height - 1, height, device=device, dtype=dtype)
    if normalized_coordinates:
        xs = (xs / (width - 1) - 0.5) * 2
        ys = (ys / (height - 1) - 0.5) * 2
    # generate grid by stacking coordinates
    base_grid: torch.Tensor = torch.stack(torch.meshgrid([xs, ys])).transpose(1, 2)  # 2xHxW
    return torch.unsqueeze(base_grid, dim=0).permute(0, 2, 3, 1)  # 1xHxWx2

def calculate_normals(pts):
    """
    pts: [B, H, W, 3]
    """
    # compute the pointcloud spatial gradients
    pts = pts.permute(0, 3, 1, 2)
    gradients: torch.Tensor = kornia.filters.spatial_gradient(pts)  # Bx3x2xHxW

    # compute normals
    a, b = gradients[:, :, 0], gradients[:, :, 1]  # Bx3xHxW

    normals: torch.Tensor = torch.cross(a, b, dim=1)  # Bx3xHxW
    return F.normalize(normals, dim=1, p=2)

def calculate_normals_finite_difference(pts, func):
    """
    compute the normal of an implicit reprentation using finite difference

    pts: (B, N, 3)

    """

    batch_size, n_pts = pts.shape[:2]
    epsilon = np.abs(np.random.normal(scale=0.01, size=(batch_size, n_pts)))
    epsilon = torch.from_numpy(epsilon).float().to(pts.device)
    epsilon = torch.clamp(epsilon, min=1e-3, max=1)
    epsilon_x = torch.stack(
        [epsilon, torch.zeros_like(epsilon), torch.zeros_like(epsilon)],
        dim=-1
    ).detach()
    epsilon_y = torch.stack(
        [torch.zeros_like(epsilon), epsilon, torch.zeros_like(epsilon)],
        dim=-1
    ).detach()
    epsilon_z = torch.stack(
        [torch.zeros_like(epsilon), torch.zeros_like(epsilon), epsilon],
        dim=-1
    ).detach()

    pts_x0 = pts - epsilon_x
    pts_x1 = pts + epsilon_x

    pts_y0 = pts - epsilon_y
    pts_y1 = pts + epsilon_y

    pts_z0 = pts - epsilon_z
    pts_z1 = pts + epsilon_z
    f_out = func(
        torch.cat([pts_x0, pts_x1, pts_y0, pts_y1, pts_z0, pts_z1], dim=0)
    )
    f_x0, f_x1, f_y0, f_y1, f_z0, f_z1 = torch.split(f_out, batch_size, dim=0)

    g_x = (f_x1 - f_x0) / (2 * epsilon)
    g_y = (f_y1 - f_y0) / (2 * epsilon)
    g_z = (f_z1 - f_z0) / (2 * epsilon)

    normals = torch.stack([g_x, g_y, g_z], dim=-1)
    normals = normals / (torch.norm(normals, dim=-1, keepdim=True) + 1e-5)
    return normals

def recenter(vertices):
    min_ = np.min(vertices, axis=0)
    max_ = np.max(vertices, axis=0)
    center = (max_ + min_) / 2.
    vertices = vertices - center[None, :]
    return vertices, center

reason: import numpy,from scipy


text:
#!/usr/local/bin/python
"""
Pindrop interview question.

Sinan,

Here is the question, it was nice talking to you:

You are a magician, and as such, love to hide balls under cups.  You have
exactly 10,000 balls and 10,000 cups, each with a different diameter drawn from
 a random distribution. (You may choose any distribution of diameters, as long
 as it is the same distribution for both balls and cups).

At each iteration, you randomly pair up all balls with cups.
(A pair consists of a single ball and single cup, and each ball/cup can only
belong to a single pair). For each pair, if the diameter of the ball is less
than the diameter of the cup, the ball gets hidden and you set the pair aside
for your show later. Repeat this process until it is impossible for any more
balls to be hidden under cups.

On average, what percentage of balls/cups are used?

Please send code with comments describing your thought process, what your code
is doing, etc. If you get stuck just explain why, or what you intended to do
next.

Thanks,
John

Short Answer: around 0.6 - 1.0%
"""
import numpy as np
import logging

DEBUG = False

logger = logging.getLogger('pindrop')
logging.basicConfig()
LOG_LEVEL = logging.DEBUG if DEBUG else logging.INFO
logger.setLevel(LOG_LEVEL)

NUMBER_OF_BALLS = 10000
# number_of_threads = 3
MAX_DIAMETER = 10

SMALLEST_CUP_SIZE = 0
SMALLEST_BALL_SIZE = 0

# The Ball/Cup pairs are always in the same structure in the list
BALL_AXIS = 0
CUP_AXIS = 1

ORDER_BY_CUPS_DESC = ['f{0}'.format(BALL_AXIS)]
ORDER_BY_BALLS_DESC = ['f{0}'.format(CUP_AXIS)]

def get_random_diameter():
    """Return a random diameter value."""
    return np.random.random_integers(1, MAX_DIAMETER)

def calculate_usage_percentage(used_pairs, unused_pairs):
    """Calculate how many of the balls are paired in an iteration."""
    total = len(used_pairs) + len(unused_pairs)
    logger.debug('used: {0}\ntotal :{1}'.format(len(used_pairs), total))
    percentage = float(len(used_pairs)) * 100.0 / float(total)
    logger.info('percentage: {0}'.format(percentage))
    return percentage

def calculate_impossibility(pairs):
    """
    Calculate if its impossible to still match balls with cups.

    Parameters
    ----------
    pairs : List. A list of (Ball, Cup) tuples.

    Returns
    -------
    Boolean
        Impossibility value.

    """
    if not pairs or len(pairs) < 2:
        return True
    _np_pairs_array = np.array(pairs)

    # Sort the pairs with respect to the ball/cup colon
    # By sorting the arrays we can get the smallest and largest diameter
    # We need to look at the largest and smallest one to make a conclusion
    min_ball_diameter_array = np.sort(
        _np_pairs_array.view('i8,i8'), order=ORDER_BY_BALLS_DESC, axis=0
    ).view(np.int)
    max_cup_diameter_array = np.sort(
        _np_pairs_array.view('i8,i8'), order=ORDER_BY_CUPS_DESC, axis=0
    ).view(np.int)

    # first element of the sorted ball/cup array
    min_ball_size = min_ball_diameter_array[0][BALL_AXIS]
    max_cup_size = max_cup_diameter_array[-1][CUP_AXIS]

    impossibility_reached = False
    if max_cup_size <= min_ball_size:
        impossibility_reached = True
    return impossibility_reached

def randomly_match_cup_ball_pairs(pairs, max_range=None):
    """
    Randomly match the cups and balls.

    If max_range is given, create the cups and balls.

    Parameters
    ----------
    pairs : List. A list of (Ball, Cup) tuples.

    max_range : Integer. The maximum number of ball/cup objects to create.

    Returns
    -------
    dict
        the used and unused pair lists with the matching metadata.

    """
    unused_pairs = []
    used_pairs = []

    count = max_range if max_range is not None else len(pairs)
    # This loop can be divided into smaller chunks to be threaded or
    # made into async jobs for both creation and sorting.
    for i in range(0, count):
        if max_range:
            ball = cup = get_random_diameter()
        else:
            if len(pairs) <= 2:
                # Edge case, when left with 2 pairs we need to be careful
                impossibility_reached = calculate_impossibility(pairs)
                if impossibility_reached:
                    unused_pairs.extend(pairs)
                    return {
                        'used_pairs': used_pairs,
                        'unused_pairs': unused_pairs
                    }
            # take a random pair from each batch
            random_index = np.random.random_integers(0, len(pairs) - 1)
            pair = pairs.pop(random_index)
            ball = pair[BALL_AXIS]

            random_index = np.random.random_integers(0, len(pairs) - 1)
            pair = pairs.pop(random_index)
            cup = pair[CUP_AXIS]
            # when we pop a ball and a cup (to use)there is also an unused pair
            unused_pairs.append((pair[BALL_AXIS], pair[CUP_AXIS]))

        if ball < cup:
            used_pairs.append((ball, cup))
        else:
            unused_pairs.append((ball, cup))

    # make sure unused pairs is always a list and not None
    return {
        'used_pairs': used_pairs,
        'unused_pairs': unused_pairs if unused_pairs else []
    }

def main():
    """Main program."""
    percent_used = []
    matched_pairs = []
    # Create the initial ball/cup sets with random diameters.
    # This function also randomly pairs balls with cups
    # From now on we will always work with pairs (balls, cup)
    result = randomly_match_cup_ball_pairs([], NUMBER_OF_BALLS)
    matched_pairs.extend(result.get('used_pairs'))
    # calculate the initial usage percentage
    percent_used.append(calculate_usage_percentage(
        result.get('used_pairs'),
        result.get('unused_pairs')
    ))
    impossibility_reached = False
    while not impossibility_reached:
        # We do a check initially for impossibilities
        impossibility_reached = calculate_impossibility(
            result.get('unused_pairs'))
        if impossibility_reached:
            break
        # we match every unused pair of ball cups until we cant anymore
        result = randomly_match_cup_ball_pairs(result.get('unused_pairs'))
        matched_pairs.extend(result.get('used_pairs'))

        # can be async
        if not impossibility_reached:
            percent_used.append(calculate_usage_percentage(
                result.get('used_pairs'),
                result.get('unused_pairs')
            ))

    logger.debug('Percentage of hidden balls at each iteration: {0}'.format(
        percent_used))
    average_ball_cup_usage = np.average(percent_used)
    logger.info('Average percentage of ball/cup usage: {0}'.format(
        average_ball_cup_usage))

if __name__ == '__main__':
    main()


text:
# --- Python 3.8 ---
"""
@File    :   test_con_derivs.py
@Time    :   2020/12/19
@Author  :   Josh Anibal
@Desc    :   Test the derivatives of the geometric constraint functions
"""

# === Standard Python modules ===
import unittest
import os

# === External Python modules ===
import numpy as np

from mpi4py import MPI

# === Extension modules ===
import openmdao.api as om
from openmdao.utils.assert_utils import assert_near_equal

import os
from mphys.multipoint import Multipoint

# for geometric DVs
from pygeo.mphys import OM_DVGEOCOMP

# only try to import this so that people can run the script w/o mdolab code
# TACS is required regardless of the structural solver used
from adflow import ADFLOW

baseDir = os.path.dirname(os.path.abspath(__file__))

class Top(om.Group):
    def setup(self):

        ################################################################################
        # MPHYS setup
        ################################################################################

        # ivc for DVs
        dvs = self.add_subsystem("dvs", om.IndepVarComp(), promotes=["*"])
        # geometry parametrization with FFD and general geometric constraints
        self.add_subsystem(
            "geo",
            OM_DVGEOCOMP(
                ffd_file=os.path.join(baseDir, "../input_files/ffd.xyz"),
            ),
        )

        # # create the multiphysics multipoint group.
        # mp = self.add_subsystem(
        #     "mp", Multipoint(aero_builder=None, struct_builder=None, xfer_builder=None)
        # )

        # mp.mphys_add_scenario("s0")

    def configure(self):

        # ################################################################################
        # # GEOMETRY SETUP
        # ################################################################################

        from adflow.mphys import ADflowBuilder

        aero_options = {
            # I/O Parameters
            "gridFile": os.path.join(baseDir, "../input_files/wing_vol_L3.cgns"),
            "outputDirectory": os.path.join(baseDir, "../output_files"),
            "monitorvariables": ["resrho", "resturb", "cl", "cd"],
            "surfacevariables": ["cp", "vx", "vy", "vz", "mach"],
            # 'isovariables': ['shock'],
            "isoSurface": {"shock": 1},  # ,'vx':-0.0001},
            "writeTecplotSurfaceSolution": False,
            "writevolumesolution": False,
            "writesurfacesolution": False,
            "liftindex": 3,
            # Physics Parameters
            "equationType": "RANS",
            # Solver Parameters
            "smoother": "DADI",
            "CFL": 1.5,
            "MGCycle": "sg",
            "MGStartLevel": -1,
            # ANK Solver Parameters
            "useANKSolver": True,
            "nsubiterturb": 5,
            "anksecondordswitchtol": 1e-4,
            "ankcoupledswitchtol": 1e-6,
            "ankinnerpreconits": 2,
            "ankouterpreconits": 1,
            "anklinresmax": 0.1,
            "infchangecorrection": True,
            "ankcfllimit": 1e4,
            # NK Solver Parameters
            # 'useNKSolver':True,
            "nkswitchtol": 1e-4,
            # Termination Criteria
            "L2Convergence": 1e-15,
            "L2ConvergenceCoarse": 1e-2,
            "L2ConvergenceRel": 1e-14,
            "nCycles": 10000,
            "adjointl2convergencerel": 1e-14,
            "adjointl2convergence": 1e-14
            # force integration
        }

        aero_builder = ADflowBuilder(aero_options)

        aero_builder.initialize(self.comm)
        aero_mesh = aero_builder.get_mesh_coordinate_subsystem()

        aero_mesh.aero_solver = aero_builder.get_solver()

        tri_points = aero_mesh._getTriangulatedMeshSurface()

        points = {}
        points["aero_points"] = aero_mesh.aero_solver.getSurfaceCoordinates(includeZipper=False).flatten(order="C")

        # add these points to the geometry object
        self.geo.nom_add_point_dict(points)
        # create constraint DV setup
        # tri_points = self.mp.mphys_get_triangulated_surface()
        # self.geo.nom_setConstraintSurface(tri_points)

        self.geo.nom_setConstraintSurface(tri_points)

        # add geometric DVs
        nRefAxPts = self.geo.nom_addRefAxis(name="wing", xFraction=0.25, alignIndex="k")
        nTwist = nRefAxPts - 1

        # Set up global design variables
        def twist(val, geo):
            for i in range(1, nRefAxPts):
                geo.rot_y["wing"].coef[i] = val[i - 1]

        self.geo.nom_addGeoDVGlobal(dvName="twist", value=np.zeros(nTwist), func=twist)
        self.dvs.add_output("twist", val=np.array([0] * nTwist))
        self.connect("twist", "geo.twist")

        # # shape DVs
        n_dv = self.geo.nom_addGeoDVLocal(dvName="local", axis="z")
        self.dvs.add_output("local", val=np.array([0] * n_dv))
        self.connect("local", "geo.local")

        # if we have shape DVs, we also need to have geometric constraints

        # first set the triangulated surface for projections
        # tri_points = self.mp.mphys_get_triangulated_surface()
        # self.geo.nom_setConstraintSurface(tri_points)

        # Le/Te constraints
        self.geo.nom_add_LETEConstraint("lecon", 0, "iLow")
        self.geo.nom_add_LETEConstraint("tecon", 0, "iHigh")

        # Volume constraints
        leList = [[0.1, 0.001, 0], [0.1 + 7.5, 14, 0]]
        teList = [[4.2, 0.001, 0], [8.5, 14, 0]]
        self.geo.nom_addVolumeConstraint("volcon", leList, teList, 20, 20)

        # # Thickness constraints
        self.geo.nom_addThicknessConstraints2D("2Dvolcon", leList, teList, 10, 10)

class TestDVCon(unittest.TestCase):
    def setUp(self):
        ################################################################################
        # OpenMDAO setup
        ################################################################################
        prob = om.Problem()
        prob.model = Top()

        prob.model.add_design_var("twist")
        prob.model.add_design_var("local")

        # objectives and nonlinear constraints
        prob.model.add_constraint("geo.aero_points", indices=[0, 10, 100])
        prob.model.add_constraint("geo.lecon")
        prob.model.add_constraint("geo.tecon")
        prob.model.add_constraint("geo.volcon")
        prob.model.add_constraint("geo.2Dvolcon")

        prob.setup(mode="rev")

        self.prob = prob
        # om.n2(
        #     prob,
        #     show_browser=False,
        #     outfile="test_con_derivs.html",
        # )

    def test_run_model(self):
        self.prob.run_model()

    def test_derivatives(self):
        self.prob.run_model()
        print("----------------strating check totals--------------")
        # data = self.prob.check_totals(step=1e-7, form="forward")
        data = self.prob.check_totals(step=1e-8, compact_print=True)  # out_stream=None
        for var, err in data.items():

            rel_err = err[
                "abs error"
            ]  # abs error used here because for some values should be zero so the ref err is undefined
            assert_near_equal(rel_err.forward, 0.0, 1e-6)

if __name__ == "__main__":
    unittest.main()

max_stars_repo_name: timryanb/mphys

max_issues_repo_name: timryanb/mphys

max_forks_repo_name: timryanb/mphys


text:
import os
import gzip
import numpy as np

import paddle.v2 as paddle

from utils import logger, load_dict
from beam_search import BeamSearch
import config as conf
from network_conf import rnn_lm

def rnn_generate(gen_input_file, model_path, max_gen_len, beam_size,
                 word_dict_file):
    """
    use RNN model to generate sequences.

    :param word_id_dict: vocab.
    :type word_id_dict: dictionary with content of "{word, id}",
                        "word" is string type , "id" is int type.
    :param num_words: the number of the words to generate.
    :type num_words: int
    :param beam_size: beam width.
    :type beam_size: int
    :return: save prediction results to output_file
    """

    assert os.path.exists(gen_input_file), "test file does not exist!"
    assert os.path.exists(model_path), "trained model does not exist!"
    assert os.path.exists(
        word_dict_file), "word dictionary file does not exist!"

    # load word dictionary
    word_2_ids = load_dict(word_dict_file)
    try:
        UNK_ID = word_2_ids["<unk>"]
    except KeyError:
        logger.fatal("the word dictionary must contain a <unk> token!")
        sys.exit(-1)

    # initialize paddle
    paddle.init(use_gpu=conf.use_gpu, trainer_count=conf.trainer_count)

    # load the trained model
    pred_words = rnn_lm(
        len(word_2_ids),
        conf.emb_dim,
        conf.hidden_size,
        conf.stacked_rnn_num,
        conf.rnn_type,
        is_infer=True)

    parameters = paddle.parameters.Parameters.from_tar(
        gzip.open(model_path, "r"))

    inferer = paddle.inference.Inference(
        output_layer=pred_words, parameters=parameters)

    generator = BeamSearch(inferer, word_dict_file, beam_size, max_gen_len)
    # generate text
    with open(conf.gen_file, "r") as fin, open(conf.gen_result, "w") as fout:
        for idx, line in enumerate(fin):
            fout.write("%d\t%s" % (idx, line))
            for gen_res in generator.gen_a_sentence([
                    word_2_ids.get(w, UNK_ID)
                    for w in line.lower().strip().split()
            ]):
                fout.write("%s\n" % gen_res)
            fout.write("\n")

if __name__ == "__main__":
    rnn_generate(conf.gen_file, conf.model_path, conf.max_gen_len,
                 conf.beam_size, conf.vocab_file)


text:
import sys

import cv2
import numpy as np

import obj_loader
from utils import (do_everything, compare_markers, four_point_transform, convert_to_grayscale, draw_rectangle, find_homographies, get_matrix, get_homography_from_corners,
                   load_ref_images, calculate_dist_matches, calculate_dist_corners, render, display_image_with_matched_keypoints)

if __name__ == "__main__":
    OBJ_PATH = sys.argv[1]

    OBJ = obj_loader.OBJ(OBJ_PATH, swapyz=True)
    REF_IMAGES, REF_DSC = load_ref_images()
    VID_FEED = cv2.VideoCapture(-1)

    REACHED_X, REACHED_Y = 0, 0

    # iterate = 0

    # flag_g = True
    # flag_b = True
    methods = [
    ("THRESH_BINARY", cv2.THRESH_BINARY),
    ("THRESH_BINARY_INV", cv2.THRESH_BINARY_INV),
    ("THRESH_TRUNC", cv2.THRESH_TRUNC),
    ("THRESH_TOZERO", cv2.THRESH_TOZERO),
    ("THRESH_TOZERO_INV", cv2.THRESH_TOZERO_INV)]

    while True:
        # box = None
        # box_g = None
        RET, FRAMED = VID_FEED.read()
        if not RET:
            print("Unable to capture video")
            sys.exit()

        SIZE = FRAMED.shape
        FOC_LEN = SIZE[1]
        CENTER = (SIZE[1]/2, SIZE[0]/2)

        CAM_MAT = np.array(
            [[FOC_LEN, 0, CENTER[0]],
             [0, FOC_LEN, CENTER[1]],
             [0, 0, 1]], dtype="double")

        REACHED_X, REACHED_Y = do_everything(FRAMED, REF_IMAGES, OBJ, CAM_MAT, REACHED_X, REACHED_Y)
        if cv2.waitKey(1) & 0xFF == ord("q"):
            break

        # FRAMED = render(FRAMED, OBJ, PROJ_MAT, REF_IMAGE1, False)
        cv2.imshow("frame", FRAMED)
        # cv2.imshow("frame", thresh)
        if cv2.waitKey(1) & 0xFF == ord("q"):
            break

    VID_FEED.release()
    cv2.destroyAllWindows()


text:
import collections
import json

import jsonpickle
import numpy as np
import tensorflow as tf
from tensorflow.python.framework import ops
from tensorflow.python.tpu import tpu_feed

from .. import tf_wrapper as tfw
from ..dataclass import ModelParameter

tf1 = tf.compat.v1
Dataset = tf1.data.Dataset

def place_dataloader(params: ModelParameter, input_fn):
    num_cores = params.mesh_impl.device_assignment.num_replicas

    ordered_ordinals = []
    ordered_hosts = []
    ordered_host_ids = []
    host_id_to_its_pnums = collections.defaultdict(list)
    d_assignment = params.mesh_impl.device_assignment

    for pnum in range(num_cores):
        physical_pnum = params.mesh_impl.l2p(pnum)

        # For MTF, there's always 1 core per replica. So logical_core=0.
        ordered_ordinals.append(d_assignment.tpu_ordinal(replica=physical_pnum, logical_core=0))
        host_device = d_assignment.host_device(replica=physical_pnum)
        host_id = int(host_device.lower().split("/task:")[1].split("/device:")[0])
        ordered_hosts.append(host_device)
        ordered_host_ids.append(host_id)
        host_id_to_its_pnums[host_id].append(pnum)

    num_hosts = len(set(ordered_hosts))

    pnum_maps = []
    macro_batching_multi = params.macro_batching if params.train else 1
    batch_size = params.input_pipeline_shape[0].to_integer_list[0] * macro_batching_multi
    for mtf_shape in params.input_pipeline_shape:
        # Make sure that the batch size is the same across all input tensors.
        assert batch_size == mtf_shape.to_integer_list[0] * macro_batching_multi

        s_shape = params.mesh_impl.slice_shape(mtf_shape)
        shape_list = [dim_size // s_dim_size for dim_size, s_dim_size in zip(mtf_shape.to_integer_list, s_shape)]

        pnum_map_shape = shape_list + [num_cores // np.prod(shape_list)]
        assert np.prod(pnum_map_shape) == num_cores

        # Initialize the pnum_map to None.
        pnum_map = np.empty(pnum_map_shape, dtype=object)
        pnum_map[:] = None

        for pnum in range(num_cores):
            s_begin = params.mesh_impl.slice_begin(mtf_shape, pnum)
            coord = [dim_size // s_dim_size for dim_size, s_dim_size in zip(s_begin, s_shape)]
            # put pnum in pnum_map[coord]
            pnum_array_ref = pnum_map[tuple(coord)]
            for idx, value in enumerate(pnum_array_ref):
                if value is None:
                    pnum_array_ref[idx] = pnum
                    break

        pnum_maps.append(pnum_map)

    # For each sub-batch, we need to know which host should read it.
    if params.train:

        # This records how many datasets (ds) are already stored on each host.
        num_dss_per_host = [0] * num_hosts

        # A list of host_ids that holds datasets (ds).
        hosts_to_hold_ds = []

        for sub_batch_pnum_map in pnum_maps[0]:

            num_pnums_per_host = [0] * num_hosts
            for pnum in sub_batch_pnum_map.flatten():
                num_pnums_per_host[ordered_host_ids[pnum]] += 1

            host_metrics = [(host_id, num_pnums_per_host[host_id], num_dss_per_host[host_id]) for host_id in
                            range(num_hosts)]
            host_id, _, _ = max(host_metrics, key=lambda keys: (keys[1], -keys[2]))

            num_dss_per_host[host_id] += 1
            hosts_to_hold_ds.append(host_id)

    else:
        # There should be just one dataset-holding host. Make the last host do it.
        hosts_to_hold_ds = [num_hosts - 1]

    sub_batch_size = batch_size // len(hosts_to_hold_ds)
    tf1.logging.info("MTF sub_batch_size: {}".format(sub_batch_size))
    assert sub_batch_size * len(hosts_to_hold_ds) == batch_size

    # Slots for all laidout tensors.
    all_laidout_tensors = [[None] * len(params.input_pipeline_shape) for _ in range(num_cores)]

    log_path = params.model_path + "/DataLog.log"
    _run_log = []
    run_log = None

    if params.use_checkpointing:
        if tf.io.gfile.exists(log_path):
            _run_log = json.load(tf.io.gfile.GFile(log_path, 'r'))

        curran_stats = {'steps': params.current_step, 'ctx': params.sequence_length,
                        'slice_count': len(hosts_to_hold_ds),
                        'interleave_size': params.interleaved_datasets,
                        'batch_size': params.train_batch_size,
                        'grad_accumulation': params.grad_accumulation,
                        'token_patch_size': params.token_patch_size
                        }

        size_dump = jsonpickle.dumps(_run_log + [curran_stats], indent=4)
        with tf.io.gfile.GFile(f"{params.model_path}/model_size.info", 'w') as f:
            f.write(size_dump)

        if len(_run_log) > 0 and not params.use_random_dataloader:
            _run_log = [r for r in _run_log if r['steps'] != params.current_step]
            if len(_run_log) > 0:
                run_log = [_run_log.pop(-1)]
                for r in _run_log[::-1]:
                    if run_log[-1]['steps'] != r['steps'] and r['steps'] != params.current_step:
                        run_log.append(r)
                run_log = run_log[::-1]

                for run_idx in range(len(run_log) - 1):
                    run_log[run_idx]['steps'] = run_log[run_idx + 1]['steps'] - run_log[run_idx]['steps']

                run_log[-1]['steps'] = params.current_step - run_log[-1]['steps']

                if run_log[-1]['steps'] <= 0:
                    run_log = None

    ds_iterator = []
    # For each sub-batch, create a SubBatchSlicer object.
    for sub_batch_i, host_id in enumerate(hosts_to_hold_ds):
        # Get the list of pnums for each input.
        if params.train:

            all_sub_batch_pnums = []
            for pnum_map in pnum_maps:
                sub_batch_pnums = pnum_map[sub_batch_i, ...].flatten().tolist()
                all_sub_batch_pnums.append(sub_batch_pnums)

        else:

            all_sub_batch_pnums = [pnum_map.flatten().tolist() for pnum_map in pnum_maps]

        with ops.device(f"/job:worker/task:{host_id}/device:CPU:0"):
            dataset = input_fn(params, sub_batch_size, sub_batch_i, len(hosts_to_hold_ds), run_log)
            if not params.use_random_dataloader and params.train and params.use_video:
                dataset = dataset.skip(params.current_step // params.macro_batching)
            dataset = dataset.prefetch(params.buffer_size)
            options = tf.data.Options()
            options.autotune.enabled = True
            options.deterministic = not params.train
            options.experimental_distribute.auto_shard_policy = tf.data.experimental.AutoShardPolicy.AUTO
            options.experimental_optimization.filter_fusion = True
            options.experimental_optimization.apply_default_optimizations = True
            options.experimental_optimization.map_and_batch_fusion = True
            options.experimental_optimization.map_and_filter_fusion = True
            options.experimental_optimization.map_fusion = True
            options.experimental_optimization.map_parallelization = True
            options.experimental_optimization.noop_elimination = True
            options.experimental_optimization.parallel_batch = True
            options.experimental_optimization.shuffle_and_repeat_fusion = True
            options.experimental_slack = True
            options.threading.private_threadpool_size = 96
            options.threading.max_intra_op_parallelism = 1
            dataset: Dataset = dataset.with_options(options)
            _ds_iterator = tf1.data.make_initializable_iterator(dataset)
            ds_iterator.append(_ds_iterator)
            all_input_tensors = _ds_iterator.get_next()

            if isinstance(all_input_tensors, tf.Tensor):
                all_input_tensors = [all_input_tensors]
            assert len(all_input_tensors) == len(all_sub_batch_pnums)

            for input_i in range(len(all_input_tensors)):
                input_tensor = all_input_tensors[input_i]
                sub_batch_pnums = all_sub_batch_pnums[input_i]
                mtf_input_shape = params.input_pipeline_shape[input_i]

                # Initialize the cache for each input_i
                _slice_dict = collections.defaultdict(list)

                for idx, pnum in enumerate(sub_batch_pnums):

                    s_begin = params.mesh_impl.slice_begin(mtf_input_shape, pnum)
                    if not not params.train:
                        # Always slice from 0 in the first dimension (batch dimension), since
                        # input_tensor a sub-batch tensor.
                        s_begin[0] = 0
                    if tuple(s_begin) in _slice_dict:
                        input_slice = _slice_dict[tuple(s_begin)]
                    else:
                        s_shape = params.mesh_impl.slice_shape(mtf_input_shape)
                        s_shape[0] = s_shape[0] * macro_batching_multi
                        input_slice = tfw.slice(input_tensor, s_begin, s_shape)

                    all_laidout_tensors[pnum][input_i] = input_slice

    # Make sure that there are no Nones in all_laidout_tensors.
    for laidout_tensors in all_laidout_tensors:
        assert None not in laidout_tensors

    with ops.device(f"/job:worker/task:{hosts_to_hold_ds[0]}/device:CPU:0"):

        def _tpu_ordinal_function_impl(pnum):
            return ordered_ordinals[pnum]

        def _placement_function_impl(pnum):
            return ordered_hosts[pnum]

        laidout_tensors0 = all_laidout_tensors[0]
        infeed_queue = tpu_feed.InfeedQueue(
            number_of_tuple_elements=len(laidout_tensors0),
            tuple_types=[x.dtype for x in laidout_tensors0],
            tuple_shapes=[x.shape for x in laidout_tensors0])
        enqueue_ops = infeed_queue.generate_enqueue_ops(
            all_laidout_tensors,
            tpu_ordinal_function=_tpu_ordinal_function_impl,
            placement_function=_placement_function_impl)

    input_initializers = [ds.initializer for ds in ds_iterator]

    return input_initializers, enqueue_ops, infeed_queue

def infeed_from_session(params: ModelParameter):
    num_cores = params.mesh_impl.device_assignment.num_replicas
    d_assignment = params.mesh_impl.device_assignment
    ordered_ordinals = []
    ordered_hosts = []

    for pnum in range(num_cores):
        physical_pnum = params.mesh_impl.l2p(pnum)
        host_device = d_assignment.host_device(replica=physical_pnum)
        ordered_hosts.append(host_device)

        # For MTF, there's always 1 core per replica. So logical_core=0.
        ordered_ordinals.append(d_assignment.tpu_ordinal(replica=physical_pnum, logical_core=0))

    def _tpu_ordinal_function_impl(pnum):
        return ordered_ordinals[pnum]

    def _placement_function_impl(pnum):
        return ordered_hosts[pnum]

    prompt = tf1.placeholder(dtype=tf.int32, shape=[t.size for t in params.token_dim_shape])
    iter_pos = tf1.placeholder(dtype=tf.int32, shape=[1])
    samp_temp = tf1.placeholder(dtype=tf.float32, shape=[1])
    end_iter = tf1.placeholder(dtype=tf.int32, shape=[1])

    all_laidout_tensors = [[prompt, prompt, iter_pos, samp_temp, end_iter] for _ in range(params.num_cores)]

    laidout_tensors0 = all_laidout_tensors[0]
    infeed_queue = tpu_feed.InfeedQueue(
        number_of_tuple_elements=len(laidout_tensors0),
        tuple_types=[x.dtype for x in laidout_tensors0],
        tuple_shapes=[x.shape for x in laidout_tensors0])
    enqueue_ops = infeed_queue.generate_enqueue_ops(all_laidout_tensors,
                                                    tpu_ordinal_function=_tpu_ordinal_function_impl,
                                                    placement_function=_placement_function_impl)

    place_holders = [prompt, iter_pos, samp_temp, end_iter]
    return enqueue_ops, infeed_queue, place_holders

max_stars_repo_name: Ziems/OBST

max_issues_repo_name: Ziems/OBST

max_forks_repo_name: Ziems/OBST