text:
#!/usr/bin/env python

"""

The rendering widget

It renderes a projection via the OpenCL (defined in volumerender.py)
into a texture which is drawn by simple OpenGL calls onto the canvas.

It should handle all user interaction via a transformation model.

author: Martin Weigert
email: mweigert@mpi-cbg.de

understanding glBlendFunc:

first color:     d
second color:    s
resulting color: c

c = s*S + d*D

where S and D are set with glBlendFunc(S,D)

e.g. glBlendFunc(GL_SRC_ALPHA,GL_ONE_MINUS_SRC_ALPHA)

c = s*s.w + d*(1-s.w)

"""

from __future__ import absolute_import, print_function, unicode_literals, division

import logging

logger = logging.getLogger(__name__)
# logger.setLevel(logging.DEBUG)

import sys
import os
from PyQt5 import QtCore
from PyQt5 import QtGui, QtWidgets
from PyQt5 import QtOpenGL
from PyQt5.QtGui import QOpenGLShaderProgram, QOpenGLShader
from OpenGL.GL import *
import OpenGL.arrays.vbo as glvbo
import spimagine
from spimagine.volumerender.volumerender import VolumeRenderer
from spimagine.utils.transform_matrices import *
from spimagine.models.transform_model import TransformModel
from spimagine.models.data_model import DataModel
from spimagine.gui.mesh import Mesh, SphericalMesh, EllipsoidMesh
import numpy as np
from spimagine.gui.gui_utils import *

# on windows numpy.linalg.inv crashes without notice, so we have to import scipy.linalg
if os.name == "nt":
    from scipy import linalg
else:
    from numpy import linalg

import pyopencl.array as cl_array
from gputools import OCLArray

import time
from spimagine.utils.quaternion import Quaternion

# from spimagine.shaders import vertShaderTex, fragShaderTex, vertShaderCube, fragShaderCube

# logger.setLevel(logging.DEBUG)

def _next_golden(n):
    res = round((np.sqrt(5) - 1.) / 2. * n)
    return int(round((np.sqrt(5) - 1.) / 2. * n))

def absPath(myPath):
    """ Get absolute path to resource, works for dev and for PyInstaller """
    try:
        # PyInstaller creates a temp folder and stores path in _MEIPASS
        base_path = sys._MEIPASS
        return os.path.join(base_path, os.path.basename(myPath))
    except Exception:
        base_path = os.path.abspath(os.path.dirname(__file__))
        return os.path.join(base_path, myPath)

class GLWidget(QtOpenGL.QGLWidget):
    _dataModelChanged = QtCore.pyqtSignal()

    _BACKGROUND_BLACK = (0., 0., 0., 0.)
    _BACKGROUND_WHITE = (1., 1., 1., 0.)

    def __init__(self, parent=None, N_PREFETCH=0, interpolation="linear", **kwargs):
        logger.debug("init")
        #
        # fmt = QtOpenGL.QGLFormat(QtOpenGL.QGL.AlphaChannel)
        #
        # super(GLWidget, self).__init__(fmt,parent, **kwargs)

        super(GLWidget, self).__init__(parent, **kwargs)

        self.parent = parent
        self.texture_LUT = None

        self.setAcceptDrops(True)

        self.renderer = VolumeRenderer((spimagine.config.__DEFAULT_TEXTURE_WIDTH__,
                                        spimagine.config.__DEFAULT_TEXTURE_WIDTH__),
                                       interpolation=interpolation)

        self.renderer.set_projection(mat4_perspective(60, 1., .1, 100))
        # self.renderer.set_projection(projMatOrtho(-2,2,-2,2,-10,10))

        self.output = np.zeros([self.renderer.height, self.renderer.width], dtype=np.float32)
        self.output_alpha = np.zeros([self.renderer.height, self.renderer.width], dtype=np.float32)

        self.sliceOutput = np.zeros((100, 100), dtype=np.float32)

        self.setTransform(TransformModel())

        self.renderTimer = QtCore.QTimer(self)
        self.renderTimer.setInterval(10)
        self.renderTimer.timeout.connect(self.onRenderTimer)
        self.renderTimer.start()
        self.renderedSteps = 0

        self.N_PREFETCH = N_PREFETCH

        self.NSubrenderSteps = 1

        self.dataModel = None

        self.meshes = []

        # self.setMouseTracking(True)

        self._dataModelChanged.connect(self.dataModelChanged)

        self.refresh()

        # self.installEventFilter(self)

    def set_background_mode_black(self, mode_back=True):
        self._background_mode_black = mode_back
        self.refresh()

    def setModel(self, dataModel):
        logger.debug("setModel to %s" % dataModel)
        if self.dataModel is None or (self.dataModel != dataModel):
            self.dataModel = dataModel
            self.transform.setModel(dataModel)
            self.dataModel._dataSourceChanged.connect(self.dataSourceChanged)
            self.dataModel._dataPosChanged.connect(self.dataPosChanged)
            self._dataModelChanged.emit()

    def dragEnterEvent(self, event):
        if event.mimeData().hasUrls():
            event.accept()
        else:
            event.ignore()

    def dropEvent(self, event):
        def _url_to_path(url):
            path = url.toLocalFile()
            if spimagine.config.__SYSTEM_DARWIN__:
                path = spimagine.config._parseFileNameFix(path)
            return path

        self.setCursor(QtCore.Qt.BusyCursor)
        urls = event.mimeData().urls()

        if len(urls) == 0:
            return
        elif len(urls) == 1:
            path = _url_to_path(urls[0])
        elif len(urls) > 1:
            path = tuple(_url_to_path(url) for url in urls)

        try:
            if self.dataModel:
                self.dataModel.loadFromPath(path, prefetchSize=self.N_PREFETCH)
            else:
                self.setModel(DataModel.fromPath(path, prefetchSize=self.N_PREFETCH))

            self.setCursor(QtCore.Qt.ArrowCursor)
        except Exception as e:
            QtWidgets.QMessageBox.warning(self, "", "Error loading Data:\n %s" % str(e))

    def set_colormap(self, name):
        """name should be either jet, hot, gray, coolwarm"""

        try:
            arr = spimagine.config.__COLORMAPDICT__[name]
            self._set_colormap_array(arr)
        except KeyError:
            print("could not load colormap '%s'" % name)
            print("valid names: %s" % list(spimagine.config.__COLORMAPDICT__.keys()))

    def set_colormap_rgb(self, color=[1., 1., 1.]):
        self._set_colormap_array(np.outer(np.linspace(0, 1., 255), np.array(color)))

    def _set_colormap_array(self, arr):
        """arr should be of shape (N,3) and gives the rgb components of the colormap"""
        if not arr.ndim == 2 and arr.shape[-1] == 3:
            raise ValueError("wrong shape of color array: should be (N,3) but is %s")

        self.makeCurrent()

        self.texture_LUT = fillTexture2d(arr.reshape((1,) + arr.shape), self.texture_LUT)
        self.refresh()

    def _shader_from_file(self, fname_vert, fname_frag):
        shader = QOpenGLShaderProgram()
        shader.addShaderFromSourceFile(QOpenGLShader.Vertex, fname_vert)
        shader.addShaderFromSourceFile(QOpenGLShader.Fragment, fname_frag)
        shader.link()
        shader.bind()
        logger.debug("GLSL program log:%s", shader.log())
        return shader

    def initializeGL(self):

        self.resized = True

        logger.debug("initializeGL")

        self.programTex = self._shader_from_file(absPath("shaders/texture.vert"),
                                                 absPath("shaders/texture.frag"))

        self.programCube = self._shader_from_file(absPath("shaders/box.vert"),
                                                  absPath("shaders/box.frag"))

        self.programSlice = self._shader_from_file(absPath("shaders/slice.vert"),
                                                   absPath("shaders/slice.frag"))

        self.programMesh = self._shader_from_file(absPath("shaders/mesh.vert"),
                                                  absPath("shaders/mesh.frag"))

        self.programMeshLight = self._shader_from_file(
            absPath("shaders/mesh_light.vert"),
            absPath("shaders/mesh_light.frag"))

        self.texture = None
        self.textureAlpha = None
        self.textureSlice = None

        self.quadCoord = np.array([[-1., -1., 0.],
                                   [1., -1., 0.],
                                   [1., 1., 0.],
                                   [1., 1., 0.],
                                   [-1., 1., 0.],
                                   [-1., -1., 0.]])

        self.quadCoordTex = np.array([[0, 0],
                                      [1., 0.],
                                      [1., 1.],
                                      [1., 1.],
                                      [0, 1.],
                                      [0, 0]])

        # self.cubeCoords = create_cube_coords([-1,1,-1,1,-1,1])

        self.set_colormap(spimagine.config.__DEFAULTCOLORMAP__)

        glEnable(GL_BLEND)

        # glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA)

        # self.set_interpolation(spimagine.config.__DEFAULT_INTERP__ != "nearest")

        # glLineWidth(1.0);
        glBlendFunc(GL_ONE, GL_ONE)

        glEnable(GL_LINE_SMOOTH);
        glDisable(GL_DEPTH_TEST)

        glLineWidth(spimagine.config.__DEFAULT_BOX_LINEWIDTH__)

        # self.set_background_color(0,0,0,.0)
        self.set_background_mode_black(True)
        self.clear_canvas()

        # self.set_background_color(1,1,1,.6)

    def clear_canvas(self):
        if self._background_mode_black:
            glClearColor(*self._BACKGROUND_BLACK)
        else:
            glClearColor(*self._BACKGROUND_WHITE)

        if glCheckFramebufferStatus(GL_FRAMEBUFFER) == GL_FRAMEBUFFER_COMPLETE:
            glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)

    def set_interpolation(self, interpolate=True):
        interp = "linear" if interpolate else "nearest"
        self.renderer.rebuild_program(interpolation=interp)
        self.refresh()

    def setTransform(self, transform):
        self.transform = transform
        self.transform._transformChanged.connect(self.refresh)
        self.transform._stackUnitsChanged.connect(self.setStackUnits)
        self.transform._boundsChanged.connect(self.setBounds)

    def dataModelChanged(self):
        logger.debug("+++++++++ data model changed")

        if self.dataModel:
            logger.debug("dataModelchanged")

            self.renderer.set_data(self.dataModel[0], autoConvert=True)

            mi, ma = self._get_min_max()

            self.transform.reset(minVal=mi,
                                 maxVal=ma,
                                 stackUnits=self.dataModel.stackUnits())

            self.meshes = []
            self.refresh()

    def _get_min_max(self):
        # as amax is too slow for bug arrays, do it on the gpu

        if self.dataModel:
            try:
                im = self.renderer.dataImg
                tmp_buf = OCLArray.empty(im.shape, im.dtype)
                tmp_buf.copy_image(im)
                mi = float(cl_array.min(tmp_buf).get())
                ma = float(cl_array.max(tmp_buf).get())

            except Exception as e:
                print(e)
                mi = np.amin(self.dataModel[0])
                ma = np.amax(self.dataModel[0])
        return mi, ma

    def set_background_color(self, r, g, b, a=1.):
        self._background_color = (r, g, b, a)
        glClearColor(r, g, b, a)

    def dataSourceChanged(self):

        logger.debug("dataSourcechanged")

        self.renderer.set_data(self.dataModel[0], autoConvert=True)

        mi, ma = self._get_min_max()

        self.transform.reset(minVal=mi,
                             maxVal=ma,
                             stackUnits=self.dataModel.stackUnits())

        self.refresh()

    def setBounds(self, x1, x2, y1, y2, z1, z2):
        self.cubeCoords = create_cube_coords([x1, x2, y1, y2, z1, z2])
        self.renderer.set_box_boundaries([x1, x2, y1, y2, z1, z2])

    def setStackUnits(self, px, py, pz):
        logger.debug("setStackUnits to %s" % [px, py, pz])
        self.renderer.set_units([px, py, pz])

    def dataPosChanged(self, pos):
        self.renderer.update_data(self.dataModel[pos])
        self.refresh()

    def refresh(self):
        # if self.parentWidget() and self.dataModel:
        #     self.parentWidget().setWindowTitle("SpImagine %s"%self.dataModel.name())

        self.renderUpdate = True
        self.renderedSteps = 0

    def resizeGL(self, width, height):
        # somehow in qt5 the OpenGLWidget width/height parameters above are double the value of self.width/height
        self._viewport_width, self._viewport_height = width, height

    def add_mesh(self, mesh=SphericalMesh()):
        """
        adds a mesh with vertices and facecolor/edgecolor to be drawn

        mesh is an instance of spimagine.gui.Mesh, e.g.

        mesh = Mesh(vertices = [[0,1,0],[0,1,0],...],
                    normals = [[0,1,0],[0,1,0],...],
                    facecolor = (1.,.4,.4,.2),
                    edgecolor = None,...)

        there are some predefined meshes like
        SphericalMesh, EllipsoidMesh ...
        """

        self.meshes.append([mesh,
                            glvbo.VBO(mesh.vertices.astype(np.float32, copy=False)),
                            glvbo.VBO(np.array(mesh.normals).astype(np.float32, copy=False)),
                            glvbo.VBO(np.array(mesh.indices).astype(np.uint32, copy=False),
                                      target=GL_ELEMENT_ARRAY_BUFFER)])

        self.refresh()
        # sort according to opacity as the opaque objects should be drawn first
        # self.meshes.sort(key=lambda x: x[0].alpha, reverse=True)

    def _paintGL_render(self):
        # Draw the render texture

        self.programTex.bind()

        self.texture = fillTexture2d(self.output, self.texture)
        # self.textureAlpha = fillTexture2d(self.output_alpha, self.textureAlpha)

        glEnable(GL_BLEND)
        glEnable(GL_TEXTURE_2D)
        glDisable(GL_DEPTH_TEST)

        self.programTex.enableAttributeArray("position")
        self.programTex.enableAttributeArray("texcoord")
        self.programTex.setAttributeArray("position", self.quadCoord)
        self.programTex.setAttributeArray("texcoord", self.quadCoordTex)

        self.programTex.setUniformValue("is_mode_black", self._background_mode_black)
        glActiveTexture(GL_TEXTURE0)
        glBindTexture(GL_TEXTURE_2D, self.texture)
        self.programTex.setUniformValue("texture", 0)

        glActiveTexture(GL_TEXTURE1)
        glBindTexture(GL_TEXTURE_2D, self.textureAlpha)
        self.programTex.setUniformValue("texture_alpha", 1)

        glActiveTexture(GL_TEXTURE2)
        glBindTexture(GL_TEXTURE_2D, self.texture_LUT)
        self.programTex.setUniformValue("texture_LUT", 2)

        glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA)

        glDrawArrays(GL_TRIANGLES, 0, len(self.quadCoord))

    def _paintGL_slice(self):
        # draw the slice
        self.programSlice.bind()
        self.programSlice.setUniformValue("mvpMatrix", QtGui.QMatrix4x4(*self._mat_modelviewproject.flatten()))

        self.programSlice.setUniformValue("is_mode_black", self._background_mode_black)
        self.programSlice.enableAttributeArray("position")

        pos, dim = self.transform.slicePos, self.transform.sliceDim

        coords = slice_coords(1. * pos / self.dataModel.size()[2 - dim + 1], dim)

        texcoords = [[0., 0.], [1, 0.], [1., 1.],
                     [1., 1.], [0., 1.], [0., 0.]]

        self.programSlice.setAttributeArray("position", coords)
        self.programSlice.setAttributeArray("texcoord", texcoords)

        self.textureSlice = fillTexture2d(self.sliceOutput, self.textureSlice)

        glActiveTexture(GL_TEXTURE0)
        glBindTexture(GL_TEXTURE_2D, self.textureSlice)
        self.programSlice.setUniformValue("texture", 0)

        glActiveTexture(GL_TEXTURE1)
        glBindTexture(GL_TEXTURE_2D, self.texture_LUT)
        self.programSlice.setUniformValue("texture_LUT", 1)

        glDrawArrays(GL_TRIANGLES, 0, len(coords))

    def _paintGL_box(self):

        glEnable(GL_BLEND)
        # Draw the cube
        self.programCube.bind()
        self.programCube.setUniformValue("mvpMatrix", QtGui.QMatrix4x4(*self._mat_modelviewproject.flatten()))
        self.programCube.enableAttributeArray("position")

        if self._background_mode_black:
            self.programCube.setUniformValue("color",
                                             QtGui.QVector4D(1, 1, 1, 0.6))
        else:
            self.programCube.setUniformValue("color",
                                             QtGui.QVector4D(0, 0, 0, 0.6))

        self.programCube.setAttributeArray("position", self.cubeCoords)

        glActiveTexture(GL_TEXTURE0)
        glBindTexture(GL_TEXTURE_2D, self.textureAlpha)
        self.programCube.setUniformValue("texture_alpha", 0)

        glEnable(GL_DEPTH_TEST)
        # glBlendFunc(GL_ONE_MINUS_SRC_ALPHA, GL_SRC_ALPHA)
        glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA)

        glDrawArrays(GL_LINES, 0, len(self.cubeCoords))

        glDisable(GL_DEPTH_TEST)

    def _paintGL_mesh(self, mesh, vbo_vertices, vbo_normals, vbo_indices):
        """
        paint a mesh (which has all the coordinates and colors in it
        """
        glEnable(GL_DEPTH_TEST)
        glDisable(GL_BLEND)

        prog = self.programMeshLight
        prog.bind()
        prog.setUniformValue("mvpMatrix",
                             QtGui.QMatrix4x4(*self._mat_modelviewproject.flatten()))

        prog.setUniformValue("mvMatrix",
                             QtGui.QMatrix4x4(*self._mat_modelview.flatten()))

        prog.setUniformValue("normMatrix",
                             QtGui.QMatrix4x4(*self._mat_normal.flatten()))

        if mesh.light:
            prog.setUniformValue("light",
                                 QtGui.QVector3D(*mesh.light))
            prog.setUniformValue("light_components",
                                 QtGui.QVector3D(.2, .5, .3))
        else:
            prog.setUniformValue("light",
                                 QtGui.QVector3D(0, 0, 0))
            prog.setUniformValue("light_components",
                                 QtGui.QVector3D(1., 0, 0))

        if not mesh.facecolor is None:
            r, g, b = mesh.facecolor[:3]
            a = mesh.alpha
            prog.setUniformValue("color",
                                 QtGui.QVector4D(r, g, b, a))

            prog.enableAttributeArray("position")
            vbo_vertices.bind()

            glVertexAttribPointer(prog.attributeLocation("position"), 3, GL_FLOAT, GL_FALSE, 0, vbo_vertices)

            prog.enableAttributeArray("normal")
            vbo_normals.bind()
            glVertexAttribPointer(prog.attributeLocation("normal"), 3, GL_FLOAT, GL_FALSE, 0, vbo_normals)

            vbo_indices.bind()

            glDrawElements(GL_TRIANGLES, len(vbo_indices.data), GL_UNSIGNED_INT, None)

            vbo_indices.unbind()
            vbo_vertices.unbind()

            glDisable(GL_DEPTH_TEST)
            prog.disableAttributeArray("position")
            prog.disableAttributeArray("normal")

            #
            # if not mesh.edgecolor is None:
            #     r, g, b = mesh.edgecolor
            #     a = mesh.alpha
            #
            #     prog.enableAttributeArray("position")
            #     vbo_vertices.bind()
            #     glVertexAttribPointer(prog.attributeLocation("position"), 2, GL_FLOAT, GL_FALSE, 0, vbo_edges)
            #
            #     prog.setUniformValue("color",
            #                          QtGui.QVector4D(r, g, b, a))
            #
            #     glDrawArrays(GL_LINES, 0, len(mesh.edges))

    def paintGL(self):

        self.makeCurrent()

        if not glCheckFramebufferStatus(GL_FRAMEBUFFER) == GL_FRAMEBUFFER_COMPLETE:
            return

        w = max(self._viewport_width, self._viewport_height)
        # force viewport to always be a square
        glViewport((self._viewport_width - w) // 2, (self._viewport_height - w) // 2, w, w)

        self.clear_canvas()

        self._mat_modelview = self.transform.getModelView()
        self._mat_proj = self.transform.getProjection()
        self._mat_normal = np.linalg.inv(self._mat_modelview).T

        self._mat_modelviewproject = np.dot(self._mat_proj, self._mat_modelview)

        if self.dataModel:

            self.textureAlpha = fillTexture2d(self.output_alpha, self.textureAlpha)

            if self.transform.isBox:
                self._paintGL_box()

            if self.transform.isSlice and self.sliceOutput is not None:
                self._paintGL_slice()

            self._paintGL_render()

        for (m, vbo_verts, vbo_normals, vbo_indices) in self.meshes:
            self._paintGL_mesh(m, vbo_verts, vbo_normals, vbo_indices)

    def render(self):
        logger.debug("render")

        if self.dataModel:

            self.renderer.set_modelView(self.transform.getUnscaledModelView())
            self.renderer.set_projection(self.transform.getProjection())
            self.renderer.set_min_val(self.transform.minVal)

            self.renderer.set_max_val(self.transform.maxVal)
            self.renderer.set_gamma(self.transform.gamma)
            self.renderer.set_alpha_pow(self.transform.alphaPow)

            self.renderer.set_occ_strength(self.transform.occ_strength)
            self.renderer.set_occ_radius(self.transform.occ_radius)
            self.renderer.set_occ_n_points(self.transform.occ_n_points)

            if self.transform.isIso:
                renderMethod = "iso_surface"

            else:
                renderMethod = "max_project"

            self.renderer.render(method=renderMethod, return_alpha=True, numParts=self.NSubrenderSteps, currentPart=(
                                                                                self.renderedSteps * _next_golden(
                                                                                self.NSubrenderSteps)) % self.NSubrenderSteps)
            self.output, self.output_alpha = self.renderer.output, self.renderer.output_alpha

            if self.transform.isSlice:
                if self.transform.sliceDim == 0:
                    out = self.dataModel[self.transform.dataPos][:, :, self.transform.slicePos]
                elif self.transform.sliceDim == 1:
                    out = self.dataModel[self.transform.dataPos][:, self.transform.slicePos, :]
                elif self.transform.sliceDim == 2:
                    out = self.dataModel[self.transform.dataPos][self.transform.slicePos, :, :]

                min_out, max_out = np.amin(out), np.amax(out)
                if max_out > min_out:
                    self.sliceOutput = (1. * (out - min_out) / (max_out - min_out))
                else:
                    self.sliceOutput = np.zeros_like(out)

    # def getFrame(self):
    #     self.render()
    #     self.paintGL()
    #     glFlush()
    #     im = self.grabFrameBuffer()
    #     im = im.convertToFormat(QtGui.QImage.Format_RGB32)
    #
    #     width = im.width()
    #     height = im.height()
    #
    #     ptr = im.bits()
    #     ptr.setsize(im.byteCount())
    #     arr = np.array(ptr).reshape(height, width, 4)  # Copies the data
    #     return arr[..., [2, 1, 0, 3]].copy()

    def saveFrame(self, fName, with_alpha=False):
        """FIXME: scaling behaviour still hast to be implemented (e.g. after setGamma)"""
        logger.info("saving frame as %s", fName)

        # has to be png

        name, ext = os.path.splitext(fName)
        if ext != ".png":
            fName = name + ".png"

        self.render()
        self.paintGL()
        glFlush()
        im = self.grabFrameBuffer(withAlpha=with_alpha)
        im.save(fName)

    def onRenderTimer(self):
        # if self.renderUpdate:
        #     self.render()
        #     self.renderUpdate = False
        #     self.updateGL()
        if self.renderedSteps < self.NSubrenderSteps:
            # print ((self.renderedSteps*7)%self.NSubrenderSteps)
            s = time.time()
            self.render()
            logger.debug("time to render:  %.2f" % (1000. * (time.time() - s)))
            self.renderedSteps += 1
            self.updateGL()

    def wheelEvent(self, event):
        """ self.transform.zoom should be within [1,2]"""
        newZoom = self.transform.zoom * 1.2 ** (event.angleDelta().y() / 1000.)
        newZoom = np.clip(newZoom, .4, 3)
        self.transform.setZoom(newZoom)

        logger.debug("newZoom: %s", newZoom)
        # self.refresh()

    def posToVec3(self, x, y, r0=.8, isRot=True):
        x, y = 2. * x / self.width() - 1., 1. - 2. * y / self.width()

        r = np.sqrt(x * x + y * y)
        if r > r0 - 1.e-7:
            x, y = 1. * x * r0 / r, 1. * y * r0 / r
        z = np.sqrt(max(0, r0 ** 2 - x * x - y * y))

        if isRot:
            M = np.linalg.inv(self.transform.quatRot.toRotation3())
            x, y, z = np.dot(M, [x, y, z])

        return x, y, z

    def posToVec2(self, x, y):
        x, y = 2. * x / self.width() - 1., 1. - 2. * y / self.width()
        return x, y

    def mousePressEvent(self, event):
        super(GLWidget, self).mousePressEvent(event)

        if event.buttons() == QtCore.Qt.LeftButton:
            self._x0, self._y0, self._z0 = self.posToVec3(event.x(), event.y())

        if event.buttons() == QtCore.Qt.RightButton:
            (self._x0, self._y0), self._invRotM = self.posToVec2(event.x(), event.y()), linalg.inv(
                self.transform.quatRot.toRotation3())

            # self.setCursor(QtCore.Qt.ClosedHandCursor)

    def mouseReleaseEvent(self, event):
        super(GLWidget, self).mouseReleaseEvent(event)

        # self.setCursor(QtCore.Qt.ArrowCursor)

    def mouseMoveEvent(self, event):

        # c = append(self.cubeCoords,ones(24)[:,newaxis],axis=1)
        # cUser = dot(c,self.finalMat)
        # cUser = cUser[:,:3]/cUser[:,-1,newaxis]
        # print self.finalMat
        # print c[0], cUser[0]
        # Rotation

        if event.buttons() == QtCore.Qt.LeftButton:
            x1, y1, z1 = self.posToVec3(event.x(), event.y())
            logger.debug("mouse position: %s %s %s " % (x1, y1, z1))
            n = np.cross(np.array([self._x0, self._y0, self._z0]), np.array([x1, y1, z1]))
            nnorm = linalg.norm(n)
            if np.abs(nnorm) >= 1.:
                nnorm *= 1. / np.abs(nnorm)
            w = np.arcsin(nnorm)
            n *= 1. / (nnorm + 1.e-10)
            q = Quaternion(np.cos(.5 * w), *(np.sin(.5 * w) * n))
            self.transform.setQuaternion(self.transform.quatRot * q)

        # Translation
        if event.buttons() == QtCore.Qt.RightButton:
            x, y = self.posToVec2(event.x(), event.y())

            dx, dy, foo = np.dot(self._invRotM, [x - self._x0, y - self._y0, 0])

            self.transform.addTranslate(dx, dy, foo)
            self._x0, self._y0 = x, y

        self.refresh()

    def resizeEvent(self, event):
        # enforce each dimension to be divisable by 4 (and so the saved frames)
        super(GLWidget, self).resizeEvent(event)

        size = event.size()
        w, h = size.width(), size.height()
        if not ((w % 4 == 0) and (h % 4 == 0)):
            self.resize(QtCore.QSize((w // 4) * 4, (h // 4) * 4))

    def _enforce_resize(self):
        """ this is to enforce the resizeGL event """
        self.resize(self.width() + 1, self.height())
        self.resize(self.width() - 1, self.height())

    def onScreenNumberChange(self, evt):
        self._enforce_resize()

    def _get_screen_number(self):
        return QtGui.QGuiApplication.instance().desktop().screenNumber(QtGui.QCursor.pos())

    def moveEvent(self, evt):
        current_screen = self._get_screen_number()
        if hasattr(self, "_current_screen") and self._current_screen != current_screen:
            self.onScreenNumberChange(evt)

        self._current_screen = current_screen

def test_sphere():
    from spimagine import DataModel, NumpyData, SpimData, TiffData

    app = QtWidgets.QApplication(sys.argv)

    win = GLWidget(size=QtCore.QSize(500, 500), interpolation="nearest")

    x = np.linspace(-1, 1, 64)
    Z, Y, X = np.meshgrid(x, x, x)
    # R = sqrt(Z**2+Y**2+(X-.35)**2)
    # R2 = sqrt(Z**2+Y**2+(X+.35)**2)

    # d = 100.*exp(-10*R**2)+.0*np.random.normal(0,1.,X.shape)

    # d += 100.*exp(-10*R2**2)+.0*np.random.normal(0,1.,X.shape)

    Ns = 5
    r = .6
    phi = np.linspace(0, 2 * np.pi, Ns + 1)[:-1]
    d = np.zeros_like(X)
    for p in phi:
        d += 100. * np.exp(-10 * (Z ** 2 + (Y - r * np.sin(p)) ** 2 + (X - r * np.cos(p)) ** 2))

    win.setModel(DataModel(NumpyData(d)))

    win.transform.setValueScale(0, 40)

    win.show()

    win.raise_()

    sys.exit(app.exec_())

def test_empty():
    from spimagine import DataModel, NumpyData, SpimData, TiffData

    app = QtWidgets.QApplication(sys.argv)

    win = GLWidget(size=QtCore.QSize(1000, 1000))

    d = np.zeros((800,) * 3, np.float32)

    d[0, 0, 0] = 1.

    win.show()

    win.raise_()

    QtCore.QThread.msleep(1000)

    win.setModel(DataModel(NumpyData(d)))

    sys.exit(app.exec_())

def test_demo():
    from data_model import DataModel, DemoData, SpimData, TiffData, NumpyData

    app = QtWidgets.QApplication(sys.argv)

    win = GLWidget(size=QtCore.QSize(800, 800))

    win.setModel(DataModel(DemoData()))

    win.show()

    win.raise_()

    sys.exit(app.exec_())

def test_demo_simple():
    from spimagine import DataModel, DemoData

    app = QtWidgets.QApplication(sys.argv)

    win = GLWidget(size=QtCore.QSize(800, 800))

    win.setModel(DataModel(DemoData()))
    win.show()

    win.raise_()

    sys.exit(app.exec_())

def test_surface():
    from spimagine import DataModel, DemoData

    app = QtWidgets.QApplication(sys.argv)

    win = GLWidget(size=QtCore.QSize(800, 600))

    win.setModel(DataModel(DemoData()))

    # win.add_surface_sphere((0,0,0), 1., facecolor = (.0,.3,1.,.5),
    #                                 Nphi = 30, Ntheta=20)

    # win.add_mesh(SphericalMesh(r = .8,
    #                            facecolor = (1.,0.,0.),
    #                            #edgecolor = (1.,1.,1.),
    #                            edgecolor = None,
    #                            alpha = .3))

    # win.add_mesh(EllipsoidMesh(rs = (.3,.6,.6),
    #                             pos = (0,0,-.5),
    #
    #                            facecolor = (0.,1.,1.),
    #                            #edgecolor = (1.,1.,1.),
    #                            edgecolor = None,
    #                            alpha = .4))
    #

    win.show()

    win.raise_()

    sys.exit(app.exec_())

if __name__ == '__main__':
    # test_empty()

    test_sphere()

    # test_demo_simple()
    # test_surface()

max_stars_repo_name: maweigert/spimagine

max_issues_repo_name: maweigert/spimagine

max_forks_repo_name: maweigert/spimagine

reason: import numpy,from numpy,import scipy,from scipy


text:
"""Unit tests for utility.py"""

import unittest
import numpy as np
import xarray as xr
from littlenet import neural_net as nn
from littlenet import utility as util

class UtilityTest(unittest.TestCase):

    def test_make_onehot(self):
        int_labels = xr.DataArray(
            np.array([[0, 0], [1, 2]]), dims=('batches', nn.DIM_CASE))
        int_symbols = [0, 1]
        int_onehot = util.make_onehot(int_labels, int_symbols)
        int_expected = np.array([[[1, 0], [1, 0]], [[0, 1], [0, 0]]])
        self.assertDictEqual(dict(int_onehot.sizes), {
                             'batches': int_labels.sizes['batches'],
                             nn.DIM_CASE: int_labels.sizes[nn.DIM_CASE],
                             nn.DIM_LABEL: len(int_symbols)})
        for i, j in zip(int_onehot, int_expected):
            np.testing.assert_array_equal(i, j)

    def test_rotate_images_90_deg(self):
        images = xr.DataArray(np.arange(27).reshape((3, 3, 3)),
                              dims=(nn.DIM_CASE, nn.DIM_Y, nn.DIM_X)).stack(
                                  inputs=(nn.DIM_Y, nn.DIM_X))
        img1 = util.rotate_images_90_deg(images)
        self.assertFalse(img1.equals(images))
        np.testing.assert_allclose(img1.unstack(nn.DIM_IN).coords[nn.DIM_Y], np.arange(3))
        np.testing.assert_allclose(img1[0], [6, 3, 0, 7, 4, 1, 8, 5, 2])
        img2 = util.rotate_images_90_deg(images, num_clockwise=-3)
        np.testing.assert_allclose(img1, img2)
        img1 = util.rotate_images_90_deg(images, num_clockwise=3)
        np.testing.assert_allclose(img1[0], [2, 5, 8, 1, 4, 7, 0, 3, 6])
        self.assertFalse(img1.equals(images))
        img2 = util.rotate_images_90_deg(images, num_clockwise=-1)
        np.testing.assert_allclose(img1, img2)
        img1 = util.rotate_images_90_deg(images, num_clockwise=2)
        np.testing.assert_allclose(img1[0], [8, 7, 6, 5, 4, 3, 2, 1, 0])
        self.assertFalse(img1.equals(images))
        img2 = util.rotate_images_90_deg(images, num_clockwise=-2)
        np.testing.assert_allclose(img1, img2)

    def test_flip_images_on_angle(self):
        images = xr.DataArray(np.arange(9).reshape((3, 3)),
                              dims=(nn.DIM_Y, nn.DIM_X)).stack(inputs=(nn.DIM_Y, nn.DIM_X))
        img1 = util.flip_images_on_angle(images)
        img2 = util.flip_images_on_angle(images, topright_to_bottomleft=True)
        np.testing.assert_allclose(img1, [0, 3, 6, 1, 4, 7, 2, 5, 8])
        np.testing.assert_allclose(img2, [8, 5, 2, 7, 4, 1, 6, 3, 0])
        np.testing.assert_allclose(img1.unstack(nn.DIM_IN).coords[nn.DIM_Y], np.arange(3))

    def test_flip_images(self):
        images = xr.DataArray(np.arange(27).reshape((3, 3, 3)),
                              dims=(nn.DIM_CASE, nn.DIM_Y, nn.DIM_X)).stack(
                                  inputs=(nn.DIM_Y, nn.DIM_X))
        img1 = util.flip_images(images)
        np.testing.assert_allclose(img1.unstack(nn.DIM_IN).coords[nn.DIM_Y], np.arange(3))
        np.testing.assert_allclose(img1[0], [6, 7, 8, 3, 4, 5, 0, 1, 2])
        self.assertFalse(img1.equals(images))
        img2 = util.flip_images(img1)
        np.testing.assert_allclose(img2, images)
        img1 = util.flip_images(images, dim=nn.DIM_X)
        self.assertFalse(img1.equals(images))
        img2 = util.flip_images(img1, dim=nn.DIM_X)
        np.testing.assert_allclose(img2, images)

    def test_quarter_images(self):
        images = xr.DataArray(np.arange(16).reshape((4, 4)),
                              dims=(nn.DIM_Y, nn.DIM_X)).stack(inputs=(nn.DIM_Y, nn.DIM_X))
        img1 = util.quarter_images(images)
        np.testing.assert_allclose(img1.unstack(nn.DIM_IN).coords[nn.DIM_Y], np.arange(4))
        np.testing.assert_allclose(
            img1, [10, 11, 8, 9, 14, 15, 12, 13, 2, 3, 0, 1, 6, 7, 4, 5])

    def test_compose_params(self):
        params1 = ['a', 'b']
        params2 = [0, 1, 2]
        params3 = [30, 40, 50, 60]
        params = util.compose_params(params1, params2, params3)
        np.testing.assert_allclose(
            [len(params1) * len(params2) * len(params3), 3], params.shape)
        np.testing.assert_array_equal(params[0], ['a', 0, 30])
        np.testing.assert_array_equal(params[-1], ['b', 2, 60])

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

ext: py

max_stars_repo_name: devinkwok/littlenet

max_issues_repo_name: devinkwok/littlenet

max_forks_repo_name: devinkwok/littlenet

reason: import numpy


text:
import numpy as np
import miepy

def get_tmatrix(rad, dS, eps, eps_m, wavelength, lmax):
    _, Ntheta, Nphi = dS.shape
    theta = np.linspace(0, np.pi, Ntheta)
    phi = np.linspace(0, 2*np.pi, Nphi)
    THETA, PHI = np.meshgrid(theta, phi, indexing='ij')
    dS = miepy.coordinates.vec_cart_to_sph(dS, THETA, PHI)

    rmax = lmax*(lmax+2)
    k1 = 2*np.pi/wavelength
    k2 = 2*np.pi*np.sqrt(eps)/wavelength

    class ingoing:
        M_k2  = np.empty((rmax,3) + THETA.shape, dtype=complex)
        N_k2  = np.empty((rmax,3) + THETA.shape, dtype=complex)
        M_k1  = np.empty((rmax,3) + THETA.shape, dtype=complex)
        N_k1  = np.empty((rmax,3) + THETA.shape, dtype=complex)

    class outgoing:
        M_k2 = np.empty((rmax,3) + THETA.shape, dtype=complex)
        N_k2 = np.empty((rmax,3) + THETA.shape, dtype=complex)
        M_k1 = np.empty((rmax,3) + THETA.shape, dtype=complex)
        N_k1 = np.empty((rmax,3) + THETA.shape, dtype=complex)

    for r,n,m in miepy.mode_indices(lmax):
        Emn = miepy.vsh.Emn(m, n)

        Nfunc, Mfunc = miepy.vsh.VSH(n, m, miepy.vsh_mode.incident)
        ingoing.M_k2[r] = Emn*Mfunc(rad, THETA, PHI, k2)
        ingoing.N_k2[r] = Emn*Nfunc(rad, THETA, PHI, k2)
        ingoing.M_k1[r] = Emn*Mfunc(rad, THETA, PHI, k1)
        ingoing.N_k1[r] = Emn*Nfunc(rad, THETA, PHI, k1)

        Nfunc, Mfunc = miepy.vsh.VSH(n, m, miepy.vsh_mode.outgoing)
        outgoing.M_k2[r] = Emn*Mfunc(rad, THETA, PHI, k2)
        outgoing.N_k2[r] = Emn*Nfunc(rad, THETA, PHI, k2)
        outgoing.M_k1[r] = Emn*Mfunc(rad, THETA, PHI, k1)
        outgoing.N_k1[r] = Emn*Nfunc(rad, THETA, PHI, k1)

    def get_Q(p, q):
        Q = np.zeros([2,rmax,2,rmax], dtype=complex)
        p_modes = ingoing if p == 1 else outgoing
        q_modes = ingoing if q == 1 else outgoing

        for r1,n1,m1 in miepy.mode_indices(lmax):
            M1_p_k2 = p_modes.M_k2[r1]
            N1_p_k2 = p_modes.N_k2[r1]
            M1_p_k1 = p_modes.M_k1[r1]
            N1_p_k1 = p_modes.N_k1[r1]

            for r2,n2,m2 in miepy.mode_indices(lmax):
                M2_q_k2 = q_modes.M_k2[r2]
                N2_q_k2 = q_modes.N_k2[r2]
                M2_q_k1 = q_modes.M_k1[r2]
                N2_q_k1 = q_modes.N_k1[r2]

                # factor = (-1)**(m2 - m1)
                factor = 1

                integrand = np.sum(np.cross(dS, M2_q_k2, axis=0)*N1_p_k1, axis=0) \
                            + np.sqrt(eps/eps_m)*np.sum(np.cross(dS, N2_q_k2, axis=0)*M1_p_k1, axis=0)
                integrand *= 1j*k1**2/np.pi*factor
                Q[1,r1,1,r2] = miepy.vsh.misc.trapz_2d(theta, phi, integrand)

                integrand = np.sum(np.cross(dS, N2_q_k2, axis=0)*N1_p_k1, axis=0) \
                            + np.sqrt(eps/eps_m)*np.sum(np.cross(dS, M2_q_k2, axis=0)*M1_p_k1, axis=0)
                integrand *= 1j*k1**2/np.pi*factor
                Q[1,r1,0,r2] = miepy.vsh.misc.trapz_2d(theta, phi, integrand)

                integrand = np.sum(np.cross(dS, M2_q_k2, axis=0)*M1_p_k1, axis=0) \
                            + np.sqrt(eps/eps_m)*np.sum(np.cross(dS, N2_q_k2, axis=0)*N1_p_k1, axis=0)
                integrand *= 1j*k1**2/np.pi*factor
                Q[0,r1,1,r2] = miepy.vsh.misc.trapz_2d(theta, phi, integrand)

                integrand = np.sum(np.cross(dS, N2_q_k2, axis=0)*M1_p_k1, axis=0) \
                            + np.sqrt(eps/eps_m)*np.sum(np.cross(dS, M2_q_k2, axis=0)*N1_p_k1, axis=0)
                integrand *= 1j*k1**2/np.pi*factor
                Q[0,r1,0,r2] = miepy.vsh.misc.trapz_2d(theta, phi, integrand)

        return Q

    Q31 = -get_Q(3, 1)
    Q11 = get_Q(1, 1)

    T = -np.einsum('aibj,bjck->aick', Q11, np.linalg.tensorinv(Q31))

    return T

nm = 1e-9

lmax = 4

wavelength = 600*nm
radius = 60*nm
eps = 4

material = miepy.constant_material(eps)
sphere = miepy.sphere([0,0,0], radius, material)
T = sphere.compute_tmatrix(lmax, wavelength, 1)
# print(T[0,:,0,0])

theta = np.linspace(0, np.pi, 80)
phi = np.linspace(0, 2*np.pi, 80)
THETA, PHI = np.meshgrid(theta, phi, indexing='ij')
rhat, that, phat = miepy.coordinates.sph_basis_vectors(THETA, PHI)
dS = rhat*np.sin(THETA)*radius**2

# T = get_tmatrix(radius, dS, eps, 1, wavelength, lmax)
# print(T[0,:,0,0])

def ellipsoid_dS(a, b, c, theta, phi):
    rad = 1/np.sqrt(np.sin(theta)**2*np.cos(phi)**2/a**2 + np.sin(theta)**2*np.sin(phi)**2/b**2 + np.cos(theta)**2/c**2)
    dr_theta = -rad**3*(np.sin(theta)*np.cos(theta)*(np.cos(phi)**2/a**2 + np.sin(phi)**2/b**2 - 1/c**2))
    dr_phi = -rad**3*np.sin(phi)*np.cos(phi)*np.sin(theta)*(1/b**2 - 1/a**2)

    rhat, that, phat = miepy.coordinates.sph_basis_vectors(theta, phi)
    sigma = (rhat - 1/rad*dr_theta*that - 1/rad*dr_phi*phat)*rad**2*np.sin(theta)

    return rad, sigma

rad, dS = ellipsoid_dS(1*radius, 1.3*radius, 1.6*radius, THETA, PHI)
T = get_tmatrix(rad, dS, eps, 1, wavelength, lmax)
print(T[1,:,0,0])
T1 = np.copy(T)

spheroid = miepy.spheroid([0,0,0], radius, 1.3*radius, material)
spheroid = miepy.ellipsoid([0,0,0], radius, 1.3*radius, 1.6*radius, material)
T = spheroid.compute_tmatrix(lmax, wavelength, 1)
print(T[1,:,0,0])
T2 = np.copy(T)

import matplotlib.pyplot as plt
rmax = lmax*(lmax+2)
T1 = np.reshape(T1, [2*rmax, 2*rmax])
T2 = np.reshape(T2, [2*rmax, 2*rmax])

y = np.abs(T1-T2)
y /= np.max(y)
plt.pcolormesh(y, vmax=.1)
plt.gca().set_aspect('equal')
plt.colorbar()

plt.show()

max_stars_repo_name: johnaparker/MiePy

max_issues_repo_name: johnaparker/MiePy

max_forks_repo_name: johnaparker/MiePy


text:
"""Using Auto Encoder with Outlier Detection."""
# Author: Yue Zhao <zhaoy@cmu.edu>
# License: BSD 2 clause

from __future__ import division, print_function

import numpy as np
from sklearn.preprocessing import StandardScaler
from sklearn.utils import check_array
from sklearn.utils.validation import check_is_fitted

from ..utils.stat_models import pairwise_distances_no_broadcast
from ..utils.utility import check_parameter
from .base import BaseDetector
from .base_dl import _get_tensorflow_version

# if tensorflow 2, import from tf directly
if _get_tensorflow_version() == 1:
    from keras.models import Sequential
    from keras.layers import Dense, Dropout
    from keras.regularizers import l2
    from keras.losses import mean_squared_error
else:
    from tensorflow.keras.models import Sequential
    from tensorflow.keras.layers import Dense, Dropout
    from tensorflow.keras.regularizers import l2
    from tensorflow.keras.losses import mean_squared_error

# noinspection PyUnresolvedReferences,PyPep8Naming,PyTypeChecker
class AutoEncoder(BaseDetector):
    """Auto Encoder (AE) is a type of neural networks for learning useful data
    representations unsupervisedly. Similar to PCA, AE could be used to detect
    outlying objects in the data by calculating the reconstruction errors. See
    :cite:`aggarwal2015outlier` Chapter 3 for details.

    Parameters
    ----------
    hidden_neurons : list, optional (default=[64, 32, 32, 64])
        The number of neurons per hidden layers.

    hidden_activation : str, optional (default='relu')
        Activation function to use for hidden layers.
        All hidden layers are forced to use the same type of activation.
        See https://keras.io/activations/

    output_activation : str, optional (default='sigmoid')
        Activation function to use for output layer.
        See https://keras.io/activations/

    loss : str or obj, optional (default=keras.losses.mean_squared_error)
        String (name of objective function) or objective function.
        See https://keras.io/losses/

    optimizer : str, optional (default='adam')
        String (name of optimizer) or optimizer instance.
        See https://keras.io/optimizers/

    epochs : int, optional (default=100)
        Number of epochs to train the model.

    batch_size : int, optional (default=32)
        Number of samples per gradient update.

    dropout_rate : float in (0., 1), optional (default=0.2)
        The dropout to be used across all layers.

    l2_regularizer : float in (0., 1), optional (default=0.1)
        The regularization strength of activity_regularizer
        applied on each layer. By default, l2 regularizer is used. See
        https://keras.io/regularizers/

    validation_size : float in (0., 1), optional (default=0.1)
        The percentage of data to be used for validation.

    preprocessing : bool, optional (default=True)
        If True, apply standardization on the data.

    verbose : int, optional (default=1)
        Verbosity mode.

        - 0 = silent
        - 1 = progress bar
        - 2 = one line per epoch.

        For verbose >= 1, model summary may be printed.

    random_state : random_state: int, RandomState instance or None, optional
        (default=None)
        If int, random_state is the seed used by the random
        number generator; If RandomState instance, random_state is the random
        number generator; If None, the random number generator is the
        RandomState instance used by `np.random`.

    contamination : float in (0., 0.5), optional (default=0.1)
        The amount of contamination of the data set, i.e.
        the proportion of outliers in the data set. When fitting this is used
        to define the threshold on the decision function.

    Attributes
    ----------
    encoding_dim_ : int
        The number of neurons in the encoding layer.

    compression_rate_ : float
        The ratio between the original feature and
        the number of neurons in the encoding layer.

    model_ : Keras Object
        The underlying AutoEncoder in Keras.

    history_: Keras Object
        The AutoEncoder training history.

    decision_scores_ : numpy array of shape (n_samples,)
        The outlier scores of the training data.
        The higher, the more abnormal. Outliers tend to have higher
        scores. This value is available once the detector is
        fitted.

    threshold_ : float
        The threshold is based on ``contamination``. It is the
        ``n_samples * contamination`` most abnormal samples in
        ``decision_scores_``. The threshold is calculated for generating
        binary outlier labels.

    labels_ : int, either 0 or 1
        The binary labels of the training data. 0 stands for inliers
        and 1 for outliers/anomalies. It is generated by applying
        ``threshold_`` on ``decision_scores_``.
    """

    def __init__(self,
                 hidden_neurons=None,
                 hidden_activation='relu',
                 output_activation='sigmoid',
                 loss=mean_squared_error,
                 optimizer='adam',
                 epochs=100,
                 batch_size=32,
                 dropout_rate=0.2,
                 l2_regularizer=0.1,
                 validation_size=0.1,
                 preprocessing=True,
                 verbose=1,
                 random_state=None,
                 contamination=0.1):
        super(AutoEncoder, self).__init__(contamination=contamination)
        self.hidden_neurons = hidden_neurons
        self.hidden_activation = hidden_activation
        self.output_activation = output_activation
        self.loss = loss
        self.optimizer = optimizer
        self.epochs = epochs
        self.batch_size = batch_size
        self.dropout_rate = dropout_rate
        self.l2_regularizer = l2_regularizer
        self.validation_size = validation_size
        self.preprocessing = preprocessing
        self.verbose = verbose
        self.random_state = random_state

        # default values
        if self.hidden_neurons is None:
            self.hidden_neurons = [64, 32, 32, 64]

        # Verify the network design is valid
        if not self.hidden_neurons == self.hidden_neurons[::-1]:
            print(self.hidden_neurons)
            raise ValueError('Hidden units should be symmetric')

        self.hidden_neurons_ = self.hidden_neurons

        check_parameter(
            dropout_rate, 0, 1, param_name='dropout_rate', include_left=True)

    def _build_model(self):
        model = Sequential()
        # Input layer
        model.add(
            Dense(
                self.hidden_neurons_[0],
                activation=self.hidden_activation,
                input_shape=(self.n_features_, ),
                activity_regularizer=l2(self.l2_regularizer)))
        model.add(Dropout(self.dropout_rate))

        # Additional layers
        for i, hidden_neurons in enumerate(self.hidden_neurons_, 1):
            model.add(
                Dense(
                    hidden_neurons,
                    activation=self.hidden_activation,
                    activity_regularizer=l2(self.l2_regularizer)))
            model.add(Dropout(self.dropout_rate))

        # Output layers
        model.add(
            Dense(
                self.n_features_,
                activation=self.output_activation,
                activity_regularizer=l2(self.l2_regularizer)))

        # Compile model
        model.compile(loss=self.loss, optimizer=self.optimizer)
        if self.verbose >= 1:
            print(model.summary())
        return model

    # noinspection PyUnresolvedReferences
    def fit(self, X, y=None):
        """Fit detector. y is ignored in unsupervised methods.

        Parameters
        ----------
        X : numpy array of shape (n_samples, n_features)
            The input samples.

        y : Ignored
            Not used, present for API consistency by convention.

        Returns
        -------
        self : object
            Fitted estimator.
        """
        # validate inputs X and y (optional)
        X = check_array(X)
        self._set_n_classes(y)

        # Verify and construct the hidden units
        self.n_samples_, self.n_features_ = X.shape[0], X.shape[1]

        # Standardize data for better performance
        if self.preprocessing:
            self.scaler_ = StandardScaler()
            X_norm = self.scaler_.fit_transform(X)
        else:
            X_norm = np.copy(X)

        # Shuffle the data for validation as Keras do not shuffling for
        # Validation Split
        np.random.shuffle(X_norm)

        # Validate and complete the number of hidden neurons
        if np.min(self.hidden_neurons) > self.n_features_:
            raise ValueError('The number of neurons should not exceed '
                             'the number of features')
        self.hidden_neurons_.insert(0, self.n_features_)

        # Calculate the dimension of the encoding layer & compression rate
        self.encoding_dim_ = np.median(self.hidden_neurons)
        self.compression_rate_ = self.n_features_ // self.encoding_dim_

        # Build AE model & fit with X
        self.model_ = self._build_model()
        self.history_ = self.model_.fit(
            X_norm,
            X_norm,
            epochs=self.epochs,
            batch_size=self.batch_size,
            shuffle=True,
            validation_split=self.validation_size,
            verbose=self.verbose).history
        # Reverse the operation for consistency
        self.hidden_neurons_.pop(0)
        # Predict on X itself and calculate the reconstruction error as
        # the outlier scores. Noted X_norm was shuffled has to recreate
        if self.preprocessing:
            X_norm = self.scaler_.transform(X)
        else:
            X_norm = np.copy(X)

        pred_scores = self.model_.predict(X_norm)
        self.decision_scores_ = pairwise_distances_no_broadcast(
            X_norm, pred_scores)
        self._process_decision_scores()
        return self

    def decision_function(self, X):
        """Predict raw anomaly score of X using the fitted detector.

        The anomaly score of an input sample is computed based on different
        detector algorithms. For consistency, outliers are assigned with
        larger anomaly scores.

        Parameters
        ----------
        X : numpy array of shape (n_samples, n_features)
            The training input samples. Sparse matrices are accepted only
            if they are supported by the base estimator.

        Returns
        -------
        anomaly_scores : numpy array of shape (n_samples,)
            The anomaly score of the input samples.
        """
        check_is_fitted(self, ['model_', 'history_'])
        X = check_array(X)

        if self.preprocessing:
            X_norm = self.scaler_.transform(X)
        else:
            X_norm = np.copy(X)

        # Predict on X and return the reconstruction errors
        pred_scores = self.model_.predict(X_norm)
        return pairwise_distances_no_broadcast(X_norm, pred_scores)

max_stars_repo_name: jianzhnie/AutoTabular

max_issues_repo_name: Fanxingye/Autotabular

max_forks_repo_name: Fanxingye/Autotabular


text:
'''
VBLearnAlg.py

Implementation of both EM and VB for bnpy models

Notes
-------
Essentially, EM and VB are the same iterative *algorithm*,
repeating the steps of a monotonic increasing objective function until convergence.

EM recovers the parameters for a *point-estimate* of quantities of interest
while VB learns the parameters of an approximate *distribution* over quantities of interest

For more info, see the documentation [TODO]
'''
import numpy as np
from collections import defaultdict
from LearnAlg import LearnAlg

class VBLearnAlg( LearnAlg ):

  def __init__( self, **kwargs ):
    ''' Create VBLearnAlg, subtype of generic LearnAlg
    '''
    super(type(self), self).__init__( **kwargs )
    self.BirthLog = list()

  def fit(self, hmodel, Data):
    ''' Run EM/VB learning algorithm, fit global parameters of hmodel to Data
        Returns
        --------
        LP : local params from final pass of Data
        Info : dict of run information, with fields
              evBound : final ELBO evidence bound
              status : str message indicating reason for termination
                        {'converged', 'max passes exceeded'}
    '''
    prevBound = -np.inf
    LP = None
    mergeFlags = dict(doPrecompEntropy=True, doPrecompMergeEntropy=True)
    self.set_start_time_now()
    for iterid in xrange(self.algParams['nLap'] + 1):
      lap = self.algParams['startLap'] + iterid
      self.set_random_seed_at_lap(lap)

      # M step
      if iterid > 0:
        hmodel.update_global_params(SS)

      if self.hasMove('birth') and iterid > 1:
        hmodel, LP = self.run_birth_move(hmodel, Data, SS, LP, iterid)

      # E step
      LP = hmodel.calc_local_params(Data, LP, **self.algParamsLP)

      # Suff Stat step
      if self.hasMove('merge'):
        SS = hmodel.get_global_suff_stats(Data, LP, **mergeFlags)
      else:
        SS = hmodel.get_global_suff_stats(Data, LP)

      # ELBO calculation
      evBound = hmodel.calc_evidence(Data, SS, LP)

      # Attempt merge move
      if self.hasMove('merge'):
        hmodel, SS, LP, evBound = self.run_merge_move(
                                          hmodel, Data, SS, LP, evBound)

      # Save and display progress
      self.add_nObs(Data.nObs)
      self.save_state(hmodel, iterid, lap, evBound)
      self.print_state(hmodel, iterid, lap, evBound)
      self.eval_custom_func(hmodel, iterid, float(iterid))
      # Check for Convergence!
      #  report warning if bound isn't increasing monotonically
      isConverged = self.verify_evidence( evBound, prevBound )
      if isConverged:
        break
      prevBound = evBound

    #Finally, save, print and exit
    if isConverged:
      status = "converged."
    else:
      status = "max passes thru data exceeded."
    self.save_state(hmodel,iterid, lap, evBound, doFinal=True)
    self.print_state(hmodel,iterid, lap, evBound, doFinal=True, status=status)

    return LP, self.buildRunInfo(evBound, status, nLap=lap)

  ########################################################### Birth Move
  ###########################################################
  def run_birth_move(self, hmodel, Data, SS, LP, lap):
    ''' Run birth move on hmodel
    '''
    import BirthMove # avoid circular import
    self.BirthLog = list()
    if not self.do_birth_at_lap(lap):
      return hmodel, LP

    kbirth = BirthMove.select_birth_component(SS,
                          randstate=self.PRNG,
                          **self.algParams['birth'])

    TargetData = BirthMove.subsample_data(Data, LP, kbirth,
                          randstate=self.PRNG,
                          **self.algParams['birth'])

    hmodel, SS, MoveInfo = BirthMove.run_birth_move(
                 hmodel, TargetData, SS, ktarget=kbirth, randstate=self.PRNG,
                 **self.algParams['birth'])
    self.print_msg(MoveInfo['msg'])
    self.BirthLog.extend(MoveInfo['birthCompIDs'])
    LP = None
    return hmodel, LP

  ########################################################### Merge Move
  ###########################################################
  def run_merge_move(self, hmodel, Data, SS, LP, evBound):
    ''' Run merge move on hmodel
    '''
    import MergeMove
    excludeList = list()
    excludePairs = defaultdict(lambda:set())
    nMergeAttempts = self.algParams['merge']['mergePerLap']
    trialID = 0
    while trialID < nMergeAttempts:

      # Synchronize contents of the excludeList and excludePairs
      # So that comp excluded in excludeList (due to accepted merge)
      #  is automatically contained in the set of excluded pairs
      for kx in excludeList:
        for kk in excludePairs:
          excludePairs[kk].add(kx)
          excludePairs[kx].add(kk)

      for kk in excludePairs:
        if len(excludePairs[kk]) > hmodel.obsModel.K - 2:
          if kk not in excludeList:
            excludeList.append(kk)

      if len(excludeList) > hmodel.obsModel.K - 2:
        break # when we don't have any more comps to merge

      if len(self.BirthLog) > 0:
        kA = self.BirthLog.pop()
        if kA in excludeList:
          continue
      else:
        kA = None

      oldEv = hmodel.calc_evidence(SS=SS)
      hmodel, SS, evBound, MoveInfo = MergeMove.run_merge_move(
                 hmodel, Data, SS, evBound, kA=kA, randstate=self.PRNG,
                 excludeList=excludeList, excludePairs=excludePairs,
                  **self.algParams['merge'])
      newEv = hmodel.calc_evidence(SS=SS)

      trialID += 1
      self.print_msg(MoveInfo['msg'])
      if 'kA' in MoveInfo and 'kB' in MoveInfo:
        kA = MoveInfo['kA']
        kB = MoveInfo['kB']
        excludePairs[kA].add(kB)
        excludePairs[kB].add(kA)

      if MoveInfo['didAccept']:
        assert newEv > oldEv
        kA = MoveInfo['kA']
        kB = MoveInfo['kB']
        # Adjust excludeList since components kB+1, kB+2, ... K
        #  have been shifted down by one due to removal of kB
        for kk in range(len(excludeList)):
          if excludeList[kk] > kB:
            excludeList[kk] -= 1

        # Exclude new merged component kA from future attempts
        #  since precomputed entropy terms involving kA aren't good
        excludeList.append(kA)

        # Adjust excluded pairs to remove kB and shift down kB+1, ... K
        newExcludePairs = defaultdict(lambda:set())
        for kk in excludePairs.keys():
          ksarr = np.asarray(list(excludePairs[kk]))
          ksarr[ksarr > kB] -= 1
          if kk > kB:
            newExcludePairs[kk-1] = set(ksarr)
          elif kk < kB:
            newExcludePairs[kk] = set(ksarr)
        excludePairs = newExcludePairs

        # Update LP to reflect this merge!
        LPkeys = LP.keys()
        keepLPkeys = hmodel.allocModel.get_keys_for_memoized_local_params()

        for key in LPkeys:
          if key in keepLPkeys:
            LP[key][:, kA] = LP[key][:, kA] + LP[key][:, kB]
            LP[key] = np.delete(LP[key], kB, axis=1)
    return hmodel, SS, LP, evBound


text:
'''
Copyright (c) 2020, NVIDIA CORPORATION. All rights reserved.

NVIDIA CORPORATION and its licensors retain all intellectual property
and proprietary rights in and to this software, related documentation
and any modifications thereto. Any use, reproduction, disclosure or
distribution of this software and related documentation without an express
license agreement from NVIDIA CORPORATION is strictly prohibited.
'''
import argparse
import numpy as np
import cv2    # for 2d minimum rotated rectangle fitting and similarity color assignment in LAB
import quaternion    # for convert cv rotation matrix to quaternion

from isaac import Application, Codelet
'''
Contains PyCodelets and helper functions to compute 3d pose of colored cubes from RGB-D image
from a near-top-down camera. It uses a geometry method based on surflets generated by the
superpixel algorithm to compute 3d poses.

This is intended as an example for fast-prototyping using Isaac's python API.
'''

ELEMENT_TYPE_INT = 'int32'
COLOR_TEMPLATES = {
    "red": [128, 0, 0],
    "green": [0, 128, 0],
    "blue": [0, 0, 128],
    "yellow": [128, 128, 0]
}
CUBE_LENGTH = 0.048
PLANE_DISTANCE_THRESHOLD = 0.005

def fit_plane_to_points(data):
    '''
    Fits a plane normal to a set of 3d points using SVD, and return the normal and offset of the
    plane model.

    Input:
      data: float32 2D (Nx3) numpy array of the input 3D point cloud

    Output:
      normal, offset of the plane model
    '''
    centroid = np.mean(data, axis=0)
    X = data - centroid
    u, s, vh = np.linalg.svd(X, full_matrices=True)
    normal = vh[-1, :]
    offset = -np.dot(normal, centroid)
    return normal, offset

def detect_plane(data, distance_threshold=0.005, max_iterations=10):
    '''
    Given a point cloud, return a mask for inliers and the plane model. This method iteractively
    fits plane to inliers points until the set of inliers are unchanged

    Input:
      data: float32 2D (Nx3) numpy array of the input 3D point cloud
      distance_threshold: threshold for point to plane distance for inliers

    Output:
      inliers: boolean numpy array of size N, mask for inliers
      model: plane model of (normal, offset)
    '''
    new_inliers = np.ones((data.shape[0]), dtype=bool)
    for i in range(max_iterations):
        inliers = new_inliers
        if np.sum(inliers) < 3:
            break
        normal, offset = fit_plane_to_points(data[inliers, :])
        dist = np.abs(np.matmul(data, normal) + offset)
        new_inliers = dist <= distance_threshold
        if np.all(inliers == new_inliers):
            break
    return inliers, np.append(normal, offset)

def estimate_pose(points,
                  block_size,
                  distance_threshold,
                  fit_rotation=True,
                  model=None,
                  min_size_ratio=0.5,
                  max_size_ratio=1.5):
    '''
    Estimates the pose of a block from a near-top-down viewpoint by first fitting the top plane
    of the block, then project points to the top plane and use in-plane minimum rectangle to
    estimate centroid position and in-plane rotation, and finally combine the top plane model
    with in-plane pose to output 3d pose

    Inputs:
      points: Nx3 numpy array of point clouds that belongs to the block
      block_size: edge length of the block
      distance_threshold: point-to-plane distance threshold for fitting the top surface plane
      fit_rotation: if True, use minimal rectangle to estimate in-plane rotation. Otherwise
                    default to 0 rotatino angle
      model: 3d plane model for the top surface. if None, estimate the plane model by fitting to
             points
      min(max)_size_ratio: if using the minimal rectangle, the ratio of the fitted bounding box
                           size to the block_size must fall between this two

    Output:
      centroid: 3d position of the block's center
      Q: rotation quaternion
    '''
    # get plane model for top plane of the block cluster
    if not model:
        _, model = detect_plane(points, distance_threshold)
        if np.count_nonzero(np.isnan(model)) > 0:
            raise ValueError("Unable to fit plane model")

    # get 3d axis of the block, making sure the offset is positive, e.g., the plane normal is
    # pointing towards the camera
    if model[3] < 0:
        model = [-x for x in model]
    az = model[:3]
    az /= np.linalg.norm(az)
    ax = np.cross(np.array([-1.0, 0.0, 0.0]), az)
    ax /= np.linalg.norm(ax)
    ay = np.cross(az, ax)

    # project to inplane (2d) to get in-plane centroid, extend and rotation
    projected = np.vstack([points.dot(ax), points.dot(ay)]).transpose().astype(np.float32)
    if fit_rotation:
        box_center, box_size, box_angle = cv2.minAreaRect(projected)
        # limit in-plane rotation to [-45, 45]
        while box_angle > 45.0:
            box_angle -= 90.0
        while box_angle < -45.0:
            box_angle += 90.0
        box_angle = np.deg2rad(box_angle)
    else:
        box_center = (np.max(projected, axis=0) + np.min(projected, axis=0)) * 0.5
        box_size = np.max(projected, axis=0) - np.min(projected, axis=0)

    # check box size is reasonable
    max_size = max_size_ratio * block_size
    min_size = min_size_ratio * block_size
    if not (min_size <= box_size[0] <= max_size) or not (min_size <= box_size[1] <= max_size):
        raise ValueError("Box size unexpected")

    # get 3d centroid of box
    projected_z = np.max(points.dot(az))
    centroid = box_center[0] * ax + box_center[1] * ay + (projected_z - block_size * 0.5) * az
    # align the in-plane axis and get the 3d rotation as quaternion
    if fit_rotation:
        ax, ay = ax * np.cos(box_angle) + ay * np.sin(box_angle), ay * np.cos(box_angle) - ax * \
                 np.sin(box_angle)
    M = np.vstack([ax, ay, az]).transpose()
    Q = quaternion.from_rotation_matrix(M)
    return centroid, Q

class ClusterPoseEstimation(Codelet):
    '''
    Estimates the 3d pose of blocks given their clusters of surflets, using estimation of the top
    plane of the block and in-plane projections of the ponits.
    '''

    def start(self):
        # Receives surflets
        self.rx_surflets = self.isaac_proto_rx("CompositeProto", "surflets")
        # Receives cluster assignment
        self.rx_assignment = self.isaac_proto_rx("TensorProto", "assignment")
        # Receives cluster assignment
        self.rx_plane_model = self.isaac_proto_rx("PlaneProto", "plane_model")
        # Publishes detections3 (prediction + pose)
        self.tx_output_poses = self.isaac_proto_tx("Detections3Proto", "output_poses")
        # Publishes filtered assignment
        self.tx_output_assignment = self.isaac_proto_tx("TensorProto", "output_assignment")

        self.tick_on_message(self.rx_assignment)
        self.synchronize(self.rx_assignment, self.rx_plane_model, self.rx_surflets)

    def tick(self):
        # parse surflets
        acqtime = self.rx_surflets.message.acqtime
        sfl = self.rx_surflets.message.tensor
        if sfl is None:
            self.log_info("Surflect tensor is None")
            return
        points = sfl.reshape((-1, 11))[:, 2:5].astype(np.float32)

        # get cluster assignment
        clusters = self.rx_assignment.message.tensor
        if clusters is None:
            self.log_info("Assignment tensor is None")
            return
        num_cluster = clusters.shape[0]
        assert points.shape[0] == clusters.shape[1], "Size of surflets do not match size of " \
                                                    "assignment"
        # get plane model
        if self.config.fit_plane:
            plane_normal = self.rx_plane_model.message.proto.normal
            plane_model = [
                plane_normal.x, plane_normal.y, plane_normal.z,
                self.rx_plane_model.message.proto.offset
            ]
        else:
            plane_model = None

        # get pose estimation for each cluster
        valid_clusters = []
        results = []
        for i in range(num_cluster):
            try:
                t, q = estimate_pose(points[clusters[i] > 0], self.config.block_length,
                                     self.config.distance_threhold, self.config.fit_rotation,
                                     plane_model)
                valid_clusters.append(i)
                results.append((t, q))
            except ValueError:
                pass

        num_cluster = len(valid_clusters)
        self.show("cluster", num_cluster)

        # publish updated cluster assignment
        tx_message = self.tx_output_assignment.init()
        tx_message.proto.elementType = ELEMENT_TYPE_INT
        tensor_size = tx_message.proto.init('sizes', 2)
        tensor_size[0] = num_cluster
        tensor_size[1] = clusters.shape[1]
        tx_message.proto.dataBufferIndex = 0
        tx_message.buffers = [clusters[valid_clusters]]
        tx_message.acqtime = acqtime
        self.tx_output_assignment.publish()

        # publish detection3 (output_poses)
        tx_message = self.tx_output_poses.init()
        predictions = tx_message.proto.init("predictions", num_cluster)
        poses = tx_message.proto.init("poses", num_cluster)
        for i in range(num_cluster):
            predictions[i].label = self.config.label
            # estimate poses for each cluster
            t, q = results[i]
            predictions[i].confidence = 1.0
            poses[i].translation.x = t[0].item()
            poses[i].translation.y = t[1].item()
            poses[i].translation.z = t[2].item()
            poses[i].rotation.q.w = q.w
            poses[i].rotation.q.x = q.x
            poses[i].rotation.q.y = q.y
            poses[i].rotation.q.z = q.z
        tx_message.acqtime = acqtime
        self.tx_output_poses.publish()

def rgb2lab(color):
    '''
    Converts array of rgb color (in uint8) to lab color

    Input:
      color: input numpy array of rgb colors in [0, 255], show be shape of [N, 3]
    Output:
      color converted to LAB colorspace
    '''
    num_pixel = int(color.size / 3)
    # resize to rank-3
    lab_cv = cv2.cvtColor(color.reshape((1, num_pixel, 3)), cv2.COLOR_RGB2LAB).reshape(num_pixel, 3)
    # cv2.cvtColor rescales LAB to 255, scale it back to L [0,100[, a/b [-128, 128[
    return np.vstack((lab_cv[:, 0] / 255.0 * 100.0, lab_cv[:, 1] - 128.0, lab_cv[:, 2] - 128.0)).T

def deltaE_cie76(colors_1, colors_2):
    '''
    Computes the CIE76 distance of two set of LAB colors

    Input:
      colors_1: a list of LAB colors
      colors_2: another list of LAB colors
    Output:
      numpy array of CIE76 distance between each pair of colors in colors_1 and colors_2
    '''
    colors_1 = colors_1.reshape(-1, 3)
    colors_2 = colors_2.reshape(-1, 3)
    row, col = colors_1.shape[0], colors_2.shape[0]
    result = np.zeros((row, col))
    for i in range(row):
        x = colors_1[i, :]
        for j in range(col):
            y = colors_2[j, :]
            result[i, j] = np.linalg.norm(x-y)
    # reshape output to lower rank if input colors_1 or colors_2 has only a single color
    if col == 1 and row == 1:
        return result[0, 0]
    if col == 1:
        return result[:, 0]
    if row == 1:
        return result[0, :]
    return result

class ColorMatcher(object):
    '''
    Helper class to query the closest color in a fixed set of color templates to any given input.
    Color distance is computed in the lab color space using deltaE_cie76
    '''

    def __init__(self, color_templates, threshold=60):
        '''
        Inputs:
          color_templates: a dictionary of color name (str) to rgb value (list of three ints)
          threshold: distance between a query color and template color to be considered a match
        '''
        self._color_labels = list(color_templates.keys())
        color_values = [color_templates[c] for c in self._color_labels]
        self._color_values = rgb2lab(np.asarray(color_values, dtype=np.uint8))
        self._threshold = threshold

    def match(self, color):
        '''
        Input:
          color: rgb value of the query color

        Output:
          Name (str) of the matched template color. None if the distance to all template colors are
          greater than the threshold
        '''
        query = rgb2lab(np.asarray(color, dtype=np.uint8))
        delta = deltaE_cie76(self._color_values, query)
        idx = np.argmin(delta)
        if delta[idx] < self._threshold:
            return self._color_labels[idx]
        return None

class ClusterColorAssignment(Codelet):
    '''
    This codelet subfixes detections3 labels with color label by matching the color for the
    corresponding surflet clusters to a set of color templates.
    '''

    def start(self):
        # Receives surflets
        self.rx_surflets = self.isaac_proto_rx("CompositeProto", "surflets")
        # Receives cluster assignment
        self.rx_assignment = self.isaac_proto_rx("TensorProto", "assignment")
        # Receives detections3 (prediction + pose)
        self.rx_input_poses = self.isaac_proto_rx("Detections3Proto", "input_poses")
        # Publishes detections3 (prediction + pose) with relabeld color
        self.tx_output_poses = self.isaac_proto_tx("Detections3Proto", "output_poses")

        self.tick_on_message(self.rx_input_poses)
        self.synchronize(self.rx_input_poses, self.rx_assignment, self.rx_surflets)

        # initialize the color matcher with color templates from config
        self._color_matcher = ColorMatcher(self.config.color_templates)

    def tick(self):
        # parse surflets
        acqtime = self.rx_surflets.message.acqtime
        sfl = self.rx_surflets.message.tensor
        if sfl is None:
            self.log_info("Surflect tensor is None")
            return
        sfl_colors = sfl.reshape((-1, 11))[:, 8:11]

        # get cluster assignment
        clusters = self.rx_assignment.message.tensor
        num_cluster = clusters.shape[0]
        # assign color
        colors = [np.median(sfl_colors[clusters[i] > 0], axis=0) for i in range(num_cluster)]
        color_assignment = [self._color_matcher.match(c) for c in colors]
        indices = [i for i in range(num_cluster) if color_assignment[i] is not None]
        num_valid = len(indices)

        # publish detection3 (output_poses)
        tx_message = self.tx_output_poses.init()
        predictions = tx_message.proto.init("predictions", num_valid)
        poses = tx_message.proto.init("poses", num_valid)
        input_proto = self.rx_input_poses.message.proto
        for i in range(num_valid):
            j = indices[i]
            # assign detection color as subfix to label
            predictions[i].label = "{0}_{1}".format(input_proto.predictions[j].label,
                                                    color_assignment[j])
            predictions[i].confidence = 1.0
            # estimate poses for each cluster
            poses[i] = input_proto.poses[j]
        tx_message.acqtime = acqtime
        self.tx_output_poses.publish()

def create_block_pose_estimation(app,
                                 block_length=CUBE_LENGTH,
                                 distance_threhold=PLANE_DISTANCE_THRESHOLD,
                                 use_refinement=False):
    '''
    Creates the block pose estimation subgraph, adds ClusterPoseEstimation and
    ClusterColorAssignment PyCodelets and connects the appropriate edges with or without refinement
    '''
    app.load(
        filename="apps/samples/pick_and_place/block_pose_estimation.subgraph.json",
        prefix="detection_pose_estimation")
    perception_interface = app.nodes["detection_pose_estimation.interface"]["Subgraph"]

    cluster_pose = app.nodes["detection_pose_estimation.block_detection"].add(
        ClusterPoseEstimation, "ClusterPoseEstimation")
    cluster_pose.config.label = "block"
    cluster_pose.config.block_length = block_length
    cluster_pose.config.fit_rotation = not use_refinement
    cluster_pose.config.distance_threhold = distance_threhold

    color_assignment = app.nodes["detection_pose_estimation.block_detection"].add(
        ClusterColorAssignment, "ClusterColorAssignment")
    color_assignment.config.color_templates = COLOR_TEMPLATES
    if use_refinement:
        pose_refinement_interface = app.nodes[
            'detection_pose_estimation.object_pose_refinement.subgraph']['interface']
        app.connect(pose_refinement_interface, "output_poses", color_assignment, "input_poses")
    else:
        app.nodes['detection_pose_estimation.object_pose_refinement' \
                  '.object_pose_refinement'].config['disable_automatic_start'] = True
        app.connect(cluster_pose, "output_poses", color_assignment, "input_poses")

    return perception_interface

def main(args):
    '''
    Creates and runs block pose estimation app on RGBD image. Input image may come from sim (
    default), cask recording or realsense camera

    To run with sim (sim needs to be running first):
      bazel run apps/samples/pick_and_place:block_pose_estimation
    To see all options:
      bazel run apps/samples/pick_and_place:block_pose_estimation -- --help
    '''
    app = Application(name="block_pose_estimation")
    perception_interface = create_block_pose_estimation(app, use_refinement=args.refinement)

    if args.mode == 'cask':
        assert args.cask is not None, "cask path not specified. use --cask to set."
        # Load replay subgraph and configure interface node
        app.load("packages/cask/apps/replay.subgraph.json", prefix="replay")
        source = app.nodes["replay.interface"].components["output"]
        source.config.cask_directory = args.cask
        source.config.loop = True
    elif args.mode == "realsense":
        # Create realsense camera codelet
        app.load_module("realsense")
        source = app.add("camera").add(app.registry.isaac.RealsenseCamera)
        source.config.rows = 720
        source.config.cols = 1280
        source.config.color_framerate = 15
        source.config.depth_framerate = 15
    elif args.mode == "sim":
        app.load(filename='packages/navsim/apps/navsim_tcp.subgraph.json', prefix='simulation')
        source = app.nodes["simulation.interface"]["output"]

    app.connect(source, "color", perception_interface, "color")
    app.connect(source, "depth", perception_interface, "depth")
    app.connect(source, "color_intrinsics", perception_interface, "intrinsics")

    app.load_module("utils")
    pose_node = app.add("pose")
    pose_injector = pose_node.add(app.registry.isaac.utils.DetectionsToPoseTree)
    pose_injector.config.detection_frame = "camera"
    app.connect(perception_interface, "output_poses", pose_injector, "detections")

    app.run()

if __name__ == '__main__':
    parser = argparse.ArgumentParser()
    parser.add_argument(
        "--mode",
        help="Mode for data source",
        type=str,
        choices=["sim", "realsense", "cask"],
        default="sim")
    parser.add_argument("--cask", help="Name of the cask recording base directory", type=str)
    parser.add_argument("--refinement", help="In franka, use pose refinement", action='store_true')
    args = parser.parse_args()
    main(args)


text:
# License: MIT

import sys
import traceback
import copy
import numpy as np
from openbox.acquisition_function.acquisition import EI, LPEI
from openbox.surrogate.base.rf_with_instances import RandomForestWithInstances
from openbox.acq_maximizer.ei_optimization import InterleavedLocalAndRandomSearch, RandomSearch
from openbox.acq_maximizer.random_configuration_chooser import ChooserProb
from openbox.utils.util_funcs import get_types, get_rng
from openbox.utils.config_space.util import convert_configurations_to_array
from openbox.utils.constants import MAXINT, SUCCESS, FAILED, TIMEOUT
from openbox.utils.limit import time_limit, TimeoutException
from openbox.facade.bo_facade import BaseFacade

class BatchBayesianOptimization(BaseFacade):
    def __init__(self, objective_function, config_space,
                 sample_strategy='bo',
                 time_limit_per_trial=180,
                 max_runs=200,
                 logging_dir='logs',
                 initial_configurations=None,
                 initial_batch=1,
                 batch_size=3,
                 task_id='default_task_id',
                 rng=None):
        super().__init__(config_space, task_id, output_dir=logging_dir)
        self.logger = super()._get_logger(self.__class__.__name__)
        if rng is None:
            run_id, rng = get_rng()

        self.batch_size = batch_size
        self.init_batch = initial_batch
        self.max_iterations = max_runs
        self.iteration_id = 0
        self.sls_max_steps = None
        self.n_sls_iterations = 5
        self.sls_n_steps_plateau_walk = 10
        self.time_limit_per_trial = time_limit_per_trial
        self.default_obj_value = MAXINT
        self.sample_strategy = sample_strategy

        self.configurations = list()
        self.failed_configurations = list()
        self.perfs = list()

        # Initialize the basic component in BO.
        self.config_space.seed(rng.randint(MAXINT))
        self.objective_function = objective_function
        types, bounds = get_types(config_space)
        # TODO: what is the feature array.
        self.model = RandomForestWithInstances(types=types, bounds=bounds, seed=rng.randint(MAXINT))
        if self.sample_strategy == 'local_penalization':
            self.acquisition_function = LPEI(self.model)
        else:
            self.acquisition_function = EI(self.model)
        self.optimizer = InterleavedLocalAndRandomSearch(
            acquisition_function=self.acquisition_function,
            config_space=self.config_space,
            rng=np.random.RandomState(seed=rng.randint(MAXINT)),
            max_steps=self.sls_max_steps,
            n_steps_plateau_walk=self.sls_n_steps_plateau_walk,
            n_sls_iterations=self.n_sls_iterations
        )
        self._random_search = RandomSearch(
            self.acquisition_function, self.config_space, rng
        )
        self.random_configuration_chooser = ChooserProb(prob=0.25, rng=rng)

    def run(self):
        while self.iteration_id < self.max_iterations:
            self.iterate()

    def iterate(self):
        if len(self.configurations) == 0:
            X = np.array([])
        else:
            failed_configs = list() if self.max_y is None else self.failed_configurations.copy()
            X = convert_configurations_to_array(self.configurations + failed_configs)

        failed_perfs = list() if self.max_y is None else [self.max_y] * len(self.failed_configurations)
        Y = np.array(self.perfs + failed_perfs, dtype=np.float64)

        config_list = self.choose_next(X, Y)
        trial_state_list = list()
        trial_info_list = list()
        perf_list = list()

        for i, config in enumerate(config_list):
            trial_state_list.append(SUCCESS)
            trial_info_list.append(None)
            perf_list.append(None)

            if config not in (self.configurations + self.failed_configurations):
                # Evaluate this configuration.
                try:
                    args, kwargs = (config,), dict()
                    timeout_status, _result = time_limit(self.objective_function, self.time_limit_per_trial,
                                                         args=args, kwargs=kwargs)
                    if timeout_status:
                        raise TimeoutException(
                            'Timeout: time limit for this evaluation is %.1fs' % self.time_limit_per_trial)
                    else:
                        perf_list[i] = _result
                except Exception as e:
                    if isinstance(e, TimeoutException):
                        trial_state_list[i] = TIMEOUT
                    else:
                        traceback.print_exc(file=sys.stdout)
                        trial_state_list[i] = FAILED
                    perf_list[i] = MAXINT
                    trial_info_list[i] = str(e)
                    self.logger.error(trial_info_list[i])

                if trial_state_list[i] == SUCCESS and perf_list[i] < MAXINT:
                    if len(self.configurations) == 0:
                        self.default_obj_value = perf_list[i]

                    self.configurations.append(config)
                    self.perfs.append(perf_list[i])
                    self.history_container.add(config, perf_list[i])

                    self.perc = np.percentile(self.perfs, self.scale_perc)
                    self.min_y = np.min(self.perfs)
                    self.max_y = np.max(self.perfs)
                else:
                    self.failed_configurations.append(config)
            else:
                self.logger.debug('This configuration has been evaluated! Skip it.')
                if config in self.configurations:
                    config_idx = self.configurations.index(config)
                    trial_state_list[i], perf_list[i] = SUCCESS, self.perfs[config_idx]
                else:
                    trial_state_list[i], perf_list[i] = FAILED, MAXINT

        self.iteration_id += 1
        self.logger.info(
            'Iteration-%d, objective improvement: %.4f' % (
                self.iteration_id, max(0, self.default_obj_value - min(perf_list))))
        return config_list, trial_state_list, perf_list, trial_info_list

    def choose_next(self, X: np.ndarray, Y: np.ndarray):
        # Select a batch of configs to evaluate next.
        _config_num = X.shape[0]
        batch_configs_list = list()

        if _config_num < self.init_batch * self.batch_size or self.sample_strategy == 'random':
            for i in range(self.batch_size):
                batch_configs_list.append(self.sample_config())
            return batch_configs_list

        if self.sample_strategy == 'median_imputation':
            estimated_y = np.mean(Y)
            batch_history_container = copy.deepcopy(self.history_container)
            for i in range(self.batch_size):
                self.model.train(X, Y)

                incumbent_value = batch_history_container.get_incumbents()[0][1]

                self.acquisition_function.update(model=self.model, eta=incumbent_value,
                                                 num_data=len(batch_history_container.data))

                challengers = self.optimizer.maximize(
                    runhistory=batch_history_container,
                    num_points=5000,
                    random_configuration_chooser=self.random_configuration_chooser
                )

                is_repeated_config = True
                repeated_time = 0
                curr_batch_config = None
                while is_repeated_config:
                    try:
                        curr_batch_config = challengers.challengers[repeated_time]
                        batch_history_container.add(curr_batch_config, estimated_y)
                    except ValueError:
                        is_repeated_config = True
                        repeated_time += 1
                    else:
                        is_repeated_config = False

                batch_configs_list.append(curr_batch_config)
                X = np.append(X, curr_batch_config.get_array().reshape(1, -1), axis=0)
                Y = np.append(Y, estimated_y)
                estimated_y = np.mean(Y)

        elif self.sample_strategy == 'local_penalization':
            self.model.train(X, Y)
            incumbent_value = self.history_container.get_incumbents()[0][1]
            # L = self.estimate_L(X)
            for i in range(self.batch_size):
                self.acquisition_function.update(model=self.model, eta=incumbent_value,
                                                 num_data=len(self.history_container.data),
                                                 batch_configs=batch_configs_list)

                challengers = self.optimizer.maximize(
                    runhistory=self.history_container,
                    num_points=5000,
                    random_configuration_chooser=self.random_configuration_chooser
                )
                batch_configs_list.append(challengers.challengers[0])
        else:
            raise ValueError('Invalid sampling strategy - %s.' % self.sample_strategy)

        return batch_configs_list

    def sample_config(self):
        config = None
        _sample_cnt, _sample_limit = 0, 10000
        while True:
            _sample_cnt += 1
            config = self.config_space.sample_configuration()
            if config not in (self.configurations + self.failed_configurations):
                break
            if _sample_cnt >= _sample_limit:
                config = self.config_space.sample_configuration()
                break
        return config

    # def estimate_L(self, X):  # X: N*D
    #     """
    #         Estimate the Lipschitz constant of f by taking maximum norm of the gradient of f.
    #         Calculate numerical gradient = f(x+h) - f(x-h) / 2h , where h is very small
    #     """
    #     eps = 1e-3
    #     for i in range(10000):
    #         X = np.vstack((X, self.sample_config().get_array()))
    #
    #     # print("=" * 20, "X", X)
    #     dX = np.zeros_like(X)
    #     for j in range(X.shape[1]):
    #         h = np.zeros_like(X)
    #         h[:, j] = eps
    #         # print("=" * 20, "X+h", X + h)
    #         f_pos, _ = self.surrogate.predict_marginalized_over_instances(X + h)
    #         f_neg, _, = self.surrogate.predict_marginalized_over_instances(X - h)
    #         # print("=" * 20, "f_pos", dX)
    #         # print("=" * 20, "f_neg", dX)
    #         f_pos = f_pos.reshape(-1)
    #         f_neg = f_neg.reshape(-1)
    #         dX[:, j] = (f_pos - f_neg) / (2 * eps)
    #     # print("="*20,"dX",dX)
    #
    #     norm_dX = np.linalg.norm(dX, axis=1)  # shape=(N,)
    #
    #     return np.max(norm_dX)

max_stars_repo_name: Dee-Why/lite-bo

max_issues_repo_name: ZongWei-HUST/open-box

max_forks_repo_name: ZongWei-HUST/open-box


text:
import numpy as np

def loader_np():
    with open('../.input/day04', 'r') as f:
        data = f.read().splitlines()

    vector = np.array([int(x) for x in data[0].split(",")])
    tensor = np.array([
        np.array(list(map(np.int32, map(str.split, board))))
        for board
        in [data[i:i+5] for i in range(2, len(data), 6)]
    ])
    return vector, tensor

def solve1() -> int:
    vector, tensor = loader_np()
    for scale in range(1, len(vector)+1):
        for matrix in tensor:
            mask = np.isin(matrix, vector[:scale])
            if mask.all(axis=0).any() or mask.all(axis=1).any():
                return (matrix * ~mask).sum() * vector[scale-1]

def solve2() -> int:
    vector, tensor = loader_np()
    for scale in range(len(vector)+1, 0, -1):
        for matrix in tensor:
            mask = np.isin(matrix, vector[:scale])
            if not mask.all(axis=0).any() and not mask.all(axis=1).any():
                prev_mask = np.isin(matrix, vector[:scale+1])
                return (matrix * ~prev_mask).sum() * vector[scale]

if __name__ == '__main__':
    print(solve1())  # 67716
    print(solve2())  # 1830

max_stars_repo_name: kamoshi/Advent-of-Code

max_issues_repo_name: kamoshi/advent-of-code

max_forks_repo_name: kamoshi/advent-of-code


text:
# -*- coding: utf-8 -*-
"""320 HS_.ipynb

Automatically generated by Colaboratory.

"""

# Libraries
import pandas as pd
import re
from sklearn.feature_extraction.text import TfidfVectorizer
from sklearn.feature_extraction.text import ENGLISH_STOP_WORDS
from sklearn.linear_model import LogisticRegression
from sklearn.multiclass import OneVsRestClassifier
from sklearn.model_selection import train_test_split, GridSearchCV
from sklearn import linear_model
import numpy as np
from sklearn.metrics import classification_report, confusion_matrix
import seaborn as sn
import matplotlib.pyplot as plt

# Data loading
print(data)

#data["PRODUCT_DESCRIPTION"] = data["PRODUCT_DESCRIPTION"].str.replace('[\d{1,2}.\d{1,2}.\d{2,4}]',"")

data["PRODUCT_DESCRIPTION"] = data["PRODUCT_DESCRIPTION"].str.replace('\(|\)'," ")

data["PRODUCT_DESCRIPTION"] = data["PRODUCT_DESCRIPTION"].str.replace('[^a-zA-Z\d\s+]',"")

data["PRODUCT_DESCRIPTION"] = data["PRODUCT_DESCRIPTION"].apply(lambda x: re.sub('\d\d\d\d\d+', '',x))

data["PRODUCT_DESCRIPTION"]

data['Product_ID'] = data.Product_ID.astype(str).str.replace('\[|\]|\'', '')

data["Product_ID"]

data.head()

fe = TfidfVectorizer(
    min_df=1,
    max_df=1.0,
    ngram_range=(1, 4),
    max_features=5000,
    stop_words=ENGLISH_STOP_WORDS,
    decode_error='ignore'

)

x_data = fe.fit_transform(data["PRODUCT_DESCRIPTION"])

X_train, X_test, y_train, y_test = train_test_split(x_data, data['Product_ID'], test_size=0.2, random_state=1)

X_train.shape, X_test.shape,y_train.shape, y_test.shape

"""logistic = linear_model.LogisticRegression()
penalty = ['l1', 'l2']
C = np.logspace(0, 4, 10)
hyperparameters = dict(C=C, penalty=penalty)
clf1 = GridSearchCV(logistic, hyperparameters, cv=5, verbose=0)

best_model = clf1.fit(X_train, y_train)
best_model.best_estimator_.get_params()"""

model = OneVsRestClassifier(estimator=LogisticRegression(C=1, class_weight=None,
                                                 dual=False, fit_intercept=True,
                                                 intercept_scaling=1,
                                                 l1_ratio=None, max_iter=100,
                                                 multi_class='warn',
                                                 n_jobs=None, penalty='l1',
                                                 random_state=None,
                                                 solver='warn', tol=0.0001,
                                                 verbose=0, warm_start=False))

model.fit(X_train, y_train)

predicted = model.predict(X_test)

print(classification_report(y_test, predicted))

print(confusion_matrix(y_test, predicted))

test_data=pd.read_excel("sampe_test_.xlsx")

test_data["PRODUCT_DESCRIPTION"] = test_data["PRODUCT_DESCRIPTION"].str.replace('[^a-zA-Z\d\s+]',"")

test_data['Product_ID'] = test_data.Product_ID.astype(str).str.replace('\[|\]|\'', '')

test_data["Product_ID"]

x_test_data = fe.transform(test_data["PRODUCT_DESCRIPTION"])

predictedtest = model.predict(x_test_data)

test_data['predicted'] = predictedtest

max_stars_repo_name: AdiptaBiswas/company-name-mapping

max_issues_repo_name: AdiptaBiswas/company-name-mapping

max_forks_repo_name: AdiptaBiswas/company-name-mapping


text:
import time
from numba import jit

@jit
def loop(repeats=10000000):
    x = 0
    for i in range(repeats):
        x += 1
    return x

for i in range(3):
    start_time = time.time()
    loop()
    print(time.time() - start_time)

max_stars_repo_name: butla/experiments

max_issues_repo_name: butla/experiments

max_forks_repo_name: butla/experiments

reason: from numba


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

# 合并evo和vgg两个特征值

import numpy as np
from vggface import verify as verify1
from face_evoLVe import verify as verify2

# 定位人脸，然后人脸的特征值列表，可能不止一个脸
def get_features(filename):
    encoding_list1, face_boxes1 = verify1.get_features(filename)
    encoding_list2, face_boxes2 = verify2.get_features(filename)

    encoding_list = []
    for i in range(len(face_boxes1)):
        encoding_list.append(np.concatenate((encoding_list1[i], encoding_list2[i]), axis=0))

    return encoding_list, face_boxes1


text:
import matplotlib.pyplot as plt
import numpy as np
fig, ax = plt.subplots()
x = np.linspace(0, 10, 200)
y = np.sin(x)
ax.plot(x, y, 'r-', lw=2, label=r'$y=\sin(x)$', alpha=0.6)
ax.legend(loc='upper center')
plt.show()

max_stars_repo_name: sglyon/quant-econ


text:
from __future__ import division, print_function, absolute_import

from os.path import splitext

from ..utils.six import string_types

import numpy as np

def read_bvals_bvecs(fbvals, fbvecs):
    """
    Read b-values and b-vectors from the disk

    Parameters
    ----------
    fbvals : str
             path of file with b-values, or None if you don't want to read bvals
    fbvecs : str
             path of file with b-vectorsl, or None if you don't want to read bvecs

    Returns
    -------
    bvals : array, (N,) or None
    bvecs : array, (N, 3) or None

    Notes
    -----
    Files can be either '.bvals'/'.bvecs' or '.txt' or '.npy' (containing arrays
    stored with the appropriate values).

    """

    # Loop over the provided inputs, reading each one in turn and adding them
    # to this list:
    vals = []
    for this_fname in [fbvals, fbvecs]:
        # If the input was None, we don't read anything and move on:
        if this_fname is None:
            vals.append(None)
        else:
            if isinstance(this_fname, string_types):
                base, ext = splitext(this_fname)
                if ext in ['.bvals', '.bval', '.bvecs', '.bvec', '.txt', '']:
                    vals.append(np.squeeze(np.loadtxt(this_fname)))
                elif ext == '.npy':
                    vals.append(np.squeeze(np.load(this_fname)))
                else:
                    e_s = "File type %s is not recognized"%ext
                    raise ValueError(e_s)
            else:
                raise ValueError('String with full path to file is required')

    # Once out of the loop, unpack them:
    bvals, bvecs = vals[0], vals[1]

    # If bvecs is None, you can just return now w/o making more checks:
    if bvecs is None:
        return bvals, bvecs

    if bvecs.shape[1] > bvecs.shape[0]:
        bvecs = bvecs.T
    if min(bvecs.shape) != 3:
        raise IOError('bvec file should have three rows')
    if bvecs.ndim != 2:
        raise IOError('bvec file should be saved as a two dimensional array')

    # If bvals is None, you don't need to check that they have the same shape:
    if bvals is None:
        return bvals, bvecs

    if max(bvals.shape) != max(bvecs.shape):
            raise IOError('b-values and b-vectors shapes do not correspond')

    return bvals, bvecs

max_stars_repo_name: oesteban/dipy

max_issues_repo_name: jyeatman/dipy

max_forks_repo_name: jyeatman/dipy