text:
import pickle
import codecs
import tables
import numpy as np
from tqdm import tqdm
import logging
logger = logging.getLogger(__name__)
logging.basicConfig(level=logging.INFO, format="%(asctime)s: %(name)s: %(levelname)s: %(message)s")

def load_pickle(filename):
    return pickle.load(open(filename, 'rb'))

##### Data Set #####
def load_codebase(path, chunk_size):
    """load codebase
    codefile: h5 file that stores raw code
    """
    logger.info('Loading codebase (chunk size={})..'.format(chunk_size))
    codebase=[]
    #codes=codecs.open(self.path+self.data_params['use_codebase']).readlines()
    codes=codecs.open(path, encoding='utf8',errors='replace').readlines()
        #use codecs to read in case of encoding problem
    for i in tqdm(range(0,len(codes), chunk_size)):
        codebase.append(codes[i:i+chunk_size])
    return codebase

### Results Data ###
def load_code_reprs(path, chunk_size):
    logger.debug(f'Loading code vectors (chunk size={chunk_size})..')
    """reads vectors (2D numpy array) from a hdf5 file"""
    codereprs=[]
    h5f = tables.open_file(path)
    vecs = h5f.root.vecs
    for i in range(0, len(vecs), chunk_size):
        codereprs.append(vecs[i: i+ chunk_size])
    h5f.close()
    return codereprs

def save_code_reprs(vecs, path):
    npvecs=np.array(vecs)
    fvec = tables.open_file(path, 'w')
    atom = tables.Atom.from_dtype(npvecs.dtype)
    filters = tables.Filters(complib='blosc', complevel=5)
    ds = fvec.create_carray(fvec.root, 'vecs', atom, npvecs.shape,filters=filters)
    ds[:] = npvecs
    fvec.close()

def load_hdf5(vecfile, start_offset, chunk_size):
    """reads training sentences(list of int array) from a hdf5 file"""
    table = tables.open_file(vecfile)
    data = table.get_node('/phrases')[:].astype(np.int)
    index = table.get_node('/indices')[:]
    data_len = index.shape[0]
    if chunk_size==-1:#if chunk_size is set to -1, then, load all data
        chunk_size=data_len
    start_offset = start_offset%data_len
    logger.debug("{} entries".format(data_len))
    logger.debug("starting from offset {} to {}".format(start_offset,start_offset+chunk_size))
    sents = []
    for offset in tqdm(range(start_offset, start_offset+chunk_size)):
        offset = offset%data_len
        len, pos = index[offset]['length'], index[offset]['pos']
        sents.append(data[pos:pos + len])
    table.close()
    return sents


text:
"""
Library Features:

Name:          lib_jupyter_data_io_ascii
Author(s):     Fabio Delogu (fabio.delogu@cimafoundation.org)
Date:          '20210113'
Version:       '1.0.0'
"""
#######################################################################################
# Library
import logging
import os
import re

import xarray as xr
import pandas as pd
import numpy as np

from library.jupyter_generic.lib_jupyter_utils_system import get_dict_values, fill_tags2string, make_folder
from library.jupyter_generic.lib_jupyter_info_args import logger_name

# Logging
log_stream = logging.getLogger(logger_name)
#######################################################################################

# -------------------------------------------------------------------------------------
# Method to read ascii hydrograph file
def read_file_ascii_hydrograph(file_name, file_header=None, file_columns=None,
                               file_terminator='\n', file_separator='\t', file_names=None):

    if os.path.exists(file_name):

        if file_names is None:
            file_names = ['data']

        # Read ts
        file_data = pd.read_table(file_name, header=file_header, lineterminator=file_terminator, names=['data'])

        file_time_ws = []
        file_data_ws = []
        for id, ts_obj in enumerate(file_data['data']):

            ts_raw = ts_obj.split(file_separator)
            ts_parts = ts_raw[0].split()

            ts_time = pd.Timestamp(ts_parts[0])
            ts_values = ts_parts[1:]
            ts_values = [float(i) for i in ts_values]

            file_time_ws.append(ts_time)
            file_data_ws.append(ts_values)

        file_time_idx = pd.DatetimeIndex(file_time_ws)
        file_data_n = file_data_ws[0].__len__()

        if file_columns is None:
            file_columns = ["section_{:}".format(idx) for idx in range(1, file_data_n + 1)]

        file_dframe = pd.DataFrame(file_data_ws, columns=file_columns, index=file_time_idx)

        return file_dframe

    else:
        logging.warning(' ===> File ' + file_name + ' not found!')

# -------------------------------------------------------------------------------------

ext: py

max_stars_repo_name: c-hydro/fp-labs-hmc

max_issues_repo_name: c-hydro/fp-labs-hmc

max_forks_repo_name: c-hydro/fp-labs-hmc

reason: import numpy


text:
import pandas as pd
import json
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
from pycocotools.coco import  COCO
import numpy as np
import skimage.io as io
import cv2
import os
import random
from tensorflow.keras.utils import to_categorical

def read_annotation_json(train_json_path = None):
    # make a file object for reading json path
    file_obj = open(train_json_path)

    # load the json data
    json_data = json.load(file_obj)

    # get all the keys
    keys = json_data.keys()

    usable_df_list = []

    for key in keys:
        temp_dict ={}
        try:
            df = pd.DataFrame(json_data[key])
            temp_dict['key']=key
            temp_dict['df'] = df

            print('-----------------------------------------')
            print(key)
            print(df.columns)
            print(df.head(5))

            usable_df_list.append(temp_dict)
        except:
            continue

    # close the file object
    file_obj.close()

    return usable_df_list

class DataManipulationCoco:
    def __init__(self, annotation_file=None, batch_size=None, training=False):
        self.coco = COCO(annotation_file)
        self.batch_size = batch_size
        self.num_class = self.get_number_of_class()

        if training:
            self.img_folder = r'../../data/trainingset/'
        else:
            self.img_folder=r'../../data/testset/'

    def get_categories(self):
        category_ids = self.coco.getCatIds()
        self.categories = self.coco.loadCats(category_ids)
        return self.categories

    def get_number_of_class(self):
        classes = [category['name'] for category in self.get_categories()]
        return len(classes)

    def get_class_for_category_id(self, category_id):
        categories = self.get_categories()
        for i in range(categories):
            if categories[i]['id'] == category_id:
                return categories[i]['name']
        return None

    def get_all_combination_subset(self, filter_classes=None):
        images = []
        if filter_classes!=None:
            for f_class in filter_classes:
                category_id = self.coco.getCatIds(catNms=f_class)
                image_ids = self.coco.getImgIds(catIds=category_id)
                images += self.coco.loadImgs(image_ids)
        else:
            image_ids = self.coco.getImgIds()
            images = self.coco.loadImgs(image_ids)

        unique_images = []
        for i in range(len(images)):
            if images[i] not in unique_images:
                unique_images.append(images[i])
        print("Number of unique images ", len(unique_images))

        random.shuffle(unique_images)

        return unique_images

    def getClassName(self, classID, cats):
        for i in range(len(cats)):
            if cats[i]['id'] == classID:
                return cats[i]['name']
        return None

    def getImage(self, imageObj, img_folder, input_image_size):
        # Read and normalize an image
        train_img = io.imread(img_folder + '/' + imageObj['file_name']) / 255.0
        # Resize
        train_img = cv2.resize(train_img, input_image_size)
        if (len(train_img.shape) == 3 and train_img.shape[2] == 3):  # If it is a RGB 3 channel image
            return train_img
        else:  # To handle a black and white image, increase dimensions to 3
            stacked_img = np.stack((train_img,) * 3, axis=-1)
            return stacked_img

    def getNormalMask(self, imageObj, classes,  catIds, input_image_size):
        annIds = self.coco.getAnnIds(imageObj['id'], catIds=catIds, iscrowd=None)
        anns = self.coco.loadAnns(annIds)
        cats = self.coco.loadCats(catIds)
        train_mask = np.zeros(input_image_size)
        classes = [category['name'] for category in self.get_categories()]
        for a in range(len(anns)):
            className = self.getClassName(anns[a]['category_id'], cats)
            pixel_value = classes.index(className)
            new_mask = cv2.resize(self.coco.annToMask(anns[a]) * pixel_value, input_image_size)
            train_mask = np.maximum(new_mask, train_mask)

        # Add extra dimension for parity with train_img size [X * X * 3]
        train_mask = train_mask.reshape(input_image_size[0], input_image_size[1], 1)
        return train_mask

    def getBinaryMask(self, imageObj, catIds, input_image_size):
        annIds = self.coco.getAnnIds(imageObj['id'], catIds=catIds, iscrowd=None)
        anns = self.coco.loadAnns(annIds)
        train_mask = np.zeros(input_image_size)
        for a in range(len(anns)):
            new_mask = cv2.resize(self.coco.annToMask(anns[a]), input_image_size)

            # Threshold because resizing may cause extraneous values
            new_mask[new_mask >= 0.5] = 1
            new_mask[new_mask < 0.5] = 0

            train_mask = np.maximum(new_mask, train_mask)

        # Add extra dimension for parity with train_img size [X * X * 3]
        train_mask = train_mask.reshape(input_image_size[0], input_image_size[1], 1)
        return train_mask

    def dataGeneratorCoco(self, images, classes=None,
                          input_image_size=(224, 224), batch_size=None,  mask_type='normal'):

        img_folder = self.img_folder
        dataset_size = len(images)
        if classes is None:
            classes = [category['name'] for category in self.get_categories()]
        catIds = self.coco.getCatIds(catNms=classes)
        batch_size=self.batch_size

        c = 0
        while (True):
            img = np.zeros((batch_size, input_image_size[0], input_image_size[1], 3)).astype('float')
            mask = np.zeros((batch_size, input_image_size[0], input_image_size[1], 1)).astype('float')

            for i in range(c, c + batch_size):  # initially from 0 to batch_size, when c = 0
                imageObj = images[i]

                ### Retrieve Image ###
                train_img = self.getImage(imageObj, img_folder, input_image_size)

                ### Create Mask ###
                if mask_type == "binary":
                    train_mask = self.getBinaryMask(imageObj, catIds, input_image_size)

                elif mask_type == "normal":
                    train_mask = self.getNormalMask(imageObj, classes, catIds, input_image_size)

                    # Add to respective batch sized arrays
                img[i - c] = train_img
                mask[i - c] = train_mask

            c += batch_size
            if (c + batch_size >= dataset_size):
                c = 0
                random.shuffle(images)
            mask_categorical = to_categorical(mask, num_classes=self.num_class)
            yield img, mask_categorical

    def get_image_subset(self, filter_class=None):
        # Get all category class in case no class has been defined
        if filter_class==None:
            filter_class = [category['name'] for category in self.get_categories()]

        # get all the category ids
        category_ids = self.coco.getCatIds(catNms=filter_class)
        print(category_ids)

        self.get_all_combination_subset(filter_classes=filter_class)

        # get all images containing the above category ids
        image_ids = self.coco.getImgIds(catIds=category_ids)
        print(image_ids)

        # print number of images containing all classes
        print("Number of images containing all classes ", len(image_ids))

    def visualizeGenerator(self, gen):
        img, mask = next(gen)

        mask = np.argmax(mask, axis=3)

        fig = plt.figure(figsize=(20, 10))
        outerGrid = gridspec.GridSpec(1, 2, wspace=0.1, hspace=0.1)

        for i in range(2):
            innerGrid = gridspec.GridSpecFromSubplotSpec(2, 2,
                                                         subplot_spec=outerGrid[i], wspace=0.05, hspace=0.05)

            for j in range(4):
                ax = plt.Subplot(fig, innerGrid[j])
                if (i == 1):
                    ax.imshow(img[j])
                else:
                    ax.imshow(mask[j][:, :, 0])

                ax.axis('off')
                fig.add_subplot(ax)
        plt.savefig('graph.png')

if __name__ == '__main__':
    # usable_data_list = read_annotation_json()
    # for dict_obj in usable_data_list:
    #     if dict_obj['key']=='annotations':
    #         print(dict_obj['df']['category_id'].value_counts().plot.bar())
    #         print(dict_obj['df']['segmentation'])
    # plt.show()
    data_manipulation = DataManipulationCoco()
    images = data_manipulation.get_all_combination_subset()
    train_gen = data_manipulation.dataGeneratorCoco(images, None)
    data_manipulation.visualizeGenerator(train_gen)


text:
import numpy as np
import math
import cv2
from skimage.color import rgb2hsv,hsv2rgb

e = np.e
esp = 2.2204e-16

# 用于排序时存储原来像素点位置的数据结构
class Node(object):
	def __init__(self,x,y,value):
		self.x = x
		self.y = y
		self.value = value
	def printInfo(self):
		print(self.x,self.y,self.value)

def rayleighStrUpper(nodes, height, width,lower_Position):
    allSize = height * width
    alpha = 0.5
    selectedRange = [0, 255]
    NumPixel = np.zeros(256)
    temp = np.zeros(256)
    for i in range(lower_Position, allSize):
            NumPixel[nodes[i].value] = NumPixel[nodes[i].value] + 1
    ProbPixel = NumPixel / (allSize-lower_Position)
    CumuPixel = np.cumsum(ProbPixel)
    # print('UpperCumuPixel', CumuPixel)
    valSpread = selectedRange[1] - selectedRange[0]
    hconst = 2 * alpha ** 2
    vmax = 1 - e ** (-1 / hconst)
    val = vmax * (CumuPixel)
    val = np.array(val)

    for i in range(256):
        if (val[i] >= 1):
            val[i] = val[i] - esp
    for i in range(256):
        temp[i] = np.sqrt(-hconst * math.log((1 - val[i]), e))
        normalization = temp[i] * valSpread
        if(normalization > 255):
            CumuPixel[i] = 250
        else:
            CumuPixel[i] = normalization
    # print('UpperCumuPixel',CumuPixel)
    for i in range(lower_Position, allSize):
        nodes[i].value = CumuPixel[nodes[i].value]
    return nodes

def uperLower(r, height, width):
    allSize = height * width
    R_max = np.max(r)
    R_min = np.min(r)
    R__middle = (R_max - R_min) / 2 + R_min
    # R__middle = np.mean(r)
    nodes = []
    for i in range(0, height):
        for j in range(0, width):
            oneNode = Node(i, j, r[i, j])
            nodes.append(oneNode)
    nodes = sorted(nodes, key=lambda node: node.value, reverse=False)
    # print('R__middle',R__middle)

    for i in range(240000):
        if (nodes[i].value > R__middle):
            # print('nodes[i].value',nodes[i].value)
            middle_Position = i
            break
    lower_Position = middle_Position
    # print('lower_Position',lower_Position)

    for i in range(240000):
        nodes[i].value = np.int(nodes[i].value)
    # print('nodes[0].value)', nodes[0].value)
    # print('nodes[lower_Position + 2].value)', nodes[lower_Position + 2].value)

    nodesUpper  = rayleighStrUpper(nodes, height, width,lower_Position)

    array_upper_histogram_stretching = np.zeros((height, width))

    # print('lower_Position', lower_Position)
    # print('allSize', allSize)
    for i in range(0, allSize):
        if(i > lower_Position):
            array_upper_histogram_stretching[nodes[i].x, nodes[i].y] = nodesUpper[i].value
        else:
            array_upper_histogram_stretching[nodes[i].x, nodes[i].y] = 5

    # print('np.mean(array_upper_histogram_stretching))',np.mean(array_upper_histogram_stretching))

    return array_upper_histogram_stretching

def rayleighStretching_Upper(sceneRadiance, height, width):
    img_hsv = rgb2hsv(sceneRadiance)
    h, s, v = cv2.split(img_hsv)

    s = s * 255
    v = v * 255

    v_array_upper_histogram_stretching = uperLower(v, height, width)/255
    s_array_upper_histogram_stretching = uperLower(s, height, width)/255
    h_array_upper_histogram_stretching = h
    # print('R_array_upper_histogram_stretching',R_array_upper_histogram_stretching)

    sceneRadiance_Upper = np.zeros((height, width, 3))
    sceneRadiance_Upper[:, :, 0] = h_array_upper_histogram_stretching
    sceneRadiance_Upper[:, :, 1] = s_array_upper_histogram_stretching
    sceneRadiance_Upper[:, :, 2] = v_array_upper_histogram_stretching
    img_rgb = hsv2rgb(sceneRadiance_Upper) * 255
    for i in range(0, 3):
        for j in range(0, img_rgb.shape[0]):
            for k in range(0, img_rgb.shape[1]):
                if img_rgb[j, k, i] > 250:
                    img_rgb[j, k, i] = 250
                if img_rgb[j, k, i] < 5:
                    img_rgb[j, k, i] = 5

    sceneRadiance_Upper = np.uint8(img_rgb)

    return sceneRadiance_Upper


text:
import gym
import gym.envs.toy_text.frozen_lake as fl
import numpy as np
import pdb

SEED = 12345

np.set_printoptions(precision=3)

# Create environment
env = fl.FrozenLakeEnv()  # gym.make('FrozenLake-v1')
env.seed(SEED)

# It should print a rendering of the FrozenLake environment
env.render()
pdb.set_trace(header="\n*** *** ***\nHello, you entered the interactive source code debugger.\nHere you can "
                     "experiment with elements "
                     "that have been defined. For instance you can call again env.render().\nTo continue, just type "
                     "continue.\nTo quit the program, just type quit.\n*** *** ***\n")

# Of course, you should RTFM (please do a quick peek into the doc)
help(fl)
# ... and again
help(gym.envs.toy_text.discrete)
pdb.set_trace(header="\nTake some time to read the doc.\n")

# Some printing to understand the 'env' class
print("Action Space - ", env.action_space, env.nA)
print("State Space - ", env.observation_space, env.nS)
# isd is the initial state distribution
print("\nInitial State Distribution  - ", env.isd, np.sum(env.isd) == 1)
print("\nTransitions from STATE 1 when going LEFT\n", env.P[1][0])
pdb.set_trace(header="\nTake some time to understand the few previous lines.\nDon't forget that you can play a little "
                     "with the variables (but be careful to not modify them).\n")

# If you understand those two functions
# you understand how states are linked to env.render()
def to_s(row, col):
    return row * env.ncol + col

def to_row_col(s):
    col = s % env.ncol
    row = s // env.ncol
    return row, col

# Documentation for interacting with a gym environment
help(gym.Env)
pdb.set_trace(header="\nRemember to read the doc if you don't understand/recall something.\n")


text:
import copy
import json
import unittest
import numpy as np
from numpy import testing
try:
    import tinyfk
except:
    import _tinyfk as tinyfk
from math import *

urdf_model_path = tinyfk.pr2_urdfpath()
test_data_path = tinyfk._test_data_urdfpath()

with open(test_data_path, 'r') as f:
    test_data = json.load(f)
joint_names = test_data['joint_names']
link_names = test_data['link_names']
angle_vector = test_data['angle_vector']
gt_pose_list_ = np.array(test_data['pose_list'])
gt_pose_list = copy.copy(gt_pose_list_)
gt_pose_list[:, 3] = gt_pose_list_[:, 5] # note in skrobot rpy order is z-y-x
gt_pose_list[:, 5] = gt_pose_list_[:, 3] # so, the two lines are swapped.

# common setups
fksolver = tinyfk.RobotModel(urdf_model_path)

# adding new link `mylink` to `r_upper_arm_link`
parent_id = fksolver.get_link_ids(["r_upper_arm_link"])[0]
fksolver.add_new_link('mylink', parent_id, [0.1, 0.1, 0.1])

# To interact with fksolver, we must get the correspoinding link_ids and joint_ids.
link_ids = fksolver.get_link_ids(link_names)
joint_ids = fksolver.get_joint_ids(joint_names)

def test_fksovler():
    # check P (array of poses [pos, rpy] of each link) coincides with the ground truth
    use_rotation = True # If true P[i, :] has 6 dim, otherwise has 3 dim.
    use_base = True # If true, assumes that angle_vector takes the form of [q_joints, q_base (3dof)]
    with_jacobian = True # If true, jacobian is computed, otherewise returns J = None.
    P, _ = fksolver.solve_forward_kinematics(
            [angle_vector], link_ids, joint_ids, use_rotation, use_base, with_jacobian)
    testing.assert_almost_equal(P, gt_pose_list)

    # The following test assumes: that the above test without jacobian succeeded.
    # check resulting jacbian J_analytical coincides with J_numerical witch is
    # computed via numerical differentiation.
    for link_id, link_name in zip(link_ids, link_names):
        P_tmp, J_analytical = fksolver.solve_forward_kinematics(
                [angle_vector], [link_id], joint_ids, True, True, True)
        P0, _ = fksolver.solve_forward_kinematics(
                [angle_vector], [link_id], joint_ids, True, True, False)
        testing.assert_almost_equal(P_tmp, P0) # P computed with and without jacobian must match

        eps = 1e-7
        J_numerical = np.zeros(J_analytical.shape)
        for i in range(len(angle_vector)):
            angle_vector_p = copy.copy(angle_vector)
            angle_vector_p[i] += eps
            P1, _ = fksolver.solve_forward_kinematics(
                    [angle_vector_p], [link_id], joint_ids, True, True, False)
            P_diff = (P1 - P0)/eps
            J_numerical[:, i] = P_diff.flatten()

        rpy = P0[0][3:]
        # test position jacobian
        testing.assert_almost_equal(J_numerical[:3, :], J_analytical[:3, :])
        # test rpy jacobian
        testing.assert_almost_equal(J_numerical[3:, :], J_analytical[3:, :])

max_stars_repo_name: HiroIshida/tinyfk

max_issues_repo_name: HiroIshida/tinyfk

max_forks_repo_name: HiroIshida/tinyfk

reason: import numpy,from numpy


text:
# Copyright (c) 2016, The Bifrost Authors. All rights reserved.
# Copyright (c) 2016, NVIDIA CORPORATION. All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions
# are met:
# * Redistributions of source code must retain the above copyright
#   notice, this list of conditions and the following disclaimer.
# * Redistributions in binary form must reproduce the above copyright
#   notice, this list of conditions and the following disclaimer in the
#   documentation and/or other materials provided with the distribution.
# * Neither the name of The Bifrost Authors nor the names of its
#   contributors may be used to endorse or promote products derived
#   from this software without specific prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS ``AS IS'' AND ANY
# EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
# PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR
# CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
# EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
# PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
# PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY
# OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
# (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

import ctypes
import numpy as np
from bifrost.memory import raw_malloc, raw_free, memset, memcpy, memcpy2D
from bifrost.libbifrost import _bf, _check, _string2space
from bifrost.array import _array2bifrost # This doesn't exist!

class GPUArray(object):
    def __init__(self, shape, dtype, buffer=None, offset=0, strides=None):
        itemsize = dtype().itemsize
        shape = tuple(np.array(shape).ravel().astype(np.uint64))
        if strides is None:
            # This magic came from http://stackoverflow.com/a/32874295
            strides = itemsize * np.r_[1, np.cumprod(shape[::-1][:-1],
                                                     dtype=np.int64)][::-1]
        self.shape   = shape
        self.dtype   = dtype
        self.buffer  = buffer
        self.offset  = offset
        self.strides = strides
        self.base    = None
        self.flags   = {'WRITEABLE':    True,
                        'ALIGNED':      buffer % (itemsize == 0
                                                  if buffer is not None
                                                  else True),
                        'OWNDATA':      False,
                        'UPDATEIFCOPY': False,
                        'C_CONTIGUOUS': self.nbytes == strides[0] * shape[0],
                        'F_CONTIGUOUS': False,
                        'SPACE':        'cuda'}
        class CTypes(object):
            def __init__(self, parent):
                self.parent = parent
            @property
            def data(self):
                return self.parent.data
        self.ctypes = CTypes(self)
        if self.buffer is None:
            self.buffer = raw_malloc(self.nbytes, space='cuda')
            self.flags['OWNDATA'] = True
            self.flags['ALIGNED'] = True
            memset(self, 0)
        else:
            self.buffer += offset
    def __del__(self):
        if self.flags['OWNDATA']:
            raw_free(self.buffer, self.flags['SPACE'])
    @property
    def data(self):
        return self.buffer
    # def reshape(self, shape):
    #    # TODO: How to deal with strides?
    #    #         May be non-contiguous but the reshape still works
    #    #           E.g., splitting dims
    #    return GPUArray(shape, self.dtype,
    #                    buffer=self.buffer,
    #                    offset=self.offset,
    #                    strides=self.strides)
    @property
    def size(self):
        return int(np.prod(self.shape))
    @property
    def itemsize(self):
        return self.dtype().itemsize
    @property
    def nbytes(self):
        return self.size * self.itemsize
    @property
    def ndim(self):
        return len(self.shape)
    def get(self, dst=None):
        hdata = dst if dst is not None else np.empty(self.shape, self.dtype)
        # hdata = dst if dst is not None else np.zeros(self.shape, self.dtype)
        assert(hdata.shape == self.shape)
        assert(hdata.dtype == self.dtype)
        if self.flags['C_CONTIGUOUS'] and hdata.flags['C_CONTIGUOUS']:
            memcpy(hdata, self)
        elif self.ndim == 2:
            memcpy2D(hdata, self)
        else:
            raise RuntimeError("Copying with this data layout is unsupported")
        return hdata
    def set(self, hdata):
        assert(hdata.shape == self.shape)
        hdata = hdata.astype(self.dtype)
        if self.flags['C_CONTIGUOUS'] and hdata.flags['C_CONTIGUOUS']:
            memcpy(self, hdata)
        elif self.ndim == 2:
            memcpy2D(self, hdata)
        else:
            raise RuntimeError("Copying with this data layout is unsupported")
        return self
    def as_BFarray(self):
        return _array2bifrost(self)

max_stars_repo_name: ajameson/bifrost

max_issues_repo_name: ajameson/bifrost

max_forks_repo_name: ajameson/bifrost


text:
from __future__ import print_function

import numpy as np

from keras.layers import Layer, Input, Dense, Lambda, Flatten, Reshape, Conv2D, Conv2DTranspose, Subtract
from keras.models import Model
from keras import backend as K
from keras.losses import mse, binary_crossentropy

def sampling(args):
    """Reparameterization trick by sampling from an isotropic unit Gaussian.

    # Arguments:
        args (tensor): mean and log of variance of Q(z|X)

    # Returns:
        z (tensor): sampled latent vector
    """
    z_mean, z_log_var, scale_var = args
    batch = K.shape(z_mean)[0]
    dim = K.int_shape(z_mean)[1]
    epsilon = K.random_normal(shape=(batch, dim))
    return z_mean + K.exp(0.5 * z_log_var) * epsilon * scale_var

def build_vae_encoder_model(filters      = 32,
                            latent_dims  = 2,
                            z_activation = 'relu'):
    input_shape = (28, 28, 1)
    input_img = Input(shape=input_shape) # the input image
    input_scale_var = Input(shape=(1,))        # used to scale z_log_var
    x         = Conv2D(filters, (3,3), activation='relu', strides=(2,2), padding='same', name='conv0_encode')(input_img)
    x         = Conv2D(filters*2, (3,3), activation='relu', strides=(2,2), padding='same', name='conv1_encode')(x)
    x         = Flatten(name='flatten0_encode')(x)
    x         = Dense(32, activation='relu', name='dense0_encode')(x)
    z_mean    = Dense(latent_dims, name='z_mean_encode')(x)
    z_log_var = Dense(latent_dims, name='z_log_var_encode')(x)
    z         = Lambda(sampling, output_shape=(latent_dims,), name='z_encode')([z_mean, z_log_var, input_scale_var])
    encoder   = Model([input_img, input_scale_var], [z, z_mean, z_log_var], name='encoder')
    encoder.summary()
    return encoder

def build_vae_decoder_model(filters     = 32,
                            latent_dims = 2,
                            output_cnls = 1):
    z_input = Input(shape=(latent_dims,))
    x = Dense(7*7*filters*2, activation='relu')(z_input)
    x = Reshape((7, 7, filters*2))(x)
    x = Conv2DTranspose(filters*2, (3,3), activation='relu', strides=(2,2), padding='same', name='conv0_decode')(x)
    x = Conv2DTranspose(filters, (3,3), activation='relu', strides=(2,2), padding='same', name='conv1_decode')(x)
    decoded = Conv2DTranspose(1, (3,3), activation='sigmoid', padding='same', name='output')(x)
    decoder = Model(z_input, decoded, name='decoder')
    decoder.summary()
    return decoder

def build_vae_conv_model():
    """ build a convolutional variational autoencoder. """
    inputs_img   = Input(shape=(28, 28, 1))
    inputs_scale_var = Input(shape=(1,))
    encoder  = build_vae_encoder_model()
    decoder  = build_vae_decoder_model()
    z, z_mean, z_log_var  = encoder([inputs_img, inputs_scale_var])
    outputs  = decoder(z)
    vae = Model([inputs_img, inputs_scale_var], [outputs], name='vae')
    # define custom loss
    reconstruction_loss = mse(K.flatten(inputs_img), K.flatten(outputs))
    reconstruction_loss *= (28*28)
    kl_loss = 1 + z_log_var - K.square(z_mean) - K.exp(z_log_var)
    kl_loss = K.sum(kl_loss, axis=-1)
    kl_loss *= -0.5
    vae_loss = K.mean(reconstruction_loss + kl_loss)
    vae.add_loss(vae_loss)
    vae.compile(optimizer='adam')
    vae.summary()
    return vae, encoder, decoder

max_stars_repo_name: RyanDellana/autoencoder-attractors

max_issues_repo_name: RyanDellana/autoencoder-attractors

max_forks_repo_name: RyanDellana/autoencoder-attractors


text:
# -----------------------------------------------------------------------------
# From Numpy to Python
# Copyright (2017) Nicolas P. Rougier - BSD license
# More information at https://github.com/rougier/numpy-book
# -----------------------------------------------------------------------------

import numpy as np
import numba as nb

@nb.jit(nopython=True, parallel=False, fastmath=True)
def linspace(start, stop, num, dtype):
    X = np.empty((num, ), dtype=dtype)
    dist = (stop - start) / (num - 1)
    for i in range(num):
        X[i] = start + i * dist
    return X

@nb.jit(nopython=True, parallel=False, fastmath=True)
def mandelbrot(xmin, xmax, ymin, ymax, xn, yn, maxiter, horizon=2.0):
    # Adapted from https://www.ibm.com/developerworks/community/blogs/jfp/...
    #              .../entry/How_To_Compute_Mandelbrodt_Set_Quickly?lang=en
    # X = np.linspace(xmin, xmax, xn, dtype=np.float64)
    # Y = np.linspace(ymin, ymax, yn, dtype=np.float64)
    X = linspace(xmin, xmax, xn, dtype=np.float64)
    Y = linspace(ymin, ymax, yn, dtype=np.float64)
    # C = X + Y[:,None]*1j
    C = X + np.reshape(Y, (yn, 1)) * 1j
    N = np.zeros(C.shape, dtype=np.int64)
    Z = np.zeros(C.shape, dtype=np.complex128)
    for n in range(maxiter):
        # I = np.less(abs(Z), horizon)
        I = np.less(np.absolute(Z), horizon)
        # N[I] = n
        for j in range(C.shape[0]):
            for k in range(C.shape[1]):
                if I[j, k]:
                    N[j, k] = n
        # Z[I] = Z[I]**2 + C[I]
        for j in range(C.shape[0]):
            for k in range(C.shape[1]):
                if I[j, k]:
                    Z[j, k] = Z[j, k]**2 + C[j, k]
    # N[N == maxiter-1] = 0
    for j in range(C.shape[0]):
        for k in range(C.shape[1]):
            if N[j, k] == maxiter - 1:
                N[j, k] = 0
    return Z, N

max_stars_repo_name: frahlg/npbench

max_issues_repo_name: frahlg/npbench

max_forks_repo_name: frahlg/npbench

reason: import numpy,import numba


text:
from __future__ import absolute_import
from __future__ import division
from __future__ import print_function

import argparse
from functools import partial
import json
import traceback

import imlib as im
import numpy as np
import pylib
import tensorflow as tf
import tflib as tl

import data
import models

import os

# ==============================================================================
# =                                    param                                   =
# ==============================================================================

parser = argparse.ArgumentParser()
parser.add_argument('--experiment_name', help='experiment_name')
parser.add_argument('--gpu', type=str, default='all', help='gpu')
parser.add_argument('--dataroot', type=str, default='/data/Datasets/CelebA/Img')
# if assigned, only given images will be tested.
parser.add_argument('--img', type=int, nargs='+', default=None, help='e.g., --img 182638 202599')
# for multiple attributes
parser.add_argument('--test_att', type=str, default='Bangs')
args_ = parser.parse_args()
with open('./output/%s/setting.txt' % args_.experiment_name) as f:
    args = json.load(f)

# model
atts = args['atts']
img_size = args['img_size']

dataroot = args_.dataroot
img = args_.img
print('Using selected images:', img)

gpu = args_.gpu
if gpu != 'all':
    os.environ['CUDA_VISIBLE_DEVICES'] = gpu

# others
use_cropped_img = args['use_cropped_img']
experiment_name = args_.experiment_name

# ==============================================================================
# =                                   graphs                                   =
# ==============================================================================

# data
sess = tl.session()
te_data = data.Celeba(dataroot, atts, img_size, 1, part='test', sess=sess, crop=not use_cropped_img, im_no=img)
# models
G = partial(models.G)
D = partial(models.D)

# inputs
xa_sample = tf.placeholder(tf.float32, shape=[1, img_size, img_size, 3])
_b_sample = tf.placeholder(tf.float32, shape=[1, ])
raw_b_sample = tf.placeholder(tf.float32, shape=[1, ])

# sample
test_label = _b_sample
x_sample, mask_sample = G(xa_sample,  test_label, is_training=False)

# ==============================================================================
# =                                    test                                    =
# ==============================================================================

# initialization
ckpt_dir = './output/%s/checkpoints' % experiment_name
tl.load_checkpoint(ckpt_dir, sess)

# test
try:
    for idx, batch in enumerate(te_data):
        xa_sample_ipt = batch[0]
        a_sample_ipt = batch[1]
        b_sample_ipt_list = [a_sample_ipt.copy() for _ in range(1)]
        # test_single_attributes
        tmp = np.array(a_sample_ipt, copy=True)
        tmp[:] = 1 - tmp[:]   # inverse attribute
        b_sample_ipt_list.append(tmp)

        x_sample_opt_list = [xa_sample_ipt, np.full((1, img_size, img_size // 10, 3), -1.0)]
        mask_sample_opt_list = [np.full((1, img_size, img_size // 10, 1), -1.0)]
        for i, b_sample_ipt in enumerate(b_sample_ipt_list):
            x_sample_opt_list.append(sess.run(x_sample, feed_dict={xa_sample: xa_sample_ipt,
                                                                   _b_sample: b_sample_ipt,}))
            mask_sample_opt_list.append(sess.run(mask_sample, feed_dict={xa_sample: xa_sample_ipt,
                                                                         _b_sample: b_sample_ipt}))
        sample = np.concatenate(x_sample_opt_list, 2)
        masks = np.concatenate(mask_sample_opt_list, 2)

        save_folder = 'sample_testing'
        save_dir = './output/%s/%s' % (experiment_name, save_folder)
        pylib.mkdir(save_dir)
        im.imwrite(sample.squeeze(0), '%s/%06d%s.png' % (save_dir,
                                                         idx + 182638 if img is None else img[idx],
                                                         '_%s' % ''))
        im.imwrite(masks.squeeze(0), '%s/Mask%06d%s.png' % (save_dir,
                                                            idx + 182638 if img is None else img[idx],
                                                            '_%s' % ''))

        print('%06d.png done!' % (idx + 182638 if img is None else img[idx]))
except:
    traceback.print_exc()
finally:
    sess.close()

max_stars_repo_name: DateBro/SaGAN-TensorFlow

max_issues_repo_name: DateBro/SaGAN-TensorFlow

max_forks_repo_name: DateBro/SaGAN-TensorFlow


text:
"""
Copyright (c) 2017, Florian Engels
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
    * Redistributions of source code must retain the above copyright
      notice, this list of conditions and the following disclaimer.
    * Redistributions in binary form must reproduce the above copyright
      notice, this list of conditions and the following disclaimer in the
      documentation and/or other materials provided with the distribution.
    * Neither the name "Framework on High-Resolution Frequency Estimation"
      nor the names of its contributors may be used to endorse or promote
      products derived from this software without specific prior written
      permission.

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

class FrameworkOnHighResolutionFrequencyEstimation:
    """
    A framework on computationally efficient, high-resolution frequency estimation

    Implementation of the framework proposed in

    [1] F. Engels, P. Heidenreich, A. Zoubir, F. Jondral, and M. Wintermantel,
        "Advances in Automotive Radar: A framework on computationally efficient
         high-resolution frequency estimation",
        IEEE Signal Processing Magazine, Vol. 34, No. 2, Pp. 36-46.

    extended to process a fourth frequency dimension corresponding to elevation angle.

    Instance Variables

    w_lambda,          -- window functions in range, velocity, and angular (azimuth/elevation) dimension, respectively
    w_mu,
    w_nu,
    w_omega
    Lf,Mf,Nf,Of        -- support of Fourier transform in range, velocity, and angular (azimuth/elevation) dimension, respectively
    L,M,N,O            -- support of peak neighborhood in range, velocity, and angular (azimuth/elevation) dimension, respectively
    T_db               -- threshold for multiple target indication based on the single-target fitting error
    gamma_db           -- threshold for the generalized likelihood ratio test

    Ls,Ms,Ns,Os        -- sampling support in range, velocity, and angular (azimuth/elevation) dimension, respectively

    W_lambda,          -- model matrix in range, velocity, and angular (azimuth/elevation) dimension, respectively
    W_mu,
    W_nu,
    W_omega

    dfe                -- decoupled frequency estimation comprising K-target model fit in the resolution
                          dimension and subsequent single-target model fits in the remaining dimensions

    """
    def __init__(self,
                 w_lambda,w_mu,w_nu,
                 Lf,Mf,Nf,L,M,N,
                 T_db,gamma_db,
                 w_omega=None,Of=1,O=1,
                 average=True):

        self.w_lambda,self.w_mu,self.w_nu = w_lambda,w_mu,w_nu
        self.Lf,self.Mf,self.Nf = Lf,Mf,Nf
        self.L,self.M,self.N = L,M,N
        self.T_db = T_db
        self.gamma_db = gamma_db
        self.average = average

        self.Ls,self.Ms,self.Ns = w_lambda.shape[0],w_mu.shape[0],w_nu.shape[0]
        self.L_nls,self.M_nls,self.N_nls = 8,8,8
        self.grid_lambda = arange(self.L)*2*pi/Lf
        self.grid_mu = arange(self.M)*2*pi/Mf
        self.grid_nu = arange(self.N)*2*pi/Nf
        self.grid_lambda_nls = linspace(0,L*2*pi/Lf,self.L_nls)
        self.grid_mu_nls = linspace(0,M*2*pi/Mf,self.M_nls)
        self.grid_nu_nls = linspace(0,self.N*2*pi/Nf,self.N_nls)

        self.W_lambda = ModelMatrixFourierDomain(self.w_lambda,self.grid_lambda)
        self.W_mu =     ModelMatrixFourierDomain(self.w_mu,self.grid_mu)
        self.W_nu =     ModelMatrixFourierDomain(self.w_nu,self.grid_nu)

        if (w_omega is not None) and (Of is not 1) and (O is not 1):
            self.w_omega,self.Of,self.O = w_omega,Of,O

            self.Os = w_omega.shape[0]
            self.os = arange(self.Os)
            self.op = (self.O-1)//2
            self.O_nls = 8

            self.grid_omega = arange(self.O)*2*pi/Of
            self.grid_omega_nls = linspace(0,self.O*2*pi/Of,self.O_nls)
            self.W_omega = ModelMatrixFourierDomain(self.w_omega,self.grid_omega)
            self.dfe = DecoupledFrequencyEstimation((self.L,self.M,self.N,self.O),
                                                    [
                                                     NonlinearLeastSquares(self.Ls,self.L,self.grid_lambda_nls,self.W_lambda),
                                                     NonlinearLeastSquares(self.Ms,self.M,self.grid_mu_nls,self.W_mu),
                                                     NonlinearLeastSquares(self.Ns,self.N,self.grid_nu_nls,self.W_nu),
                                                     NonlinearLeastSquares(self.Os,self.O,self.grid_omega_nls,self.W_omega)
                                                     ],
                                                    [
                                                     FrequencyEstimationPeakBased(self.w_lambda,self.Lf,self.grid_lambda),
                                                     FrequencyEstimationPeakBased(self.w_mu,self.Mf,self.grid_mu),
                                                     FrequencyEstimationPeakBased(self.w_nu,self.Nf,self.grid_nu),
                                                     FrequencyEstimationPeakBased(self.w_omega,self.Of,self.grid_omega)
                                                    ],
                                                    [self.W_lambda,self.W_mu,self.W_nu,self.W_omega])
            self.Md = 4
            self.psi_hat = zeros((4,2))
            self.peak_selector = tuple((array([self.L,self.M,self.N,self.O])-1)//2)

        else:

            self.dfe = DecoupledFrequencyEstimation((self.L,self.M,self.N),
                                                    [
                                                     NonlinearLeastSquares(self.Ls,self.L,self.grid_lambda_nls,self.W_lambda),
                                                     NonlinearLeastSquares(self.Ms,self.M,self.grid_mu_nls,self.W_mu),
                                                     NonlinearLeastSquares(self.Ns,self.N,self.grid_nu_nls,self.W_nu)],
                                                    [
                                                     FrequencyEstimationPeakBased(self.w_lambda,self.Lf,self.grid_lambda),
                                                     FrequencyEstimationPeakBased(self.w_mu,self.Mf,self.grid_mu),
                                                     FrequencyEstimationPeakBased(self.w_nu,self.Nf,self.grid_nu)
                                                    ],
                                                    [self.W_lambda,self.W_mu,self.W_nu])
            self.Md = 3
            self.psi_hat = zeros((3,2))
            self.peak_selector = tuple((array([self.L,self.M,self.N])-1)//2)

    def fit_single_target_model(self,X,dim):
        """
        Fits a single-target model in the peak neighborhood for the specified dimension.

        X       -- 3D or 4D peak neighborhood in the Fourier domain, with
                   dimensions corresponding to range, velocity, and
                   angle (azimuth and/or elevation)
        dim     -- Dimension in which the model is fitted.

        mse     -- Mean squared error of the model fit.
        """

        s=X.shape
        sel_p = tuple((asarray(s)-1)//2)
        sl = slice_1d(self.peak_selector,dim)
        P = 20*log10(abs(X))

        R = self.dfe.covariance_matrix(X,dim)

        fe_peak_based = self.dfe.frequency_estimation_peak_based[dim]
        W =  self.dfe.model_matrices[dim]
        psi_hat_0 = fe_peak_based(self.peak_selector[dim],P[sl])
        self.psi_hat[dim,:] = psi_hat_0
        w_vec_psi = W(psi_hat_0)
        P_W  = dot(w_vec_psi,w_vec_psi.conj().T)/dot(w_vec_psi.conj().T,w_vec_psi)
        mse = real(trace(dot((eye(s[dim]) - P_W),R)))/s[dim]

        return mse

    def __call__(self,X,psi_a):
        """
        Processing steps as shown in Figure 2 in [1] extended to process a
        fourth frequency dimension corresponding to elevation angle.

        X       -- 3D or 4D peak neighborhood in the Fourier domain, with
                   dimensions corresponding to range, velocity, and
                   angle (azimuth and/or elevation)
        psi_a   -- lower border of the 3D or 4D frequency grid in the peak
                   neighborhood, with dimensions corresponding to range, velocity, and
                   angle (azimuth and/or elevation)

        psi_hat -- 3D or 4D frequency estimates, with dimensions corresponding to
                   range, velocity, and angle (azimuth and/or elevation)

        Note, that two frequencies are returned even if the single-target model
        is accepted. In that case both estimates are identical.
        """
        self.mse_2,self.mse_1,self.psi_a = 100.,zeros(0),asarray(psi_a)
        self.dfe.Ta = 20*log10(abs(X[self.peak_selector])) - 6

        for dim in arange(self.Md):
            self.mse_1 = append(self.mse_1,self.fit_single_target_model(X,dim))

        res_dim = self.mse_1.argmax()

        if 10*log10(self.mse_1[res_dim]) > self.T_db:
            self.mse_2 = self.dfe.fit_k_target_model(X,res_dim)

            if 10*log10(self.mse_1[res_dim]/self.mse_2) > self.gamma_db:
                self.psi_hat = self.dfe.remaining_dimensions()

        for dim in arange(self.Md):
            self.psi_hat[dim,:] += self.psi_a[dim]

        return self.psi_hat

class ModelMatrix:
    """
    Model matrix in the original domain

    Instance Variables

    Ns -- sample support

    """
    def __init__(self,Ns):
        self.Ns = Ns
        self.ns = arange(self.Ns)

    def __call__(self,psi_k,A=None):
        if A is None:
            A = zeros((self.Ns,psi_k.shape[0]),complex_)
        exp(1j*outer(self.ns,psi_k),A)
        return A/sqrt(self.Ns)

    def derivative(self,psi_k,Ad=None):
        if Ad is None:
            Ad = zeros((psi_k.shape[0],self.Ns,psi_k.shape[0]),complex_)
        else:
            Ad *= 0
        for k,psi_k_c in enumerate(psi_k):
            Ad[k,:,k] = 1j*self.ns*exp(1j*self.ns*psi_k_c)/sqrt(self.Ns)
        return Ad

class ModelMatrixFourierDomain(ModelMatrix):
    """
    Model matrix in the Fourier domain

    Instance Variables

    w    -- window function
    N    -- support of Fourier transform
    grid -- frequency grid
    WF_H -- conjugate transpose of Fourier transform matrix

    """

    def __init__(self,w,grid):
        ModelMatrix.__init__(self,w.shape[0])
        self.w = w/sum(w)
        self.N = grid.shape[0]
        self.grid = grid
        self.WF_H = dot(diag(self.w),exp(1j*outer(self.ns,self.grid))).conj().T

    def __call__(self,psi_k,W=None):
        if W is None:
            W = zeros((self.N,psi_k.shape[0]),complex_)
        dot(self.WF_H, exp(1j*outer(self.ns,psi_k)),W)
        return W

    def derivative(self,psi_k,Wd=None):
        if Wd is None:
            Wd = zeros((psi_k.shape[0],self.N,psi_k.shape[0]),complex_)
        else:
            Wd *= 0
        for k in arange(psi_k.shape[0]):
            Wd[k,:,k] = dot(self.WF_H, 1j*self.ns*exp(1j*self.ns*psi_k[k]))
        return Wd

class NonlinearLeastSquares:
    """
    Nonlinear least squares (NLS) method for frequency estimation

    Instance Variables

    Ns   -- sample support
    N    -- support of Fourier transform
    grid -- frequency grid for evaluating the NLS criterion function
    K    -- Number of targets
    P    -- projection matrix calculated on the frequency grid

    """

    def __init__(self,Ns,N,grid,W_ci,K=2):
        self.Ns,self.N,self.grid,self.K,self.N_grid = Ns,N,grid,K,grid.shape[0]
        self.W,self.W_dot = W_ci,W_ci.derivative
        self.P = zeros((self.N_grid**self.K,self.N,self.N),complex_)

        for i0 in arange(self.N_grid**self.K):
            im = asarray(unravel_index(i0,(self.N_grid,)*self.K))
            W = self.W(self.grid[im])
            self.P[i0,:,:] = dot(W,linalg.pinv(W))
        self.c = zeros(self.N_grid**self.K)

    def __call__(self,R):
        """
        Perform NLS optimization by evaluating the criterion function on a
        coarse grid, followed by Gauss-Newton iterations.

        R      -- covariance matrix

        psi_hat -- estimated frequencies
        mse    -- mean squared error of the K-target model fit

        """
        psi_hat_coarse = self._evaluate_on_grid(R)
        psi_hat = self._gauss_newton(R,psi_hat_coarse)

        return psi_hat % (2*pi)

    def _evaluate_on_grid(self,R):
        """
        Perform NLS optimization by evaluating the criterion function on a
        coarse grid.

        R      -- covariance matrix

        psi_hat -- estimated frequencies

        """
        self.c = real(trace(dot(self.P,R),axis1=1,axis2=2))
        im = asarray(unravel_index(self.c.argmax(),(self.N_grid,)*self.K))
        return self.grid[im]

    def _gauss_newton(self,R,psi_hat_init,N_it=4):
        """
        Perform NLS optimization by Gauss-Newton iterations.

        R           -- covariance matrix
        psi_hat_init -- starting point for iterations
        N_it        -- maximal number of iterations (default four)

        psi_hat      -- estimated frequencies

        """
        psi_hat = psi_hat_init.copy()
        g,H = zeros(self.K),diag(ones(self.K))
        W = zeros((self.N,self.K),complex_)
        Wd = zeros((self.K,self.N,self.K),complex_)

        for i in arange(N_it):
            W,Wd = self.W(psi_hat,W),self.W_dot(psi_hat,Wd)
            Wp = linalg.pinv(W)
            self._gradient(W,Wp,Wd,R,g)
            self._hessian(W,Wp,Wd,R,H)
            Hig = dot(linalg.pinv(H),g)
            sl = self._step_length(Hig,psi_hat,R)
            psi_hat -= sl*Hig

            if sqrt(inner(Hig,Hig)) < 1e-5*2*pi/self.Ns:
                break
        return psi_hat

    def _gradient(self,W,Wp,Wd,R,g):
        """
        Calculate the gradient of the NLS criterion function

        W  -- model matrix for current iteration
        Wp -- pseudo inverse of model matrix for current iteration
        Wd -- derivative of model matrix for current iteration
        R  -- covariance matrix

        g  -- gradient

        """
        P = eye(self.N) - dot(W,Wp)
        PR = dot(P,R)
        for i in arange(self.K):
            g[i] = -2*real(trace(dot(Wp.conj().T,dot(Wd[i,:,:].conj().T,PR))))

    def _hessian(self,W,Wp,Wd,R,H):
        """
        Calculate approximate Hessian matrix of the NLS criterion function

        W  -- model matrix for current iteration
        Wp -- pseudo inverse of model matrix for current iteration
        Wd -- derivative of model matrix for current iteration
        R  -- covariance matrix

        H  -- approximated Hessian matrix

        """

        P = eye(self.N) - dot(W,Wp)
        p_hat = dot(Wp,dot(R,Wp.conj().T))
        for i in arange(self.K):
            for j in arange(self.K):
                H[i,j] = 2.0*real(trace(dot(Wd[i,:,:].conj().T,dot(P,dot(Wd[j,:,:],p_hat)))))

    def _step_length(self,Hig,psi_hat,R):
        """
        Determine step length for update in current iteration

        Hig    -- product of the inverse Hessian matrix and the gradient for current iteration
        psi_hat -- estimated frequencies for current iteration
        R      -- covariance matrix

        """
        c_last,sl_last = 1000.,0.
        for sl in 0.5**arange(10):
            W = self.W(psi_hat - sl*Hig)
            P =  eye(self.N) - dot(W,linalg.pinv(W))
            c = real(trace(dot(P,R)))
            if c_last <= c:
                break
            else:
                c_last = c
                sl_last = sl
        return sl_last

class FrequencyEstimationPeakBased:
    """
    Frequency estimation via peak position refinement.

    Sub-grid accuracy is achieved by a look-up table (LUT) approach.

    Instance Variables

    w      -- window function
    Nf     -- support of Fourier transform
    grid   -- local frequency grid in the vicinity of the periodogram peak
    lut    -- LUT for sub-grid shift over logarithmic peak power ratio
    sz_lut -- step size of logarithmic power ratio values

    """
    def __init__(self,w,Nf,grid,N_lut=100):
        self.w = w
        self.Nf = Nf
        self.grid = grid
        self.lut = None
        self.sz_lut = None
        self.calc_lut(N_lut)

    def calc_lut(self,N_lut=100):
        """
        Calculate LUT for sub-grid shift over logarithmic peak power ratio

        N_lut -- number of entries in the LUT

        """
        N = self.w.shape[0]
        psi_s = 2*pi*linspace(0.5,0.0,N_lut)/self.Nf
        x = exp(1j*outer(psi_s,arange(N)))
        P = 20.*log10(abs(fft.fft(self.w*x,int(self.Nf),axis=1)))
        ratio = P[:,0] - P[:,1]
        if ratio[-1] == inf:
            ratio[-1] = -psi_s[-3]*(ratio[-2]-ratio[-3])/(psi_s[-2]-psi_s[-3]) + ratio[-3]
        sz_lut = (ratio.max()-ratio.min())/(N_lut-1)
        lut = zeros(N_lut+1)
        ratio /= sz_lut
        for i in arange(N_lut):
            i_min = abs(i-ratio).argmin()
            if i_min == N_lut-1:
                i_min = N_lut-2
            lut[i] = psi_s[i_min] + (psi_s[i_min+1]-psi_s[i_min])*(i-ratio[i_min])/(ratio[i_min+1]-ratio[i_min])
        self.lut = lut
        self.sz_lut = sz_lut

    def __call__(self,maximizer,P,grid=None):
        """
        Peak-based frequency estimation with sub-grid accuracy

        maximizer -- peak frequency on grid
        P         -- logarithmic periodogram, i.e power spectrum
        grid      -- local frequency grid in the vicinity of the periodogram peak

        psi_hat   -- frequency, i.e. sub-grid peak position
        """

        if grid is None:
            grid = self.grid
        N = P.shape[0]
        N_lut = self.lut.size -1
        ratio_l = array([P[maximizer] - P[(maximizer - 1) % N]])
        ratio_r = array([P[maximizer] - P[(maximizer + 1) % N]])
        ratio = minimum(ratio_l,ratio_r)
        i_lut = ratio/self.sz_lut
        i_lut[i_lut > N_lut-1] = N_lut-1
        i_lut_g = floor(i_lut).astype('int')
        psi_s_lut = self.lut[i_lut_g] + (self.lut[i_lut_g+1]-self.lut[i_lut_g])*(i_lut-i_lut_g)
        psi_hat = (grid[maximizer] + sign(ratio_l-ratio_r)*psi_s_lut) % (2*pi)
        return psi_hat

class DecoupledFrequencyEstimation:
    """
    Decoupled frequency estimation for an arbitrary number of data dimensions.

    s                                    -- Shape of the data
    frequency_estimation_high_resolution -- List of high-resolution estimators for
                                            model fitting in the resolution dimension.
    frequency_estimation_peak_based      -- List of peak-based estimators for model fitting
                                            in the remaining dimensions.
    model_matrices                       -- List of model matrices per dimension.
    K                                    -- Model order, i.e. number of targets. Note
                                            that the current implementation is restricted to K=2.
    Ta                                   -- Threshold for covariance matrix averaging over the
                                            remaining dimensions .

    psi_hat                              -- Frequency estimates
    """

    def __init__(self,s,
                 frequency_estimation_high_resolution,
                 frequency_estimation_peak_based,
                 model_matrices,K=2,
                 Ta=-100):
        self.frequency_estimation_high_resolution = frequency_estimation_high_resolution
        self.frequency_estimation_peak_based = frequency_estimation_peak_based
        self.model_matrices = model_matrices
        self.s,self.K = s,K
        self.Ta = Ta
        self.psi_hat = zeros((len(s),self.K))

    def _swap_shape(self,s,dim):
        s_swap = asarray(s)
        s_swap[0],s_swap[dim] = s_swap[dim],s_swap[0]
        s_swap = tuple(s_swap)
        return s_swap

    def _swap(self,X,s=None,dim=None):
        if (s is not None) and (dim is not None):
            s_swap = asarray(self._swap_shape(s,dim))
        else:
            s_swap=asarray(self.s_swap)
            dim = self.res_dim
        s_flat = (s_swap[0],prod(s_swap[1:]))
        return reshape(swapaxes(X,0,dim),s_flat)

    def _set_resolution_dimension(self,res_dim):
        self.res_dim = res_dim
        self.s_swap = self._swap_shape(self.s,res_dim)
        self.idx_o = arange(len(self.s))
        self.idx_o[0],self.idx_o[self.res_dim] = self.res_dim,0
        self.W = self.model_matrices[self.res_dim]

    def covariance_matrix(self,X,dim=None,Ta=None):
        """
        Covariance matrix estimation in the peak neighborhood.

        X        -- Multidimensional data
        dim      -- Data dimension; if not specified the first dimension is used
        Ta       -- Threshold for averaging in the remaining dimensions

        returns  -- Covariance matrix in given dimension
        """
        if Ta is None:
            Ta = self.Ta
        if dim is not None:
            X = self._swap(X,X.shape,dim)
        P = 20*log10(abs(X))
        sel = P[P.argmax(axis=0),arange(X.shape[1])] >= Ta
        Xs = X[:,sel]
        return dot(Xs,Xs.conj().T)/Xs.shape[1]

    def fit_k_target_model(self,X,res_dim):
        """
        Frequency estimation in the resolution dimension by fitting a K-target model.

        The resulting frequencies are stored internally and are subsequently used to
        obtain corresponding frequencies in the remaining dimensions.

        X        -- Multidimensional data
        res_dim  -- Resolution dimension

        mse      -- Mean squared error of the model fit.
        """
        self._set_resolution_dimension(res_dim)
        self.X = self._swap(X)

        R = self.covariance_matrix(self.X)
        psi_hat = self.frequency_estimation_high_resolution[self.res_dim](R)
        Wk  = self.W(psi_hat)
        P   = eye(self.s_swap[0]) - dot(Wk,linalg.pinv(Wk))
        mse = real(trace(dot(P,R)))/self.s_swap[0]
        self.psi_hat[self.res_dim,:] = psi_hat

        return mse

    def remaining_dimensions(self,X=None,psi_res_dim=None):
        """
        Frequency estimation in remaining dimension based on resolved frequencies
        in the resolution dimension.

        """
        if psi_res_dim is None:
            psi_res_dim = self.psi_hat[self.res_dim,:]
        if X is None:
            X=self.X

        P = 20*log10(abs(dot(linalg.pinv(self.W(psi_res_dim)),X)))
        P = reshape(P,(self.K,) + self.s_swap[1:])

        for k in arange(self.K):
            Pk = P[k,:]
            im = unravel_index(Pk.argmax(),self.s_swap[1:])
            for d in arange(len(self.s)-1):
                sl = slice_1d(im,d)
                self.psi_hat[self.idx_o[d+1],k] = self.frequency_estimation_peak_based[self.idx_o[d+1]](im[d],Pk[sl])

        return self.psi_hat

def slice_1d(i,dim):
    sl = list(i)
    sl[dim] = slice(None)
    return sl

reason: from numpy


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

from __future__ import division, print_function, absolute_import, unicode_literals

import datetime as dt

import numpy as np
import pytest

import aacgm2 as aacgmv2
from aacgm2._aacgmv2 import A2G, G2A, TRACE, BADIDEA, ALLOWTRACE, GEOCENTRIC

date = (2015, 1, 1, 0, 0, 0)
dtObj = dt.datetime(*date)

def test_module_structure():
    assert aacgmv2
    assert aacgmv2.convert
    assert aacgmv2.convert_mlt
    assert aacgmv2.subsol
    assert aacgmv2.set_coeff_path
    assert aacgmv2.AACGM_v2_DAT_PREFIX
    assert aacgmv2.IGRF_12_COEFFS

def test_output_type():
    lat, lon, alt = aacgmv2.convert(60, 0, 300, dtObj)
    print(type(lat))
    print(lat.shape)
    print(lat.size)
    assert isinstance(lat, np.ndarray)
    assert isinstance(lon, np.ndarray)
    assert isinstance(alt, np.ndarray)

    lat, lon, alt = aacgmv2.convert([60], [0], [300], dtObj)
    assert isinstance(lat, np.ndarray)
    assert isinstance(lon, np.ndarray)
    assert isinstance(alt, np.ndarray)

    lat, lon, alt = aacgmv2.convert([60, 61], [0, 0], [300, 300], dtObj)
    assert isinstance(lat, np.ndarray)
    assert isinstance(lon, np.ndarray)
    assert isinstance(alt, np.ndarray)

    lat, lon, alt = aacgmv2.convert([60, 61, 62], 0, 300, dtObj)
    assert isinstance(lat, np.ndarray)
    assert isinstance(lon, np.ndarray)
    assert isinstance(alt, np.ndarray)

    lat, lon, alt = aacgmv2.convert(np.array([60, 61, 62]), 0, 300, dtObj)
    assert isinstance(lat, np.ndarray)
    assert isinstance(lon, np.ndarray)
    assert isinstance(alt, np.ndarray)

    lat, lon, alt = aacgmv2.convert(np.array([[60, 61, 62], [63, 64, 65]]), 0, 300, dtObj)
    assert isinstance(lat, np.ndarray)
    assert isinstance(lon, np.ndarray)
    assert isinstance(alt, np.ndarray)

def test_output_shape_size():
    lat, lon, alt = aacgmv2.convert(60, 0, 300, dtObj)
    assert lat.shape == tuple()
    assert lon.shape == tuple()
    assert alt.shape == tuple()
    assert lat.size == 1
    assert lon.size == 1
    assert alt.size == 1

    lat, lon, alt = aacgmv2.convert([60], [0], [300], dtObj)
    assert lat.shape == (1,)
    assert lon.shape == (1,)
    assert alt.shape == (1,)
    assert lat.size == 1
    assert lon.size == 1
    assert alt.size == 1

    lat, lon, alt = aacgmv2.convert([60, 61], [0, 0], [300, 300], dtObj)
    assert lat.shape == (2,)
    assert lon.shape == (2,)
    assert alt.shape == (2,)
    assert lat.size == 2
    assert lon.size == 2
    assert alt.size == 2

    lat, lon, alt = aacgmv2.convert([60, 61, 62], 0, 300, dtObj)
    assert lat.shape == (3,)
    assert lon.shape == (3,)
    assert alt.shape == (3,)
    assert lat.size == 3
    assert lon.size == 3
    assert alt.size == 3

    lat, lon, alt = aacgmv2.convert(np.array([60, 61, 62]), 0, 300, dtObj)
    assert lat.shape == (3,)
    assert lon.shape == (3,)
    assert alt.shape == (3,)
    assert lat.size == 3
    assert lon.size == 3
    assert alt.size == 3

    lat, lon, alt = aacgmv2.convert(np.array([[60, 61, 62],
                                         [63, 64, 65]]),
                               0, 300, dtObj)
    assert lat.shape == (2, 3)
    assert lon.shape == (2, 3)
    assert alt.shape == (2, 3)
    assert lat.size == 6
    assert lon.size == 6
    assert alt.size == 6

def test_convert_result_values_shape():
    lat, lon, alt = aacgmv2.convert(np.array([[60, 61, 62],
                                         [63, 64, 65]]),
                               0, 300, dtObj)
    aacgmv2._aacgmv2.setDateTime(*date)
    np.testing.assert_allclose((lat[0, 0], lon[0, 0], alt[0, 0]),
            aacgmv2._aacgmv2.aacgmConvert(60, 0, 300, G2A), rtol=1e-5)
    np.testing.assert_allclose((lat[0, 1], lon[0, 1], alt[0, 1]),
            aacgmv2._aacgmv2.aacgmConvert(61, 0, 300, G2A), rtol=1e-5)
    np.testing.assert_allclose((lat[0, 2], lon[0, 2], alt[0, 2]),
            aacgmv2._aacgmv2.aacgmConvert(62, 0, 300, G2A), rtol=1e-5)
    np.testing.assert_allclose((lat[1, 0], lon[1, 0], alt[1, 0]),
            aacgmv2._aacgmv2.aacgmConvert(63, 0, 300, G2A), rtol=1e-5)
    np.testing.assert_allclose((lat[1, 1], lon[1, 1], alt[1, 1]),
            aacgmv2._aacgmv2.aacgmConvert(64, 0, 300, G2A), rtol=1e-5)
    np.testing.assert_allclose((lat[1, 2], lon[1, 2], alt[1, 2]),
            aacgmv2._aacgmv2.aacgmConvert(65, 0, 300, G2A), rtol=1e-5)

def test_convert_datetime_date():
    lat_1, lon_1, alt_1 = aacgmv2.convert(60, 0, 300, dt.date(2013, 12, 1))
    lat_2, lon_2, alt_2 = aacgmv2.convert(60, 0, 300, dt.datetime(2013, 12, 1, 0, 0, 0))
    assert lat_1 == lat_2
    assert lon_1 == lon_2
    assert alt_1 == alt_2

def test_convert_result_values_G2A_coeff():
    lat_p, lon_p, alt_p = aacgmv2.convert(60, 0, 300, dtObj)
    aacgmv2._aacgmv2.setDateTime(*date)
    lat_c, lon_c, alt_c = aacgmv2._aacgmv2.aacgmConvert(60, 0, 300, G2A)
    assert lat_p == lat_c
    assert lon_p == lon_c
    assert alt_p == alt_c

def test_convert_result_values_A2G_coeff():
    lat_p, lon_p, alt_p = aacgmv2.convert(60, 0, 300, dtObj, a2g=True)
    aacgmv2._aacgmv2.setDateTime(*date)
    lat_c, lon_c, alt_c = aacgmv2._aacgmv2.aacgmConvert(60, 0, 300, A2G)
    assert lat_p == lat_c
    assert lon_p == lon_c
    assert alt_p == alt_c

def test_convert_result_values_G2A_trace():
    lat_p, lon_p, alt_p = aacgmv2.convert(60, 0, 300, dtObj, trace=True)
    aacgmv2._aacgmv2.setDateTime(*date)
    lat_c, lon_c, alt_c = aacgmv2._aacgmv2.aacgmConvert(60, 0, 300, G2A | TRACE)
    assert lat_p == lat_c
    assert lon_p == lon_c
    assert alt_p == alt_c

def test_convert_result_values_A2G_trace():
    lat_p, lon_p, alt_p = aacgmv2.convert(60, 0, 300, dtObj, a2g=True, trace=True)
    aacgmv2._aacgmv2.setDateTime(*date)
    lat_c, lon_c, alt_c = aacgmv2._aacgmv2.aacgmConvert(60, 0, 300, A2G | TRACE)
    assert lat_p == lat_c
    assert lon_p == lon_c
    assert alt_p == alt_c

def test_convert_result_values_allowtrace():
    lat, lon, alt = aacgmv2.convert(60, 0, [300, 5000], dtObj, allowtrace=True)
    aacgmv2._aacgmv2.setDateTime(*date)
    np.testing.assert_allclose((lat[0], lon[0], alt[0]),
            aacgmv2._aacgmv2.aacgmConvert(60, 0, 300, ALLOWTRACE), rtol=1e-5)
    np.testing.assert_allclose((lat[1], lon[1], alt[1]),
            aacgmv2._aacgmv2.aacgmConvert(60, 0, 5000, ALLOWTRACE), rtol=1e-5)

def test_convert_result_values_badidea():
    lat, lon, alt = aacgmv2.convert(60, 0, [300, 5000], dtObj, badidea=True)
    aacgmv2._aacgmv2.setDateTime(*date)
    np.testing.assert_allclose((lat[0], lon[0], alt[0]),
            aacgmv2._aacgmv2.aacgmConvert(60, 0, 300, BADIDEA), rtol=1e-5)
    np.testing.assert_allclose((lat[1], lon[1], alt[1]),
            aacgmv2._aacgmv2.aacgmConvert(60, 0, 5000, BADIDEA), rtol=1e-5)

def test_convert_result_values_geocentric():
    lat_p, lon_p, alt_p = aacgmv2.convert(60, 0, 300, dtObj, geocentric=True)
    aacgmv2._aacgmv2.setDateTime(*date)
    lat_c, lon_c, alt_c = aacgmv2._aacgmv2.aacgmConvert(60, 0, 300, GEOCENTRIC)
    assert lat_p == lat_c
    assert lon_p == lon_c
    assert alt_p == alt_c

def test_warning_below_ground():
    with pytest.warns(UserWarning):
        aacgmv2.convert(60, 0, -1, dtObj)
    with pytest.warns(UserWarning):
        aacgmv2.convert(60, 0, [300, -1], dtObj)

def test_exception_maxalt():
    with pytest.raises(ValueError):
        aacgmv2.convert(60, 0, 2001, dtObj)
    with pytest.raises(ValueError):
        aacgmv2.convert(60, 0, [300, 2001], dtObj)

    # the following should not raise exceptions
    aacgmv2.convert(60, 0, 2001, dtObj, trace=True)
    aacgmv2.convert(60, 0, 2001, dtObj, allowtrace=True)
    aacgmv2.convert(60, 0, 2001, dtObj, badidea=True)

def test_exception_lat90():
    with pytest.raises(ValueError):
        aacgmv2.convert(91, 0, 300, dtObj)
    with pytest.raises(ValueError):
        aacgmv2.convert(-91, 0, 300, dtObj)
    with pytest.raises(ValueError):
        aacgmv2.convert([60, 91], 0, 300, dtObj)
    with pytest.raises(ValueError):
        aacgmv2.convert([60, -91], 0, 300, dtObj)

    # the following should not raise exceptions
    aacgmv2.convert(90, 0, 300, dtObj)
    aacgmv2.convert(-90, 0, 300, dtObj)

def test_forbidden():
    np.testing.assert_equal(aacgmv2.convert(7, 0, 0), (np.nan, np.nan, np.nan))

def test_MLT_forward_backward():
    mlon = aacgmv2.convert_mlt(12, dtObj, m2a=True)
    mlt = aacgmv2.convert_mlt(mlon, dtObj)
    np.testing.assert_allclose(mlt, 12)

def test_MLT_a2m():
    mlt = aacgmv2.convert_mlt([1, 12, 23], dt.datetime(2015, 2, 24, 14, 0, 15))
    np.testing.assert_allclose(mlt, [9.058207,  9.791541, 10.524874], rtol=1e-6)

def test_MLT_m2a():
    mlon = aacgmv2.convert_mlt([1, 12, 23], dt.datetime(2015, 2, 24, 14, 0, 15), m2a=True)
    np.testing.assert_allclose(mlon, [240.12689,  45.12689, 210.12689], rtol=1e-6)


text:
from numpy import concatenate,zeros,mean,std,uint16
#from CAServer import casput
from serial import Serial
from time import time, sleep,gmtime, strftime
from sys import stdout
import os.path
import matplotlib.pyplot as plt
from pdb import pm
from struct import unpack as struct_unpack
import logging
from logging import error,warn,info,debug
import circular_buffer_LL

__version__ = '1.0.5' #cleaned code, added comments
__date__ = "11-10-17"

class DI245_DL(object):
    def __init__(self,SN = 1):
        #this will create an instance
        self.SN = SN
        print('Initializing DL...')
        self.init_device(SN = self.SN)
        print('initializing buffer...')
        self.buffer = circular_buffer_LL.server(size = (4,4320000), var_type = 'int16')#4320000

    def init_device(self,SN = 1, scan_lst = ['0','1','2','3'], phys_ch_lst = ['0','1','2','3'], gain_lst = ['0.25','T-thrmc','0.25','T-thrmc']):#4320000
        self.SN = SN
        flag = True
        time_out = 0
        self._use_com_port('COM37')
        flag = False
        while flag:
            try:
                self._find_com_port()
                flag = False
            except:
                time_out += 1
                if time_out >10:
                    flag = False

        #create com port connection and flush it. It doesn't seem to be import to flush but why not.
        #self.ser = Serial(port, baudrate=115200, rtscts=True, timeout=0.1)
        self.ser.flushInput()
        self.ser.flushOutput()

        #Values
        self.cjc_value = -2.0 #this number needs to be tested with ice or another well calibrated thermocouple

        self.calib = [0.05,0,0.05,0,0] # [scale_up, offset_up, scale_down, offset_down, time] 50 mV per atm factory setting

        self.scan_lst = scan_lst
        self.phys_ch_lst = phys_ch_lst
        self.gain_lst = gain_lst

        self.channels_to_read = len(self.scan_lst)
        #self.broadcast = 0

        #dictionary with device parameters such as S\N, device name, ect
        self.dict_DI245 = {}
        self.dict_DI245['Device name'] = self._inquire('A1')[0][2:]
        self.dict_DI245['Firmware version'] = int(self._inquire('A2')[0][2:],16)/100.
        self.dict_DI245['Last Calibration date in hex'] = self._inquire('A7')[0][2:]
        self.dict_DI245['Serial Number'] = self._inquire('NZ')[0][2:]

        #Useful variables
        self.tau = 0.0001 #wait time after write (before) read. This seems to be not necessary.

        for i in self.dict_DI245.keys():
            print i, self.dict_DI245[i]
        print 'Complete: Initialization of the DI-245 with S\N', self.dict_DI245['Serial Number']

        self._config_channels()

    def _use_com_port(self,port):
        self.ser = Serial(port, baudrate=115200, rtscts=True, timeout=0.1)
        self.stop_scan()

    def _find_com_port(self):
        #this function will scan com ports and find DI-245 devices by sending command A1 and receiving a word 2450 back
        import serial.tools.list_ports
        list = serial.tools.list_ports.comports()
        connected = []

        for element in list:
            print element.device
            try:
                self.ser = Serial(element.device, baudrate=115200, rtscts=True, timeout=0.1)
                self.stop_scan()
                try:
                    print 'DI device is avaiable:',self._inquire('A1')[0][2:]
                    if self._inquire('A1')[0][2:] == '2450':
                        print("the requested device is connected")
                    else:
                        print("DI-245 with requested serial number is not connceted")
                        self.ser.close()
                        print("closing com port")
                except:
                    self.ser._close_port()
            except:
                pass

    def _read(self):
        sleep(self.tau)
        lines =  self.ser.readlines() # Friedrich doesn't like command readline, since "The method readline() reads one entire line from the file.
        #However, my tests have showed if you know what you expect to get it always works. Just make sure you have "
        stdout.flush()
        return lines[0] #I am not sure why do I have to specify index 0 in this list

    def _write(self,command):
        self.ser.flushOutput()
        self.ser.write(command)

    def _inquire(self,command):
        self.ser.flushInput()
        self.ser.flushOutput()
        self.ser.write(command)
        return self.ser.readlines()

    def _close_port(self):
        self.ser.close()

    #Scanning\data acquisition section
    #an example: chn(0x20)member(0x20)value(0x0D)
    #2-byte value needs to be converted from binary to int. The binary 2 byte start counting from right.
    #The values in the function definition are default values in case user does not specify them.
    def _config_channels(self):
        #self.scan_lst = ['0','1','2','3']
        #self.phys_ch_lst = ['0','1','2','3']
        #self.gain_lst = ['2.5','2.5','2.5','2.5']
        #This method needs to be rewritten to take a list instead of a set of parameters. This will help with scalability.

        _config_dict_gain = {}
        _config_dict_gain['0.010'] = '00101'
        _config_dict_gain['0.025'] = '00100'
        _config_dict_gain['0.05'] = '00011'
        _config_dict_gain['0.1'] = '00010'
        _config_dict_gain['0.25'] = '00001'
        _config_dict_gain['0.5'] = '00000'
        _config_dict_gain['1'] = '01101'
        _config_dict_gain['2.5'] = '01100'
        _config_dict_gain['5'] = '01011'
        _config_dict_gain['10'] = '01010'
        _config_dict_gain['25'] = '01001'
        _config_dict_gain['50'] = '01000'
        _config_dict_gain['B-thrmc'] = '10000'
        _config_dict_gain['E-thrmc'] = '10001'
        _config_dict_gain['J-thrmc'] = '10010'
        _config_dict_gain['K-thrmc'] = '10011'
        _config_dict_gain['N-thrmc'] = '10100'
        _config_dict_gain['R-thrmc'] = '10101'
        _config_dict_gain['S-thrmc'] = '10110'
        _config_dict_gain['T-thrmc'] = '10111'

        for i in range(len(self.scan_lst)):
            config_byte = str(int('000'+_config_dict_gain[self.gain_lst[i]]+'0000' + bin(int(self.phys_ch_lst[i]))[2:].zfill(4),2))
            ch_config_command = 'chn '+self.scan_lst[i]+' '+config_byte+' \x0D'

            self.ser.flushInput()
            self.ser.flushOutput()
            result1 = self._inquire(ch_config_command)

        xrate_config_command = 'xrate 4099 2000 \x0D'
        #xrate_config_command = 'xrate 62 127 \x0D' #good slow
        #xrate_config_command = 'xrate 1871 10 \x0D'
        #xrate_config_command = 'xrate 2051 2000\x0D'
        result2 = self._inquire(xrate_config_command) #config xrate for the data acquisition

        print 'Complete: Channel(s) and xrate configuration'
        return result1, result2

    #this method sends a proper command to start the scan.
    def start_scan(self):
        read_byte_temp = ""
        while read_byte_temp != 'S1':
            self._write('S1')
            read_byte_temp = self.ser.read(2) #read 2-byte echo respons
        print 'The configured measurement(s) has(have) started'

    #This method reads one
    def read_number(self):
        value_array = zeros((self.channels_to_read))
        for j in range(self.channels_to_read):
            #value_array[2*j] = time.time()
            tempt_t = time()
            read_byte = self.ser.read(2)
            read_byte = bin(struct_unpack("H", read_byte)[0])[2:].zfill(16)
            read_byte_lst = list(read_byte)
            sync_byte = read_byte_lst[15] #this is the byte 0 that is issued in DI-245 for sync. 0 stand for the beginning of channel(s) data stream. Hence, every set of readouts  starts with 0.
            del(read_byte_lst[15]) #This is so called sync bit, it is 0 in the beginning of a new sequence of data
            del(read_byte_lst[7])
            read_byte = ""
            for i in read_byte_lst:
                read_byte += str(i)
            int_val = int(read_byte,2)
            value_array[j] = int_val
        return value_array

    def _waiting(self):
        logging.debug('In ' + str(self.ser.in_waiting) + '  Out '+ str(self.ser.out_waiting)) #these commands tell how many bytes are in RS-232 buffer in and out

    def stop_scan(self):
        self._inquire('S0')

    def full_stop(self):
        self.stop_scan()
        self._close_port()

    def _int_to_units(self, array_in, channel):
        if channel == 'all':
            array_out = zeros(4)
            array_out[0]  = (((array_in[0]-8192.)/8192.)*float(self.gain_lst[0]))/self.calib[0]
            array_out[1]  = ((array_in[1]-8192.)*0.036621+100. - self.calib[1])
            array_out[2]  = (((array_in[2]-8192.)/8192.)*float(self.gain_lst[2]))/self.calib[2]
            array_out[3]  = ((array_in[3]-8192.)*0.036621+100. - self.calib[3])
        elif channel == 0:
            array_out = (((array_in-8192.)/8192.)*float(self.gain_lst[channel]))/self.calib[0]
        elif channel == 1:
            array_out  = ((array_in-8192.)*0.036621+100. + self.cjc_value)
        elif channel == 2:
            array_out = (((array_in-8192.)/8192.)*float(self.gain_lst[channel]))/self.calib[2]
        elif channel == 3:
            array_out  = ((array_in-8192.)*0.036621+100. + self.cjc_value)
        elif channel == 'rms_ch0':
            array_out = ((float(self.gain_lst[0]))/self.calib[0])/8192.
        elif channel == 'rms_ch1':
            array_out = (float(self.gain_lst[1]))*0.036621
        elif channel == 'rms_ch2':
            array_out = ((float(self.gain_lst[2]))/self.calib[2])/8192.
        elif channel == 'rms_ch3':
            array_out = (float(self.gain_lst[3]))*0.036621
        return array_out

    def set_calib(self):
        print '_set_calibration'
        #self.calib = [0.05,0,0.05,0,0] # [scale_up, offset_up, scale_down, offset_down, time] 50 mV per atm factory setting
        #while self.RingBufferPointer < 50*1:
        #    time.sleep(1)
        a = mean(((self.buffer.get_last_N(N = 300)[0]-8192.)/8192.)*float(self.gain_lst[0]))
        b = mean(self.buffer.get_last_N(N = 300)[1]-8192.)*0.036621+100.-self.calib[1]
        c = mean(((self.buffer.get_last_N(N = 300)[2]-8192.)/8192.)*float(self.gain_lst[2]))
        d = mean(self.buffer.get_last_N(N = 300)[3]-8192.)*0.036621+100.-self.calib[3]
        self.calib = [a,b ,c ,d , time()]
        print self.calib
        print "-------"

    def save_to_a_file(self):
        debug('save to a file pressed %r' % time())
        pass

    def broadcast(self, method = "on demand"):

        if method == "on demand":
            debug('Broadcasting on demand')
            casput("NIH:DI245."+str(self.dict_DI245['Serial Number'])+".CH1.pressure",self._int_to_units(int(mean(self.buffer.get_last_N(N = 50)[0])),0))
            casput("NIH:DI245."+str(self.dict_DI245['Serial Number'])+".CH2.temperature",self._int_to_units(int(mean(self.buffer.get_last_N(N = 50)[1])),1))
            casput("NIH:DI245."+str(self.dict_DI245['Serial Number'])+".CH3.pressure",self._int_to_units(int(mean(self.buffer.get_last_N(N = 50)[2])),2))
            casput("NIH:DI245."+str(self.dict_DI245['Serial Number'])+".CH4.temperature",self._int_to_units(int(mean(self.buffer.get_last_N(N = 50)[3])),3))
            print("NIH:DI245."+str(self.dict_DI245['Serial Number'])+".CH4.temperature",self._int_to_units(int(mean(self.buffer.get_last_N(N = 50)[3])),3))

        else:
            pass#debug("unknown method(%r) requested" % method)

if __name__ == "__main__":
    logging.basicConfig(filename='di_245_DL.log', level=logging.DEBUG, format="%(asctime)s %(levelname)s: %(message)s")
    dev = DI245_DL()
    print("----------------------")
    print("DI-245 Device Level library")
    print("----------------------")
    print("-----Basic functions-----")
    print("RING BUFFER commands need to be added")
    print("dev.init_device() #-> will use preconfigured sequence")
    print("dev._use_com_port(port)")
    print("dev._find_com_port()")
    print("dev._read()")
    print("dev._write(self,command)")
    print("dev._inquire(command)")
    print("dev._waiting()")
    print("dev._close_port()")
    print("-----Advance functions-----")
    print("dev.config_channels()")
    print("dev.start_scan()")
    print("dev.read_number()")
    print("dev.stop_scan()")
    print("dev.full_stop()")
    print("dev._int_to_units()")
    print("dev.set_calib()")
    print("dev.save_to_a_file()")
    print("dev.broadcast()")

max_stars_repo_name: vstadnytskyi/drivers

max_issues_repo_name: vstadnytskyi/drivers

max_forks_repo_name: vstadnytskyi/drivers


text:
''' module face_detection_engine.py

    Purpose: implement a face detecting engine
'''

from enum import Enum

from numpy import asarray
from PIL import Image

from face_embedding_engine import FaceEmbeddingModelEnum

# We use descriptive variable and function names so
# disable the pylint warning for long lines
# pylint: disable=line-too-long

SSD_MOBILENET_V2_FACE_MODEL = 'models/ssd_mobilenet_v2_face_quant_postprocess_edgetpu.tflite' # trained with Celebrity imageset

class FaceDetectionMethodEnum(Enum):
    ''' enum DetectionMethod

        Enumerates all methods supported for detecting faces
    '''
    MTCNN = 1 # currently only supported on Ubuntu
    SSD_MOBILENET_V2 = 2 # currently only supported on Coral dev board

class FaceDetectionEngine:
    ''' class FaceDetectionEngine

        Purpose: detect faces in an image
    '''

    def __init__(self, detection_method):
        ''' function constructor

        Constructor for FaceDetectionEngine

        Args:
            detection_method (DetectionMethod): Method to use for detection

        Returns:
            None
        '''

        # We only want to import these modules at run-time since
        # they will only be installed on certain platforms.
        # pylint: disable=import-outside-toplevel, import-error

        self.detection_method = detection_method
        if self.detection_method == FaceDetectionMethodEnum.MTCNN:
            # create the MTCNN detector, using default weights
            print("Using MTCNN for face detection")
            from mtcnn.mtcnn import MTCNN
            self.face_detection_engine = MTCNN()
        elif self.detection_method == FaceDetectionMethodEnum.SSD_MOBILENET_V2:
            # load the MobileNet V2 SSD Face model
            print("Using SSD MobileNet V2 for face detection")
            from edgetpu.detection.engine import DetectionEngine
            self.face_detection_engine = DetectionEngine(SSD_MOBILENET_V2_FACE_MODEL)
        else:
            raise Exception("Invalid detection method: {}".format(detection_method))

    # detect faces
    def detect_faces(self, rgb_array):
        ''' function detect_faces

        Detect any faces that are present in the given image.

        Args:
            rgb_array (numpy.ndarray): An image that may or may not contain faces

        Returns:
            An array of bounding boxes (top_left_x, top_left_y, width, height)
            for each face detected in the given image
        '''

        results = [] # assume no faces are detected

        # detect faces in the image
        if self.detection_method == FaceDetectionMethodEnum.MTCNN:
            detected_faces = self.face_detection_engine.detect_faces(rgb_array)

            # extract the bounding box from the first face
            if len(detected_faces) == 0:
                return results

            for detected_face in detected_faces:
                # note the bounding box is in the format we want
                results.append(tuple(detected_face['box']))

        else: # DetectionMethod.SSD_MOBILENET_V2
            frame_as_image = Image.fromarray(rgb_array)
            detected_faces = self.face_detection_engine.detect_with_image(
                frame_as_image,
                threshold=0.5,
                keep_aspect_ratio=True,
                relative_coord=False,
                top_k=5,
                resample=Image.BOX)

            if len(detected_faces) == 0:
                return results

            # extract the bounding box from the first face
            for detected_face in detected_faces:
                # convert the bounding box to the format we want
                x_1, y_1, x_2, y_2 = detected_face.bounding_box.flatten().astype("int")
                width = abs(x_2 - x_1)
                height = abs(y_2 - y_1)
                result = (x_1, y_1, width, height)
                results.append(result)

        return results

    def extract_face(self, rgb_array, embedding_model):
        ''' function extract_face

        Extract a single face from the given frame

        Args:
            rgb_array (numpy.ndarray): The image that may or may not contain
                            one or more faces
            embedding_model (FaceEmbeddingModelEnum): The model being used for generating
                            embeddings for face images

        Returns:
            If a face is detected, returns an RGB numpy.ndarray of the face extracted from
            the given frame of the dimensions required for the given embedding model.
            Otherwise it returns an empty array.
        '''

        detected_faces = self.detect_faces(rgb_array)
        if len(detected_faces) == 0:
            return []

        if detected_faces[0][2] == 0:
            return []

        x_1, y_1, width, height = tuple(detected_faces[0])
        x_1, y_1 = abs(x_1), abs(y_1)
        x_2, y_2 = x_1 + width, y_1 + height

        # extract a cropped image of the detected face
        face = rgb_array[y_1:y_2, x_1:x_2]

        # resize pixels to the dimension required for the specified embedding model
        image = Image.fromarray(face)
        image = image.resize((160, 160))

        # convert image to numpy array
        face_rgb_array = asarray(image)
        return face_rgb_array


text:
import os
import pandas
import shutil
import numpy as np

raw_images_path = '/home/public/car_exam/raw_images'
output_path = '/home/wang/vin'
ipl_par = os.listdir(raw_images_path)
GROUND_TRUTH_CSV = './ground_truth.csv'

df = pandas.read_csv(GROUND_TRUTH_CSV, encoding='utf-8')
no = np.asarray(df['作业流水号']).astype(str)
vin_no = np.asarray(df['车架号']).astype(str)

lp_dict = {}
for i, n in enumerate(no):
    lp_dict[n] = vin_no[i]

for ip_par in ipl_par:
    if os.path.isdir(raw_images_path + '/' + ip_par):
        for ip in os.listdir(raw_images_path + '/' + ip_par):
            if ip.endswith('.jpg') and '0113' in ip:
                print(ip_par)
                try:
                    shutil.copy(raw_images_path + '/' + ip_par + '/' + ip,
                            output_path + '/{}_{}_'.format(lp_dict[ip_par], ip_par) + ip)

                except KeyError:
                    shutil.copy(raw_images_path + '/' + ip_par + '/' + ip,
                                output_path + '/{}_{}_'.format('UNKNOWN', ip_par) + ip)
                except FileNotFoundError:
                    pass


text:
import numpy as np

def ellipsoid_gravity_field(x, a, b , c, rho):

    G = 6.6695e-11                     # Gravitational constant [kg m**3/s**2]
    aa = a**2
    bb = b**2
    cc = c**2
    errtol = 1e-10

    # find the greatest root LAM
    B = -(x[0]**2 + x[1]**2 + x[2]**2 - aa - bb - cc)
    C = -cc*(x[0]**2 + x[1]**2) + bb*(cc - x[0]**2 - x[2]**2) + aa*(bb +cc - x[1]**2 - x[2]**2)
    D = -bb*cc*x[0]**2 + aa*(-cc*x[1]**2 + bb*(c - x[2])*(c + x[2]))
    poly = np.asarray([1, B, C, D])
    LAM = np.max(np.roots(poly))
    # find elliptic integrals F1 and E1
    phi = np.arcsin(np.sqrt((aa-cc)/(LAM + aa)))# Argument phi, s.t. 0 < phi <= np.pi/2
    k = np.sqrt((aa-bb)/(aa-cc))            # Modulus k, s.t. 0 < k < 1
    F1 = lellipf(phi, k, errtol)
    E1 = lellipe(phi, k, errtol)
    fac1 = 4*np.pi*G*rho*a*b*c/np.sqrt(aa-cc)
    fac2 = np.sqrt((aa-cc)/((aa+LAM)*(bb+LAM)*(cc+LAM)))

    # Components of attraction.
    X = fac1/(aa-bb)*(E1-F1)
    Y = fac1*((-aa+cc)*E1/((aa-bb)*(bb-cc)) + F1/(aa-bb)+ (cc+LAM)*fac2/(bb-cc))

    Z = fac1*(-E1 + (bb+LAM)*fac2)/(cc-bb)
    return np.asarray([X,Y,Z]), LAM

#
# lellipe(phi, k, errtol)
#
# Inputs:
#
#   phi     Input angle vector size 1xN.
#   k       Input parameter vector size 1 or 1xN.
#   errtol  Error tolerance for Carlson's algorithms.
#
# Matlab function to compute Legendre's (incomplete) elliptic integral
# E(phi, k).  Uses a vectorized implementation of Carlson's Duplication Algorithms
# for symmetric elliptic integrals as found in "Computing Elliptic
# Integrals by Duplication," by B. C. Carlson, Numer. Math. 33, 1-16 (1979)
# and also found in ACM TOMS Algorithm 577.  Section 4 in the paper cited
# here describes how to convert between the symmetric elliptic integrals
# and Legendre's elliptic integrals.
#
# Returns NaN's for any argument values outside input range.
#

def lellipe(phi, k, errtol):

    # Argument checking for vectorization:

    phivec = phi
    kvec = k

    snphi = np.sin(phivec)
    csphi = np.cos(phivec)
    snphi2 = snphi**2
    csphi2 = csphi**2
    k2 = kvec**2
    y = 1.0 - k2*snphi2
    onesvec = 1
    f = snphi * rf(csphi2,  y, onesvec, errtol) - k2 * snphi * snphi2 * rd(csphi2, y, onesvec, errtol)/3.0

    return f

#
# lellipf(phi, k, errtol)
#
# Inputs:
#
#   phi     Input angle vector size 1 or 1xN.
#   k       Input parameter vector size 1 or 1xN.
#   errtol  Error tolerance for Carlson's algorithms.
#
# Matlab function to compute Legendre's (incomplete) elliptic integral
# F(phi, k).  Uses a vectorized implementation of Carlson's Duplication Algorithms
# for symmetric elliptic integrals as found in "Computing Elliptic
# Integrals by Duplication," by B. C. Carlson, Numer. Math. 33, 1-16 (1979)
# and also found in ACM TOMS Algorithm 577.  Section 4 in the paper cited
# here describes how to convert between the symmetric elliptic integrals
# and Legendre's elliptic integrals.
#
# Returns NaN's for any argument values outside input range.
#

def lellipf(phi, k, errtol):

    phivec = phi
    kvec = k

    snphi = np.sin(phivec)
    csphi = np.cos(phivec)
    csphi2 = csphi * csphi
    onesvec = 1
    y = onesvec - kvec*kvec * snphi*snphi
    f = snphi * rf(csphi2,  y, onesvec, errtol)
    return f

# Elliptic Integrals by Duplication," by B. C. Carlson, Numer. Math.
# 33, 1-16 (1979).
#
# Returns NaN's for any argument values outside input range.
#
# Algorithm is also from Carlson's ACM TOMS Algorithm 577.
#
# This code is a complete rewrite of the algorithm in vectorized form.
# It was not produced by running a FORTRAN to Matlab converter.
#
# The following text is copied from ACM TOMS Algorithm 577 FORTRAN code:
#
#   X AND Y ARE THE VARIABLES IN THE INTEGRAL RC(X,Y).
#
#   ERRTOL IS SET TO THE DESIRED ERROR TOLERANCE.
#   RELATIVE ERROR DUE TO TRUNCATION IS LESS THAN
#   16 * ERRTOL ** 6 / (1 - 2 * ERRTOL).
#
#   SAMPLE CHOICES:  ERRTOL     RELATIVE TRUNCATION
#                               ERROR LESS THAN
#                    1.D-3      3.D-19
#                    3.D-3      2.D-16
#                    1.D-2      3.D-13
#                    3.D-2      2.D-10
#                    1.D-1      3.D-7
#
# Note by TRH:
#
#   Absolute truncation error when the integrals are order 1 quantities
#   is closer to errtol, so be careful if you want high absolute precision.
#
# Thomas R. Hoffend Jr., Ph.D.
# 3M Company
# 3M Center Bldg. 236-GC-26
# St. Paul, MN 55144
# trhoffendjr@mmm.com
#

def rf(x, y, z, errtol):
    realmin = 1e-100
    realmax = 1e100
    # Argument limits as set by Carlson:
    LoLim = 5.0 * realmin
    UpLim = 5.0 * realmax

    # Define internally acceptable variable ranges for iterations:
    Xi = x
    Yi = y
    Zi = z

    # Carlson's duplication algorithm for Rf:
    Xn = Xi
    Yn = Yi
    Zn = Zi
    Mu = (Xn + Yn + Zn) / 3.0
    Xndev = 2.0 - (Mu + Xn) / Mu
    Yndev = 2.0 - (Mu + Yn) / Mu
    Zndev = 2.0 - (Mu + Zn) / Mu
    epslon = np.max( np.abs([Xndev, Yndev, Zndev]) )
    while epslon >= errtol:
        Xnroot = np.sqrt(Xn)
        Ynroot = np.sqrt(Yn)
        Znroot = np.sqrt(Zn)
        lambda1 = Xnroot * (Ynroot + Znroot) + Ynroot * Znroot
        Xn = 0.25 * (Xn + lambda1)
        Yn = 0.25 * (Yn + lambda1)
        Zn = 0.25 * (Zn + lambda1)
        Mu = (Xn + Yn + Zn) / 3.0
        Xndev = 2.0 - (Mu + Xn) / Mu
        Yndev = 2.0 - (Mu + Yn) / Mu
        Zndev = 2.0 - (Mu + Zn) / Mu
        epslon = np.max( np.abs([Xndev , Yndev , Zndev]) )

    C1 = 1.0 / 24.0
    C2 = 3.0 / 44.0
    C3 = 1.0 / 14.0
    E2 = Xndev * Yndev - Zndev * Zndev
    E3 = Xndev * Yndev * Zndev
    S = 1.0 + (C1 * E2 - 0.1 - C2 * E3) * E2 + C3 * E3
    f = S / np.sqrt(Mu)

    # Return NaN's where input argument was out of range:
    return f

#
# rd(x, y, z, errtol)
#
# Inputs:
#
#   x       Input vector size 1xN.
#   y       Input vector size 1xN.
#   z       Input vector size 1xN.
#   errtol  Error tolerance.
#
# Matlab function to compute Carlson's symmetric elliptic integral Rd.
# Implementation of Carlson's Duplication Algorithm 4 in "Computing
# Elliptic Integrals by Duplication," by B. C. Carlson, Numer. Math.
# 33, 1-16 (1979).
#
# Returns NaN's for any argument values outside input range.
#
# Algorithm is also from Carlson's ACM TOMS Algorithm 577.
#
# This code is a complete rewrite of the algorithm in vectorized form.
# It was not produced by running a FORTRAN to Matlab converter.
#
# The following text is copied from ACM TOMS Algorithm 577 FORTRAN code:
#
#   X AND Y ARE THE VARIABLES IN THE INTEGRAL RC(X,Y).
#
#   ERRTOL IS SET TO THE DESIRED ERROR TOLERANCE.
#   RELATIVE ERROR DUE TO TRUNCATION IS LESS THAN
#   16 * ERRTOL ** 6 / (1 - 2 * ERRTOL).
#
#   SAMPLE CHOICES:  ERRTOL     RELATIVE TRUNCATION
#                               ERROR LESS THAN
#                    1.D-3      3.D-19
#                    3.D-3      2.D-16
#                    1.D-2      3.D-13
#                    3.D-2      2.D-10
#                    1.D-1      3.D-7
#
# Note by TRH:
#
#   Absolute truncation error when the integrals are order 1 quantities
#   is closer to errtol, so be careful if you want high absolute precision.
#
# Thomas R. Hoffend Jr., Ph.D.
# 3M Company
# 3M Center Bldg. 236-GC-26
# St. Paul, MN 55144
# trhoffendjr@mmm.com
#

def rd(x, y, z, errtol ):

    realmin = 1e-100
    realmax = 1e100

    # Argument limits as set by Carlson:
    LoLim = 5.0 * realmin
    UpLim = 5.0 * realmax

    # Define internally acceptable variable ranges for iterations:
    Xi = x
    Yi = y
    Zi = z

    # Carlson's duplication algorithm for Rf:
    Xn = Xi
    Yn = Yi
    Zn = Zi
    sigma = 0.0
    power4 = 1.0

    Mu = (Xn + Yn + 3.0 * Zn) * 0.2
    Xndev = (Mu - Xn) / Mu
    Yndev = (Mu - Yn) / Mu
    Zndev = (Mu - Zn) / Mu
    epslon = np.max( np.abs([Xndev, Yndev, Zndev]) )
    while epslon >= errtol:
        Xnroot = np.sqrt(Xn)
        Ynroot = np.sqrt(Yn)
        Znroot = np.sqrt(Zn)
        lambda1 = Xnroot * (Ynroot + Znroot) + Ynroot * Znroot
        sigma = sigma + power4 / (Znroot * (Zn + lambda1))
        power4 = 0.25 * power4
        Xn = 0.25 * (Xn + lambda1)
        Yn = 0.25 * (Yn + lambda1)
        Zn = 0.25 * (Zn + lambda1)
        Mu = (Xn + Yn + 3.0 * Zn) * 0.2
        Xndev = (Mu - Xn) / Mu
        Yndev = (Mu - Yn) / Mu
        Zndev = (Mu - Zn) / Mu
        epslon = np.max( np.abs([Xndev, Yndev, Zndev]) )

    C1 = 3.0 / 14.0
    C2 = 1.0 / 6.0
    C3 = 9.0 / 22.0
    C4 = 3.0 / 26.0
    EA = Xndev * Yndev
    EB = Zndev * Zndev
    EC = EA - EB
    ED = EA - 6.0 * EB
    EF = ED + EC + EC
    S1 = ED * (-C1 + 0.25 * C3 * ED - 1.50 * C4 * Zndev * EF)
    S2 = Zndev * (C2 * EF + Zndev * (-C3 * EC + Zndev * C4 * EA))
    f = 3.0 * sigma + power4 * (1.0 + S1 + S2) / (Mu * np.sqrt(Mu))

    return f

max_stars_repo_name: Aerospace-AI/Hovering-with-Altimetry

max_issues_repo_name: Aerospace-AI/Hovering-with-Altimetry

max_forks_repo_name: Aerospace-AI/Hovering-with-Altimetry