text:
################################################################################
##
## This library is free software; you can redistribute it and/or
## modify it under the terms of the GNU Lesser General Public
## License as published by the Free Software Foundation; either
## version 2.1 of the License, or (at your option) any later version.
##
## This library is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
## Lesser General Public License for more details.
##
## You should have received a copy of the GNU Lesser General Public
## License along with this library; if not, write to the Free Software
## Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301 USA
##
## (C) Copyrights Dr. Michel F. Sanner and TSRI 2016
##
################################################################################

"""
#// $Id: AtomAndBondShader.py,v 1.6.4.1 2017/07/13 22:27:37 annao Exp $
#// ****************************************************************** #//
#//                                                                    #//
#// Copyright (C) 2010-2011 by                                         #//
#// Laboratoire de Biochimie Theorique (CNRS),                         #//
#// Laboratoire d'Informatique Fondamentale d'Orleans (Universite      #//
#// d'Orleans),                                                        #//
#// (INRIA) and                                                        #//
#// Departement des Sciences de la Simulation et de l'Information      #//
#// (CEA).                                                             #//
#// ALL RIGHTS RESERVED.                                               #//
#//                                                                    #//
#// contributors :                                                     #//
#// Matthieu Chavent,                                                  #//
#// Antoine Vanel,                                                     #//
#// Alex Tek,                                                          #//
#// Marc Piuzzi,                                                       #//
#// Jean-Denis Lesage,                                                 #//
#// Bruno Levy,                                                        #//
#// Sophie Robert,                                                     #//
#// Sebastien Limet,                                                   #//
#// Bruno Raffin and                                                   #//
#// Marc Baaden                                                        #//
#// Adapted to Python by Ludovic Autin                                 #//
#// May 2011                                                           #//
#//                                                                    #//
#// Contact: Marc Baaden                                               #//
#// E-mail: baaden@smplinux.de                                         #//
#// Webpage: http:#//hyperballs.sourceforge.net                         #//
#//                                                                    #//
#// This software is a computer program whose purpose is to visualize  #//
#// molecular structures. The source code is part of FlowVRNano, a     #//
#// general purpose library and toolbox for interactive simulations.   #//
#//                                                                    #//
#// This software is governed by the CeCILL-C license under French law #//
#// and abiding by the rules of distribution of free software. You can #//
#// use, modify and/or redistribute the software under the terms of    #//
#// the CeCILL-C license as circulated by CEA, CNRS and INRIA at the   #//
#// following URL "http:#//www.cecill.info".                            #//
#//                                                                    #//
#// As a counterpart to the access to the source code and  rights to   #//
#// copy, modify and redistribute granted by the license, users are    #//
#// provided only with a limited warranty and the software's author,   #//
#// the holder of the economic rights, and the successive licensors    #//
#// have only limited liability.                                       #//
#//                                                                    #//
#// In this respect, the user's attention is drawn to the risks        #//
#// associated with loading, using, modifying and/or developing or     #//
#// reproducing the software by the user in light of its specific      #//
#// status of free software, that may mean  that it is complicated to  #//
#// manipulate, and that also therefore means that it is reserved for  #//
#// developers and experienced professionals having in-depth computer  #//
#// knowledge. Users are therefore encouraged to load and test the     #//
#// software's suitability as regards their requirements in conditions #//
#// enabling the security of their systems and/or data to be ensured   #//
#// and, more generally, to use and operate it in the same conditions  #//
#// as regards security.                                               #//
#//                                                                    #//
#// The fact that you are presently reading this means that you have   #//
#// had knowledge of the CeCILL-C license and that you accept its      #//
#// terms.                                                             #//
#// ****************************************************************** #//
"""
import numpy as np
from math import *
useopengltk=False
try :
    import opengltk
    useopengltk=True
except:
    useopengltk=False
if useopengltk :
    from opengltk.OpenGL import GL
    from opengltk.OpenGL import GLU
    from opengltk.extent import _gllib
    from opengltk.extent import glextlib
    GL.GL_STREAM_DRAW = glextlib.GL_STREAM_DRAW_ARB
    GL.GL_ELEMENT_ARRAY_BUFFER = glextlib.GL_ELEMENT_ARRAY_BUFFER # 0x8893
    GL.GL_ARRAY_BUFFER = glextlib.GL_ARRAY_BUFFER_ARB #0x8892
    glBindBuffer = glextlib.glBindBufferARB#GL.glBindBuffer
    glBufferData = glextlib.glBufferDataARB
    GL.glGenBuffers = glextlib.glGenBuffersARB
    GL.glTexImage2D = _gllib.glTexImage2D
    for k in glextlib.__dict__:
        if not hasattr(GL,k):
            setattr(GL,k, glextlib.__dict__[k])
else :
    from OpenGL import GL
    from OpenGL import GLU

GL_TEXTURE_RECTANGLE_NV = 0x84F5
GL_RGB32F_ARB = 0x8815
GL_RGBA32F_ARB = 0x8814#from OpenGL import GL as oGL
import sys

#from OpenGL.arrays import vbo

sizeof = sys.getsizeof

#in case
#35040
#// Identifier of the different buffers  used by the shaders
BUF_VERTICE=0
BUF_INDICE=1
BUF_TCOORD=2
BUF_TCOORD0=3
BUF_TCOORD1=4
BUF_TCOORD2=5
BUF_STICKV=2
BUF_STICKI=3

#// idenfifier of the different textures used by the shaders
TEX_POSITIONS=0
TEX_COLORS=1
TEX_SIZES=2
TEX_ASCALES=3
TEX_BSCALES=4
TEX_SHRINKS=5
#    	#// vertex, indices, and 4 for coordinates
NB_BUFFERS=4#7+2;
#    	#// atom position, atom color, atom size, atom scale, bond scale and bond shrink,framebuffer
NB_TEXTURES=7;

CUBE_VERTICES =  np.array([
         [-1,  -1, -1],
         [1,  -1, -1],
         [1,  -1, 1],
         [-1,  -1, 1],
         [-1,  1, -1],
         [1,  1, -1],
         [1,  1, 1],
         [-1,  1, 1],
    ],'f')

#6 quad face
CUBE_INDICES = np.array([
         [0,  1, 2],
         [0,  2, 3],
         [1,  5, 6],
         [1,  6, 2],
         [4,  6, 5],
         [4,  7, 6],
         [0,  7, 4],
         [0,  3, 7],
         [0,  5, 1],
         [0,  4, 5],
         [3,  2, 6],
         [3,  6, 7],    ],np.uint32)

import array

v=np.array([[ 1.,  1.,  0.],
       [ 0.,  1.,  0.],
       [ 0.,  0.,  0.],
       [ 1.,  0.,  0.],
       [ 1.,  1.,  1.],
       [ 0.,  1.,  1.],
       [ 0.,  0.,  1.],
       [ 1.,  0.,  1.]])

verts = np.array( [
-0.285437,-0.744976,-0.471429,
-0.285437,-0.744976,-2.471429,
1.714563,-0.744976,-2.471429,
1.714563,-0.744976,-0.471429,
-0.285437,1.255024,-0.471429,
-0.285437,1.255024,-2.471429,
1.714563,1.255024,-2.471429,
1.714563,1.255024,-0.471429] )

faces = np.array( [
4,5,1,
5,6,2,
6,7,3,
4,0,7,
0,1,2,
7,6,5,
0,4,1,
1,5,2,
2,6,3,
7,0,3,
3,0,2,
4,7,5] )

def hexahedron(radius):
        """
        Create the mesh data of a hexahedron of a given radius

        @type  radius: float
        @param radius: radius of the embeding sphere

        @rtype:   array
        @return:  vertex,face, face normal of the hexahedron
        """
        faceIndices = ( (3,2,0), (1,2,3), #bot
                        (4,6,7), (7,6,5), #top
                        (0,2,6), (6,4,0), #front
                        (2,1,5), (5,6,2), #right
                        (1,3,7), (7,5,1), #rear
                        (3,0,4), (4,7,3), #left
                        )
        a = 1.0
        faceNormals = None
        diameter = radius * 2
        width  = diameter
        height = diameter
        depth  = radius*sin(radians(45.))
        boundingRadius = radius
        surfaceArea = (sqrt(3) / 2) * (radius ** 2)
        volume = (4/3) * (radius ** 3)
        _corners = (
            (-radius, 0.0, -depth), (radius, 0.0, -depth),
            (0.0, -radius, -depth), (0.0, radius, -depth),
            (-radius, 0.0, depth), (radius, 0.0, depth),
            (0.0, -radius, depth), (0.0, radius, depth),)

        return np.array(_corners,'f'),np.array(faceIndices,np.uint32),np.array(faceNormals,float)

#verts,faces,fn= hexahedron(1.0)
verts=CUBE_VERTICES
faces=CUBE_INDICES
#verts = np.array([[-1,0,0],[1,0,0],[1,1,0],[-1,1,0]],'f')*5.0#+np.array([9.2688, 9.7873, 7.9671],'f')
#faces = np.array([ [0,1,2],[0,2,3]],np.uint32)
#from OpenGLContext.scenegraph.basenodes import Sphere
#v_vbo,i_vbo,count = Sphere(  radius = 1  ).compile()

##//Return power of 2 superior than nb
##//Used for setting size of textures
def nextPower2(nb):
    power = 2;
    while(power <= nb):
        power = power << 1;
    return power;

class AtomAndBondShaders:
    def __init__(self,) :
#    	#// GL variables
        self.buffer_id =np.ones(NB_BUFFERS,int);
        self.texture_id= np.ones(NB_TEXTURES,int);#GLuint
        self.vertice=None;
        self.vertice_sticks=None
        self.texturecoord=None;
        self.texturecoord0=None;
        self.texturecoord1=None;
        self.texturecoord2=None;
        self.vertice_offset=None
        self.vertice_offset1=None
        self.vertice_offset2=None
        self.indice=None;
        self.indice_sticks=None
        self.colors=None;
        self.positions=None;
        self.atomScales=None;
        self.bondScales=None;
        self.sizes=None;
        self.shrinks=None;

        self.aScales=0.05;
        self.bScales=0.05;
        self.bshrinks=0.1;

        self.nbAtoms=nbBonds=0;
        self.textureSize=0
        self.textureSizeBonds=0;
        #// Parameter passed to the shaders??
        self.normMax=0.0;
        self.useVBO=True
        self.links = None

        self.atomVertexShaderProgramName=""
        self.bondVertexShaderProgramName=""
        self.atomFragmentShaderProgramName=""
        self.bondFragmentShaderProgramName=""

        self.atomVertexShaderProgramSource=""
        self.bondVertexShaderProgramSource=""
        self.atomFragmentShaderProgramSource=""
        self.bondFragmentShaderProgramSource=""

        self.ivert = np.reshape(np.mgrid[-1:2:2,-1:2:2,-1:2:2].T, (8,3))#one cube coordinate
        #one cube faces indices
        self.ifaces= np.array([[0,1,2],[0,2,3],
        [1,5,6],[1,6,2],
        [4,6,5],[4,7,6],
        [0,7,4],[0,3,7],
        [0,5,1],[0,4,5],
        [3,2,6],[3,6,7],
        ],int )

##//Destructor
#AtomAndBondShaders::~AtomAndBondShaders() {
#	if(buffer_id) 		delete [] buffer_id;
#	if(texture_id) 		delete [] texture_id;
#	if(vertice) 		delete [] vertice;
#	if(texturecoord) 	delete [] texturecoord;
#	if(texturecoord0)	delete [] texturecoord0;
#	if(texturecoord1) 	delete [] texturecoord1;
#	if(texturecoord2) 	delete [] texturecoord2;
#	if(indice)			delete [] indice;
#	if(colors) 			delete [] colors;
#	if(positions) 		delete [] positions;
#	if(atomScales) 		delete [] atomScales;
#	if(bondScales) 		delete [] bondScales;
#	if(sizes) 			delete [] sizes;
#	if(shrinks) 		delete [] shrinks;
#}

#pragma mark -
#pragma mark Setters
##//-----------------------------------------------------------------------------
    def setColors(self,colors):#GLfloat
        self.colors = colors#memcpy((GLfloat*)self.colors,(GLfloat*) colors,4*nbAtoms*sizeof(GLfloat));

##//-----------------------------------------------------------------------------
    def setPositions(self,positions):#GLfloat
        self.positions = positions;

##//-----------------------------------------------------------------------------
    def setAtomScales(self,scales):
        self.atomScales = scales #   memcpy((GLfloat*)self.atomScales,(GLfloat*) scales,nbAtoms*sizeof(GLfloat));

##//-----------------------------------------------------------------------------
    def setBondScales(self,scales):
        self.bondScales=scales#memcpy((GLfloat*)self.bondScales,(GLfloat*) scales,nbBonds*sizeof(GLfloat));

##//-----------------------------------------------------------------------------
    def setSizes(self,sizes):
        self.sizes = sizes#memcpy((GLfloat*)self.sizes,(GLfloat*) sizes,nbAtoms*sizeof(GLfloat));

##//-----------------------------------------------------------------------------
    def setShrinks(self, shrinks):
        self.shrinks=shrinks#memcpy((GLfloat*)self.shrinks,(GLfloat*) shrinks,nbAtoms*sizeof(GLfloat));

##//-----------------------------------------------------------------------------
    def setAtomFragmentShaderProgramName(self,fragmentShaderProgramName):
        self.atomFragmentShaderProgramName = fragmentShaderProgramName

##//-----------------------------------------------------------------------------
    def setBondFragmentShaderProgramName(self, fragmentShaderProgramName):
        self.bondFragmentShaderProgramName = fragmentShaderProgramName;

##//-----------------------------------------------------------------------------
    def setAtomVertexShaderProgramName(self, vertexShaderProgramName):
        self.atomVertexShaderProgramName = vertexShaderProgramName;

##//-----------------------------------------------------------------------------
    def setBondVertexShaderProgramName(self,vertexShaderProgramName):
        self.bondVertexShaderProgramName = vertexShaderProgramName;

    #//-----------------------------------------------------------------------------
    def setAtomFragmentShaderProgramSource(self,src):
    	self.atomFragmentShaderProgramSource = str(src);

    #//-----------------------------------------------------------------------------
    def setAtomVertexShaderProgramSource(self,src):
    	self.atomVertexShaderProgramSource = str(src);

    #//-----------------------------------------------------------------------------
    def setBondFragmentShaderProgramSource(self,src):
    	self.bondFragmentShaderProgramSource = str(src);

    #//-----------------------------------------------------------------------------
    def setBondVertexShaderProgramSource(self,src):
    	self.bondVertexShaderProgramSource = str(src);

    #pragma mark -
    #pragma mark Initialization Methods

    #//-----------------------------------------------------------------------------
    #//
    #// Init atoms parameters
    #//
    def initAtomTextures(self,colors,radius):
        #need to check shape of colors thatshould be (nbAtoms, 4)
        self.colors = np.array(colors,'f')#.flatten()
        if self.colors.shape[1] != 4 :
            one = np.ones( (self.colors.shape[0], 1), self.colors.dtype.char )
            self.colors = np.concatenate( (self.colors, one), 1 )
        self.sizes=np.array(radius,'f')
#        print self.sizes
        self.atomScales=np.ones(self.nbAtoms,'f')*self.aScales
#        print self.atomScales
#        print self.aScales

    def initAtomTextures_mol(self,dataAtoms):
        self.colors = np.array(map(lambda x: [x.colors['cpk'][0],x.colors['cpk'][1],x.colors['cpk'][2],1.0], dataAtoms),float)#.flatten()
        self.sizes=np.array(dataAtoms.vdwRadius,'f')
        self.atomScales=np.ones(self.nbAtoms,'f')*0.05
    #//-----------------------------------------------------------------------------
    #//
    #// Init bond parameters
    #//
    def initBondTextures(self,):
        self.shrinks=np.ones(self.nbBonds,'f')*self.bshrinks
        self.bondScales=np.ones(self.nbBonds,'f')*self.bScales

    def initBondTextures_mol(self, dataBonds):
        self.shrinks=np.ones(self.nbBonds,'f')*0.1
        self.bondScales=np.ones(self.nbBonds,'f')*0.05
#        for i in range(self.nbBonds):
##    	for(int i = 0; i < nbBonds; ++i){
#    		self.shrinks[i] = dataBonds[i].shrink;
#    		self.bondScales[i] = dataBonds[i].scale;

    #//-----------------------------------------------------------------------------
    #//
    #// Init vertex positions and storage in textures
    #//
    def initVertTCoordIndice(self,):
        #cube vertices and faces
        self.ivert,self.ifaces,fn= hexahedron(1.0)
        #Ball data
        nvert = self.nbAtoms# if  self.nbAtoms > self.nbBonds else self.nbBonds
        if nvert > 0:
            incr = np.repeat(np.arange(nvert)*len(self.ivert),len(self.ifaces))
            idata = np.tile(self.ifaces,(nvert,1))
            self.indice=idata+incr.reshape((len(incr),1))
            self.vertice = np.array(np.tile(self.ivert,(nvert,1)),'f')#added dirrectly the pshere position?
            self.vertice_offset = np.array(np.repeat(self.positions,len(self.ivert),axis=0),'f')
            self.texturecoord = np.vstack([np.repeat(np.arange(self.nbAtoms,dtype='f') % self.textureSize,len(self.ivert)),
                                           np.repeat(np.arange(self.nbAtoms,dtype='f') / self.textureSize,len(self.ivert))]).transpose().astype('f')#/self.textureSize
            self.texturecoord = np.floor(self.texturecoord)#/(self.textureSize-1)
            self.indice=np.array(self.indice,np.uint32)
        #add color, cant get it work with Texture !
#        self.vcolors = np.array(np.repeat(self.colors,len(self.ivert),axis=0),'f')
        #sticks data
        nvert = self.nbBonds# if  self.nbAtoms > self.nbBonds else self.nbBonds
        if nvert > 0:
            incr = np.repeat(np.arange(nvert)*len(self.ivert),len(self.ifaces))
            idata = np.tile(self.ifaces,(nvert,1))
            self.indice_sticks=idata+incr.reshape((len(incr),1))
            self.vertice_sticks = np.array(np.tile(self.ivert,(nvert,1)),'f')#added dirrectly the pshere position?
    #        self.vertice_sticks_offset = np.array(np.repeat(self.positions,len(self.ivert),axis=0),'f')
            vtext0 = np.array(self.links,int).transpose()
            #to fix
            self.texturecoord0=np.vstack([np.repeat(vtext0[0] % self.textureSize,len(self.ivert)),
                                          np.repeat(vtext0[0] / self.textureSize,len(self.ivert))]).transpose().astype('f')#/self.textureSize
            self.texturecoord1=np.vstack([np.repeat(vtext0[1] % self.textureSize,len(self.ivert)),
                                          np.repeat(vtext0[1] / self.textureSize,len(self.ivert))]).transpose().astype('f')#/self.textureSize
            self.texturecoord2 = np.vstack([np.repeat(np.arange(self.nbBonds,dtype='f') % self.textureSizeBonds,len(self.ivert)),
                                            np.repeat(np.arange(self.nbBonds,dtype='f') / self.textureSizeBonds,len(self.ivert))]).transpose().astype('f')#/self.textureSizeBonds
            self.indice_sticks=np.array(self.indice_sticks,np.uint32)
            self.vertice_offset1 = np.array(np.repeat(np.take(self.positions,vtext0[0].astype(int),axis=0),len(self.ivert),axis=0),'f')
            self.vertice_offset2 = np.array(np.repeat(np.take(self.positions,vtext0[1].astype(int),axis=0),len(self.ivert),axis=0),'f')
            self.texturecoord0 = np.floor(self.texturecoord0)#/(self.textureSize-1)
            self.texturecoord1 = np.floor(self.texturecoord1)#/(self.textureSize-1)
            self.texturecoord2 = np.floor(self.texturecoord2)#/(self.textureSize-1)

#        self.vcolors1 = np.array(np.repeat(np.take(self.colors,vtext0[0].astype(int),axis=0),len(self.ivert),axis=0),'f')
#        self.vcolors2 = np.array(np.repeat(np.take(self.colors,vtext0[1].astype(int),axis=0),len(self.ivert),axis=0),'f')
    #//-----------------------------------------------------------------------------
    #//
    #// Init arrays used to fill the textures
    #//
    def setupBuffersAndTextures(self, nbAtoms, nbBonds, coords, bonds, colors, radii):
    	#//initialize nbAtoms
    	self.nbAtoms = nbAtoms;
    	self.nbBonds = nbBonds;
    	self.positions = np.array(coords,'f')
    	#print self.positions
    	nbMax=self.nbAtoms if self.nbAtoms>self.nbBonds else self.nbBonds

    	#// textureSize contains the side size of the  textures for atoms
    	self.textureSize=int(ceil(sqrt(self.nbAtoms)));
    	#// textureSizeBonds contains the side size of the  textures for bonds
    	self.textureSizeBonds=int(ceil(sqrt(self.nbBonds)));

    	self.links=np.array(bonds,np.uint32)#new LinkAtom[nbBonds];
    	#for(int z=0; z<nbBonds; z++){
    	#	links[z].atom1 = dataBonds[z].idsrc;
    	#	links[z].atom2 = dataBonds[z].iddest;
    	#}
    	print "setupBuffersAndTextures"
    	self.initAtomTextures(colors,radii);
    	self.initBondTextures();
    	self.initVertTCoordIndice();
    	print "OK setupBuffersAndTextures"

    def setupBuffersAndTextures_mol(self, nbAtoms, nbBonds, dataAtoms, dataBonds):
    	#//initialize nbAtoms
    	self.nbAtoms = nbAtoms;
    	self.nbBonds = nbBonds;
    	self.positions = dataAtoms.coords
    	print self.positions
    	nbMax=self.nbAtoms if self.nbAtoms>self.nbBonds else self.nbBonds

    	#// textureSize contains the side size of the  textures for atoms
    	self.textureSize=int(ceil(sqrt(self.nbAtoms)));
    	#// textureSizeBonds contains the side size of the  textures for bonds
    	self.textureSizeBonds=int(ceil(sqrt(self.nbBonds)));

    	self.links=dataBonds#new LinkAtom[nbBonds];
    	#for(int z=0; z<nbBonds; z++){
    	#	links[z].atom1 = dataBonds[z].idsrc;
    	#	links[z].atom2 = dataBonds[z].iddest;
    	#}
    	self.initAtomTextures_mol(dataAtoms);
    	self.initBondTextures_mol(dataBonds);
    	self.initVertTCoordIndice();

    #pragma mark -
    #pragma mark Update Parameters Methods

    def updateBuffersAndTextures(self, nbAtoms, nbBonds, coords, bonds, colors, radii):
        #//initialize nbAtoms
        self.nbAtoms = nbAtoms;
        self.nbBonds = nbBonds;
        self.positions = np.array(coords,'f')
        #print self.positions
        nbMax=self.nbAtoms if self.nbAtoms>self.nbBonds else self.nbBonds
        #// textureSize contains the side size of the  textures for atoms
        self.textureSize=int(ceil(sqrt(self.nbAtoms)));
        #// textureSizeBonds contains the side size of the  textures for bonds
        self.textureSizeBonds=int(ceil(sqrt(self.nbBonds)));
        self.links=np.array(bonds,np.uint32)#new LinkAtom[nbBonds];
        self.initAtomTextures(colors,radii);
        self.initBondTextures();
        self.initVertTCoordIndice();
        self.activateAllTextures(update=True)
        self.updateVBO()

    def updateTextureSize(self):
        nbMax=self.nbAtoms if self.nbAtoms>self.nbBonds else self.nbBonds
        #// textureSize contains the side size of the  textures for atoms
        self.textureSize=int(ceil(sqrt(self.nbAtoms)));
        #// textureSizeBonds contains the side size of the  textures for bonds
        self.textureSizeBonds=int(ceil(sqrt(self.nbBonds)));

    #//-----------------------------------------------------------------------------
    def updatePositions(self, positions):
        #update textureSize
        self.activeTexture(GL.GL_TEXTURE0,self.texture_id[TEX_POSITIONS],np.array(positions,'f'),self.textureSize,size=3)
        #radius ?

    def updateSizes(self, sizes):
        self.sizes = np.array(sizes,'f')
        self.activeTexture(GL.GL_TEXTURE2,self.texture_id[TEX_SIZES],self.sizes,self.textureSize)

    #//-----------------------------------------------------------------------------
    def updateColors(self, colors):
        self.colors = np.array(colors,'f')#.flatten()
        if self.colors.shape[1] != 4 :
            one = np.ones( (self.colors.shape[0], 1), self.colors.dtype.char )
            self.colors = np.concatenate( (self.colors, one), 1 )
        self.activeTexture(GL.GL_TEXTURE1,self.texture_id[TEX_COLORS],self.colors,self.textureSize,size=4)
        print  self.texture_id,self.texture_id[TEX_COLORS]

    #//-----------------------------------------------------------------------------
    def updateAtomScales(self, newScale=1,scales=None):
        if scales is None :
            self.atomScales=np.ones(self.nbAtoms,'f')*newScale
        else :
            self.atomScales=np.array(scales[:],'f')
#    	self.activeTexture(GL.GL_TEXTURE2,self.texture_id[TEX_SIZES],self.sizes,self.textureSize)
    	self.activeTexture(GL.GL_TEXTURE3,self.texture_id[TEX_ASCALES],self.atomScales,self.textureSize)

    #//-----------------------------------------------------------------------------
    def updateBondScales(self, newScale=1,scales=None):
        if scales is None :
            self.bondScales=np.ones(self.nbBonds,'f')*newScale
        else :
            self.bondScales=np.array(scales[:],'f')
    	self.activeTexture(GL.GL_TEXTURE4,self.texture_id[TEX_BSCALES],self.bondScales,self.textureSizeBonds)

    #//-----------------------------------------------------------------------------
    def updateShrinks(self, newShrink, shrinks=None):
        if shrinks is None :
            self.shrinks=np.ones(self.nbBonds,'f')*newShrink
        else :
            self.shrinks=np.array(shrinks[:],'f')
    	self.activeTexture(GL.GL_TEXTURE5,self.texture_id[TEX_SHRINKS],self.shrinks,self.textureSizeBonds)

    #pragma mark -
    #pragma mark Opengl Drawing Methods

    #//-----------------------------------------------------------------------------
    def drawAtoms(self):
        pass

    #//-----------------------------------------------------------------------------
    def drawBonds(self):
        pass

    def generateVBOsimple(self):
        #Create the Vertex Array Object
#        v_vbo,i_vbo,count
        self.vertice_vbo=v_vbo#vbo.VBO( verts )#self.vertice
        self.indice_vbo=i_vbo#vbo.VBO( faces, target=GL.GL_ELEMENT_ARRAY_BUFFER )#self.indice
#        self.texturecoord_vbo=vbo.VBO( self.texturecoord,usage = GL.GL_STREAM_DRAW)

    def makeVBO(self,np_data,buff_indice,target=GL.GL_ARRAY_BUFFER,usage=GL.GL_STREAM_DRAW):
        buffer_id=int(self.buffer_id[buff_indice])#GL.glGenBuffers()
        glBindBuffer(target,int(buffer_id));
        glBufferData(target,np_data.nbytes,np_data.flatten(),usage)#,GL.GL_STREAM_DRAW);
        return buffer_id
        #4*len(vertices)*len(vertices[0])

    def bindVBO(self,buffer_id,np_data,target=GL.GL_ARRAY_BUFFER,usage=GL.GL_STREAM_DRAW):
        glBindBuffer(target,int(buffer_id))
        glBufferData(target,np_data.nbytes,np_data.flatten(),usage)#,GL.GL_STREAM_DRAW);

    def generateVBO(self):
        #Create the Vertex Array Object thatcarry everthing coords, uv, normal if any
        #ball_vbo
        self.ball_data = np.hstack([self.vertice,self.texturecoord])#,self.vertice_offset])#,self.vcolors])
        self.vertice_vbo=self.makeVBO(self.ball_data,BUF_VERTICE) #vbo.VBO( self.ball_data, usage = GL.GL_STREAM_DRAW)#self.vertice
        self.indice_vbo=self.makeVBO( self.indice,BUF_INDICE,usage = GL.GL_STREAM_DRAW, target=GL.GL_ELEMENT_ARRAY_BUFFER )#self.indice
        #stick_vbo
        self.stick_data = np.hstack([self.vertice_sticks,self.texturecoord0,self.texturecoord1,self.texturecoord2])#,
#                                     self.vertice_offset1,self.vertice_offset2,])
#                                     self.vcolors1,self.vcolors2])
        self.stick_vertice_vbo=self.makeVBO( self.stick_data,BUF_STICKV)# usage = GL.GL_STREAM_DRAW)#self.vertice
        self.stick_indice_vbo=self.makeVBO( self.indice_sticks,BUF_STICKI)#usage = GL.GL_STREAM_DRAW, target=GL.GL_ELEMENT_ARRAY_BUFFER )#self.indice

    def updateVBO(self):
        #Create the Vertex Array Object thatcarry everthing coords, uv, normal if any
        #ball_vbo
        self.ball_data = np.hstack([self.vertice,self.texturecoord])#,self.vertice_offset])#,self.vcolors])
#        self.vertice_vbo.set_array(self.ball_data)#  =vbo.VBO( self.ball_data, usage = GL.GL_STREAM_DRAW)#self.vertice
#        self.indice_vbo.set_array(self.indice)#=vbo.VBO( self.indice,usage = GL.GL_STREAM_DRAW, target=GL.GL_ELEMENT_ARRAY_BUFFER )#self.indice
        #stick_vbo
        self.stick_data = np.hstack([self.vertice_sticks,self.texturecoord0,self.texturecoord1,self.texturecoord2])#,
#                                     self.vertice_offset1,self.vertice_offset2,])
#                                     self.vcolors1,self.vcolors2])
#        self.stick_vertice_vbo.set_array(self.stick_data)#=vbo.VBO( self.stick_data, usage = GL.GL_STREAM_DRAW)#self.vertice
#        self.stick_indice_vbo.set_array(self.indice_sticks)#=vbo.VBO( self.indice_sticks,usage = GL.GL_STREAM_DRAW, target=GL.GL_ELEMENT_ARRAY_BUFFER )#self.indice

    def generateVBO_VBO(self):
        #Create the Vertex Array Object thatcarry everthing coords, uv, normal if any
        #ball_vbo
        self.ball_data = np.hstack([self.vertice,self.texturecoord])#,self.vertice_offset])#,self.vcolors])
        self.vertice_vbo=vbo.VBO( self.ball_data, usage = GL.GL_STREAM_DRAW)#self.vertice
        self.indice_vbo=vbo.VBO( self.indice,usage = GL.GL_STREAM_DRAW, target=GL.GL_ELEMENT_ARRAY_BUFFER )#self.indice
        #stick_vbo
        self.stick_data = np.hstack([self.vertice_sticks,self.texturecoord0,self.texturecoord1,self.texturecoord2])#,
#                                     self.vertice_offset1,self.vertice_offset2,])
#                                     self.vcolors1,self.vcolors2])
        self.stick_vertice_vbo=vbo.VBO( self.stick_data, usage = GL.GL_STREAM_DRAW)#self.vertice
        self.stick_indice_vbo=vbo.VBO( self.indice_sticks,usage = GL.GL_STREAM_DRAW, target=GL.GL_ELEMENT_ARRAY_BUFFER )#self.indice

    def updateVBO_VBO(self):
        #Create the Vertex Array Object thatcarry everthing coords, uv, normal if any
        #ball_vbo
        self.ball_data = np.hstack([self.vertice,self.texturecoord,self.vertice_offset])#,self.vcolors])
        self.vertice_vbo.set_array(self.ball_data)#  =vbo.VBO( self.ball_data, usage = GL.GL_STREAM_DRAW)#self.vertice
        self.indice_vbo.set_array(self.indice)#=vbo.VBO( self.indice,usage = GL.GL_STREAM_DRAW, target=GL.GL_ELEMENT_ARRAY_BUFFER )#self.indice
        #stick_vbo
        self.stick_data = np.hstack([self.vertice_sticks,self.texturecoord0,self.texturecoord1,self.texturecoord2,
                                     self.vertice_offset1,self.vertice_offset2,])
#                                     self.vcolors1,self.vcolors2])
        self.stick_vertice_vbo.set_array(self.stick_data)#=vbo.VBO( self.stick_data, usage = GL.GL_STREAM_DRAW)#self.vertice
        self.stick_indice_vbo.set_array(self.indice_sticks)#=vbo.VBO( self.indice_sticks,usage = GL.GL_STREAM_DRAW, target=GL.GL_ELEMENT_ARRAY_BUFFER )#self.indice

    def useTexture(self,gl_id,tex_id):
        GL.glActiveTexture(gl_id);
#        GL.glBindTexture(GL_TEXTURE_RECTANGLE_NV,int(tex_id));#GL.GL_TEXTURE_2D
        GL.glBindTexture(GL.GL_TEXTURE_2D,int(tex_id));#GL.GL_TEXTURE_2D

    def activeTexture(self,gl_id,tex_id,data,texSize,size=1,update=False):
    	GL.glActiveTexture(gl_id);
#    	ttype = GL_TEXTURE_RECTANGLE_NV
    	ttype = GL.GL_TEXTURE_2D#or GL_TEXTURE_RECTANGLE_NV
    	print "texSize ",texSize,ttype,int(tex_id),data.shape#,data
    	#flatten the data ?
    	if size == 1 :
            GL.glEnable(ttype);
            GL.glBindTexture(ttype,int(tex_id))#int(self.texture_id[tex_id]));
            GL.glTexImage2D (ttype, 0, GL.GL_INTENSITY ,
                                  texSize, texSize, 0, GL.GL_LUMINANCE, GL.GL_FLOAT,
                      data.flatten());
    	elif size == 4 :
#            ttype = GL_TEXTURE_RECTANGLE_NV#or GL_TEXTURE_RECTANGLE_NV
            GL.glEnable(ttype);
            GL.glBindTexture(ttype,int(tex_id))#int(self.texture_id[tex_id]));
            GL.glTexImage2D (ttype, 0, GL_RGBA32F_ARB ,
                      texSize, texSize, 0, GL.GL_RGBA, GL.GL_FLOAT,
                      data.flatten());
    	else :#3
            try :
                GL.glEnable(ttype);
                GL.glBindTexture(ttype,int(tex_id))#int(self.texture_id[tex_id]));
                #print "try glTexImage2D",GL_RGB32F_ARB,GL.GL_RGB,data.flatten()
                GL.glTexImage2D (ttype, 0, GL_RGB32F_ARB ,
                      texSize, texSize, 0, GL.GL_RGB, GL.GL_FLOAT,
                      data.flatten());
            except :
                print "something rong"
    	print "ok teximage"
    	if not update :
         GL.glTexParameteri(ttype,GL.GL_TEXTURE_MIN_FILTER,GL.GL_NEAREST);
         GL.glTexParameteri(ttype,GL.GL_TEXTURE_MAG_FILTER,GL.GL_NEAREST);
    	GL.glDisable(ttype);

    def useAllTextures(self):
    	#// atom positions
    	self.useTexture(GL.GL_TEXTURE0,self.texture_id[TEX_POSITIONS])#,self.positions,size=3)
    	#// atom colors
    	self.useTexture(GL.GL_TEXTURE1,self.texture_id[TEX_COLORS])#,self.colors,size=4)
    	#// atom sizes
    	self.useTexture(GL.GL_TEXTURE2,self.texture_id[TEX_SIZES])#,self.sizes)
    	#// atom scales
    	self.useTexture(GL.GL_TEXTURE3,self.texture_id[TEX_ASCALES])#,self.atomScales)
    	#//   bond scales
    	self.useTexture(GL.GL_TEXTURE4,self.texture_id[TEX_BSCALES])#,self.bondScales)
    	#//   bond shrinks
    	self.useTexture(GL.GL_TEXTURE5,self.texture_id[TEX_SHRINKS])#,self.shrinks)

    def activateAllTextures(self,update=False):
    	if not update :
            self.texture_id=GL.glGenTextures(NB_TEXTURES).tolist()
    	print "activateAllTextures",self.texture_id
    	print  self.texture_id,self.texture_id[TEX_COLORS]
    	print "#// atom positions"
    	self.activeTexture(GL.GL_TEXTURE0,self.texture_id[TEX_POSITIONS],self.positions,self.textureSize,size=3)
    	print "#// atom colors"
    	self.activeTexture(GL.GL_TEXTURE1,self.texture_id[TEX_COLORS],self.colors,self.textureSize,size=4)
    	print "#// atom sizes"
    	self.activeTexture(GL.GL_TEXTURE2,self.texture_id[TEX_SIZES],self.sizes,self.textureSize)
    	print "#// atom scales"
    	self.activeTexture(GL.GL_TEXTURE3,self.texture_id[TEX_ASCALES],self.atomScales,self.textureSize)
    	print "#//   bond scales"
    	self.activeTexture(GL.GL_TEXTURE4,self.texture_id[TEX_BSCALES],self.bondScales,self.textureSizeBonds)
    	print "#//   bond shrinks"
    	self.activeTexture(GL.GL_TEXTURE5,self.texture_id[TEX_SHRINKS],self.shrinks,self.textureSizeBonds)

    def bindAllBuffer(self):
    	nbMax=self.nbAtoms if self.nbAtoms>self.nbBonds else self.nbBonds;
    	print "nmax is ",  nbMax

#    	print BUF_VERTICE,int(self.buffer_id[BUF_VERTICE]),nbMax*3*8*sizeof(float)
#    	glBindBuffer(GL.GL_ARRAY_BUFFER,int(self.buffer_id[BUF_VERTICE]));
#    	glBufferData(GL.GL_ARRAY_BUFFER,self.vertice.nbytes,self.vertice,GL.GL_STREAM_DRAW)#,GL.GL_STREAM_DRAW);
#    	print "bin buffer   BUF_VERTICE",self.vertice.flatten(),self.buffer_id[BUF_INDICE],BUF_INDICE,self.indice.flatten()
#
#    	glBindBuffer(GL.GL_ELEMENT_ARRAY_BUFFER,int(self.buffer_id[BUF_INDICE]));
#    	glBufferData(GL.GL_ELEMENT_ARRAY_BUFFER,self.indice.nbytes,self.indice,GL.GL_STREAM_DRAW);
#    	print "bin buffer   BUF_INDICE"

    	glBindBuffer(GL.GL_ARRAY_BUFFER,int(self.buffer_id[BUF_TCOORD]));
    	glBufferData(GL.GL_ARRAY_BUFFER,self.texturecoord.nbytes,self.texturecoord,GL.GL_STREAM_DRAW);
    	print "bin buffer   BUF_TCOORD"

    	glBindBuffer(GL.GL_ARRAY_BUFFER,int(self.buffer_id[BUF_TCOORD0]));
    	glBufferData(GL.GL_ARRAY_BUFFER,self.texturecoord0.nbytes,self.texturecoord0,GL.GL_STREAM_DRAW);

    	print "bin buffer   BUF_TCOORD0"
    	glBindBuffer(GL.GL_ARRAY_BUFFER,int(self.buffer_id[BUF_TCOORD1]));
    	glBufferData(GL.GL_ARRAY_BUFFER,self.texturecoord1.nbytes,self.texturecoord1,GL.GL_STREAM_DRAW);
         #GL_TEXTURE_BUFFER??
    	print "bin buffer   BUF_TCOORD1"
    	glBindBuffer(GL.GL_ARRAY_BUFFER,int(self.buffer_id[BUF_TCOORD2]));
    	glBufferData(GL.GL_ARRAY_BUFFER,self.texturecoord2.nbytes,self.texturecoord2,GL.GL_STREAM_DRAW);
    	print "bin buffer   BUF_TCOORD2"

    def setupInstanceMatrice(self,):
        self.matrix_loc=pos=GL.glGetAttribLocation(self.program, "transformmatrix");
        pos1 = pos + 0;
        pos2 = pos + 1;
        pos3 = pos + 2;
        pos4 = pos + 3;
#        GL.glEnableVertexAttribArray(pos1);
#        GL.glEnableVertexAttribArray(pos2);
#        GL.glEnableVertexAttribArray(pos3);
#        GL.glEnableVertexAttribArray(pos4);
#
#        glBindBuffer(GL_ARRAY_BUFFER, VBO_containing_matrices);
#        glVertexAttribPointer(pos1, 4, GL_FLOAT, GL_FALSE, sizeof(GLfloat) * 4 * 4, (void*)(0));
#        glVertexAttribPointer(pos2, 4, GL_FLOAT, GL_FALSE, sizeof(GLfloat) * 4 * 4, (void*)(sizeof(float) * 4));
#        glVertexAttribPointer(pos3, 4, GL_FLOAT, GL_FALSE, sizeof(GLfloat) * 4 * 4, (void*)(sizeof(float) * 8));
#        glVertexAttribPointer(pos4, 4, GL_FLOAT, GL_FALSE, sizeof(GLfloat) * 4 * 4, (void*)(sizeof(float) * 12));
#        glVertexAttribDivisor(pos1, 1);
#        glVertexAttribDivisor(pos2, 1);
#        glVertexAttribDivisor(pos3, 1);
#        glVertexAttribDivisor(pos4, 1);
    #//-----------------------------------------------------------------------------
    def _initGL(self) :
    	print "#// Generation of the buffers and textures ID"
    	self.buffer_id=GL.glGenBuffers(NB_BUFFERS)#GL.glGenBuffers(NB_BUFFERS,self.buffer_id);
    	self.texture_id=GL.glGenTextures(NB_TEXTURES).tolist()
    	print self.buffer_id
    	print self.texture_id
    	if self.useVBO :
         self.generateVBO()#simple()
#         self.activateAllTextures()
         return
    	self.buffer_id=GL.glGenBuffers(NB_BUFFERS);

    	print "#// loading the buffers ",self.buffer_id
    	nbMax=self.nbAtoms if self.nbAtoms>self.nbBonds else self.nbBonds;
    	self.bindAllBuffer()
    	print "#//Generation of the textures used by the shaders"
    	self.texture_id=GL.glGenTextures(NB_TEXTURES);
    	print self.texture_id
    	#// atom positions
    	self.activeTexture(GL.GL_TEXTURE0,TEX_POSITIONS,self.positions)

    	#// atom colors
    	self.activeTexture(GL.GL_TEXTURE1,TEX_COLORS,self.colors)

    	#// atom sizes
    	self.activeTexture(GL.GL_TEXTURE2,TEX_SIZES,self.sizes)

    	#// atom scales
    	self.activeTexture(GL.GL_TEXTURE3,TEX_ASCALES,self.atomScales)

    	#//   bond scales
    	self.activeTexture(GL.GL_TEXTURE4,TEX_BSCALES,self.bondScales)

    	#//   bond shrinks
    	self.activeTexture(GL.GL_TEXTURE5,TEX_SHRINKS,self.shrinks)
    	print "all textures activated "


text:
import os
from PIL import Image
import numpy

"""
    Script that creates result images that puts the label results over input images, with some opacity.

"""

# image_path = "../../IR_images/combined_dataset/"
# label_path = "../result_imgs/for_combined_images/"
# image_path = "../../aerial_datasets/IR_RGB_0.1res/IR_images/combined_dataset/test_images/images" #Input images
image_path = "../../aerial_datasets/IR_RGB_0.1res/IR_images/combined_dataset/test_images/images/" #Ground truth
label_path = "../result_imgs_IR/"
outpath = "../result_imgs_combined/"

images = sorted(os.listdir(image_path))
labels = sorted(os.listdir(label_path))

#IR test = 309
#IR labels = 309

#First rename all image-labels based on what number in sorted list they are? Smallest number = 0, next smallest number = 1 and so on?

processed = 0

for image_name in images:
    for label_name in labels:
        labelnr = label_name.split("testing_image")[1].split(".png")[0]

        if(int(labelnr) == int(processed)):
            label = Image.open(os.path.join(label_path, label_name))
            print("same number!")
            image = Image.open(os.path.join(image_path, image_name))
            # image = Image.open(os.path.join(image_path, image_name)).convert('RGB')
            print(label)
            print(image)
            blend_img = Image.blend(label, image, .7)
            blend_img.save(os.path.join(outpath, image_name))
            # Image.blend(label, image, .7).save(os.path.join(outpath, label_name+"_label"))
    processed +=1

# for label_name in labels:
#     labelnr = label_name.split("testing_image")[1].split(".png")[0]
#     label = Image.open(os.path.join(label_path, label_name))
#     if(len(labelnr) != -1):
#         # labelnr = "400"+labelnr
#         # print(int(float(labelnr)))
#         labelnr = int(float(labelnr)) + 4271
#         print("new label number is: ")
#         print(labelnr)
#     # elif(len(labelnr) == 2):
#     #     labelnr = "40"+labelnr
#     # elif(len(labelnr) == 4):
#     #     labelnr = "4"+labelnr
#     # print(labelnr)
#     for image_name in images:
#         imagenr = image_name.split("map_color_")[1].split(".png")[0]
#         print('imagenr is:')
#         print(imagenr)
#         if(int(labelnr) == int(imagenr)):
#             print("same number!")
#             image = Image.open(os.path.join(image_path, image_name))
#             Image.blend(label, image, .7).save(os.path.join(outpath, label_name))

ext: py

max_stars_repo_name: mathildor/TF-SegNet

reason: import numpy


text:
import argparse

import pegasusio as io
import numpy as np
import pandas as pd

from collections import namedtuple
from typing import List, Dict, Tuple

demux_type_dict = {'SNG': 'singlet', 'DBL': 'doublet', 'AMB': 'unknown'}
SNP = namedtuple('SNP', ['CHROM', 'POS', 'REF', 'ALT'])

def check_colnames(fields: List[str]) -> bool:
    template = ['CHROM', 'POS', 'ID', 'REF', 'ALT', 'QUAL', 'FILTER', 'INFO', 'FORMAT']
    for i, key in enumerate(template):
        if fields[i] != key:
            return False
    return True

def parse_denovo_vcf(input_vcf: str) -> Tuple[List[str], Dict[object, object]]:
    sample_names = None
    snp2geno = {}

    file_ext = input_vcf.split('.')[-1]
    if file_ext == 'gz':
        import gzip
        fin = gzip.open(input_vcf, 'rt')
    else:
        fin = open(input_vcf, 'rt')

    for line in fin:
        if line.startswith('##'):
            continue
        if line.startswith('#'):
            fields = line.strip()[1:].split('\t')
            assert check_colnames(fields)
            sample_names = fields[9:]
        else:
            fields = line.strip().split('\t')
            assert fields[8].startswith('GT')
            snp = SNP(fields[0], fields[1], fields[3], fields[4])
            snp2geno[snp] = [x.split(':')[0] for x in fields[9:]]
    fin.close()
    return sample_names, snp2geno

def calc_matching(denovo_geno: List[str], ref_geno: List[str], mmat: np.array) -> None:
    for i, r_geno in enumerate(ref_geno):
        for j, d_geno in enumerate(denovo_geno):
            mmat[i, j] += (r_geno == d_geno)

def parse_reference_vcf(reference_vcf: str, snp2geno: dict, sample_names: List[str]) -> Tuple[List[str], np.array]:
    nsample = len(sample_names)

    cnt = 0
    nbingo = 0
    ref_names = None
    mmat = None # mmat: matching matrix

    file_ext = reference_vcf.split('.')[-1]
    if file_ext == 'gz':
        import gzip
        fin = gzip.open(reference_vcf, 'rt')
    else:
        fin = open(reference_vcf, 'rt')

    for line in fin:
        if line.startswith('##'):
            continue
        if line.startswith('#'):
            fields = line.strip()[1:].split('\t')
            if 'FORMAT' in fields:
                assert check_colnames(fields)
            else:
                return None, None, False
            ref_names = fields[9:]
            mmat = np.zeros((len(ref_names), nsample), dtype = int)
        else:
            fields = line.strip().split('\t')
            snp = SNP(fields[0], fields[1], fields[3], fields[4])
            d_geno = snp2geno.get(snp, None)
            if d_geno is not None:
                assert fields[8].startswith('GT')
                r_geno = [x.split(':')[0] for x in fields[9:]]
                calc_matching(d_geno, r_geno, mmat)
                nbingo += 1
            cnt += 1
            if cnt % 100000 == 0:
                print("Parsed {0} variants, matched {1} variants.".format(cnt, nbingo))

    fin.close()

    print("\n{0} variants are parsed and {1} SNPs are matched.".format(cnt, nbingo))

    return ref_names, mmat, True

def replace_ref_names(ref_str: str, ref_names: List[str], has_ref_genotypes: bool) -> List[str]:
    if ref_str is not None:
        res_arr = ref_str.split(',')
        assert len(res_arr) == len(ref_names)
        ref_names = res_arr
    elif not has_ref_genotypes:  # Rename freemuxlet's cluster names from "CLUST<num>" to "Donor<num>" format.
        ref_names = ['Donor' + str(int(x[5:]) + 1) for x in ref_names]
    ref_names = ['_ref_' + x for x in ref_names]
    return ref_names

def find_max_matching(ref_names: List[str], sample_names: List[str], mmat: np.array) -> dict:
    import itertools
    import networkx as nx
    from networkx.algorithms import bipartite

    nref = len(ref_names)
    nsample = len(sample_names)

    G = nx.Graph()
    G.add_nodes_from(ref_names, bipartite = 0)
    G.add_nodes_from(sample_names, bipartite = 1)
    G.add_weighted_edges_from([(ref_names[x], sample_names[y], -mmat[x, y]) for x, y in itertools.product(range(nref), range(nsample))])
    result = bipartite.matching.minimum_weight_full_matching(G)

    for j, sample_name in enumerate(sample_names):
        if sample_name not in result:
            i = np.where(mmat[:, j] == mmat[:, j].max())[0][0]
            result[sample_name] = ref_names[i]
            if isinstance(result[ref_names[i]], str):
                result[ref_names[i]] = [result[ref_names[i]]]
            result[ref_names[i]].append(sample_name)

    ref_n2i = {}
    for i, ref_name in enumerate(ref_names):
        ref_n2i[ref_name] = i

    for j, sample_name in enumerate(sample_names):
        i = ref_n2i[result[sample_name]]
        if mmat[i, j] != mmat[:, j].max():
            k = np.where(mmat[:, j] == mmat[:, j].max())[0][0]
            print("Warning: popscle donor {} shares most SNPs with ref donor {}, but matches to ref donor {}!".format(sample_name, ref_names[k][5:], ref_names[i][5:]))

    print()
    for sample_name in sample_names:
        print("Popscle donor {} matches reference donor {}.".format(sample_name, result[sample_name][5:]))
    print()

    return result

def translate_donor_name(inp_str: str, matching: dict) -> str:
    if inp_str == '':
        return inp_str

    res_str = []
    for donor_id in inp_str.split(','):
        res_str.append(matching[f"CLUST{donor_id}"][5:])
    return ','.join(res_str)

def gen_summary_of_status(df):
    type_counts = df['DROPLET.TYPE'].value_counts()
    for type in demux_type_dict.keys():
        if type not in type_counts.index:
            type_counts[type] = 0
    return type_counts

def set_matching_no_reference(sample_numbers: List[str], ref_names: List[str]) -> dict:
    assert len(sample_numbers) == len(ref_names)
    matching = dict()
    pairs = list(zip(sample_numbers, ref_names))
    for (sample_number, ref_name) in pairs:
        matching[sample_number] = ref_name
        matching[ref_name] = sample_number

    return matching

def write_output(assignment_file: str, input_mat_file: str, output_zarr_file: str, matching: dict = None) -> None:
    algorithm = 'demuxlet' if matching is None else 'freemuxlet'

    df = pd.read_csv(assignment_file, sep='\t', header=0, index_col='BARCODE')
    df.index = pd.Index([x.split('-')[0] for x in df.index])

    type_counts = gen_summary_of_status(df)
    print("\nSinglets = {}, doublets = {}, unknown = {}.".format(type_counts['SNG'], type_counts['DBL'], type_counts['AMB']))

    df['demux_type'] = df['DROPLET.TYPE'].apply(lambda s: demux_type_dict[s])
    df['assignment'] = ''
    df.loc[df['demux_type']=='singlet', 'assignment'] = df.loc[df['demux_type']=='singlet', 'SNG.BEST.GUESS']
    if algorithm == 'demuxlet':
        df.loc[df['demux_type']=='doublet', 'assignment'] = df.loc[df['demux_type']=='doublet', 'DBL.BEST.GUESS'].apply(lambda s: ','.join(s.split(',')[:-1]))
    else:
        df.loc[df['demux_type']=='doublet', 'assignment'] = df.loc[df['demux_type']=='doublet', 'DBL.BEST.GUESS']
        f = np.vectorize(translate_donor_name)
        df['assignment'] = f(df['assignment'].astype('str').values, matching)

    data = io.read_input(input_mat_file)
    data.obs['demux_type'] = ''
    data.obs['assignment'] = ''

    idx = data.obs_names.isin(df.index)
    barcodes = data.obs_names[idx]
    df_valid = df.loc[barcodes, ['demux_type', 'assignment']]
    data.obs.loc[idx, 'demux_type'] = df_valid['demux_type'].values
    data.obs.loc[idx, 'assignment'] = df_valid['assignment'].values

    io.write_output(data, output_zarr_file)

if __name__ == '__main__':
    parser = argparse.ArgumentParser(description='Merge popscle result with gene-count matrix, matching donor names if needed.')
    parser.add_argument('demux_res', metavar='(demux_result.best | demux_result.clust1.samples.gz)', help='Demuxlet or Freemuxlet demultiplexing results.')
    parser.add_argument('raw_mat', metavar='raw_feature_bc_matrix.h5', help='Raw gene count matrix in 10x format.')
    parser.add_argument('out_file', metavar='output_result.zarr.zip', help='Output zarr file.')
    parser.add_argument('--cluster-genotypes', metavar = '(demux_result.clust1.vcf.gz | demux_result.clust1.vcf)', dest='cluster_genotypes', help = 'Genotypes detected by freemuxlet from RNA-seq reads.')
    parser.add_argument('--ref-genotypes', metavar = '(reference_genotypes.vcf.gz | reference_genotypes.vcf)', dest = 'ref_genotypes', help = 'Reference genotypes called from exome or genome sequencing data.')
    parser.add_argument('--donor-names', dest = 'ref_names', help = 'A comma-separated list containing donor names that are used to replace the ones in reference_genotypes.vcf.gz. Must match the order in the .vcf.gz file.')
    args = parser.parse_args()

    if args.cluster_genotypes is not None:  # freemuxlet
        sample_names, snp2geno = parse_denovo_vcf(args.cluster_genotypes)
        has_format = False
        if args.ref_genotypes is not None:
            ref_names, mmat, has_format = parse_reference_vcf(args.ref_genotypes, snp2geno, sample_names)

        if has_format:
            ref_names = replace_ref_names(args.ref_names, ref_names, has_ref_genotypes = True)
            matching = find_max_matching(ref_names, sample_names, mmat)
        else:
            ref_names = replace_ref_names(args.ref_names, sample_names, has_ref_genotypes = False)
            matching = set_matching_no_reference(sample_names, ref_names)
        write_output(args.demux_res, args.raw_mat, args.out_file, matching)
    else:  # demuxlet
        write_output(args.demux_res, args.raw_mat, args.out_file)

max_stars_repo_name: Richard-Hansen/cumulus

max_issues_repo_name: Richard-Hansen/cumulus

max_forks_repo_name: Richard-Hansen/cumulus

reason: import numpy,import networkx,from networkx


text:
import copy
import logging
import string

import numpy as np
import matplotlib as mpl
import tensorflow as tf

from IPython.display import HTML
from tensorflow.keras.models import Model
from typing import Callable, TYPE_CHECKING, Union, List, Tuple

if TYPE_CHECKING:  # pragma: no cover
    import keras  # noqa

logger = logging.getLogger(__name__)

def hlstr(string, color="white"):
    """
    Return HTML highlighting text with given color.

    Args:
      string (string): The string to render
      color (string): HTML color for background of the string
    """
    return f"<mark style=background-color:{color}>{string} </mark>"

def colorize(attrs, cmap="PiYG"):
    """
    Compute hex colors based on the attributions for a single instance.

    Args:
      attrs (np.ndarray): Attributions for an instance
      cmap (string): Type of color map to use
    """
    cmap_bound = np.abs(attrs).max()
    norm = mpl.colors.Normalize(vmin=-cmap_bound, vmax=cmap_bound)
    cmap = mpl.cm.get_cmap(cmap)

    colors = list(map(lambda x: mpl.colors.rgb2hex(cmap(norm(x))), attrs))
    return colors

def select_target(ps, ts):
    """
    Select the target

    Args:
        ps (tf.Tensor): Predictions from the model
        ts (tf.Tensor): Target tensor
    """
    ps = tf.linalg.diag_part(tf.gather(ps, ts, axis=1))
    return ps

def gradients_input(model, x, target):
    """
    Calculates the gradients of the target class output with respect
    to each input feature.

    Args:
        model (Model): Tensorflow model.
        x (tf.Tensor): Input data
        target (tf.Tensor or None): Target for which the gradients are calculated
    """
    with tf.GradientTape() as tape:
        tape.watch(x)
        preds = model(x)
        if len(model.output_shape) > 1 and model.output_shape[1] > 1:
            preds = select_target(preds, target)

    grads = tape.gradient(preds, x)

    return grads

def _gradients_layer(
    model,
    layer,
    orig_call,
    x,
    target,
):
    """
    Calculates the gradients of the target class output (or the output if the output dimension is equal to 1)
    with respect to each element of `layer`.

    Args:
        model (keras.models.Model): Tensorflow model.
        layer (keras.layers.Layer): Layer wrt which the gradients are calculated.
        orig_call (Callable): `call` method of the layer.
        x (tf.Tensor): Input data point.
        target (None or tf.Tensor): Target for which the gradients are calculated.
    """

    def watch_layer(layer, tape):
        """
        Make an intermediate hidden `layer` watchable by the `tape`.

        Args:
            layer (keras.layers.Layer): Intermediate layer.
            tape (tf.GradientTape): Tape for computing gradients.
        """

        def dec(func):
            def wrap(*args, **kwargs):
                layer.result = func(*args, **kwargs)
                tape.watch(layer.result)
                return layer.result

            return wrap

        layer.call = dec(layer.call)
        return layer

    with tf.GradientTape() as tape:
        watch_layer(layer, tape)
        preds = model(x)
        if len(model.output_shape) > 1 and model.output_shape[1] > 1:
            preds = select_target(preds, target)

    grads = tape.gradient(preds, layer.result)
    delattr(layer, "result")
    layer.call = orig_call
    return grads

def _sum_integral_terms(step_sizes, grads):
    """
    Sums the terms in the integral path.

    Args:
    step_sizes (list): Weights in the path integral sum.
    grads (tf.Tensor or np.ndarray): Gradients to sum for
                                     each feature.
    """
    input_str = string.ascii_lowercase[1 : len(grads.shape)]
    if isinstance(grads, tf.Tensor):
        step_sizes = tf.convert_to_tensor(step_sizes)
    einstr = "a,a{}->{}".format(input_str, input_str)
    sums = (
        tf.einsum(einstr, step_sizes, grads)
        if isinstance(grads, tf.Tensor)
        else np.einsum(einstr, step_sizes, grads)
    )
    return sums.numpy() if isinstance(grads, tf.Tensor) else sums

def _format_input_baseline(X, baselines):
    """
    Formats baselines to return a numpy array.

    Args:
        X (np.adarray):  Input data points.
        baselines (None or np.ndarray): Baselines.
    """
    if baselines is None:
        bls = np.zeros(X.shape).astype(X.dtype)
    else:
        bls = baselines.astype(X.dtype)

    return bls

def _format_target(target, nb_samples):
    """
    Formats target to return a list.

    Args:
        target (None or np.ndarray): Original target.
        nb_samples (int): Number of samples in the batch.
    """
    return [t.astype(int) for t in target] if target is not None else None

class IntegratedGradientsText(object):
    def __init__(
        self,
        model,
        layer,
        n_steps,
        internal_batch_size,
    ):
        """
        An implementation of Integrated gradients method for Tensorflow/Keras text models.
        For details of the method see the original paper:
        https://arxiv.org/abs/1703.01365 .

        Args:
            model (keras.Model): Tensorflow model.
            layer (keras.layers.Layer): Layer with respect to which the gradients are calculated.
                                        If not provided, the gradients are calculated with respect to the input.
            n_steps (int): Number of step in the path integral approximation from the
                           baseline to the input instance.
            internal_batch_size (None or int): Batch size for the internal batching.
        """
        self.model = model
        self.input_dtype = self.model.input.dtype
        self.layer = layer
        self.n_steps = n_steps
        self.internal_batch_size = internal_batch_size

    def step_sizes(self, n):
        """
        Compute step sizes.

        Args:
            n (int): The number of integration steps
        """
        return list(0.5 * np.polynomial.legendre.leggauss(n)[1])

    def alphas(self, n):
        """
        Compute alpha coefficients.

        Args:
            n (int): The number of integration steps
        """
        return list(0.5 * (1 + np.polynomial.legendre.leggauss(n)[0]))

    def explain_instance(
        self,
        X,
        baselines,
        target,
    ):
        """
        Compute attributions for each input feature.

        Args:
            X (np.ndarray): Instance for which integrated gradients attribution are computed.
            baselines (None or np.ndarray): Baselines for each instance.
            target (None or np.ndarray):  Defines which element of the model output is considered to compute the gradients.
        """

        if (
            len(self.model.output_shape) == 1 or self.model.output_shape[1] == 1
        ) and target is None:
            logger.warning(
                "It looks like you are passing a model with a scalar output and target is set to `None`."
                "If your model is a regression model this will produce correct attributions. If your model "
                "is a classification model, targets for each datapoint must be defined. "
                "Not defining the target may lead to incorrect values for the attributions."
                "Targets can be either the true classes or the classes predicted by the model."
            )

        nb_samples = len(X)

        # format and check inputs and targets
        baselines = _format_input_baseline(X, baselines)
        target = _format_target(target, nb_samples)

        # defining integral method
        step_sizes_func, alphas_func = self.step_sizes, self.alphas
        step_sizes, alphas = step_sizes_func(self.n_steps), alphas_func(self.n_steps)

        # construct paths and prepare batches
        paths = np.concatenate(
            [baselines + alphas[i] * (X - baselines) for i in range(self.n_steps)],
            axis=0,
        )
        if target is not None:
            target_paths = np.concatenate([target for _ in range(self.n_steps)], axis=0)
            paths_ds = tf.data.Dataset.from_tensor_slices((paths, target_paths)).batch(
                self.internal_batch_size
            )
        else:
            paths_ds = tf.data.Dataset.from_tensor_slices(paths).batch(
                self.internal_batch_size
            )
        paths_ds.prefetch(tf.data.experimental.AUTOTUNE)

        # fix orginal call method for layer
        if self.layer is not None:
            orig_call = self.layer.call
        else:
            orig_call = None

        # calculate gradients for batches
        batches = []
        for path in paths_ds:

            if target is not None:
                paths_b, target_b = path
            else:
                paths_b, target_b = path, None

            if self.layer is not None:
                grads_b = _gradients_layer(
                    self.model,
                    self.layer,
                    orig_call,
                    tf.dtypes.cast(paths_b, self.input_dtype),
                    target_b,
                )
            else:
                grads_b = gradients_input(
                    self.model, tf.dtypes.cast(paths_b, self.input_dtype), target_b
                )

            batches.append(grads_b)

        # tf concatatation
        grads = tf.concat(batches, 0)
        shape = grads.shape[1:]
        if isinstance(shape, tf.TensorShape):
            shape = tuple(shape.as_list())

        # invert sign of gradients for target 0 examples if classifier returns only positive class probability
        if (
            len(self.model.output_shape) == 1 or self.model.output_shape[1] == 1
        ) and target is not None:
            sign = 2 * target_paths - 1
            grads = np.array([s * g for s, g in zip(sign, grads)])

        grads = tf.reshape(grads, (self.n_steps, nb_samples) + shape)

        # sum integral terms and scale attributions
        sum_int = _sum_integral_terms(step_sizes, grads.numpy())
        if self.layer is not None:
            layer_output = self.layer.output
            model_layer = Model(self.model.input, outputs=layer_output)
            norm = (model_layer(X) - model_layer(baselines)).numpy()
        else:
            norm = X - baselines
        attributions = norm * sum_int

        attributions = attributions.sum(axis=2)

        return attributions

    def visualize(self, attrs_i, words):
        """
        Visualize the attributions as computed by integrated gradients.

        Args:
            attrs_i (np.ndarray): Attributions of each feature as computed by integrated gradients.
            words (list): List of words in the text corresponding to each attribution.
        """
        colors = colorize(attrs_i)
        return HTML("".join(list(map(hlstr, words, colors))))


text:
"""
File: examples/model/conditional_fit_model.py
Author: Keith Tauscher
Date: 4 Oct 2020

Description: Example script showing how objects of the
             SingleConditionalFitModel class relate to the model given at
             initialization.
"""
import numpy as np
import matplotlib.pyplot as pl
from distpy import ChiSquaredDistribution, GaussianDistribution,\
    DistributionSet
from pylinex import LegendreBasis, SumModel, ProductModel, LorentzianModel,\
    BasisModel, GaussianModel, SingleConditionalFitModel

fontsize = 24
triangle_plot_fontsize = 10
include_noise = True
first_seed = None
second_seed = None

half_num_channels = 50
num_channels = (2 * half_num_channels) + 1
x_values = np.linspace(-1, 1, num_channels)

gaussian_model = GaussianModel(x_values)
lorentzian_model = LorentzianModel(x_values)
num_basis_vectors = 10
basis = LegendreBasis(num_channels, num_basis_vectors - 1)
basis_model = BasisModel(basis)
sum_model = SumModel(['basis', 'lorentzian'], [basis_model, lorentzian_model])

gaussian_parameters = np.array([1, 0.2, 1])
lorentzian_parameters = np.array([15, 0, 0.5])
np.random.seed(first_seed)
basis_parameters = np.random.normal(0, 100, size=num_basis_vectors)

noise_level = 1
error = np.ones((num_channels,)) * noise_level
np.random.seed(second_seed)
if include_noise:
    noise = np.random.normal(0, 1, size=error.shape) * error
else:
    noise = np.zeros_like(error)

full_model = ProductModel(['gaussian', 'sum'], [gaussian_model, sum_model])
full_parameters = np.concatenate([gaussian_parameters, basis_parameters,\
    lorentzian_parameters])

data = full_model(full_parameters) + noise
unknown_name_chain = ['sum', 'basis']

model = SingleConditionalFitModel(full_model, data, error, unknown_name_chain)
parameters =\
    np.concatenate([gaussian_parameters, lorentzian_parameters,\
    [num_basis_vectors]])
print("model.parameters={}".format(model.parameters))

degrees_of_freedom = num_channels - num_basis_vectors
chi_squared = np.sum(np.power((data - model(parameters)) / error, 2)) /\
    degrees_of_freedom
chi_squared_distribution =\
    ChiSquaredDistribution(degrees_of_freedom, reduced=True)
probability_level = 0.95
threshold =\
    chi_squared_distribution.left_confidence_interval(probability_level)[1]

try:
    assert(chi_squared < threshold)
except:
    print(("Chi squared, {0:.4g}, was above threshold for probability " +\
        "level {1:.2f}, {2:.4g}.").format(chi_squared, probability_level,\
        threshold))

(recreation, conditional_mean, conditional_covariance) = model(parameters,\
    return_conditional_mean=True, return_conditional_covariance=True)
conditional_distribution =\
    GaussianDistribution(conditional_mean, conditional_covariance)

fig = pl.figure(figsize=(12,9))
ax = fig.add_subplot(211)
ax.scatter(x_values, data, color='k', label='full model + noise')
ax.plot(x_values, recreation, color='r', label='SingleConditionalFitModel')
ax.set_xlabel('x', size=fontsize)
ax.set_ylabel('Total y', size=fontsize)
ax.set_title('SingleConditionalFitModel test')
ax.legend(fontsize=fontsize)
ax = fig.add_subplot(212)
ax.scatter(x_values, data - model(parameters), color='k',\
    label='$\chi^2={:.4f}$'.format(chi_squared))
ax.legend(fontsize=fontsize)
ax.set_xlabel('x', size=fontsize)
ax.set_ylabel('Residual y', size=fontsize)

num_samples = int(1e6)
marginalized_parameters = [parameter for parameter in full_model.parameters\
    if parameter not in model.parameters]
conditional_distribution_set =\
    DistributionSet([(conditional_distribution, marginalized_parameters)])

fig = conditional_distribution_set.triangle_plot(num_samples,\
    reference_value_mean=basis_parameters, nbins=200,\
    fontsize=triangle_plot_fontsize,\
    contour_confidence_levels=[0.68, 0.95, 0.997], plot_type='contour',\
    parameter_renamer=(lambda parameter: '_'.join(parameter.split('_')[-2:])))
conditional_distribution_set.triangle_plot(num_samples, fig=fig,\
    reference_value_mean=basis_parameters, nbins=200,\
    fontsize=triangle_plot_fontsize, plot_type='histogram',\
    parameter_renamer=(lambda parameter: '_'.join(parameter.split('_')[-2:])))

pl.show()

max_stars_repo_name: CU-NESS/pylinex

max_issues_repo_name: CU-NESS/pylinex

max_forks_repo_name: CU-NESS/pylinex


text:
#! /usr/bin/env python3
"""
Focal loss: a case study with credit card fraud detection.
Credits to:
    - https://www.kaggle.com/mlg-ulb/creditcardfraud/code
    - https://www.dlology.com/blog/multi-class-classification-with-focal-loss-for-imbalanced-datasets/
    - https://www.tensorflow.org/tutorials/structured_data/imbalanced_data
"""
from tensorflow.keras.models import Sequential
from tensorflow.keras.layers import Dense, Dropout
from tensorflow.keras.optimizers import Adam
from tensorflow.keras.losses import BinaryCrossentropy
from tensorflow.keras.metrics import TruePositives, FalseNegatives, Recall
from tensorflow_addons.losses import SigmoidFocalCrossEntropy
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from prettytable import PrettyTable
import matplotlib as mpl
import matplotlib.pyplot as plt
import seaborn as sns
import pandas as pd
import numpy as np

#
# PARAMETERS
#
np.random.seed(42)
#mpl.rcParams['figure.figsize'] = (12, 10)
plt.rcParams['figure.dpi'] = 300
plt.rcParams['savefig.dpi'] = 300
colors = plt.rcParams['axes.prop_cycle'].by_key()['color']

# Set seaborn theme.
sns.set_theme(style='darkgrid')

EPOCHS = 5
BATCH_SIZE = 2048

#
# DATA LOADING/CLEANSING
#

# Load data from CSV file.
df_raw = pd.read_csv('creditcardfraud/creditcard.csv')
n_samples = len(df_raw)
print(f'Num. of samples: {n_samples}.')

# Check size of samples.
df_pos = df_raw[df_raw['Class'] == 1]
n_pos_samples = len(df_pos)
pos_ratio = 100 * n_pos_samples / n_samples
print(f'Num. of positive samples: {n_pos_samples} ({pos_ratio:.2f}% of total).')

# Drop useless data and convert amount to log space.
df_cleaned = df_raw.copy()
df_cleaned.pop('Time')
df_cleaned['log-amount'] = np.log(df_cleaned.pop('Amount') + 0.001)

# Double train/test split for testing and validation data.
df_train, df_test = train_test_split(df_cleaned, test_size = 0.2, shuffle = True)
df_train, df_valid = train_test_split(df_train, test_size = 0.2, shuffle = True)

print(f'Size of training data: {len(df_train)}.')
print(f'Size of validation data: {len(df_valid)}.')
print(f'Size of test data: {len(df_test)}.')

# Extract labels and features from data.
labels_train = np.array(df_train.pop('Class'))
labels_valid = np.array(df_valid.pop('Class'))
labels_test = np.array(df_test.pop('Class'))
features_train = np.array(df_train)
features_valid = np.array(df_valid)
features_test = np.array(df_test)

# Normalize data.
scaler = StandardScaler()
features_train = scaler.fit_transform(features_train)
features_valid = scaler.transform(features_valid)
features_test = scaler.transform(features_test)

# Enforce lower/upper bounds.
features_train = np.clip(features_train, -5, 5)
features_valid = np.clip(features_valid, -5, 5)
features_test = np.clip(features_test, -5, 5)
n_features = features_train.shape[-1]

#
# MODEL TRAINING
#

# Model parameters.
opt = Adam(learning_rate = 1e-3)

metrics = [
    TruePositives(name = 'tp'),
    FalseNegatives(name = 'fn'),
    Recall(name = 'recall')
]

losses = [
    BinaryCrossentropy(),
    SigmoidFocalCrossEntropy(gamma = 2, alpha = 4)
]

loss_names = [
    'binary cross-entropy',
    'focal loss'
]

logs_loss = []
logs_recall = []

for loss in losses:

    # Setup/compile the model.
    model = Sequential()
    model.add(Dense(16, input_dim = n_features, activation = 'relu',
        kernel_initializer = 'he_uniform'))
    model.add(Dropout(0.5))
    model.add(Dense(1, activation = 'sigmoid'))
    model.compile(optimizer = opt, loss = loss, metrics = metrics)

    # Fit the model.
    logs = model.fit(features_train, labels_train, validation_data = (features_valid,
        labels_valid), epochs = EPOCHS, verbose = 0)

    logs_loss.append(logs.history['loss'])
    logs_recall.append(logs.history['recall'])

    # Evaluate the model.
    eval_train = model.evaluate(features_train, labels_train, verbose = 0)
    eval_test = model.evaluate(features_valid, labels_valid, verbose = 0)

    table = PrettyTable()
    table.field_names = ['Data', 'Loss', 'TruePositives', 'FalseNegatives', 'Recall']

    for stage, eval_info in zip(('training', 'test'), (eval_train, eval_test)):

        row = [stage]

        for ii, lbl in enumerate(model.metrics_names):
            row.append(f'{eval_info[ii]:.3f}')

        table.add_row(row)

    print('\n')
    print(table)

#
# Plotting
#

fig, ax = plt.subplots()
for log_loss, name in zip(logs_loss, loss_names):
    ax.plot(log_loss, label = name)
plt.xlabel('Epoch')
plt.ylabel('Loss')
plt.legend()
plt.title('Loss values with different loss functions')
fig.savefig(f'losses.png')

fig, ax = plt.subplots()
for log_recall, name in zip(logs_recall, loss_names):
    ax.plot(log_recall, label = name)
plt.xlabel('Epoch')
plt.ylabel('Recall')
plt.legend()
plt.title('Evaluation metric with different loss functions')
fig.savefig(f'accuracies.png')

max_stars_repo_name: andrea-ci/misc-stuff

max_issues_repo_name: andrea-ci/misc-stuff

max_forks_repo_name: andrea-ci/misc-stuff


text:
import numpy as np
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process import kernels

class Gaussian_Process_Regression(GaussianProcessRegressor):
    def xsample2meanstd(self, _xsample, _ysample, _x, normalization = True):
        if (normalization):
            xmin = np.amin(_x)
            xmax = np.amax(_x)
            xmid = 0.5*(xmin + xmax)
            _x = (_x - xmid)/(xmax - xmin)
            _xsample = (_xsample - xmid)/(xmax - xmin)
            ymin = np.amin(_ysample)
            ymax = np.amax(_ysample)
            ymid = 0.5*(ymin + ymax)
            _ysample = (_ysample - ymid)/(ymax - ymin)
        _Xsample = _xsample.reshape(-1,1)
        self.fit(_Xsample, _ysample)
        _X = _x.reshape(-1,1)
        mean, std = self.predict(_X, return_std=True) #Get the maen and standard error for the first trial set
        if (normalization):
            mean = mean*(ymax - ymin) + ymid
            std = std*(ymax -ymin)
        return mean, std


text:
import torch
import pickle
import numpy as np
from PIL import Image
import random
from dataset.cifar10 import TransformTwice

class ImdbDatasetForGen(torch.utils.data.Dataset):
    """Imdb dataloader, output image and gender label"""

    def __init__(self, lines, transform=None):
        pathlist, labellist = [], []
        for line in lines:
            path, _, label = line.split()
            pathlist.append(path)
            labellist.append(label)
        self.pathlist = pathlist
        self.labellist = labellist
        self.transform = transform

    def __getitem__(self, index):
        path = self.pathlist[index]
        img = Image.open(path).convert("RGB")
        label = int(self.labellist[index])
        if self.transform is not None:
            img = self.transform(img)
        return img, label

    def __len__(self):
        return len(self.labellist)

class ImdbDatasetwoLabel(torch.utils.data.Dataset):
    """Imdb dataloader, output image and gender label"""

    def __init__(self, lines, transform=None):
        pathlist = []
        for line in lines:
            path, _ , _ = line.split()
            pathlist.append(path)
        self.pathlist = pathlist
        self.transform = transform

    def __getitem__(self, index):
        path = self.pathlist[index]
        img = Image.open(path).convert("RGB")
        if self.transform is not None:
            img = self.transform(img)
        return img

    def __len__(self):
        return len(self.pathlist)

def get_imdb(tr_data_path, te_data_path, n_labeled, transform_train, transform_val):
    with open(tr_data_path, "r") as f:
        lines = f.readlines()
    random.shuffle(lines)
    lines_labeled = lines[:n_labeled]
    lines_unlabeled = lines[n_labeled:]
    with open(te_data_path, "r") as f:
        te_lines = f.readlines()
    labeled_ds = ImdbDatasetForGen(lines_labeled, transform=transform_train)
    unlabeled_ds = ImdbDatasetwoLabel(lines_labeled, transform=TransformTwice(transform_train))
    test_ds = ImdbDatasetForGen(te_lines, transform=transform_val)

    print (f"#Labeled: {len(lines_labeled)} #Unlabeled: {len(lines_unlabeled)} #Val: {len(te_lines)}")

    return labeled_ds, unlabeled_ds, test_ds, test_ds


text:
import numpy as np, itertools
from .type import is_dict, is_seq_of

def concat_list_of_array(x, axis=0):
    assert is_seq_of(x)
    if len(x) == 0:
        return None
    if isinstance(x[0], dict):
        ret = {}
        for k in x[0].keys():
            ret[k] = concat_list_of_array([_[k] for _ in x], axis)
        return ret
    elif isinstance(x[0], (list, tuple)):
        ret = []
        for i in range(len(x[0])):
            ret.append(concat_list_of_array([_[i] for _ in x]))
        return type(x[0])(ret)
    elif isinstance(x[0], np.ndarray):
        if x[0].ndim == 0:
            return np.array(x)
        else:
            return np.concatenate(x, axis=axis)
    elif np.isscalar(x[0]):
        return np.array(x)
    else:
        import torch
        if is_seq_of(x, torch.Tensor):
            return torch.cat(x, dim=axis)
        else:
            raise NotImplementedError(x)

def concat_dict_of_list_array(x, axis=0):
    assert isinstance(x, dict)
    return {k: concat_list_of_array(x[k], axis) for k in x}

def stack_list_of_array(x, axis=0):
    assert is_seq_of(x), type(x)
    if len(x) == 0:
        return None
    if isinstance(x[0], dict):
        ret = {}
        for k in x[0].keys():
            ret[k] = stack_list_of_array([_[k] for _ in x], axis)
        return ret
    elif isinstance(x[0], np.ndarray):
        if x[0].ndim == 0:
            return np.array(x)
        else:
            return np.stack(x, axis=axis)
    elif np.isscalar(x[0]):
        if hasattr(x, 'dtype'):
            return np.array(x, dtype=x.dtype)
        else:
            return np.array(x)
    else:
        import torch
        if is_seq_of(x, torch.Tensor):
            return torch.stack(x, dim=axis)
        else:
            raise NotImplementedError(x)

def stack_dict_of_list_array(x, axis=0):
    assert isinstance(x, dict)
    return {k: stack_list_of_array(x[k], axis) for k in x}

def concat_seq(in_list, dtype=list):
    assert dtype in [list, tuple]
    return dtype(itertools.chain(*in_list))

def concat_list(in_list):
    return concat_seq(in_list, list)

def repeat_interleave(x, n, axis):
    if isinstance(x, np.ndarray):
        return np.repeat(x, n, axis=axis)
    elif isinstance(x, dict):
        return {key: repeat_interleave(x[key], n, axis=axis) for key in x}
    elif isinstance(x, (list, tuple)):
        return type(x)([repeat_interleave(_, n, axis=axis) for _ in x])
    else:
        import torch
        if isinstance(x, torch.Tensor):
            return torch.repeat_interleave(x, n, dim=axis)
        else:
            raise NotImplementedError()

max_stars_repo_name: art-e-fact/ManiSkill-Learn

max_issues_repo_name: art-e-fact/ManiSkill-Learn

max_forks_repo_name: art-e-fact/ManiSkill-Learn


text:
import numpy as np
from mpi4py import MPI
from keras.models import Model
from keras.layers import Input
from create_neural_network_models import create_model_max, create_discriminator
from autoencoder_training import train_autoencoder, train_AAE_surr

comm = MPI.COMM_WORLD
rank = comm.Get_rank()
size = comm.Get_size()

#possibly also doing pretraining or not

################################################################
####                 Load data                              ####
################################################################
filename_simp='Sample_data/X_rand.npy'
filename_C='Sample_data/C_rand.npy'

if rank==0:
    X_simp=np.load(filename_simp)

    C=np.load(filename_C)

    Index=np.arange(C.size)
    np.random.shuffle(Index)

    X_simp=X_simp[Index[:4000],:,:]
    C=C[Index[:4000]]

    volfrac=np.mean(X_simp)
else:
    X_simp=None
    C=None
    Index=None
    volfrac=None

comm.Barrier()

X_simp=comm.bcast(X_simp,root=0)
C=comm.bcast(C,root=0)
volfrac=comm.bcast(volfrac,root=0)

(_,nely,nelx)=X_simp.shape

n=len(C)
perrank=int(n/size)

X_simp_rank=X_simp[rank*perrank:(rank+1)*perrank,:,:]
C_rank=C[rank*perrank:(rank+1)*perrank,np.newaxis]

W_simp_rank=np.ones(X_simp_rank.shape)
W_simp_rank=W_simp_rank+1-2*np.abs(X_simp_rank-0.5)

perrank_train=int(0.8*perrank)
perrank_val=perrank-perrank_train

X_simp_rank_train=X_simp_rank[:perrank_train,:,:]
X_simp_rank_val=X_simp_rank[perrank_train:,:,:]

W_simp_rank_train=W_simp_rank[:perrank_train,:,:]
W_simp_rank_val=W_simp_rank[perrank_train:,:,:]

C_rank_train=C_rank[:perrank_train,:]

edim=100
pre=1
dis=1
surr=1

n_epochs=250
n_epochs_pre=25
typ='max'

if rank==0:
    print('edim={}, pre={}, AAE={}, surr={}, typ='.format(edim,pre,dis,surr)+typ)
comm.Barrier()
################################################################################################
####                                   Training                                             ####
################################################################################################
Eo, Ei, Di, Do=create_model_max(nelx,nely,edim)
X_rank_train=X_simp_rank_train
X_rank_val=X_simp_rank_val
W_rank_train=W_simp_rank_train
W_rank_val=W_simp_rank_val

Design_space=Input((nely,nelx))
Latent_space=Input((edim,))
Middle_layer=Input((Ei.input.shape[1],))
AEo=Model(Design_space, Do(Eo(Design_space)))
AEi=Model(Middle_layer, Di(Ei(Middle_layer)))
# Train the outer layers
AEo_weights=train_autoencoder(AEo,X_rank_train,np.ones(X_rank_train.shape), comm, rank, size, perrank_train,n_epochs_pre)
AEo.set_weights(AEo_weights)
Xm_rank=Eo.predict(X_rank_train)
# Train the inner layers
AEi_weights=train_autoencoder(AEi,Xm_rank,np.ones(Xm_rank.shape), comm, rank, size, perrank_train,n_epochs_pre)
AEi.set_weights(AEi_weights)
# Put together networks
Encoder=Model(Design_space,Ei(Eo(Design_space)))
Decoder=Model(Latent_space,Do(Di(Latent_space)))
Autoencoder=Model(Design_space, Decoder(Encoder(Design_space)))
# Train main networks
Discriminator=create_discriminator(edim)
Surrogate=create_discriminator(edim)
AE_weights=train_AAE_surr(Autoencoder,Encoder,Discriminator,Surrogate,X_rank_train, W_rank_train,C_rank_train,  comm,rank, size, perrank_train,n_epochs)
Autoencoder.set_weights(AE_weights)
#Save results
if rank==0:
    Decoder.save('Sample_data/Decoder_random.h5')

max_stars_repo_name: julianschumann/ae-opt

max_issues_repo_name: julianschumann/ae-opt

max_forks_repo_name: julianschumann/ae-opt