text:
# coding=utf-8
# #!/usr/bin/env python

# default library
import os
import sys
import struct
import copy

import yo_network_info
sys.path.append(os.path.join(yo_network_info.PATH_BASE, 'lib'))

# sys.path.append('/usr/lib/python2.7/dist-packages')

# Faster-RCNN_TF
from networks.factory import get_network
from fast_rcnn.config import cfg
from fast_rcnn.test import im_detect
from fast_rcnn.nms_wrapper import nms
from utils.timer import Timer

# ModMan module & Pose estimation
from PoseEst.Function_Pose_v5 import *
import math

import tensorflow as tf
import cv2
import matplotlib.pyplot as plt
import numpy as np
from xml.etree.ElementTree import Element, dump, SubElement, ElementTree
from datetime import datetime

# socket
import socket
import errno

# debugging
import logging
logging.basicConfig(level = logging.INFO)

import select
from time import sleep

# realsense
# import pyrealsense as pyrs

CLASSES = yo_network_info.CLASSES
Candidate_CLASSES = yo_network_info.Candidate_CLASSES
NUM_CLASSES = yo_network_info.NUM_CLASSES

CONF_THRESH = yo_network_info.DETECTION_TH
NMS_THRESH = 0.3

global_do_send_msg = False
global_msg = {'msg': [], 'who': ''}

# Checks if a matrix is a valid rotation matrix.
def isRotationMatrix(R):
    Rt = np.transpose(R)
    shouldBeIdentity = np.dot(Rt, R)
    I = np.identity(3, dtype=R.dtype)
    n = np.linalg.norm(I - shouldBeIdentity)
    return n < 1e-6

# Calculates rotation matrix to euler angles
# The result is the same as MATLAB except the order
# of the euler angles ( x and z are swapped ).
def rotationMatrixToEulerAngles(R):
    assert (isRotationMatrix(R))

    sy = math.sqrt(R[0, 0] * R[0, 0] + R[1, 0] * R[1, 0])

    singular = sy < 1e-6

    if not singular:
        x = math.atan2(R[2, 1], R[2, 2])
        y = math.atan2(-R[2, 0], sy)
        z = math.atan2(R[1, 0], R[0, 0])
    else:
        x = math.atan2(-R[1, 2], R[1, 1])
        y = math.atan2(-R[2, 0], sy)
        z = 0

    return np.array([x, y, z])

def list_files(path, ext):
    filelist = []
    for name in os.listdir(path):
        if os.path.isfile(os.path.join(path, name)):
            if name.endswith(ext):
                filelist.append(name)
    return filelist

def vis_detections(im, class_name, dets, thresh=0.5):
    """Draw detected bounding boxes."""
    inds = np.where(dets[:, -1] >= thresh)[0]
    if len(inds) == 0:
        return

    im = im[:, :, (2, 1, 0)]
    fig, ax = plt.subplots(figsize=(12, 12))
    ax.imshow(im, aspect='equal')

    for i in inds:
        bbox = dets[i, :4]
        score = dets[i, -1]

        ax.add_patch(
            plt.Rectangle((bbox[0], bbox[1]),       # (x, y)
                          bbox[2] - bbox[0],        # width
                          bbox[3] - bbox[1],        # height
                          fill=False,
                          edgecolor='red', linewidth=3.5)
            )
        ax.text(bbox[0], bbox[1] - 2,
                '{:s} {:.3f}'.format(class_name, score),
                bbox=dict(facecolor='blue', alpha=0.5),
                fontsize=14, color='white')

        # print('{:s} {:.3f}'.format(class_name, score))

    ax.set_title(('{} detections with '
                  'p({} | box) >= {:.1f}').format(class_name, class_name,
                                                  thresh),
                  fontsize=14)
    plt.axis('off')
    plt.tight_layout()
    plt.draw()

def demo(sess, net, im):
    """Detect object classes in an image using pre-computed object proposals."""

    # Detect all object classes and regress object bounds
    timer = Timer()
    timer.tic()
    scores, boxes = im_detect(sess, net, im)
    timer.toc()
    # print ('Detection took {:.3f}s for '
    #        '{:d} object proposals').format(timer.total_time, boxes.shape[0])

    # Visualize detections for each class
    for cls_ind, cls in enumerate(CLASSES[1:]):
        cls_ind += 1 # because we skipped background
        cls_boxes = boxes[:, 4*cls_ind:4*(cls_ind + 1)]
        cls_scores = scores[:, cls_ind]
        dets = np.hstack((cls_boxes,
                          cls_scores[:, np.newaxis])).astype(np.float32)
        keep = nms(dets, NMS_THRESH)
        dets = dets[keep, :]
        vis_detections(im, cls, dets, thresh=CONF_THRESH)

def myimshow(strname, image):
    if image.shape[0] > 600:
        im_resized = cv2.resize(image, None, fx=0.5, fy=0.5)
        cv2.imshow(strname, im_resized)
        cv2.waitKey(10)
    else:
        cv2.imshow(strname, image)
        cv2.waitKey(10)

def demo_all(sess, snet, im_org, strEstPathname, extMat=None, FeatureDB=None, CoorDB=None, GeoDB=None):
    # scalefactor = 300. / float(min(im.shape[0], im.shape[1]))
    # tw = int(im.shape[1] * scalefactor)
    # th = int(im.shape[0] * scalefactor)
    # im = cv2.resize(im, (tw, th))

    ret_list_forKIST = []
    ret_list_BB = []

    # Detect all object classes and regress object bounds
    timer = Timer()
    timer.tic()
    scores, boxes = im_detect(sess, net, im_org)
    timer.toc()

    fontFace = cv2.FONT_HERSHEY_PLAIN
    fontScale = 2
    fontThickness = 2

    if len(strEstPathname) > 0:
        tag_anno = Element('annotation')

    im = im_org.copy()

    # Visualize detections for each class
    for cls_ind, class_name in enumerate(CLASSES[1:]):
        cls_ind += 1  # because we skipped background
        cls_boxes = boxes[:, 4 * cls_ind:4 * (cls_ind + 1)]
        cls_scores = scores[:, cls_ind]
        dets = np.hstack((cls_boxes,
                          cls_scores[:, np.newaxis])).astype(np.float32)
        keep = nms(dets, NMS_THRESH)
        dets = dets[keep, :]

        inds = np.where(dets[:, -1] >= CONF_THRESH)[0]

        if len(inds) > 0:
            for i in inds:
                bbox = dets[i, :4]      # [xmin, ymin, xmax, ymax]
                score = dets[i, -1]

                if class_name in Candidate_CLASSES:
                    if score > 0.8:
                        fontColor = (255,0,0)
                        cv2.rectangle(im, (int(bbox[0]), int(bbox[1])), (int(bbox[2]), int(bbox[3])), fontColor, fontThickness)
                    elif score > 0.6:
                        fontColor = (0, 255, 0)
                        cv2.rectangle(im, (int(bbox[0]), int(bbox[1])), (int(bbox[2]), int(bbox[3])), fontColor, fontThickness)
                    else:
                        fontColor = (255, 255, 255)
                        cv2.rectangle(im, (int(bbox[0]), int(bbox[1])), (int(bbox[2]), int(bbox[3])), fontColor, fontThickness)

                    cv2.putText(im, '{:s} {:.3f}'.format(class_name, score), (int(bbox[0]), int(bbox[1] - 2)), fontFace, fontScale, fontColor, thickness = fontThickness)

                    ret_list_BB.append({'bbox': bbox, 'score': score, 'name': class_name})

                    # print('{:s} {:.3f} {:d}'.format(class_name, score, cls_ind))

                    if extMat is not None:
                        FeatureDB2 = FeatureDB[Candidate_CLASSES.index(class_name)]
                        CoorDB2 = CoorDB[Candidate_CLASSES.index(class_name)]
                        GeoDB2 = GeoDB[Candidate_CLASSES.index(class_name)]

                        width = np.min([int(bbox[2]) - int(bbox[0]),10*int(bbox[0])])
                        height = np.min([int(bbox[3]) - int(bbox[1]),10*int(bbox[1])])

                        cropbox_ly = int(bbox[1]-height*0.1)
                        cropbox_ry = int(bbox[3]+height*0.1)
                        cropbox_lx = int(bbox[0]-width*0.1)
                        cropbox_rx = int(bbox[2]+width*0.1)

                        cropimg = im_org[cropbox_ly:cropbox_ry, cropbox_lx:cropbox_rx, :]

                        print('bbox:')
                        print(bbox)

                        init_coord = np.array([cropbox_lx,  cropbox_ly, 0])    # init_coord[x, y, -], lefttop_point
                        rmat, tvec = PoseEstimate(cropimg, FeatureDB2, CoorDB2, extMat, init_coord)

                        print(rmat)
                        print(tvec)
                        if rmat.sum() == 0 or np.isnan(rmat).any() or np.isnan(tvec).any() == True:
                            print('cannot find the pose information and fill with dummy values with all zeros')
                            # return for KIST
                            obj_info = {'object': class_name,
                                        'score': score,
                                        'RMat': rmat,
                                        'TVec': tvec * 0.0,
                                        'x_center': (bbox[0] + bbox[2]) / 2,
                                        'y_center': (bbox[1] + bbox[3]) / 2,
                                        'left': bbox[0],
                                        'top': bbox[1],
                                        'right': bbox[2],
                                        'bottom': bbox[3]
                                        }
                            ret_list_forKIST.append(obj_info)

                        else:
                            init_coord = np.array([0, 0, 0])
                            Result = cornerpointsTransform2(GeoDB2, rmat, tvec, extMat, init_coord)

                            # return for KIST
                            obj_info = {'object': class_name,
                                        'score': score,
                                        'RMat': rmat,
                                        'TVec': tvec,
                                        'x_center': (bbox[0] + bbox[2]) / 2,
                                        'y_center': (bbox[1] + bbox[3]) / 2,
                                        'left': bbox[0],
                                        'top': bbox[1],
                                        'right': bbox[2],
                                        'bottom': bbox[3]
                                        }
                            ret_list_forKIST.append(obj_info)

                            print('\tRot info: ')
                            print(rmat)
                            # print('\tTrn info:\n\t\tx: %d\n\t\ty: %d\n\t\tz: %d' % (tvec[1] * 0.1, -tvec[0] * 0.1, tvec[2] * 0.1))      # *0.1 --> mm to cm
                            print('\tTvec info: ')
                            print(tvec)
                            #print('\tTrn info:\n\t\tx: %d\n\t\ty: %d\n\t\tz: %d' % (tvec[1]/tvec[0], -tvec[0]//tvec[0], tvec[2]/tvec[0]))

                            # draw axis
                            drawLineIndeces = ((0, 1), (0, 2), (0, 3))
                            colorList = ((255, 0, 0), (0, 255, 0), (0, 0, 255))  # z, y, x
                            for (idxStart, idxEnd), color in zip(drawLineIndeces, colorList):
                                cv2.line(im, (int(Result[idxStart][0]), int(Result[idxStart][1])),
                                         (int(Result[idxEnd][0]), int(Result[idxEnd][1])), color, thickness=4)

                            # draw point
                            for ptDisp in Result:
                                cv2.circle(im, (int(ptDisp[0]), int(ptDisp[1])), 5, (255,255,255,0), -1)

                            # draw center position from a camera (mm -> cm by x 0.1)
                            cv2.putText(im, '(%d, %d, %d)'%(tvec[1] * 0.1, -tvec[0] * 0.1, tvec[2] * 0.1), (int(Result[0][0]), int(Result[0][1])), fontFace, fontScale, fontColor, thickness=fontThickness)

                    if len(strEstPathname) > 0:
                        tag_object = Element('object')
                        SubElement(tag_object, 'name').text = class_name
                        SubElement(tag_object, 'score').text = str(score)
                        tag_bndbox = Element('bndbox')
                        SubElement(tag_bndbox, 'xmin').text = str(int(bbox[0]))
                        SubElement(tag_bndbox, 'ymin').text = str(int(bbox[1]))
                        SubElement(tag_bndbox, 'xmax').text = str(int(bbox[2]))
                        SubElement(tag_bndbox, 'ymax').text = str(int(bbox[3]))
                        tag_anno.append(tag_object)
                        tag_object.append(tag_bndbox)

                        if extMat is not None and rmat.sum() != 0 and np.isnan(rmat).any() != True:
                            SubElement(tag_object, 'rotation_matrix').text = str(rmat);
                            SubElement(tag_object, 'traslation_vector').text = str(tvec);

                            xyz = rotationMatrixToEulerAngles(rmat)
                            SubElement(tag_object, 'EulerAngles').text = str(xyz)

    myimshow('display', im)

    if len(strEstPathname) > 0:
        cv2.imwrite(strEstPathname + '_est.jpg', im)
        ElementTree(tag_anno).write(strEstPathname)

    return im, ret_list_forKIST, ret_list_BB

def read_list_linebyline(fname):
    with open(fname) as fid:
        content = fid.readlines()
    content = [item.rstrip('\n') for item in content]

    return content

def write_list_linebyline(fname, thelist):
    fid = open(fname, 'w')

    for item in thelist:
        fid.write('%s\n' % (item))

    fid.close()

def get_curtime_as_string():
    timestamp = '%d_%d_%d_%d_%d_%d_%d' % (
        datetime.now().year,
        datetime.now().month,
        datetime.now().day,
        datetime.now().hour,
        datetime.now().minute,
        datetime.now().second,
        datetime.now().microsecond
    )

    return timestamp

gmp_ix, gmp_iy = -1, -1
gmp_event = -1
def get_mouse_position(event, x, y, flags, param):
    global gmp_ix, gmp_iy, gmp_event
    if event == cv2.EVENT_LBUTTONDOWN:
        gmp_ix, gmp_iy, gmp_event = x, y, event

def getCamIntParams(nameDevice):
    extMat = np.zeros((3, 3))
    if nameDevice == 'SKKU':
        print('Set SKKU Camera intrinsic param')
        # # value from SKKU
        # extMat[0, 0] = 2398.9479815687  # fy
        # extMat[1, 1] = 2399.8176947948  # fx
        # extMat[0, 2] = 492.5551416881  # py
        # extMat[1, 2] = 622.6012000501  # px
        # extMat[2, 2] = 1

        # # SKKU Camera - Calculate at KIST by myself - 1.5m distance btw camera and obj
        # extMat[0, 0] = 2532.719  # fy
        # extMat[1, 1] = 2534.086  # fx
        # extMat[0, 2] = 493.783  # py
        # extMat[1, 2] = 627.673  # px
        # extMat[2, 2] = 1

        # SKKU Camera - Calculate at KIST by myself - 0.7m distance btw camera and obj
        extMat[0, 0] = 2559.335  # fy
        extMat[1, 1] = 2559.955  # fx
        extMat[0, 2] = 467.309  # py
        extMat[1, 2] = 605.107  # px
        extMat[2, 2] = 1

    elif nameDevice == 'C920':
        # Logitech C920
        print('Set MS C920 Webcam intrinsic param')
        extMat[0, 0] = 603.257  # fy
        extMat[1, 1] = 603.475  # fx
        extMat[0, 2] = 237.936  # py
        extMat[1, 2] = 323.122  # px
        extMat[2, 2] = 1
    elif nameDevice == 'KinectV2':
        # KinectV2
        print('Set KINECTv2 intrinsic param')
        extMat[0, 0] = 1101.658925  # fy
        extMat[1, 1] = 1102.363184  # fx
        extMat[0, 2] = 513.583814  # py
        extMat[1, 2] = 937.202226  # px
        extMat[2, 2] = 1
    elif nameDevice == 'iphone6S':
        print('Set iPhone6S intrinsic param')
        extMat[0, 0] = 1122.555  # fy
        extMat[1, 1] = 1124.040  # fx
        extMat[0, 2] = 371.990  # py
        extMat[1, 2] = 661.668  # px
        extMat[2, 2] = 1
    elif nameDevice == 'iPadPro10.5':
        print('Set iPad intrinsic param')
        extMat[0, 0] = 2398.881  # fy
        extMat[1, 1] = 2400.970  # fx
        extMat[0, 2] = 815.355  # py
        extMat[1, 2] = 1096.002  # px
        extMat[2, 2] = 1
    elif nameDevice == 'SR300':
        print('Set SR300 intrinsic param')
        extMat[0, 0] = 620.330732  # fy
        extMat[1, 1] = 617.397217  # fx
        extMat[0, 2] = 241.907642  # py
        extMat[1, 2] = 321.188038  # px
        extMat[2, 2] = 1
    elif nameDevice == 'client':
        print('Get intrinsic param from the client PC')

    return extMat

import threading

class ThreadedClient(object):
    def __init__(self):
        self.owner = '_'
        self.extMat = np.zeros((3, 3))
        self.IMG_WIDTH = 0
        self.IMG_HEIGHT = 0
    def parsing_cmd(self, cmd_raw):
        ret = 0
        buf_size = 3    # 'MME'
        buf_size_intr_params = 0
        buf_size_PadInfo = 0
        buf_size_RTInfo = 0

        cmd = '%c%c'%(struct.unpack('c', cmd_raw[0])[0], struct.unpack('c', cmd_raw[1])[0])

        size_iPad_image = 2224 * 1668 * 3
        if USE_REDUCED_IPAD_RESOL == True:
            size_iPad_image = 640 * 480 * 3
        size_sr300_image = 640 * 480 * 3

        if len(cmd) == 2:
            if cmd == 'ec': # etri + calibration
                ret = 1
                buf_size_intr_params = 4 * 4
                buf_size_PadInfo = 4 * 12
                buf_size_RTInfo = 0
                buf_size = buf_size + size_iPad_image
            elif cmd == 'er':
                ret = 1
                buf_size_intr_params = 4 * 4
                buf_size_PadInfo = 4 * 12
                buf_size_RTInfo = 0
                buf_size = buf_size + size_iPad_image
            elif cmd == 'ed':
                ret = 1
                buf_size_PadInfo = 4 * 12
                buf_size_RTInfo = 4 * 12
                buf_size = buf_size + 0
            elif cmd == 'ks':
                ret = 1
                buf_size = buf_size + size_sr300_image

        return ret, cmd, buf_size_intr_params, buf_size_PadInfo, buf_size_RTInfo, buf_size

    def listenToClient(self, clientSocket, address_tuple):
        global global_do_send_msg
        global global_msg

        address = address_tuple[0]
        print('Thread start: %s'%address)

        if address == yo_network_info.KIST_STATIC_IP:
            self.owner = 'k'
            print('client (%s) was set as KIST'%address)
        elif address == yo_network_info.ETRI_STATIC_IP:
            self.owner = 'e'
            print('client (%s) was set as ETRI'%address)

        CONNECTED = True
        while CONNECTED == True:  # do image -> get info without interception.
            '''
            Receiving part
            '''
            ing_rcv = False
            len_rcv = 0
            len_rcv_info = 0

            do_capture_success = False

            # get a one image
            NET_BUFSIZE = 3 + 2     # MMS + cmd(2 bytes)
            cmd = '  '
            info_padPosition = np.zeros((12))
            info_RT = np.zeros((12))

            print('Server (%s, %c): connected' % (address, self.owner))

            while CONNECTED == True:
                try:
                    do_read = False

                    timeout_in_seconds = 0.5
                    clientSocket.setblocking(0)

                    while True:
                        try:
                            r, _, _ = select.select([clientSocket], [], [], timeout_in_seconds)
                            do_read = bool(r)
                        except socket.Timeouterror:
                            pass

                        if do_read == True:
                            print('Server (%s, %c): try to receive'%(address, self.owner))
                            data = clientSocket.recv(NET_BUFSIZE)

                            break
                        else:
                            # check there is a data to send
                            if global_do_send_msg == True:
                                if global_msg['who'] == self.owner:
                                    print('Server (%s, %c): This is my message'%(address, self.owner))

                                    # send all in one time
                                    # clientSocket.send(global_msg['msg'])
                                    # print('Server (%s, %c): send to the client'%(address, self.owner))

                                    LEN_CHUNK = yo_network_info.LEN_CHUNK

                                    len_send = 0
                                    for ipack in xrange(0, len(global_msg['msg']), LEN_CHUNK):
                                        len_send = len_send + clientSocket.send(global_msg['msg'][ipack: ipack + LEN_CHUNK])
                                        # print('Server (%s, %c): msg sent: len %d, %d' % (
                                        # address, self.owner, len(global_msg['msg']), len_send))
                                        sleep(0.004)

                                    print('Server (%s, %c): msg sent: len %d, %d' % (address, self.owner, len(global_msg['msg']), len_send))

                                    global_msg = {'msg': [], 'who': ''}
                                    global_do_send_msg = False

                except socket.error, socket.e:
                    if isinstance(socket.e.args, tuple):
                        print("Server (%s, %c): errno is"%(address, self.owner))
                        if socket.e.errno == errno.EPIPE:
                            # remote peer disconnected
                            print("Server (%s, %c): detected remote disconnect"%(address, self.owner))
                        elif socket.e.errno == errno.ECONNRESET:
                            print('Server (%s, %c): disconnected'%(address, self.owner))
                        else:
                            # determine and handle different error
                            print("socket error - type1", socket.e)
                            print('Server (%s, %c): maybe not ready to receive.: '%(address, self.owner))
                    else:
                        print "socket error -type2", socket.e

                    CONNECTED = False
                    clientSocket.close()

                if len(data) == 0:
                    print('Server (%s, %c): we guess the client is disconnected.'%(address, self.owner))
                    CONNECTED = False
                    clientSocket.close()

                if ing_rcv == False and 'MMS' in data:  # first receive - resolve - command, params, then image
                    ing_rcv = True

                    index_start = data.index('MMS')

                    suc_cmd, cmd, buf_size_params, buf_size_PadInfo, buf_size_RTinfo, buf_size = self.parsing_cmd(data[index_start + 3:])

                    if suc_cmd == 0:
                        print('Server (%s, %c): client sends incorrect command'%(address, self.owner))
                        CONNECTED = False
                        clientSocket.close()
                    else:
                        self.owner = cmd[0]
                        NET_BUFSIZE = buf_size
                        print('Server (%s, %c): client sends command [%s]' % (address, self.owner, cmd))

                    if cmd[0] == 'e' and buf_size_params > 0: # update intrinsic params.
                        # get more 4 * 4 bytes and decode them
                        data_params = clientSocket.recv(4 * 4)

                        fx = struct.unpack('f', data_params[:4])[0]  # fx
                        fy = struct.unpack('f', data_params[4:8])[0]  # fy
                        px = struct.unpack('f', data_params[8:12])[0]  # px
                        py = struct.unpack('f', data_params[12:])[0]  # py

                        if USE_REDUCED_IPAD_RESOL == True:
                            fx = fx * 0.2878
                            fy = fy * 0.2878
                            px = px * 0.2878
                            py = py * 0.2878

                        print('Set intrinsic param: fxfy (%f, %f), pxpy (%f, %f)' % (fx, fy, px, py))

                        self.extMat[0, 0] = fy  # fy
                        self.extMat[1, 1] = fx  # fx
                        self.extMat[0, 2] = py  # py
                        self.extMat[1, 2] = px  # px
                        self.extMat[2, 2] = 1

                    elif cmd == 'ks':
                        self.extMat = getCamIntParams('SR300')

                    if cmd[0] == 'e':
                        # iPad original resolution
                        # if you change below code related to resolution,
                        # you also change the code about buf_size in self.parsing_cmd().
                        self.IMG_WIDTH = 2224
                        self.IMG_HEIGHT = 1668
                        if USE_REDUCED_IPAD_RESOL == True:
                            # iPad reduced resolution
                            self.IMG_WIDTH = 640
                            self.IMG_HEIGHT = 480
                            print('iPad reduced resolution (640x480) applied')
                    elif cmd[0] == 'k':
                        self.IMG_WIDTH = 640
                        self.IMG_HEIGHT = 480

                    if buf_size_PadInfo > 0:
                        # get more 12 * 4 bytes and decode them
                        data_dummy = clientSocket.recv(12 * 4)

                        for i in range(12):
                            info_padPosition[i] = struct.unpack('f', data_dummy[i * 4:i * 4 + 4])[0]

                        print('Pad Pos. : ', info_padPosition)

                    if buf_size_RTinfo > 0:
                        # get more 12 * 4 bytes and decode them
                        data_dummy = clientSocket.recv(12 * 4)

                        for i in range(12):
                            info_RT[i] = struct.unpack('f', data_dummy[i*4:i*4+4])[0]  # fx

                        print('Dst Pos.: ', info_RT)

                    fid_bin = open('./skku_img.bin', 'wb')

                    len_rcv = 0

                elif ing_rcv == True:  # receive again.
                    # debug for message ending.
                    #print('MME' in data)
                    #if 'MME' in data:
                    #    print(data.index('MME') + len_rcv == buf_size-3)
                    if ('MME' in data) and (data.index('MME') + len_rcv == buf_size-3):
                        ing_rcv = False
                        index_end = data.index('MME')
                        fid_bin.write(data[:index_end])
                        fid_bin.close()

                        len_rcv = len_rcv + len(data[:index_end])
                        print('final: %d == %d' % (buf_size-3, len_rcv))

                        len_rcv = 0

                        self.img = np.fromfile('./skku_img.bin', dtype='uint8')

                        if len(self.img) != self.IMG_HEIGHT * self.IMG_WIDTH * 3:
                            print('captured image size is not predefined size')
                            do_capture_success = False
                        else:
                            self.img = self.img.reshape(self.IMG_HEIGHT, self.IMG_WIDTH, 3)
                            do_capture_success = True

                        break
                    else:
                        fid_bin.write(data)
                        len_rcv = len_rcv + len(data)

                        # print(data)
                # print('intermediate: %d (exclude MMS, cmd, int/ext params=69bytes)' % (len_rcv))

            print('Server (%s, %c): received data completely'%(address, self.owner))

            # do ETRI job
            if do_capture_success is True:
                myimshow('display', self.img)

                if DO_WRITE_RESULT_IMG == True:
                    timestamp = get_curtime_as_string()
                    cv2.imwrite('../debug_dropbox_upload/%s.png' % timestamp, self.img)
                    print('captured %s image', timestamp)

                cv2.imwrite('./debug_img.png', self.img)

                if cmd == 'er' or cmd == 'ks':
                    print('Server (%s, %c): process image...' % (address, self.owner))
                    if USE_POSEESTIMATE is True:
                        result_img, list_objs_forAR, _ = demo_all(sess, net, self.img, '', self.extMat, FeatureDB, CoorDB, GeoDB)
                    else:
                        result_img, list_objs_forAR, _ = demo_all(sess, net, self.img, '')

                if DO_WRITE_RESULT_IMG == True:
                    cv2.imwrite('../debug_dropbox_upload/%s_est.png' % timestamp, result_img)

                    fid_info = open('../debug_dropbox_upload/%s_est.txt' % timestamp, 'w')
                    print >> fid_info, list_objs_forAR
                    fid_info.close()

                    print('save dummy image as a file')
            else:
                print('no frame\n')

            # send to kist
            # 4 byte: num of object (int)
            # 1 byte: obj ID (char)
            # 4 * (9 + 3) = 48 bytes = rot + trans mat
            # 4 * 2 = 8 bytes = x, y

            '''
            Sending part
            '''
            if CONNECTED == True:
                if cmd == 'ec':
                    msg = 'MMS'
                    msg = addCommand(msg, cmd)
                    msg = addPadRTInfo(msg, info_padPosition, floatType='double')
                    msg = addImage(msg, self.img)
                    msg = msg + 'MME'

                    global_msg['msg'] = copy.copy(msg)
                    global_msg['who'] = 'k'
                    global_do_send_msg = True
                    print('Server (%s, %c): set global msg _%s_'%(address, self.owner, msg[:5]))
                elif cmd == 'ed':
                    msg = 'MMS'
                    msg = addCommand(msg, cmd)
                    msg = addPadRTInfo(msg, info_padPosition, floatType='double')
                    msg = addPadRTInfo(msg, info_RT, floatType='double')
                    msg = msg + 'MME'

                    global_msg['msg'] = copy.copy(msg)
                    global_msg['who'] = 'k'
                    global_do_send_msg = True
                    print('Server (%s, %c): set global msg _%s_' % (address, self.owner, msg))
                    print('Server (%s, %c): (cont.) ' % (address, self.owner), msg[5:5 + 24 * 8])
                else:
                    # message to loopback.
                    msg = 'MMS'
                    msg = addCommand(msg, cmd)

                    if cmd == 'er':
                        msg = addListObjs(msg, list_objs_forAR, floatType='float')
                    elif cmd =='ks':
                        msg = addListObjs(msg, list_objs_forAR, floatType='double')
                    msg = msg + 'MME'

                    # for list_objs
                    clientSocket.send(msg)
                    print('Server (%s, %c): send _%s_ to the client' % (address, self.owner, msg))

                    if cmd == 'er':
                        msg = 'MMS'
                        msg = addCommand(msg, cmd)
                        # msg = addImage(msg, self.img)
                        msg = addPadRTInfo(msg, info_padPosition, floatType='double')
                        msg = addListObjs(msg, list_objs_forAR, floatType='double')
                        msg = msg + 'MME'

                        global_msg['msg'] = copy.copy(msg)
                        global_msg['who'] = 'k'
                        global_do_send_msg = True

                        print('Server (%s, %c): set global msg '%(address, self.owner))

# https://stackoverflow.com/questions/23828264/how-to-make-a-simple-multithreaded-socket-server-in-python-that-remembers-client
class ThreadedServer(object):
    def __init__(self, host, port):
        self.host = host
        self.port = port
        self.sock = socket.socket(socket.AF_INET, socket.SOCK_STREAM)
        self.sock.setsockopt(socket.SOL_SOCKET, socket.SO_REUSEADDR, 1)
        self.sock.bind((self.host, self.port))

    def listen(self):
        self.sock.listen(5)     # accept max. 5 clients
        while True:
            client, address = self.sock.accept()
            client.settimeout(30*60)               # if there is no reaction for 30*60s, it will disconnect and close.

            clientInst = ThreadedClient()

            threading.Thread(target = clientInst.listenToClient,\
                             args = (client, address)\
                             ).start()

def addPadRTInfo(msg, padinfo, floatType):
    if floatType == 'float':
        for i in range(12):
            msg = msg + struct.pack('f', padinfo[i])
    elif floatType == 'double':
        for i in range(12):
            msg = msg + struct.pack('d', padinfo[i])
    else:
        raise AssertionError('type is incorrect!!!')

    print(padinfo)

    return msg

def addImage(msg, image):
    # Image
    for i in range(image.shape[0]):
        for j in range(image.shape[1]):
            for c in range(3):
                str_c = '%c' % image[i, j, c]
                # msg = msg + struct.pack('c', str_c.encode('ascii'))
                msg = msg + struct.pack('c', str_c)

    return msg

def addListObjs(msg, listObjs, floatType):
    # num_of_objects: int x 1
    msg = msg + struct.pack('i', len(listObjs))  # int
    print('num of objs: %d' % len(listObjs))

    for obj_AR in listObjs:
        msg = msg + struct.pack('c', obj_AR['object'][0])
        print('obj name: %c' % obj_AR['object'][0])

        if floatType == 'float':
            for j_RMat in range(0, 3):
                for i_RMat in range(0, 3):
                    msg = msg + struct.pack('f', obj_AR['RMat'][i_RMat][j_RMat])
            print('RMat:')
            print(obj_AR['RMat'])

            for j_TVec in range(0, 3):
                msg = msg + struct.pack('f', obj_AR['TVec'][j_TVec][0])
            print('TVec:')
            print(obj_AR['TVec'])
        elif floatType == 'double':
            for j_RMat in range(0, 3):
                for i_RMat in range(0, 3):
                    msg = msg + struct.pack('d', obj_AR['RMat'][i_RMat][j_RMat])
            print('RMat:')
            print(obj_AR['RMat'])

            for j_TVec in range(0, 3):
                msg = msg + struct.pack('d', obj_AR['TVec'][j_TVec][0])
            print('TVec:')
            print(obj_AR['TVec'])
        else:
            raise AssertionError('type is incorrect!!!')

        msg = msg + struct.pack('i', int(obj_AR['left'])) \
              + struct.pack('i', int(obj_AR['top'])) \
              + struct.pack('i', int(obj_AR['right'])) \
              + struct.pack('i', int(obj_AR['bottom']))

        print('left: %d' % (int(obj_AR['left'])))
        print('top: %d' % (int(obj_AR['top'])))
        print('right: %d' % (int(obj_AR['right'])))
        print('bottom: %d' % (int(obj_AR['bottom'])))

    return msg

def addCommand(msg, cmd):
    # command: 2bytes
    msg = msg + struct.pack('c', cmd[0])
    msg = msg + struct.pack('c', cmd[1])

    return msg

def setTestCFG(cfg):
    cfg.NCLASSES = NUM_CLASSES
    cfg.TEST.HAS_RPN = True  # Use RPN for proposals
    cfg.TEST.BBOX_REG = True  # Use BBox regressor
    # # Scales to use during testing (can list multiple scales)
    # # Each scale is the pixel size of an image's shortest side
    cfg.TEST.SCALES = (1600,)  # only one for test  ,480, 576, 688, 864, 1200
    cfg.TEST.RPN_NMS_THRESH = 0.7

    ## Number of top scoring boxes to keep before apply NMS to RPN proposals
    # cfg.TEST.RPN_PRE_NMS_TOP_N = 6000
    cfg.TEST.RPN_PRE_NMS_TOP_N = 12000
    ## Number of top scoring boxes to keep after applying NMS to RPN proposals
    # cfg.TEST.RPN_POST_NMS_TOP_N = 300
    cfg.TEST.RPN_POST_NMS_TOP_N = 2000
    # cfg.TEST.MAX_SIZE = 2000

    return cfg

# variable to communicate between threads.
if __name__ == '__main__':
    os.environ["CUDA_DEVICE_ORDER"] = "PCI_BUS_ID"
    os.environ["CUDA_VISIBLE_DEVICES"] = "0"

    '''
    Settings
    '''
    USE_POSEESTIMATE = True
    USE_REDUCED_IPAD_RESOL = True       # True: 640x480, False: 2224x1668

    Svr_IP = yo_network_info.SERVER_IP
    Svr_PORT = yo_network_info.SERVER_PORT

    DO_WRITE_RESULT_AVI = False
    name_output_avi = 'output.avi'

    DO_WRITE_RESULT_IMG = False     # IMG name will be the current time.

    # for test
    cfg = setTestCFG(cfg)

    avgwindow = 0   # parameter to show pose estimation in stable

    '''
    Write Result
    '''
    if DO_WRITE_RESULT_AVI == True:
        # Define the codec and create VideoWriter object
        frame_width = 640
        frame_height = 480
        outavi = cv2.VideoWriter(name_output_avi, cv2.VideoWriter_fourcc('M','J','P','G'), 15.0, (frame_width, frame_height))

    # this is for pose estimation
    objNum = len(Candidate_CLASSES)
    DotLog = np.zeros((4, 3, avgwindow, objNum))
    current = np.zeros((objNum), dtype="int")
    tvecLog = np.zeros((3, 1, avgwindow, objNum))

    FeatureDB = []
    CoorDB = []
    keypointDB = []
    GeoDB = []
    if USE_POSEESTIMATE == True:
        gap = 4

        for ith, obj in enumerate(Candidate_CLASSES):
            temp_FeatureDB, temp_CoorDB, temp_keypointDB = ReadDB(obj)
            FeatureDB.append(temp_FeatureDB[::gap,:])
            CoorDB.append(temp_CoorDB[::gap,:])
            # keypointDB.extend(temp_keypointDB)

            # read plotting information of the object
            # DBv1
            filelist2 = os.listdir(os.path.join(basePath, obj.lower(), 'geom2'))
            filelist2 = np.sort(filelist2)
            strTRSet = os.path.join(os.path.join(basePath, obj.lower(), 'geom2', filelist2[0]))
            # DBv2
            # strTRSet = os.path.join(basePath, '%s-rotate-geom_00_00_50_400mat.mat'%obj)
            ftr = h5py.File(strTRSet, 'r')
            GeoDB.append(np.transpose(np.array(ftr['img']), [2, 1, 0]))

    tfArch = 'VGGnetslsv1_test'  # prototxt
    tfmodel = '../models/VGGnet_fast_rcnn_iter_70000.ckpt'

    # init session
    sess = tf.Session(config=tf.ConfigProto(allow_soft_placement=True, device_count = {'GPU': 1}))
    tf.device('')

    # load network
    net = get_network(tfArch)
    # load model
    print ('Loading network {:s}... '.format(tfArch))
    saver = tf.train.Saver()
    saver.restore(sess, tfmodel)
    print (' done.')

    # Warmup on a dummy image
    im = 128 * np.ones((300, 300, 3), dtype=np.uint8)
    for i in xrange(2):
        _, _ = im_detect(sess, net, im)

    # working as a server
    '''
    Server Info
    '''
    print('Server: waiting of client connection')
    ThreadedServer(Svr_IP, Svr_PORT).listen()


text:
import asyncio
import os
from pathlib import Path

import numpy as np
import pytest
import time

from _pytest.monkeypatch import MonkeyPatch
from _pytest.tmpdir import TempdirFactory
from unittest.mock import patch, Mock

from rasa.core.agent import Agent
from rasa.core.channels import UserMessage
from rasa.core.constants import INTENT_MESSAGE_PREFIX, DEFAULT_LOCK_LIFETIME
from rasa.core.lock import TicketLock
from rasa.core.lock_store import InMemoryLockStore, LockError, LockStore, RedisLockStore

class FakeRedisLockStore(RedisLockStore):
    """Fake `RedisLockStore` using `fakeredis` library."""

    # noinspection PyMissingConstructor
    def __init__(self):
        import fakeredis

        self.red = fakeredis.FakeStrictRedis()

        # added in redis==3.3.0, but not yet in fakeredis
        self.red.connection_pool.connection_class.health_check_interval = 0

def test_issue_ticket():
    lock = TicketLock("random id 0")

    # no lock issued
    assert lock.last_issued == -1
    assert lock.now_serving == 0

    # no one is waiting
    assert not lock.is_someone_waiting()

    # issue ticket
    ticket = lock.issue_ticket(1)
    assert ticket == 0
    assert lock.last_issued == 0
    assert lock.now_serving == 0

    # someone is waiting
    assert lock.is_someone_waiting()

def test_remove_expired_tickets():
    lock = TicketLock("random id 1")

    # issue one long- and one short-lived ticket
    _ = list(map(lock.issue_ticket, [k for k in [0.01, 10]]))

    # both tickets are there
    assert len(lock.tickets) == 2

    # sleep and only one ticket should be left
    time.sleep(0.02)
    lock.remove_expired_tickets()
    assert len(lock.tickets) == 1

@pytest.mark.parametrize("lock_store", [InMemoryLockStore(), FakeRedisLockStore()])
def test_create_lock_store(lock_store: LockStore):
    conversation_id = "my id 0"

    # create and lock
    lock = lock_store.create_lock(conversation_id)
    lock_store.save_lock(lock)
    lock = lock_store.get_lock(conversation_id)
    assert lock
    assert lock.conversation_id == conversation_id

@pytest.mark.parametrize("lock_store", [InMemoryLockStore(), FakeRedisLockStore()])
def test_serve_ticket(lock_store: LockStore):
    conversation_id = "my id 1"

    lock = lock_store.create_lock(conversation_id)
    lock_store.save_lock(lock)

    # issue ticket with long lifetime
    ticket_0 = lock_store.issue_ticket(conversation_id, 10)
    assert ticket_0 == 0

    lock = lock_store.get_lock(conversation_id)
    assert lock.last_issued == ticket_0
    assert lock.now_serving == ticket_0
    assert lock.is_someone_waiting()

    # issue another ticket
    ticket_1 = lock_store.issue_ticket(conversation_id, 10)

    # finish serving ticket_0
    lock_store.finish_serving(conversation_id, ticket_0)

    lock = lock_store.get_lock(conversation_id)

    assert lock.last_issued == ticket_1
    assert lock.now_serving == ticket_1
    assert lock.is_someone_waiting()

    # serve second ticket and no one should be waiting
    lock_store.finish_serving(conversation_id, ticket_1)

    lock = lock_store.get_lock(conversation_id)
    assert not lock.is_someone_waiting()

# noinspection PyProtectedMember
@pytest.mark.parametrize("lock_store", [InMemoryLockStore(), FakeRedisLockStore()])
def test_lock_expiration(lock_store: LockStore):
    conversation_id = "my id 2"
    lock = lock_store.create_lock(conversation_id)
    lock_store.save_lock(lock)

    # issue ticket with long lifetime
    ticket = lock.issue_ticket(10)
    assert ticket == 0
    assert not lock._ticket_for_ticket_number(ticket).has_expired()

    # issue ticket with short lifetime
    ticket = lock.issue_ticket(0.00001)
    time.sleep(0.00002)
    assert ticket == 1
    assert lock._ticket_for_ticket_number(ticket) is None

    # newly assigned ticket should get number 1 again
    assert lock.issue_ticket(10) == 1

async def test_multiple_conversation_ids(default_agent: Agent):
    text = INTENT_MESSAGE_PREFIX + 'greet{"name":"Rasa"}'

    conversation_ids = [f"conversation {i}" for i in range(2)]

    # ensure conversations are processed in order
    tasks = [default_agent.handle_text(text, sender_id=_id) for _id in conversation_ids]
    results = await asyncio.gather(*tasks)

    assert results
    processed_ids = [result[0]["recipient_id"] for result in results]
    assert processed_ids == conversation_ids

async def test_message_order(tmp_path: Path, default_agent: Agent):
    start_time = time.time()
    n_messages = 10
    lock_wait = 0.1

    # let's write the incoming order of messages and the order of results to temp files
    results_file = tmp_path / "results_file"
    incoming_order_file = tmp_path / "incoming_order_file"

    # We need to mock `Agent.handle_message()` so we can introduce an
    # artificial holdup (`wait_time_in_seconds`). In the mocked method, we'll
    # record messages as they come and and as they're processed in files so we
    # can check the order later on. We don't need the return value of this method so
    # we'll just return None.
    async def mocked_handle_message(self, message: UserMessage, wait: float) -> None:
        # write incoming message to file
        with open(str(incoming_order_file), "a+") as f_0:
            f_0.write(message.text + "\n")

        async with self.lock_store.lock(
            message.sender_id, wait_time_in_seconds=lock_wait
        ):
            # hold up the message processing after the lock has been acquired
            await asyncio.sleep(wait)

            # write message to file as it's processed
            with open(str(results_file), "a+") as f_1:
                f_1.write(message.text + "\n")

            return None

    # We'll send n_messages from the same sender_id with different blocking times
    # after the lock has been acquired.
    # We have to ensure that the messages are processed in the right order.
    with patch.object(Agent, "handle_message", mocked_handle_message):
        # use decreasing wait times so that every message after the first one
        # does not acquire its lock immediately
        wait_times = np.linspace(0.1, 0.05, n_messages)
        tasks = [
            default_agent.handle_message(
                UserMessage(f"sender {i}", sender_id="some id"), wait=k
            )
            for i, k in enumerate(wait_times)
        ]

        # execute futures
        await asyncio.gather(*(asyncio.ensure_future(t) for t in tasks))

        expected_order = [f"sender {i}" for i in range(len(wait_times))]

        # ensure order of incoming messages is as expected
        with open(str(incoming_order_file)) as f:
            incoming_order = [line for line in f.read().split("\n") if line]
            assert incoming_order == expected_order

        # ensure results are processed in expected order
        with open(str(results_file)) as f:
            results_order = [line for line in f.read().split("\n") if line]
            assert results_order == expected_order

        # Every message after the first one will wait `lock_wait` seconds to acquire its
        # lock (`wait_time_in_seconds` kwarg in `lock_store.lock()`).
        # Let's make sure that this is not blocking and test that total test
        # execution time is less than  the sum of all wait times plus
        # (n_messages - 1) * lock_wait
        time_limit = np.sum(wait_times[1:])
        time_limit += (n_messages - 1) * lock_wait
        assert time.time() - start_time < time_limit

async def test_lock_error(default_agent: Agent):
    lock_lifetime = 0.01
    wait_time_in_seconds = 0.01
    holdup = 0.1

    # Mock message handler again to add a wait time holding up the lock
    # after it's been acquired
    async def mocked_handle_message(self, message: UserMessage) -> None:
        async with self.lock_store.lock(
            message.sender_id,
            wait_time_in_seconds=wait_time_in_seconds,
            lock_lifetime=lock_lifetime,
        ):
            # hold up the message processing after the lock has been acquired
            await asyncio.sleep(holdup)

        return None

    with patch.object(Agent, "handle_message", mocked_handle_message):
        # first message blocks the lock for `holdup`,
        # meaning the second message will not be able to acquire a lock
        tasks = [
            default_agent.handle_message(
                UserMessage(f"sender {i}", sender_id="some id")
            )
            for i in range(2)
        ]

        with pytest.raises(LockError):
            await asyncio.gather(*(asyncio.ensure_future(t) for t in tasks))

async def test_lock_lifetime_environment_variable(monkeypatch: MonkeyPatch):
    import rasa.core.lock_store

    # by default lock lifetime is `DEFAULT_LOCK_LIFETIME`
    assert rasa.core.lock_store._get_lock_lifetime() == DEFAULT_LOCK_LIFETIME

    # set new lock lifetime as environment variable
    new_lock_lifetime = 123
    monkeypatch.setenv("TICKET_LOCK_LIFETIME", str(new_lock_lifetime))

    assert rasa.core.lock_store._get_lock_lifetime() == new_lock_lifetime

async def test_redis_lock_store_timeout(monkeypatch: MonkeyPatch):
    import redis.exceptions

    lock_store = FakeRedisLockStore()
    monkeypatch.setattr(
        lock_store,
        lock_store.get_or_create_lock.__name__,
        Mock(side_effect=redis.exceptions.TimeoutError),
    )

    with pytest.raises(LockError):
        async with lock_store.lock("some sender"):
            pass

ext: py

reason: import numpy


text:
# Importando pacotes
import numpy as np
import math as mt
import pygame, sys

# Inicializando pygame
pygame.init()
pygame.font.init()
font = pygame.font.SysFont('Arial', 15)
pygame.display.set_caption('Simulação de Pêndulo Simples')
screen = pygame.display.set_mode((500,500))

# Definição dos valores iniciais das variáveis
g = 9.81 # gravidade
l = 100 # comprimento da corda (em cm)
a = 100 # coeficiente da amplitude

# Definição da função slider para controle das variáveis
def slider(xmin,xmax,yposi,yposf,button,b):
    if button[0]!=0:
        pos = pygame.mouse.get_pos()
        xc = pos[0]
        yc = pos[1]
        if (xc < b-xmin) or (xc >b+xmax):
            pass
        if (yc>yposi and yc<yposf):
            if xc<xmin:
                xc = xmin
            if xc>xmax:
                xc = xmax
            return(xc)
        else:
            return(b)
    if button[0] == 0:
        xc = b
        return(xc)

# loop principal do pygame
while True:

    A = (a/100)*mt.pi # Amplitude
    screen.fill((255,255,255))
    amplitude = font.render(('Amplitude: {0:.2f}'.format(A)), True, (13,134,0))
    screen.blit(amplitude, (270,15))
    comprimento = font.render(('Comprimento: {} cm'.format(l)), True, (255,0,0))
    screen.blit(comprimento, (270,50))
    button = pygame.mouse.get_pressed()

# Sliders para L
    lpos = slider(1,230,50,75,button,l)
    l = lpos
    pygame.draw.line(screen, (235,235,235), (0,62), (230, 62))
    pygame.draw.rect(screen, (255,0,0),(l,50,8,25))

# Slider para A
    Apos = slider(0,100,10,35,button,a)
    a = Apos
    pygame.draw.line(screen, (235,235,235), (0,22), (100, 22))
    pygame.draw.rect(screen,(13,134,0),(a,10,8,25))

# Controles
    # Quit
    for event in pygame.event.get():
        if event.type == pygame.QUIT:
            pygame.quit()
            sys.exit(0)

    # Tempo e funções de posição da massa
    time = 0.02*pygame.time.get_ticks()
    theta = A*np.sin(mt.sqrt(g/l)*time)
    x = l*np.cos(theta)
    y = l*np.sin(theta)
    # Desenhando o pêndulo e atualizando o frame
    pygame.draw.line(screen, (0,0,0), (250,250), (250+int(y), 250+int(x)))
    pygame.draw.circle(screen, (0,0,255), (250+int(y),250+int(x)), 10)
    pygame.display.flip()

max_stars_repo_name: aluizgc/sim-pendsimples

max_issues_repo_name: aluizgc/simple-pendulum-sim

max_forks_repo_name: aluizgc/simple-pendulum-sim


text:
from pathlib import Path
from plyfile import PlyData
data = PlyData.read("pclouds/data/dense/fused.ply")
data.elements
data.x
data.vertex
data[vertex]
data['vertex']
data['vertex']['x']
data['vertex']['x'].shape
data['vertex']['y'].shape
xyz = np.vstack((data['vertex']['x'], data['vertex']['y'], data['vertex']['z']))
import numpy as np
xyz = np.vstack((data['vertex']['x'], data['vertex']['y'], data['vertex']['z']))
xyz.shaper
xyz.shape
xyz = np.vstack((data['vertex']['x'], data['vertex']['y'], data['vertex']['z'])).T
xyz.shape
xyz = np.hstack((data['vertex']['x'], data['vertex']['y'], data['vertex']['z'])).shape
xyz = np.hstack((data['vertex']['x'], data['vertex']['y'], data['vertex']['z'])).shape
xyz = np.hstack((data['vertex']['x'], data['vertex']['y'], data['vertex']['z']))
xyz.shape
xyz = np.dstack((data['vertex']['x'], data['vertex']['y'], data['vertex']['z']))
xyz.shape
xyz = np.vstack((data['vertex']['x'], data['vertex']['y'], data['vertex']['z'])).T
xyz.shape
normal = np.vstack((data['vertex']['nx'], ata['vertex']['ny'], ata['vertex']['nz'])).T
normal = np.vstack((data['vertex']['nx'], data['vertex']['ny'], data['vertex']['nz'])).T
normal.shape
np.savetxt?
q
np.savetxt("pclouds/data/dense/fused.xyz", xyz)
np.savetxt("pclouds/data/dense/fused.normals", normals)
np.savetxt("pclouds/data/dense/fused.normals", normal)
%hist -f ply2txt.pt

max_stars_repo_name: hexfaker/pcpnet

max_issues_repo_name: hexfaker/pcpnet

max_forks_repo_name: hexfaker/pcpnet


text:
# -*- coding: utf-8 -*-
"""
Created on Mon Feb  3 11:13:26 2020

@author: Levillain
"""

#!/usr/bin/env python
# coding: utf-8

# In[17]:

import os.path
import os
import numpy as np
import pydicom
import matplotlib.pyplot as plt
import csv
import math

from tkinter import filedialog

#select the directoy and import the MAA and Y90 RTdose files

dirname = filedialog.askdirectory(initialdir=os.getcwd(),title='Please select directory of the pacient')
print (dirname)
RefDsMAA = filedialog.askopenfilename(initialdir = "dirname",title = "Select MAA DICOM File")
RefDsMAA = pydicom.read_file(RefDsMAA)
imgMAA=RefDsMAA.pixel_array
imgMAA.setflags(write=1)
ScaleMAA=RefDsMAA.DoseGridScaling

DsY90 = filedialog.askopenfilename(initialdir = "dirname",title = "Select Y90 DICOM File")
DsY90 = pydicom.read_file(DsY90)
imgY90=DsY90.pixel_array
imgY90.setflags(write=1)
ScaleY90=DsY90.DoseGridScaling

#plot the same slice for the MAA and Y90
f, (ax1, ax2) = plt.subplots(1, 2, sharey=True)
ax1.imshow(imgMAA[45,:,:],cmap='binary')
ax1.set_title('MAA Dose')
ax2.imshow(imgY90[45,:,:],cmap='binary')
ax2.set_title('Y90 Dose')

#open the directory of the masks and print the number of files to verify if the number of lesions is correct
dirMask = filedialog.askdirectory(title='Please select dir of the masks')
nr_masks = sum(len(files) for _, _, files in os.walk(dirMask))
print("nr files:", nr_masks)

# list of the mask_names from dir without the extension .dcm
l=os.listdir(dirMask)
mask_names=[mask.split('.')[0] for mask in os.listdir(dirMask)]
print(mask_names)

#create a list with all the dicom masks, this will be used to for loop through each mask for the QVH calc
masks = []
for i in range(nr_masks):
    masks.append(pydicom.read_file(os.path.normpath(os.path.join(dirMask,l[i]))))
print('shape:',masks[0].pixel_array.shape)

class MaskResults:

    def __init__(self, imgMAA, ScaleMAA, imgY90, ScaleY90, mask, mask_name): #creating a class using as arguments all the paramaters requiered to calculate the QVH
        self.mask = mask
        self.mask_name = mask_name
        self.Qarraysorted, self.sumWF = self.qvh_calculation (imgMAA, ScaleMAA, imgY90, ScaleY90, mask) #the first and second outputs of the qvh_calculation function are saved
        self.x = [q[0] for q in self.Qarraysorted]
        self.y = [100]
        for i in range (1,len(self.Qarraysorted)):       # %volume of a voxel that is ponderated by the WF
            self.y.append(self.y[i-1]-self.Qarraysorted[i][1]/self.sumWF*100)

#to extand the lenght of each row based on the largest lenth of data, necessary to be able to export all the data in the same csv file
    def extend_xy(self, length):
        x_len = len(self.x)
        num_empty_needed = length - x_len
        self.x.extend(["" for i in range(num_empty_needed)])#create an emplty string and add it to the list
        self.y.extend(["" for i in range(num_empty_needed)])

    def qvh_calculation(self,imgMAA, ScaleMAA, imgY90, ScaleY90, mask):
 # create a initial matrix for each factor
        wf_MAA=np.zeros_like(np.copy(imgMAA),dtype=float)
        wf_Y90=np.zeros_like(np.copy(imgMAA),dtype=float)
        wf_QVH=np.zeros_like(np.copy(imgMAA),dtype=float)
        Qarray = list()

        sumWF = 0
        qFWeightedSum = 0

# loop through each voxel of the image
        for i in (range((imgMAA.shape[0]))):
            for j in (range((imgMAA.shape[1]))):
                for k in (range ((imgMAA.shape[2]))):
                    if mask.pixel_array[i][j][k] == 65535:  #Calculate the quality division factor (qF) for each voxel inside the preselected mask
                        Dt = 0.1 if (imgMAA[i,j,k]*ScaleMAA)==0 else imgMAA[i,j,k]*ScaleMAA
                        Ds = 0.1 if (imgY90[i,j,k]*ScaleY90)==0 else imgY90[i,j,k]*ScaleY90

                        qF = math.log10(Dt/Ds)

                    #Calculate the weighting factor for MAA(wf_MAA) based on the predefined dose thresholds
                        if (imgMAA[i,j,k]*ScaleMAA)<=10 or (imgMAA[i,j,k]*ScaleMAA)>=200:
                            wf_MAA=0
                        elif (imgMAA[i,j,k]*ScaleMAA)>10 and (imgMAA[i,j,k]*ScaleMAA)<40:
                            wf_MAA=((1/30)*(imgMAA[i,j,k]*ScaleMAA))-(1/3)
                        elif (imgMAA[i,j,k]*ScaleMAA)>120 and (imgMAA[i,j,k]*ScaleMAA)<200:
                            wf_MAA=((-2/160)*(imgMAA[i,j,k]*ScaleMAA))+(20/8)
                        else:
                            wf_MAA=1

                    #Calculate the weighting factor for Y90(wf_Y90) based on the predefined dose thresholds
                        if (imgY90[i,j,k]*ScaleY90)<=10 or (imgY90[i,j,k]*ScaleY90)>=200:
                            wf_Y90=0
                        elif (imgY90[i,j,k]*ScaleY90)>10 and (imgY90[i,j,k]*ScaleY90)<40:
                            wf_Y90=((1/30)*(imgY90[i,j,k]*ScaleY90))-(1/3)
                        elif (imgY90[i,j,k]*ScaleY90)>120 and (imgY90[i,j,k]*ScaleY90)<200:
                            wf_Y90=((-2/160)*(imgY90[i,j,k]*ScaleY90))+(20/8)
                        else:
                            wf_Y90=1

                     #calculate the max wf between Y90 and MAA (wf_QVH)
                        if (imgMAA[i,j,k]*ScaleMAA)<80 and (imgY90[i,j,k]*ScaleY90)>80:
                            wf_QVH = 1
                        elif (imgY90[i,j,k]*ScaleY90)<80 and (imgMAA[i,j,k]*ScaleMAA)>80:
                            wf_QVH = 1
                        else:
                            wf_QVH=max(wf_Y90,wf_MAA)

                        sumWF += wf_QVH
                        qFWeightedSum += abs(qF)*wf_QVH #at the level of the voxel, the quality factor is multiplied with the coresponding wighting factor

                        Qarray.append((qF,wf_QVH))  # crate a 2D array:first colum result of the dose division, second column the weighting factors

        Qarraysorted = sorted(Qarray) #sort the data based on the dose division (lowest to highest)

#calculate QVH quality factor
        qualityFactor= qFWeightedSum/sumWF
        print(self.mask_name,'qualityFactor:', qualityFactor)

        return Qarraysorted, sumWF

#create a dictionary with all the masks
    mask_dict = dict()

#create a list with the outputs of the class MaskResults, for each mask
res_list = [MaskResults(imgMAA, ScaleMAA, imgY90, ScaleY90, masks[i], mask_names[i]) for i in range(nr_masks)]

#plot QVH for each mask and save it in the patient directory
fig = plt.figure(figsize = (20, 10))
print(nr_masks)
for m_res_index in range(nr_masks):
    plt.subplot(nr_masks, nr_masks, m_res_index+1)
    plt.plot(res_list[m_res_index].x,res_list[m_res_index].y)
    plt.xlabel('qF ')
    plt.ylabel('Volume (%)')
    plt.xlim(0,5)
    plt.ylim(0,100)
    plt.grid()
    plt.title (res_list[m_res_index].mask_name)

FigFile=os.path.join(os.path.realpath(dirname),"Patient.png")
fig.savefig(FigFile)

plt.show()

#create a csv file that includes the qF pondereted and volume for each mask
Xs = [res.x for res in res_list]
Ys = [res.y for res in res_list]
max_length = max([len(x) for x in Xs])

print(max_length)

for res in res_list:
    res.extend_xy(max_length)

Xs = [res.x for res in res_list]
Ys = [res.y for res in res_list]

max_stars_repo_name: hlevillain/Quality-volume-histograms

max_issues_repo_name: hlevillain/Quality-volume-histograms

max_forks_repo_name: hlevillain/Quality-volume-histograms


text:
# -*- coding: utf-8 -*-
from __future__ import absolute_import
import json
import numpy as np
from benchpress.visualizer import util
import matplotlib
import matplotlib.pyplot as plt
import pylab

def value_labels(ax, rects):
    for rect in rects:
        height = rect.get_height()
        ax.text(rect.get_x() + rect.get_width()/2., 1.05*height,
                '%.2f' % float(height),
                ha='center', va='bottom')

def one_bar_per_cmd(args):
    matplotlib.rcParams['axes.labelsize'] = args.fontsize
    matplotlib.rcParams['axes.titlesize'] = args.fontsize + 2
    matplotlib.rcParams['xtick.labelsize'] = args.fontsize
    matplotlib.rcParams['ytick.labelsize'] = args.fontsize
    matplotlib.rcParams['legend.fontsize'] = args.fontsize
    matplotlib.rcParams['font.family'] = 'serif'
    matplotlib.rcParams['font.serif'] = ['Computer Modern Roman']
    matplotlib.rcParams['text.usetex'] = True
    matplotlib.rcParams['figure.max_open_warning'] = 400
    matplotlib.rcParams['image.cmap'] = "gray"
    matplotlib.rcParams['image.interpolation'] = "none"
    matplotlib.rcParams['figure.autolayout'] = True
    pylab.clf()                     # Essential! Otherwise the plots will be f**d up!
    pylab.figure(1)
    pylab.gca().yaxis.grid(True)    # Makes it easier to compare bars
    pylab.gca().xaxis.grid(False)
    pylab.gca().set_axisbelow(True)
    plt.style.use(args.pyplot_style)

    cmd_list = json.load(args.results)['cmd_list']
    cmd_list = util.filter_cmd_list(cmd_list, args.labels_to_include, args.labels_to_exclude)
    (cmd_list, label_map) = util.translate_dict(cmd_list, args.label_map)

    # For a bar-chart, we need labels, values, and errors
    labels = [util.texsafe(cmd['label']) for cmd in cmd_list]
    values = []
    std = []
    for cmd in cmd_list:
        res = util.extract_succeed_results(cmd, args.parse_regex, args.py_type)
        if len(res) == 0:
            res = [0]
        values.append(util.mean(res))
        std.append(util.standard_deviation(res))

    # Let's write the bar-chart using matplotlib
    margin = 0.05
    width = (1. - 1. * margin) / len(cmd_list)
    ind = np.arange(len(cmd_list))  # the x locations for each bar
    fig, ax = plt.subplots()
    ax.bar(ind * width, values, width, log=args.ylog, yerr=std)

    if args.ymin is not None:
        plt.ylim(ymin=float(args.ymin))
    if args.ymax is not None:
        plt.ylim(ymax=float(args.ymax))
    if args.title is not None:
        plt.title(util.texsafe(args.title))

    # Add some text for labels, title and axes ticks
    ax.set_xticks(ind*width+width/2.)
    ax.set_xticklabels(labels, rotation=args.xticklabel_rotation)

    if args.ylabel is not None:
        ax.set_ylabel(args.ylabel)

    # Now make some labels
    value_labels(ax, ax.patches)

    if args.output is not None:
        fname = args.output.name
        fformat = fname.split('.')[-1]
        print ("Writing file '%s' using format '%s'." % (fname, fformat))
        if args.mpld3:
            import mpld3
            if fformat == "html":
                mpld3.save_html(fig, args.output)
            elif fformat == "json":
                mpld3.save_json(fig, args.output)
            else:
                raise ValueError("--mpld3: The output must be either `html` or `json`")
        else:
            pylab.savefig(fname, format=fformat)
    else:
        if args.mpld3:
            import mpld3
            mpld3.show()
        else:
            pylab.show()

def main():
    parser = util.default_argparse('Plots the result of a Benchpress JSON-file (one bar per command)')
    parser.add_argument(
        '--title',
        default=None,
        help='The title of the chart.'
    )
    parser.add_argument(
        '--warmups',
        default=0,
        type=int,
        help="Specify the amount of samples from warm-up rounds."
    )
    parser.add_argument(
        '--ymax',
        default=None,
        help="Max value of the y-axis"
    )
    parser.add_argument(
        '--ymin',
        default=None,
        help="Min value of the y-axis"
    )
    parser.add_argument(
        '--ylog',
        default=False,
        action='store_true',
        help="Makes the y-axis logarithmic"
    )
    parser.add_argument(
        '--ylabel',
        default="Elapsed time in seconds",
        help="Label on the y-axis"
    )
    parser.add_argument(
        '--xticklabel-rotation',
        default=90,
        metavar='DEGREE',
        type=int,
        help="The rotation of the labels on the x-axis"
    )
    parser.add_argument(
        '--no-legend',
        default=False,
        action='store_true',
        help="Hide the legend box"
    )
    parser.add_argument(
        '--fontsize',
        default=12,
        type=int,
        help="Fontsize"
    )
    parser.add_argument(
        '--pyplot-style',
        default='seaborn-darkgrid',
        type=str,
        choices=plt.style.available,
        help="The matplotlib.pyplot.style to use"
    )
    parser.add_argument(
        '--mpld3',
        action='store_true',
        help="Generate mpld3 plots <https://mpld3.github.io>"
    )
    args = parser.parse_args()
    one_bar_per_cmd(args)

if __name__ == "__main__":
    main()


text:
import numpy as np
import string

def encrypt(text, rows):
    encrypted = []
    text = text.lower().replace(" ", "")
    text = text + "X" * (-len(text) % rows)
    text = np.array(list(text))
    text = text.reshape(-1, rows)
    text = text.T
    # print(text)
    for i in range(text.shape[0]):
        temp = list(text[i])
        encrypted += temp

    encrypted = "".join(encrypted)
    return encrypted

def decrypt(text, rows):
    decrypted = []
    text = text.lower().replace(" ", "")

    text = np.array(list(text))
    text = text.reshape(rows, -1)

    # print(text)
    for i in range(text.shape[1]):
        temp = list(text[:, i])
        decrypted += temp

    decrypted = "".join(decrypted)
    decrypted = decrypted.replace("x", "")

    return decrypted

if __name__ == "__main__":
    message = input("Enter message: ")
    rows = int(input("Enter no. of rows: "))
    # message = "test message"
    # rows = 4
    encrypted_text = encrypt(message, rows)
    decrypted_text = decrypt(encrypted_text, rows)
    print(encrypted_text)
    print(decrypted_text)

max_stars_repo_path: Semester 6/CSS/row_transposition.py

max_issues_repo_path: Semester 6/CSS/row_transposition.py

max_forks_repo_path: Semester 6/CSS/row_transposition.py


text:
# --------------

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns

# Code starts here
df = pd.read_csv(path)
df['state'] = df['state'].apply(lambda x: x.lower())
df['total'] = df['Jan'] + df['Feb'] + df['Mar']

sum_row = df[['Jan', 'Feb', 'Mar', 'total']].sum()
print(sum_row)

df_final = df.append(sum_row, ignore_index=True)
print(df_final)
# Code ends here

# --------------
import requests

# Code starts here
url = "https://en.wikipedia.org/wiki/List_of_U.S._state_abbreviations"
response = requests.get(url)
df1 = pd.read_html(response.content)[0]
df1 = df1.iloc[11:, :]
df1 = df1.rename(columns=df1.iloc[0, :]).iloc[1:, :]
df1['United States of America'] = df1['United States of America'].apply(lambda x: x.replace(" ", "")).astype(object)
# Code ends here

# --------------
df1['United States of America'] = df1['United States of America'].astype(str).apply(lambda x: x.lower())
df1['US'] = df1['US'].astype(str)

# Code starts here
mapping = df1.set_index('United States of America')['US'].to_dict()
df_final.insert(6, 'abbr', np.nan)
df_final['abbr'] = df_final['state'].map(mapping)
print(df_final.head(15))

# Code ends here

# --------------
# Code stars here

df_mississipi = df_final[df_final['state'] == 'mississipi'].replace(np.nan, 'MS')

df_final.replace(df_final.iloc[6], df_mississipi, inplace=True)

df_tenessee = df_final[df_final['state'] == 'tenessee'].replace(np.nan, 'TN')

df_final.replace(df_final.iloc[10], df_tenessee, inplace=True)

# Code ends here

# --------------
# Code starts here

df_final = df_final.groupby('abbr')
df_sub = df_final['abbr', 'Jan', 'Feb', 'Mar','total'].sum()
formatted_df = df_sub.applymap(lambda x: "${:,.0f}".format(x))
print(formatted_df)
# Code ends here

# --------------
# Code starts here
sum_row = df[["Jan", "Feb", "Mar", "total"]].sum()

df_sub_sum = pd.DataFrame(data=sum_row).T

#apply $ to the sum
df_sub_sum = df_sub_sum.applymap(lambda x: "${:,.0f}".format(x))

# append the sum
print(formatted_df)
final_table = formatted_df.append(df_sub_sum)
print(final_table)
# rename the index
final_table = final_table.rename(index={0: "Total"})
print(final_table)
# Code ends here

# --------------
# Code starts here

df_sub['total'] = df[["Jan", "Feb", "Mar"]].sum()
plot = df_sub['total'].plot.pie(y='total', figsize=(6, 3))
# Code ends here


text:
import sys
import gdal
from osgeo import osr
from ftplib import FTP
from datetime import datetime, timedelta
from argparse import ArgumentParser
from netCDF4 import Dataset
import numpy as np

def parse_arguments():
	parser = ArgumentParser(description='Fetches latest rainrader  image from neo server')
	parser.add_argument('-u', '--url', help='url to grab from', default='data.knmi.nl',required=False)
	parser.add_argument('-o', '--outdir', help='output directory', default='./',required=False)
	args = parser.parse_args()
	return args

def main():
	args = parse_arguments()
	uri = args.url
	d = datetime.today() - timedelta(days=2)
	filelist = []
	print 'downloading data'
	ftp = FTP(uri)
	ftp.login()
	ftp.cwd('download/radar_corr_accum_24h/1.0/noversion/'+str(d.year)+'/'+str(d.month)+'/'+str(d.day)+'/')
	ftp.dir('-t',filelist.append)
	for i, r in enumerate(filelist):
		filelist[i] = r.split()
	filelist.sort(key = lambda x : x[8])
	filelist.reverse()
	filename = filelist[0][8]
	outdir = args.outdir
	ftp.retrbinary("RETR " + filename, open(outdir+'/'+filename,"wb").write)
	ftp.quit()
	print 'processing'
	rootgrp = Dataset(outdir+'/'+filename, "r", format="NETCDF4")
	PV = np.array(rootgrp.groups['image1']['image_data'])
	cal = rootgrp.groups['image1']['calibration'].calibration_formulas.replace('GEO =','')
	def calibrate(PV):
		if (PV == 65535):
			return PV
		return eval(cal)
	#caldata = map(calibrate,data)
	caldata = eval(cal)

	x_pixels = 700  # number of pixels in x
	y_pixels = 765  # number of pixels in y
	PIXEL_SIZE = 1  # size of the pixel...
	x_min = 0 #575694 - (700000/2)
	y_max = -3649.979 #6819814 # + (765000/2)  # x_min & y_max are like the "top left" corner.
	proj_projection = '+proj=stere +lat_0=90 +lon_0=0 +lat_ts=60 +a=6378.14 +b=6356.75 +x_0=0 y_0=0'
	srs = osr.SpatialReference()
	srs.ImportFromProj4(proj_projection)
	wkt_projection = srs.ExportToWkt()
	#wkt_projection = 'PROJCS["unnamed",GEOGCS["unnamed ellipse",DATUM["WGS84",SPHEROID["unnamed",6378.14,298.1832632071011]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]],PROJECTION["Polar_Stereographic"],PARAMETER["latitude_of_origin",90],PARAMETER["central_meridian",0],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0]]'
	driver = gdal.GetDriverByName('GTiff')

	dataset = driver.Create(
		outdir+'/'+filename.replace('.h5','.tif'),
		x_pixels,
		y_pixels,
		1,
		gdal.GDT_Float32, )
	dataset.SetProjection(wkt_projection)
	dataset.SetGeoTransform((x_min, PIXEL_SIZE, 0, y_max, 0, -PIXEL_SIZE))
	dataset.GetRasterBand(1).WriteArray(caldata)
	dataset.GetRasterBand(1).SetNoDataValue(655.35)
	dataset.FlushCache()  # Write to disk.
	print 'saved as ' + outdir+'/'+filename.replace('.h5','.tif')
	print 'warping a copy to wgs84'
	# Open source dataset
	src_ds = gdal.Open(outdir+'/'+filename.replace('.h5','.tif'))

	# Define target SRS
	dst_srs = osr.SpatialReference()
	dst_srs.ImportFromEPSG(4326)
	dst_wkt = dst_srs.ExportToWkt()

	error_threshold = 0.125  # error threshold --> use same value as in gdalwarp
	resampling = gdal.GRA_NearestNeighbour

	# Call AutoCreateWarpedVRT() to fetch default values for target raster dimensions and geotransform
	tmp_ds = gdal.AutoCreateWarpedVRT( src_ds,
									   None, # src_wkt : left to default value --> will use the one from source
									   dst_wkt,
									   resampling,
									   error_threshold )

	# Create the final warped raster
	dst_ds = gdal.GetDriverByName('GTiff').CreateCopy(outdir+'/'+filename.replace('.h5','_wgs84.tif'), tmp_ds)
	dst_ds = None
	print 'done'

if __name__ == '__main__':
    status = main()
    sys.exit(status)


text:
# Example 2b: Snake!
# Illustrative example of FitsGeo usage
# Shows how to easily create multiple (repeating) objects with some math laws
import fitsgeo as fg
from numpy import sin, exp, linspace  # To use math functions

# Create laws for positions of snakes parts
x = linspace(0, 5, 50)  # Create 50 evenly spaced x values from 0 to 5
# Law for z coordinate
z_p = [5 * sin(3 * x_i) * 0.3 * exp(-0.4 * x_i) for x_i in x]
r = [0.02 * exp(0.2 * x_i) for x_i in x]  # Law for radius of spheres

# Define materials for snake and hat
snake_mat = fg.Material.database("MAT_SKIN_ICRP", color="green")
hat_mat = fg.Material.database("MAT_POLYETHYLENE", color="blue")

snake = []  # List with all parts of snake
for i in range(16, len(x)):  # Take x values starting from i=16
	snake.append(
		fg.SPH(xyz0=[x[i] - 3, 0, z_p[i]], r=r[i], material=snake_mat))

last_part = snake[len(snake)-1]
hat = fg.TRC(
	xyz0=[
		last_part.x0,
		last_part.y0 + last_part.r,
		last_part.z0],
	h=[0, last_part.diameter/1.5, 0],
	r_1=last_part.r/1.5,
	r_2=1e-9, material=hat_mat)

void = fg.RCC(
	xyz0=[
		snake[0].x0-0.1,
		snake[0].y0,
		snake[0].z0], h=[4, 0, 0], r=3, material=fg.MAT_VOID)

# Define cells
outer_void = fg.Cell([+void], material=fg.MAT_OUTER)
cells = []
for i in range(len(snake)):
	cells.append(
		fg.Cell([-snake[i]], name=f"Snake Part {i}", material=snake_mat))
# Add a nice hat on head =)
hat_cell = fg.Cell([-hat], name="Hat!", material=hat_mat)

text = ""
for s in snake:
	text += +s

void_cell = fg.Cell(
	[-void, " ", text + +hat], name="Vacuum", material=fg.MAT_VOID)

fg.create_scene(ax_length=5)  # Create scene before draw

# Draw snake
for s in snake:
	s.draw()
hat.draw(truncated=False)  # Difficulties with truncated cone =(
void.draw()

fg.phits_export(to_file=True, inp_name="snake")  # Export sections

max_stars_repo_name: GordoNice/FitsGeo

max_issues_repo_name: GordoNice/FitsGeo

max_forks_repo_name: GordoNice/fitsgeo

reason: from numpy


text:
import numpy as np
import scipy as sp

#import data
x = np.array(([1,3,5]) , dtype='float')
y = np.array(([0.66]) , dtype='float')

#define layers size
inputs = 3
outputs = 1

#sigmoid fnco
def sigmoid(k):
    return(1/(1+np.exp(-k)))

#derivative of sigmoid
def d_sigmoid(k):
    return(sigmoid(k) * (1-sigmoid(k)))

#weights
w1 = np.random.rand(inputs,outputs)

#define cost func
def cost(w1,x,y):

    #define forward
    z2 = x.dot(w1)
    ht = sigmoid(z2)

    #cost
    j = (0.5/len(y)) *sum(ht - y)**2

    #gradient
    d1 = (ht-y) * d_sigmoid(z2)
    grad = d1*x

    return(j,grad)

# using advanced optimisation
from scipy import optimize as opt

w1, nfeval ,un = opt.fmin_tnc(func = cost,x0=w1,args=(x,y))
print(w1)

#call answer
z2 = x.dot(w1)
ht = sigmoid(z2)
print(ht)

max_stars_repo_path: Machine Learning/3.neural network/1.to understand NN.py

reason: import numpy,import scipy,from scipy


text:
# -*- coding: utf-8 -*-
"""
author: Miroslav Jirik
"""

import logging
logger = logging.getLogger(__name__)
import numpy as np

import scipy
# from scipy.signal import butter, lfilter, freqz
import scipy.signal
import matplotlib.pyplot as plt

def butter_bandpass_1d(lowcut, highcut, fs, order=5):
    nyq = 0.5 * fs
    low = lowcut / nyq
    high = highcut / nyq
    b, a = scipy.signal.butter(order, [low, high], btype='band')
    return b, a

def butter_bandpass_filter_1d(data, lowcut, highcut, fs, order=5):
    b, a = butter_bandpass_1d(lowcut, highcut, fs, order=order)
    y = scipy.signal.lfilter(b, a, data)
    return y

def butter_bandpass_freq_filter_1d(data, lowcut, highcut, fs, order=5):
    b, a = butter_bandpass_1d(lowcut, highcut, fs, order=order)
    dataf = np.fft.fft(data)
    # w, h = freqz(b, a, worN=len(dataf), whole=True)
    w, h = scipy.signal.freqz(b, a, worN=len(dataf), whole=True)
    plt.figure(3)
    plt.subplot(211)
    nyq = (fs * 0.5 / np.pi)

    plt.plot(w*nyq, abs(h))
    # filtered_spectrum = dataf * np.real(h)
    # filtered_spectrum = dataf * np.abs(h)
    filtered_spectrum = dataf * h
    filtered = np.fft.ifft(filtered_spectrum)
    # plt.subplot(212)
    # plt.plot(abs(dataf), label='data')
    # plt.plot(abs(h), label='filter')
    # plt.legend()
    # plt.show()
    return np.real(filtered)
    # return np.abs(filtered)

def apply_filter_on_abs(shspectrum, filter):
    radii, angle = R2P(shspectrum)
    radii *= filter
    shspectrum = P2R(radii, angle)
    return shspectrum

def apply_filter(shspectrum, filter):
    ab = np.real(shspectrum)
    im = np.imag(shspectrum)
    shspectrum = ab * filter + 1j * im * filter
    # radii, angle = R2P(shspectrum)
    # radii *= filter
    # shspectrum = P2R(radii, angle)
    return shspectrum

def power_filter(shape, e, dist=None):
    # if dist is None:
    #     dist = construct_filter_dist(shape)
    output = np.power(dist, e)
    return output

def lopass_ideal_fft_mask(shape, radius, freqs=None, **kwargs):
    # if dist is None:
    #     dist = construct_filter_dist(shape)
    output = freqs <= radius
    return output

def hipass_ideal_fft_mask(shape, radius, freqs=None, **kwargs):
    # if dist is None:
    #     dist = construct_filter_dist(shape)
    output = freqs > radius
    return output

def bandpass_butter_fft_mask(shape, flow, fhigh, freqs=None, order=2):
    # TODO implement this function

    # digital
    filter_fs = np.max([fc, np.max(freqs)])
    # filter_fs = fc * 2.0
    nyq = 0.5 * filter_fs
    low_digital = 0.5 * fc / (nyq * np.pi)
    # high_digital = 0.5 * highcut / (nyq * np.pi)
    b_digital, a_digital = scipy.signal.butter(order, Wn=low_digital, btype='lowpass', analog=False)
    w_digital, h_digital = scipy.signal.freqz(b_digital, a_digital, worN=(freqs / filter_fs))

def lopass_butter_fft_mask(shape, fc, freqs=None, order=2):
    # digital
    filter_fs = np.max([fc, np.max(freqs)])
    # filter_fs = fc * 2.0
    nyq = 0.5 * filter_fs
    low_digital = 0.5 * fc / (nyq * np.pi)
    # high_digital = 0.5 * highcut / (nyq * np.pi)
    b_digital, a_digital = scipy.signal.butter(order, Wn=low_digital, btype='lowpass', analog=False)
    w_digital, h_digital = scipy.signal.freqz(b_digital, a_digital, worN=(freqs / filter_fs))

    return h_digital

def hipass_butter_fft_mask(shape, fc, freqs=None, order=2):
    # digital
    filter_fs = np.max([fc, np.max(freqs)])
    # filter_fs = fc * 2.0
    nyq = 0.5 * filter_fs
    high_digital = 0.5 * fc / (nyq * np.pi)
    b_digital, a_digital = scipy.signal.butter(order, Wn=high_digital, btype='highpass', analog=False)
    w_digital, h_digital = scipy.signal.freqz(b_digital, a_digital, worN=(freqs / filter_fs))

    return h_digital

def dist_from_center(shape, axis_size=None):
    if axis_size is None:
        axis_size = np.ones(len(shape))
    axis_size = np.asarray(axis_size, dtype=np.float)

    center = (np.asarray(shape) - 1) / 2.0

    xi = []
    for i in range(len(shape)):
        xi.append(range(shape[i]))
    yi = np.meshgrid(*xi, indexing='ij')

    dist = np.zeros(shape)
    for i in range(len(shape)):
        dist += ((yi[i] - center[i]) * axis_size[i]) ** 2
    dist = dist**0.5
    return dist

def fftfreq(shape, spacing=None):
    """
    N dimensional variant of numpy.fft.fftfreq

    :param shape:
    :param spacing: if None it is set to 1. If it is scalar it is used for every axis
    :return:
    """
    input = np.zeros(shape, dtype=np.bool)

    # set first pixel (for general dimension) to one
    # input[np.zeros([len(input.shape), 1], dtype=int).tolist()] = 1
    # shinput = np.fft.fftshift(input)
    # shdist = scipy.ndimage.morphology.distance_transform_edt(shinput, sampling=spacing)
    # dist = np.fft.ifftshift(shdist)

    if np.isscalar(shape):
        shape = [shape]

    if spacing is None:
        spacing = np.ones(len(shape))

    if np.isscalar(spacing):
        spacing = spacing * np.ones(len(shape))

    spacing = np.asarray(spacing, dtype=np.float)

    center = np.ceil((np.asarray(shape) - 1) / 2.0)

    xi = []
    for i in range(len(shape)):
        xi.append(range(shape[i]))
    yi = np.meshgrid(*xi, indexing='ij')

    shdist = np.zeros(shape)
    for i in range(len(shape)):
        shdist += (((yi[i] - center[i]) / (spacing[i] * shape[i]) ) ** 2) # * spacing[i] / shape[i]
    shdist = (shdist**0.5)

    dist = np.fft.ifftshift(shdist)
    return dist

def spectrum_analysis_3d(spectrum, fs=[1, 1, 1]):
    spectrum = np.abs(spectrum)
    voxelsize = 1.0 / np.asarray(fs)
    freqs = fftfreq(spectrum.shape, spacing=voxelsize)

    # along axis
    fr1 = freqs[:,0]
    sp1 = spectrum[:, 0]
    plt.subplot(121)
    plt.plot(fr1, sp1)
    # plt.show()

    # round
    voxelvol = np.prod(voxelsize)
    surface = 4 * np.pi * freqs

    spf = spectrum.flatten()
    freqsf = freqs.flatten()

    unq_freqs = np.unique(freqsf)
    freqs_spect = dict(zip(unq_freqs, [0.0] * len(unq_freqs)))
    for si, fi in zip(spf, freqsf):
        freqs_spect[fi] += si #* voxelvol / (4 * np.pi * fi)

    sp = freqs_spect.values()
    fr = freqs_spect.keys()
    fr_sp = sorted(zip(fr, sp))
    fr, sp = zip(*fr_sp)

    sp = list(sp)
    # sp is density function

    for i in range(1, len(unq_freqs)):
        sp[i] += sp [i - 1]

    correction = voxelvol / (4 * np.pi * (np.asarray(fr)**2 + np.mean(voxelsize)))
    correction = voxelvol / (2 * np.pi * (np.asarray(fr) + np.mean(voxelsize)))

    # sp is now distribution function

    plt.subplot(121)
    plt.suptitle(str(np.max(freqs)))
    plt.plot(fr, sp * correction)
    #plt.show()

def P2R(radii, angles):
    radii = np.asarray(radii)
    angles = np.asarray(angles)
    return radii * np.exp(1j*angles)

def R2P(x):
    return np.abs(x), np.angle(x)

def spectrum_filtration(spectrum, fs=None, exponent=0, freq0=0, freq1=None, filter_type="butter"):
    """
    Filter spectrum based on frequency
    :param spectrum:
    :param exponent:
    :param freq0:
    :param freq1:
    :return:
    """
    log = "spectrum_filtration() - begin"
    logger.debug(log)
    if freq1 < 0:
        freq1 = None

    if freq0 <= 0:
        freq0 = None

    if fs is None:
        fs = np.ones(len(spectrum.shape))

    voxelsize = 1.0 / np.asarray(fs)

    # dist = dist_from_center(spectrum.shape, axis_size=voxelsize)
    freqs = fftfreq(spectrum.shape, spacing=voxelsize)

    # spectrum = apply_filter(spectrum, pfilter)

    filt = get_freq_filter_mask(spectrum.shape, exponent, freq0, freq1, freqs,
                                filter_type
                                )

    log = "spectrum_filtration() - before apply_filter"
    logger.debug(log)
    spectrum = apply_filter(spectrum, filt)
    # spectrum *= filt
    # spectrum = np.fft.ifftshift(shspectrum)
    signal = get_real_ifft(spectrum)
    log = "spectrum_filtration() - end"
    logger.debug(log)

    return signal, filt, spectrum, freqs

def get_real_ifft(spectrum):
    return np.real(np.fft.ifftn(spectrum))

def get_freq_filter_mask(spectrum_shape, exponent, freq_start, freq_stop, freqs, filter_type):
    filter_type = "ideal"

    if filter_type is "ideal":
        lopass_fcn = lopass_ideal_fft_mask
        hipass_fcn = hipass_ideal_fft_mask
    elif filter_type is "butter":
        lopass_fcn = lopass_butter_fft_mask
        hipass_fcn = hipass_butter_fft_mask
    else:
        logger.error("Unknown filter type")

    # shspectrum = np.fft.fftshift(spectrum)
    filt = power_filter(spectrum_shape, exponent, dist=freqs)
    # not sure if the abs is correct
    if freq_start is not None:
        filt *= np.abs(hipass_fcn(spectrum_shape, freq_start, freqs=freqs))

    if freq_stop is not None:
        filt *= np.abs(lopass_fcn(spectrum_shape, freq_stop, freqs=freqs))

    return filt

def show(real_signal, filter=None, spectrum=None, freqs=None, log_view=False):
    import matplotlib.pyplot as plt

    # plt.gray()
    plt.subplot(331)
    plt.imshow(real_signal, cmap="gray")
    plt.colorbar()

    plt.title("signal")

    # shfilter = np.fft.fftshift(filter)
    plt.subplot(332)
    plt.imshow(np.abs(np.fft.fftshift(filter)))
    plt.colorbar()
    plt.title("filter")

    if freqs is not None:
        ax = plt.subplot(333)
        plt.imshow(np.abs(np.fft.fftshift(freqs)))
        plt.colorbar()
        ax.set_title("freqs")

    shspecturm = np.fft.fftshift(spectrum)
    if log_view:
        shspecturm = np.log(shspecturm)

    plt.subplot(334)
    plt.imshow(np.abs(shspecturm))
    plt.colorbar()
    plt.title("abs")

    plt.subplot(335)
    plt.imshow(np.angle(shspecturm))
    plt.colorbar()
    plt.title("angle")

    ax = plt.subplot(337)
    ax.imshow(np.real(shspecturm))
    plt.colorbar()
    ax.set_title("real")

    ax = plt.subplot(338)
    ax.imshow(np.imag(shspecturm))
    plt.colorbar()
    ax.set_title("imag")

    # plt.show()

max_stars_repo_name: mjirik/ndnoise

max_issues_repo_name: mjirik/ndnoise

max_forks_repo_name: mjirik/ndnoise


text:
# Profile Likelihood Generator
# Marisa Eisenberg (marisae@umich.edu)
# Yu-Han Kao (kaoyh@umich.edu) -7-9-17

# Input definitions
# params = starting parameters (all, including the one to be profiled)
# profparam = index within params for the parameter to be profiled
#   ---reminder to make this allow you to pass the name instead later on
# costfun = cost function for the model - should include params, times, and data as arguments.
#   Note costfun doesn't need to be specially set up for fixing the profiled parameter,
#   it's just the regular function you would use to estimate all the parameters
#   (it will get reworked to fix one of them inside ProfLike)
# times, data = data set (times & values, or whatever makes sense)
#   ---possibly change this so it's included in costfun and not a separate set of inputs? Hmm.
# perrange = the percent/fraction range to profile the parameter over (default is 0.5)
# numpoints = number of points to profile at in each direction (default is 10)

# Output
# A list with:
#   - profparvals: the values of the profiled parameter that were used
#   - fnvals: the cost function value at each profiled parameter value
#   - convergence: the convergence value at each profiled parameter value
#   - paramestvals: the estimates of the other parameters at each profiled parameter value

import numpy as np
import scipy.optimize as optimize
import copy

import sir_ode
import sir_cost

def proflike (params, profindex, cost_func, times, data, perrange = 0.5, numpoints = 10):
	profrangedown = np.linspace(params[profindex], params[profindex] * (1 - perrange), numpoints).tolist()
	profrangeup = np.linspace(params[profindex], params[profindex] * (1 + perrange), numpoints).tolist()[1:] #skip the duplicated values
	profrange = [profrangedown, profrangeup]
	currvals = []
	currparams = []
	currflags = []

	def profcost (fit_params, profparam, profindex, data, times, cost_func):
		paramstest = fit_params.tolist()
		paramstest.insert(profindex, profparam)
		return cost_func (paramstest, data, times)

	fit_params = params.tolist() #make a copy of params so we won't change the origianl list
	fit_params.pop(profindex)
	print 'Starting profile...'
	for i in range(len(profrange)):
		for j in profrange[i]:
			print i, j
			optimizer = optimize.minimize(profcost, fit_params, args=(j, profindex, data, times, cost_func), method='Nelder-Mead')
			fit_params = np.abs(optimizer.x).tolist() #save current fitted params as starting values for next round
			#print optimizer.fun
			currvals.append(optimizer.fun)
			currflags.append(optimizer.success)
			currparams.append(np.abs(optimizer.x).tolist())

	#structure the return output
	profrangedown.reverse()
	out_profparam = profrangedown+profrangeup
	temp_ind = range(len(profrangedown))
	temp_ind.reverse()
	out_params = [currparams[i] for i in temp_ind]+currparams[len(profrangedown):]
	out_fvals = [currvals[i] for i in temp_ind]+currvals[len(profrangedown):]
	out_flags = [currflags[i] for i in temp_ind]+currflags[len(profrangedown):]
	output = {'profparam': out_profparam, 'fitparam': np.array(out_params), 'fcnvals': out_fvals, 'convergence': out_flags}
	return output

reason: import numpy,import scipy