numpy.float32

Here are the examples of the python api numpy.float32 taken from open source projects. By voting up you can indicate which examples are most useful and appropriate.

200 Examples 7

Example 101

Project: fos-legacy
Source File: tree.py
View license
    def __init__(self, vertices,
                 connectivity,
                 colors = None,
                 vertices_width = None,
                 affine = None,
                 force_centering = False,
                 *args, **kwargs):
        """ A tree
        
        vertices : Nx3
            3D Coordinates x,y,z
        connectivity : Mx1
            Tree topology
        colors : Nx4 or 1x4
            Per connection color
        vertices_width : N
            Per vertex width
        affine : 4x4
            Affine transformation of the actor

        """
        super(Tree, self).__init__()
        
        self.affine = np.eye(4, dtype = np.float32)
        #self._update_glaffine()
        
        self.vertices = vertices
        if force_centering:
            self.vertices = self.vertices - np.mean(self.vertices, axis = 0)

        self.connectivity = connectivity

        # unfortunately, we need to duplicate vertices if we
        # want per line color
        self.vertices = self.vertices[self.connectivity,:]

        # we want per line color
        # duplicating the color array, we have the colors per vertex
        self.colors =  np.repeat(colors, 2, axis=0)

        # we have a simplified connectivity
        self.connectivity = np.array( range(len(self.vertices)), dtype = np.uint32 )
        

        # this coloring section is for per/vertex color
#        if colors == None:
#            # default colors
#            self.colors = np.array( [[255,255,255,255]], dtype = np.float32).repeat(len(self.vertices),axis=0) / 255.
#        else:
#            if len(colors) == 1:
#                self.colors = np.array( [colors], dtype = np.float32).repeat(len(self.vertices),axis=0) / 255.
#            else:
#                assert(len(colors) == len(self.vertices))
#                self.colors = colors


       # the sample 1d texture data array
        if not vertices_width is None:
            self.mytex = vertices_width.astype( np.float32 )
        else:
            self.mytex = np.ones( len(self.vertices), dtype = np.float32 )

        self.mytex_ptr = self.mytex.ctypes.data
#        self.make_aabb(margin = 0)
        
        # create indicies, seems to be slow with nested loops
        self.indices = self.connectivity
        self.indices_ptr = self.indices.ctypes.data
        self.indices_nr = self.indices.size
        
        # duplicate colors to make it "per vertex"
        self.colors_ptr = self.colors.ctypes.data
        
        self.vertices_ptr = self.vertices.ctypes.data
        self.mode = GL_LINES
        self.type = GL_UNSIGNED_INT
        
        # VBO related
        self.vertex_vbo = GLuint(0)        
        glGenBuffers(1, self.vertex_vbo)
        glBindBuffer(GL_ARRAY_BUFFER_ARB, self.vertex_vbo)
        # copy the vertex data to our buffer
        #glBufferData(GL_ARRAY_BUFFER, ADT.arrayByteCount(circle), ADT.voidDataPointer(circle), GL_STATIC_DRAW_ARB)
        #glBufferData(GL_ARRAY_BUFFER, 8 * sizeof(GLfloat), diamond, GL_STATIC_DRAW);
        # print "vertices size", self.vertices.size
        
        # 4 * because float32 has 4 bytes
        glBufferData(GL_ARRAY_BUFFER_ARB, 4 * self.vertices.size, self.vertices_ptr, GL_STATIC_DRAW)
        
        # /* Specify that our coordinate data is going into attribute index 0, and contains three floats per vertex */
        glVertexAttribPointer(0, 3, GL_FLOAT, GL_FALSE, 0, 0)

        # for indices
        self.indices_vbo = GLuint(0)
        glGenBuffers(1, self.indices_vbo)
        glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, self.indices_vbo)
        # uint32 has 4 bytes
        glBufferData(GL_ELEMENT_ARRAY_BUFFER, 4 * self.indices_nr, self.indices_ptr, GL_STATIC_DRAW)

        # for colors
        self.colors_vbo = GLuint(0)        
        glGenBuffers(1, self.colors_vbo)
        glBindBuffer(GL_ARRAY_BUFFER, self.colors_vbo)

        glBufferData(GL_ARRAY_BUFFER, 4 * self.colors.size, self.colors_ptr, GL_STATIC_DRAW)
        glVertexAttribPointer(1, 4, GL_FLOAT, GL_FALSE, 0, 0)

        # texture init
        self.init_texture()

        self.shader = get_shader()

Example 102

Project: fofix
Source File: Image.py
View license
def draw3Dtex(image, vertex, texcoord, coord = None, scale = None, rot = None, color = (1,1,1), multiples = False, alpha = False, depth = False, vertscale = 0):
    '''
    Simplifies tex rendering

    @param image: self.xxx - tells the system which image/resource should be mapped to the plane
    @param vertex: (Left, Top, Right, Bottom) - sets the points that define where the plane will be drawn
    @param texcoord: (Left, Top, Right, Bottom) - sets where the texture should be drawn on the plane
    @param coord: (x,y,z) - where on the screen the plane will be rendered within the 3d field
    @param scale: (x,y,z) - scales an glplane how far in each direction

    @param rot: (degrees, x-axis, y-axis, z-axis)
    a digit in the axis is how many times you want to rotate degrees around that axis

    @param color: (r,g,b) - sets the color of the image when rendered
    0 = No Color, 1 = Full color

    @param multiples: True/False
    defines whether or not there should be multiples of the plane drawn at the same time
    only really used with the rendering of the notes, keys, and flames

    @param alpha: True/False - defines whether or not the image should have black turned into transparent
    only really used with hitglows and flames

    @param depth: True/False - sets the depth by which the object is rendered
    only really used by keys and notes

    @param vertscale: # - changes the yscale when setting vertex points
    only really used by notes
    '''


    if not isinstance(image, ImgDrawing):
        return

    if alpha:
        glBlendFunc(GL_SRC_ALPHA, GL_ONE)

    if len(color) == 4:
        col_array  = np.array([[color[0],color[1],color[2], color[3]],
                            [color[0],color[1],color[2], color[3]],
                            [color[0],color[1],color[2], color[3]],
                            [color[0],color[1],color[2], color[3]]], dtype=np.float32)
    else:
        col_array  = np.array([[color[0],color[1],color[2], 1],
                            [color[0],color[1],color[2], 1],
                            [color[0],color[1],color[2], 1],
                            [color[0],color[1],color[2], 1]], dtype=np.float32)

    glEnable(GL_TEXTURE_2D)
    image.texture.bind()

    if multiples:
        glPushMatrix()

    if coord is not None:
        glTranslate(coord[0], coord[1], coord[2])
    if rot is not None:
        glRotate(rot[0], rot[1], rot[2], rot[3])
    if scale is not None:
        glScalef(scale[0], scale[1], scale[2])

    if depth:
        glDepthMask(1)

    if not isinstance(vertex, np.ndarray):
        vertex = np.array(
          [[ vertex[0],  vertscale, vertex[1]],
            [ vertex[2],  vertscale, vertex[1]],
            [ vertex[0], -vertscale, vertex[3]],
            [ vertex[2], -vertscale, vertex[3]]], dtype=np.float32)

    if not isinstance(texcoord, np.ndarray):
        texcoord = np.array(
          [[texcoord[0], texcoord[1]],
            [texcoord[2], texcoord[1]],
            [texcoord[0], texcoord[3]],
            [texcoord[2], texcoord[3]]], dtype=np.float32)

    cmgl.drawArrays(GL_TRIANGLE_STRIP, vertices=vertex, colors=col_array, texcoords=texcoord)

    if depth:
        glDepthMask(0)

    if multiples:
        glPopMatrix()

    glDisable(GL_TEXTURE_2D)

    if alpha:
        glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA)

Example 103

Project: Theano-MPI
Source File: googlenet.py
View license
    def __init__(self,config):
        ModelBase.__init__(self)
        
        self.verbose = config['verbose']   
        self.config = config
        if self.verbose: print 'GoogLeNet 7/5'
        
        batch_size = config['batch_size']
        input_width = config['input_width']
        input_height = config['input_height']
        n_softmax_out=config['n_softmax_out']
        
        
        self.name = 'googlenet'
        self.batch_size = batch_size
        self.input_width = input_width
        self.input_height = input_height
        self.n_softmax_out = n_softmax_out
        self.lrn_func = CrossChannelNormalization()


        x = T.ftensor4('x')
        y = T.lvector('y')
        lr = T.scalar('lr')
        
        self.x = x # c01b
        self.y = y
        self.lr = lr
        
        layers = []
        params = []
        weight_types = []        
        

        conv_7x7 = ConvPool_LRN(input=x,
                                image_shape=(3, 224, 224, batch_size), #c01b (3, 224, 224, batch_size)
                                filter_shape=(3, 7, 7, 64),
                                convstride=2, padsize=3,
                                poolsize=3, poolstride=2, poolpad=1,
                                W = Weight((3, 7, 7, 64), mean = 0.0, std=0.1 ),
                                b = Weight((64,), mean = 0.2 , std=0 ), 
                                lrn=True,
                                lib_conv='cudnn',
                                )
        layers.append(conv_7x7)
        params += conv_7x7.params
        weight_types += conv_7x7.weight_type                    
        # output shape = (112x112x64)
        # output shape = (56x56x64)
                                  
        conv_r3x3 = Conv(input=conv_7x7.output,
                        image_shape=(64,56,56,batch_size),
                        filter_shape=(64, 1, 1, 64),
                        convstride=1, padsize=0,
                        W = Weight((64, 1, 1, 64), mean = 0.0, std=0.1 ),
                        b = Weight((64,), mean = 0.2 , std=0 ),
                        lib_conv='cudnn',
                        )                                           

        layers.append(conv_r3x3)
        params += conv_r3x3.params
        weight_types += conv_r3x3.weight_type   
        # output shape = (56x56x64)                       

                                
                                        
        conv_3x3 = ConvPool_LRN(input=conv_r3x3.output,
                                image_shape=(64,56,56,batch_size),
                                filter_shape=(64, 3, 3, 192),
                                convstride=1, padsize=1,
                                poolsize=3, poolstride=2, poolpad=1,
                                W = Weight((64, 3, 3, 192), mean = 0.0, std=0.03 ),
                                b = Weight((192,), mean = 0.2 , std=0 ), 
                                lrn=True,
                                lib_conv='cudnn',
                                )                                           

        layers.append(conv_3x3)
        params += conv_3x3.params
        weight_types += conv_3x3.weight_type  
        # output shape = (56x56x192) 
        # output shape = (28x28x192)


        incep3a = Incept(conv_3x3.output,input_shape = (192,28,28,batch_size))
        
        layers += incep3a.layers
        params += incep3a.params
        weight_types += incep3a.weight_types        
        print 'incep3a output shape: (28x28x256)'
        # output shape = (28x28x256)
        
        incep3b = Incept(incep3a.output,input_shape = (256,28,28,batch_size),
                          n1x1=128, nr3x3=128, n3x3=192, 
                          nr5x5=32, n5x5=96, npj=64)
        
        layers += incep3b.layers
        params += incep3b.params
        weight_types += incep3b.weight_types        
        print 'incep3b output shape: (28x28x480)'
        # output shape = (28x28x480)        

#        lrn3 = self.lrn_func(incep3b.output)
#        print 'LRN(added)'
        
        pool3 = Pool(input=incep3b.output,
                      poolsize=3, poolstride=2, poolpad=1, 
                      mode = 'max' )        
        # output shape = (14x14x480)
        
        incep4a = Incept(pool3.output, input_shape = (480,14,14,batch_size), 
                          n1x1=192, nr3x3=96, n3x3=208, 
                          nr5x5=16, n5x5=48, npj=64)
        
        layers += incep4a.layers
        params += incep4a.params
        weight_types += incep4a.weight_types        
        print 'incep4a output shape: (14x14x512)'
        # output shape = (14x14x512)
        
        incep4b = Incept(incep4a.output, input_shape = (512,14,14,batch_size), 
                          n1x1=160, nr3x3=112, n3x3=224, 
                          nr5x5=24, n5x5=64, npj=64)
        
        layers += incep4b.layers
        params += incep4b.params
        weight_types += incep4b.weight_types
        print 'incep4b output shape: (14x14x512)'        
        # output shape = (14x14x512)          
        

        incep4c = Incept(incep4b.output, input_shape = (512,14,14,batch_size), 
                          n1x1=128, nr3x3=128, n3x3=256, 
                          nr5x5=24, n5x5=64, npj=64)
        
        layers += incep4c.layers
        params += incep4c.params
        weight_types += incep4c.weight_types
        print 'incep4c output shape: (14x14x512)'        
        # output shape = (14x14x512) 

        incep4d = Incept(incep4c.output, input_shape = (512,14,14,batch_size), 
                          n1x1=112, nr3x3=144, n3x3=288, 
                          nr5x5=32, n5x5=64, npj=64)
        
        layers += incep4d.layers
        params += incep4d.params
        weight_types += incep4d.weight_types
        print 'incep4d output shape: (14x14x528)'        
        # output shape = (14x14x528) 
         
        
        incep4e = Incept(incep4d.output, input_shape = (528,14,14,batch_size), 
                          n1x1=256, nr3x3=160, n3x3=320, 
                          nr5x5=32, n5x5=128, npj=128)
        
        layers += incep4e.layers
        params += incep4e.params
        weight_types += incep4e.weight_types
        print 'incep4e output shape: (14x14x832)'       
        # output shape = (14x14x832)                
        
        lrn4 = self.lrn_func(incep4e.output)  # turn on only this for 16data, 53s/5120images
        print 'LRN(added)'
        
        pool4 = Pool(input=lrn4, #incep4e.output,
                      poolsize=3, poolstride=2, poolpad=1, 
                      mode = 'max' )        
        # output shape = (7x7x832)        
        
        incep5a = Incept(pool4.output, input_shape = (832,7,7,batch_size), 
                          n1x1=256, nr3x3=160, n3x3=320, 
                          nr5x5=32, n5x5=128, npj=128)
        
        layers += incep5a.layers
        params += incep5a.params
        weight_types += incep5a.weight_types
        print 'incep5a output shape: (7x7x832)'      
        # output shape = (7x7x832)   
        
        
        incep5b = Incept(incep5a.output, input_shape = (832,7,7,batch_size), 
                          n1x1=384, nr3x3=192, n3x3=384, 
                          nr5x5=48, n5x5=128, npj=128)
        
        layers += incep5b.layers
        params += incep5b.params
        weight_types += incep5b.weight_types
        print 'incep5b output shape: (7x7x1024)' 
        # output shape = (7x7x1024)
        
#        lrn5 = self.lrn_func(incep5b.output) # turn on only this for 16data, 51s/5120images
#        print 'LRN(added)'
        
        poolx = Pool(input=incep5b.output,
                      poolsize=7, poolstride=1, poolpad=0, 
                      mode = 'average' )
        # output shape = (1x1x1024)

           
        l_flatten = T.flatten(poolx.output.dimshuffle(3, 0, 1, 2), 2)
        # output shape = (1024)                              
    
        dropout= Dropout(input=l_flatten,n_in=1024, n_out=1024, prob_drop=0.4)
        # output shape = (1024)
               
        
        softmax_layer = Softmax(input=dropout.output ,n_in=1024, n_out=n_softmax_out)
        # output shape = (n_softmax_out)       
        
        layers.append(softmax_layer)
        params += softmax_layer.params
        weight_types += softmax_layer.weight_type        
        
        # auxilary classifier
        print 'auxilary classifier 1:'
        aux1 = aux_tower(input=incep4a.output,input_shape=(512,14,14,batch_size),config=config)
        
        layers += aux1.layers
        params += aux1.params
        weight_types += aux1.weight_types
        
        print 'auxilary classifier 2:'                               
        aux2 = aux_tower(input=incep4d.output,input_shape=(528,14,14,batch_size),config=config)
        
        layers += aux2.layers
        params += aux2.params
        weight_types += aux2.weight_types 

        self.layers = layers
        self.params = params
        self.weight_types = weight_types        
        self.output = softmax_layer.p_y_given_x
        self.cost = softmax_layer.negative_log_likelihood(y)+\
                0.3*aux1.negative_log_likelihood(y)+0.3*aux2.negative_log_likelihood(y)        
        self.error = softmax_layer.errors(y)
        self.error_top_5 = softmax_layer.errors_top_x(y)
        
        # training related
        self.base_lr = np.float32(config['learning_rate'])
        self.shared_lr = theano.shared(self.base_lr)
        self.mu = config['momentum'] # def: 0.9 # momentum
        self.eta = config['weight_decay'] #0.0002 # weight decay
        
        self.shared_x = theano.shared(np.zeros((3, config['input_width'], 
                                                  config['input_height'], 
                                                  config['file_batch_size']), 
                                                  dtype=theano.config.floatX),  
                                                  borrow=True)
                                              
        self.shared_y = theano.shared(np.zeros((config['file_batch_size'],), 
                                          dtype=int),   borrow=True)
                                          
        self.grads = T.grad(self.cost,self.params)
        
        subb_ind = T.iscalar('subb')  # sub batch index
        #print self.shared_x[:,:,:,subb_ind*self.batch_size:(subb_ind+1)*self.batch_size].shape.eval()
        self.subb_ind = subb_ind
        self.shared_x_slice = self.shared_x[:,:,:,subb_ind*self.batch_size:(subb_ind+1)*self.batch_size]
        self.shared_y_slice = self.shared_y[subb_ind*self.batch_size:(subb_ind+1)*self.batch_size]

Example 104

Project: fos
Source File: primitives.py
View license
def makeCylinder( p1, p2, r1, r2, resolution = 4 ):
    """ Cylinder or Cone generation

    Notes
    -----
    http://paulbourke.net/miscellaneous/sphere_cylinder/
    """
    # TODO: implement capping

    if r1 == 0.0 and r2 == 0.0:
        print("Not both radii can be zero!")
        return None

    # Normal pointing from p1 to p2
    n = (p1 - p2).astype( np.float32 )

    if n[0] != 0.0 or n[1] != 0.0:
        A = np.array( [n[1], -n[0], 0.0], dtype = np.float32 )
    elif n[1] != 0.0 or n[2] != 0.0:
        A = np.array( [0.0, n[2], -n[1]], dtype = np.float32 )
    elif n[0] != 0.0 or n[2] != 0.0:
        A = np.array( [n[0], 0.0, -n[1]], dtype = np.float32 )
    else:
        print("Cylinder coordinates well-posed. Cannot create perpendicular vector.")

    B = np.cross( n, A )

    # Normalize
    A = A / np.linalg.norm( A )
    B = B / np.linalg.norm( B )

    n = 0
    TWOPI =  2 * np.pi  
    p = []
    f = []
    for i in range(resolution):
        theta1 = i * TWOPI / (resolution)
        sint1 = np.sin(theta1)
        cost1 = np.cos(theta1)
        theta2 = (i+1) * TWOPI / (resolution)
        sint2 = np.sin(theta2)
        cost2 = np.cos(theta2)

        face = []
        # first point, q0
        p.append( [
            p1[0] + r1 * cost1 * A[0] + r1 * sint1 * B[0],
            p1[1] + r1 * cost1 * A[1] + r1 * sint1 * B[1],
            p1[2] + r1 * cost1 * A[2] + r1 * sint1 * B[2]
        ])
        face.append( n )
        n += 1

        # second point, q1
        p.append( [
            p2[0] + r2 * cost1 * A[0] + r2 * sint1 * B[0],
            p2[1] + r2 * cost1 * A[1] + r2 * sint1 * B[1],
            p2[2] + r2 * cost1 * A[2] + r2 * sint1 * B[2]
        ])
        face.append( n )
        n += 1

        if r2 != 0.0:
            p.append( [
                p2[0] + r2 * cost2 * A[0] + r2 * sint2 * B[0],
                p2[1] + r2 * cost2 * A[1] + r2 * sint2 * B[1],
                p2[2] + r2 * cost2 * A[2] + r2 * sint2 * B[2]
            ])
            face.append( n )
            n += 1

        if r1 != 0.0:
            p.append( [
                p1[0] + r1 * cost2 * A[0] + r1 * sint2 * B[0],
                p1[1] + r1 * cost2 * A[1] + r1 * sint2 * B[1],
                p1[2] + r1 * cost2 * A[2] + r1 * sint2 * B[2]
            ])
            face.append( n )
            n += 1

        # Create triangles
        if len(face) == 4:
            f.append( [face[0], face[1], face[2]] )
            f.append( [face[0], face[2], face[3]] )
        elif len(face) == 3:
            f.append( face )
        else:
            print("Wrong number of vertices in face.")

    return np.array( p, dtype = np.float32 ), np.array( f, dtype = np.uint32 )

Example 105

Project: lowrank-gru
Source File: RNN_adding.py
View license
def main(n_iter, n_batch, n_hidden, time_steps, learning_rate, savefile, scale_penalty, use_scale,
         model, n_hidden_lstm, loss_function, n_gru_lr_proj, initial_b_u):

    #import pdb; pdb.set_trace()
 
    # --- Set optimization params --------
    gradient_clipping = np.float32(50000)

    # --- Set data params ----------------
    n_input = 2
    n_output = 1
  

    # --- Manage data --------------------
    n_train = 1e5
    n_test = 1e4
    num_batches = n_train / n_batch
    
    train_x = np.asarray(np.zeros((time_steps, n_train, 2)),
                         dtype=theano.config.floatX)
    

    train_x[:,:,0] = np.asarray(np.random.uniform(low=0.,
                                                  high=1.,
                                                  size=(time_steps, n_train)),
                                dtype=theano.config.floatX)
    
#    inds = np.asarray([np.random.choice(time_steps, 2, replace=False) for i in xrange(train_x.shape[1])])    
    inds = np.asarray(np.random.randint(time_steps/2, size=(train_x.shape[1],2)))
    inds[:, 1] += time_steps/2  
    
    for i in range(train_x.shape[1]):
        train_x[inds[i, 0], i, 1] = 1.0
        train_x[inds[i, 1], i, 1] = 1.0
 
    train_y = (train_x[:,:,0] * train_x[:,:,1]).sum(axis=0)
    train_y = np.reshape(train_y, (n_train, 1))

    test_x = np.asarray(np.zeros((time_steps, n_test, 2)),
                        dtype=theano.config.floatX)
    

    test_x[:,:,0] = np.asarray(np.random.uniform(low=0.,
                                                 high=1.,
                                                 size=(time_steps, n_test)),
                                dtype=theano.config.floatX)
    
    inds = np.asarray([np.random.choice(time_steps, 2, replace=False) for i in xrange(test_x.shape[1])])    
    for i in range(test_x.shape[1]):
        test_x[inds[i, 0], i, 1] = 1.0
        test_x[inds[i, 1], i, 1] = 1.0
 
   
    test_y = (test_x[:,:,0] * test_x[:,:,1]).sum(axis=0)
    test_y = np.reshape(test_y, (n_test, 1)) 


   #######################################################################

    gradient_clipping = np.float32(1)

    if (model == 'LSTM'):   
        #inputs, parameters, costs = LSTM(n_input, n_hidden_lstm, n_output, loss_function=loss_function)
        inputs, parameters, costs = LSTM(n_input, n_hidden_lstm, n_output, loss_function=loss_function, initial_b_f=initial_b_u)
        gradients = T.grad(costs[0], parameters)
        gradients = [T.clip(g, -gradient_clipping, gradient_clipping) for g in gradients]
    
    #by AnvaMiba
    elif (model == 'GRU'):   
        inputs, parameters, costs = GRU(n_input, n_hidden_lstm, n_output, loss_function=loss_function, initial_b_u=initial_b_u)
        gradients = T.grad(costs[0], parameters)
        gradients = [T.clip(g, -gradient_clipping, gradient_clipping) for g in gradients]
    #by AnvaMiba
    elif (model == 'GRU_LR'):   
        inputs, parameters, costs = GRU_LR(n_input, n_hidden_lstm, n_output, n_gru_lr_proj, loss_function=loss_function, initial_b_u=initial_b_u)
        gradients = T.grad(costs[0], parameters)
        gradients = [T.clip(g, -gradient_clipping, gradient_clipping) for g in gradients]

    elif (model == 'complex_RNN'):
        inputs, parameters, costs = complex_RNN(n_input, n_hidden, n_output, scale_penalty, loss_function=loss_function)
        if use_scale is False:
            parameters.pop()
        gradients = T.grad(costs[0], parameters)

    elif (model == 'complex_RNN_LSTM'):
        inputs, parameters, costs = complex_RNN_LSTM(n_input, n_hidden, n_hidden_lstm, n_output, scale_penalty, loss_function=loss_function)

    elif (model == 'IRNN'):
        inputs, parameters, costs = IRNN(n_input, n_hidden, n_output, loss_function=loss_function)
        gradients = T.grad(costs[0], parameters)
        gradients = [T.clip(g, -gradient_clipping, gradient_clipping) for g in gradients]

    elif (model == 'RNN'):
        inputs, parameters, costs = tanhRNN(n_input, n_hidden, n_output, loss_function=loss_function)
        gradients = T.grad(costs[0], parameters)
        gradients = [T.clip(g, -gradient_clipping, gradient_clipping) for g in gradients]

    else:
        print >> sys.stderr, "Unsuported model:", model
        return
 

   




    s_train_x = theano.shared(train_x)
    s_train_y = theano.shared(train_y)

    s_test_x = theano.shared(test_x)
    s_test_y = theano.shared(test_y)


    # --- Compile theano functions --------------------------------------------------

    index = T.iscalar('i')

    updates, rmsprop = rms_prop(learning_rate, parameters, gradients)

    givens = {inputs[0] : s_train_x[:, n_batch * index : n_batch * (index + 1), :],
              inputs[1] : s_train_y[n_batch * index : n_batch * (index + 1), :]}

    givens_test = {inputs[0] : s_test_x,
                   inputs[1] : s_test_y}
    
   
    
    train = theano.function([index], costs[0], givens=givens, updates=updates)
    test = theano.function([], costs[1], givens=givens_test)

    # --- Training Loop ---------------------------------------------------------------

    # f1 = file('/data/lisatmp3/shahamar/adding/complexRNN_400.pkl', 'rb')
    # data1 = cPickle.load(f1)
    # f1.close()
    # train_loss = data1['train_loss']
    # test_loss = data1['test_loss']
    # best_params = data1['best_params']
    # best_test_loss = data1['best_test_loss']

    # for i in xrange(len(parameters)):
    #     parameters[i].set_value(data1['parameters'][i])

    # for i in xrange(len(parameters)):
    #     rmsprop[i].set_value(data1['rmsprop'][i])
    
#    import pdb; pdb.set_trace()

    train_loss = []
    test_loss = []
    best_params = [p.get_value() for p in parameters]
    best_test_loss = 1e6
    for i in xrange(n_iter):
#        start_time = timeit.default_timer()


        if (n_iter % int(num_batches) == 0):
            #import pdb; pdb.set_trace()
            inds = np.random.permutation(int(n_train))
            data_x = s_train_x.get_value()
            s_train_x.set_value(data_x[:,inds,:])
            data_y = s_train_y.get_value()
            s_train_y.set_value(data_y[inds,:])


        mse = train(i % int(num_batches))
        train_loss.append(mse)
        print >> sys.stderr, "Iteration:", i
        print >> sys.stderr, "mse:", mse
        print >> sys.stderr, ''

        #if (i % 50==0):
        if (i % 500==0):
            mse = test()
            print >> sys.stderr, ''
            print >> sys.stderr, "TEST"
            print >> sys.stderr, "mse:", mse
            print >> sys.stderr, '' 
            test_loss.append(mse)

            if mse < best_test_loss:
                best_params = [p.get_value() for p in parameters]
                best_test_loss = mse

            
            save_vals = {'parameters': [p.get_value() for p in parameters],
                         'rmsprop': [r.get_value() for r in rmsprop],
                         'train_loss': train_loss,
                         'test_loss': test_loss,
                         'best_params': best_params,
                         'best_test_loss': best_test_loss,
                         'model': model,
                         'time_steps': time_steps}

            cPickle.dump(save_vals,
                         file(savefile, 'wb'),
                         cPickle.HIGHEST_PROTOCOL)

Example 106

Project: fos-legacy
Source File: network.py
View license
    def __init__(self, node_position, affine = None, force_centering = True, *args, **kwargs):
        """ Draw a network
                
        affine : (4,4)
            The affine has a translational and rotational part
        force_centering : bool
            Subtract the mean over node_position from all node positions to center the actor.
        scale_factor : float
            Scaling the actor, scales all the vertices. You might want to multiply
            your node_position array.
        global_node_size : float
            Size for all nodes, used when node_size is None. If None, defaults to 1.0
        global_edge_width : float
            The global line width. Defaults to 1.0
        
        Node related
        ------------
        node_position : (N,3)
            Node positions as ndarray

        node_size
            The size of the node

        node_label : dictionary of dictionary
            Node labels. Keys refer to the node index (starting with 0)
            
        node_shape
            cube, sphere, pyramid, electrodes (cylinders)
            
        node_color : (N,4)
            The color of the nodes and its alpha value
            Either given [0,1] or [0,255]
            (or: cmap, vmin, vmax)
        
        node_show_labels
            Show all labels on the nodes / 
            only for specified nodes.
            node_label has to be set
        
        Edge related
        ------------
        
        edge_directed : bool
            Interpret `edge_weight` as directed
        
        edge_weight : (M,1)
            The weight determines the width of the line
            
        edge_color : (N,4)
            The color of the edges
            (or cmap, vmin, vmax)
            
        edge_style
            solid, dashed, dotted, dashdot
            What does OpenGL support natively?
            
        edge_label
            The label for the edges
        
        edge_width_granularity
            Idea: Subdivide the weight histogram into different
            bins with their own line width
        
        Font related (global or per node/edge?)
        ------------
        
        font_size: int
           Font size for text labels (default=12)
    
        font_color: string
           Font color string (default='k' black)
    
        font_weight: string
           Font weight (default='normal')
    
        font_family: string
           Font family (default='sans-serif')
           
        """
        
        super(AttributeNetwork, self).__init__()        
        if not node_position is None:
            self.vertices = node_position
        else:
            raise Exception("You have to specify the node_position array")        
        self.node_glprimitive = NodePrimitive()
        self.edge_glprimitive = EdgePrimitive()        
        if affine == None:
            # create a default affine
            self.affine = np.eye(4, dtype = np.float32)
        else:
            self.affine = affine        
        self._update_glaffine()        
        self.edge_color = None
        if kwargs.has_key('edge_connectivity'):
            self.edge_connectivity = kwargs['edge_connectivity']
            nr_edges = self.edge_connectivity.shape[0]            
            if kwargs.has_key('edge_color'):
                self.edge_color = kwargs['edge_color']
            else:
                # default edge color
                self.edge_color = np.array( [[255,255,255,255]], dtype = np.ubyte).repeat(nr_edges,axis=0)
        else:
            self.edge_connectivity = None
        if kwargs.has_key('node_label'):
            self.node_label = kwargs['node_label']
        else:
            self.node_label = None
        if kwargs.has_key('edge_weight'):
            self.edge_weight = kwargs['edge_weight']        
        if kwargs.has_key('aabb'):
            self.aabb = kwargs['aabb']
        else:
            self.aabb = None            
        if kwargs.has_key('obb'):
            self.obb = kwargs['obb']
        else:
            self.obb = None            
        if kwargs.has_key('scale_factor'):
            self.scale_factor = kwargs['scale_factor']
        else:
            self.scale_factor = 1.0
        if kwargs.has_key('global_node_size'):
            self.global_node_size = kwargs['global_node_size']
        else:
            self.global_node_size = 1.0            
        if kwargs.has_key('global_edge_width'):
            self.global_edge_width = kwargs['global_edge_width']
        else:
            self.global_edge_width = 1.0            
        if not self.vertices is None:
            nr_nodes = self.vertices.shape[0]            
            if force_centering:
                self.vertices = self.vertices - np.mean(self.vertices, axis = 0)
            if kwargs.has_key('node_size'):
                self.node_size = kwargs['node_size'].ravel()
            else:
                # default size 1.0
                self.node_size = np.ones( (nr_nodes, 1), dtype = np.float32 ).ravel() * self.global_node_size
            if kwargs.has_key('node_color'):
                self.node_color = kwargs['node_color']
            else:
                # create default colors for nodes
                self.node_color = np.array( [[200,200,200,255]], dtype = np.ubyte).repeat(nr_nodes,axis=0)
            
        # default variables
        self.internal_timestamp = 0.0        
        # network creation
        ##################
        if not self.vertices is None and not self.node_size is None:
            assert self.vertices.shape[0] == self.node_size.size            
            self.node_glprimitive._make_cubes(self.vertices, self.node_size)                
            self.scale(self.scale_factor)                        
        if not self.node_color is None:
            self.node_glprimitive._make_color(self.node_color)            
        if not self.edge_connectivity is None:            
            self.edge_glprimitive._make_edges(self.vertices, self.edge_connectivity)            
        if not self.edge_color is None:
            self.edge_glprimitive._make_color(self.edge_color)        
        self.living = False        
        # create aabb (in actor)
        self.make_aabb(margin = self.node_size.max())
        if not self.node_label is None:            
            textcenter = np.zeros( (len(self.node_label), 3), dtype = np.float32)
            texts = []
            fonts = []
            font_sizes = []
            colors = []            
            for i,k in enumerate(self.node_label.keys()):
                assert k < len(self.vertices)
                if not self.node_label[k].has_key('label'):
                    print("node_label dictionary for key %s has no label key. skip.")
                    continue                
                textcenter[i,:] = self.vertices[k, :] # can shift by nodesize
                texts.append(self.node_label[k]['label'])                
                if self.node_label[k].has_key('font'):
                    fonts.append(self.node_label[k]['font'])
                else:
                    fonts.append('Arial')                    
                if self.node_label[k].has_key('size'):
                    font_sizes.append(self.node_label[k]['size'])
                else:
                    font_sizes.append(20)                    
                if self.node_label[k].has_key('color'):
                    colors.append(self.node_label[k]['color'])
                else:
                    colors.append( (0,255,0,255) )                                

Example 107

Project: Theano-MPI
Source File: vggnet_16.py
View license
    def build_model(self):
        
        if self.verbose: print 'VGGNet_16 3/20'
        
        self.name = 'vggnet'
        
        # input shape in c01b 
        self.channels = 3 # 'c' mean(R,G,B) = (103.939, 116.779, 123.68)
        self.input_width = self.config['input_width'] # '0' single scale training 224
        self.input_height = self.config['input_height'] # '1' single scale training 224
        self.batch_size = self.config['batch_size'] # 'b'
        b = self.batch_size
        
        # output dimension
        self.n_softmax_out = self.config['n_softmax_out']
        
        # start graph construction from scratch
        self.x = T.ftensor4('x')
        
        self.y = T.lvector('y')
        
        x_shuffled = self.x.dimshuffle(3, 0, 1, 2)  # c01b to bc01
        
        layers = []
        params = []
        weight_types = [] # for distinguishing w and b later      
        
        # bc01 from now on
        
        conv_3x3 = Conv(input=x_shuffled,
                        input_shape=(b,
                                    self.channels,
                                    self.input_width,
                                    self.input_height), # (b, 3, 224, 224)
                        convstride=1,
                        padsize=1,
                        W = Normal((64, self.channels, 3, 3), std=0.3), # bc01
                        b = Constant((64,), val=0.2),
                        printinfo=self.verbose
                        #output_shape = (b, 64, 224, 224)
                        )
                        
        conv_3x3 = Conv(input=conv_3x3,
                        #input_shape=pool_2x2.output_shape, # (b, 64, 112, 112)
                        convstride=1,
                        padsize=1,
                        W = Normal((64, conv_3x3.output_shape[1], 3, 3), std=0.1), # bc01
                        b = Constant((64,), val=0.1),
                        printinfo=self.verbose
                        #output_shape = (b, 128, 112, 112)
                        )

        pool_2x2 = Pool(input=conv_3x3, 
                        #input_shape=conv_3x3.output_shape, # (b, 64, 224, 224)
                        poolsize=2, 
                        poolstride=2, 
                        poolpad=0,
                        mode = 'max',
                        printinfo=self.verbose
                        #output_shape = (b, 64, 112, 112)
                        )

        conv_3x3 = Conv(input=pool_2x2,
                        #input_shape=pool_2x2.output_shape, # (b, 64, 112, 112)
                        convstride=1,
                        padsize=1,
                        W = Normal((128, pool_2x2.output_shape[1], 3, 3), std=0.1), # bc01
                        b = Constant((128,), val=0.02),
                        printinfo=self.verbose
                        #output_shape = (b, 128, 112, 112)
                        )
                        
        conv_3x3 = Conv(input=conv_3x3,
                        #input_shape=pool_2x2.output_shape, # (b, 64, 112, 112)
                        convstride=1,
                        padsize=1,
                        W = Normal((128, conv_3x3.output_shape[1], 3, 3), std=0.1), # bc01
                        b = Constant((128,), val=0.02),
                        printinfo=self.verbose
                        #output_shape = (b, 128, 112, 112)
                        )

        pool_2x2 = Pool(input=conv_3x3, 
                        #input_shape=conv_3x3.output_shape, # (b, 128, 112, 112)
                        poolsize=2, 
                        poolstride=2, 
                        poolpad=0,
                        mode = 'max',
                        printinfo=self.verbose
                        #output_shape = (b, 128, 56, 56)
                        )

        conv_3x3 = Conv(input=pool_2x2,
                        #input_shape=pool_2x2.output_shape, # (b, 128, 56, 56)
                        convstride=1,
                        padsize=1,
                        W = Normal((256, pool_2x2.output_shape[1], 3, 3), std=0.05), # bc01
                        b = Constant((256,), val=0.02),
                        printinfo=self.verbose
                        #output_shape = (b, 256, 56, 56)
                        )
        conv_3x3 = Conv(input=conv_3x3,
                        #input_shape=conv_3x3.output_shape, # (b, 256, 56, 56)
                        convstride=1,
                        padsize=1,
                        W = Normal((256, conv_3x3.output_shape[1], 3, 3), std=0.05), # bc01
                        b = Constant((256,), val=0.01),
                        printinfo=self.verbose
                        #output_shape = (b, 256, 56, 56)
                        )
        conv_3x3 = Conv(input=conv_3x3,
                        #input_shape=conv_3x3.output_shape, # (b, 256, 56, 56)
                        convstride=1,
                        padsize=1,
                        W = Normal((256, conv_3x3.output_shape[1], 3, 3), std=0.05), # bc01
                        b = Constant((256,), val=0.01),
                        printinfo=self.verbose
                        #output_shape = (b, 256, 56, 56)
                        )

        pool_2x2 = Pool(input=conv_3x3, 
                        #input_shape=conv_3x3.output_shape, # (b, 256, 56, 56)
                        poolsize=2, 
                        poolstride=2, 
                        poolpad=0,
                        mode = 'max',
                        printinfo=self.verbose
                        #output_shape = (b, 256, 28, 28)
                        )

        conv_3x3 = Conv(input=pool_2x2,
                        #input_shape=pool_2x2.output_shape, # (b, 256, 28, 28)
                        convstride=1,
                        padsize=1,
                        W = Normal((512, pool_2x2.output_shape[1], 3, 3), std=0.05), # bc01
                        b = Constant((512,), val=0.02),
                        printinfo=self.verbose
                        #output_shape = (b, 512, 28, 28)
                        )
        conv_3x3 = Conv(input=conv_3x3,
                        #input_shape=conv_3x3.output_shape, # (b, 512, 28, 28)
                        convstride=1,
                        padsize=1,
                        W = Normal((512, conv_3x3.output_shape[1], 3, 3), std=0.01), # bc01
                        b = Constant((512,), val=0.01),
                        printinfo=self.verbose
                        #output_shape = (b, 512, 28, 28)
                        )
        conv_3x3 = Conv(input=conv_3x3,
                        #input_shape=conv_3x3.output_shape, # (b, 512, 28, 28)
                        convstride=1,
                        padsize=1,
                        W = Normal((512, conv_3x3.output_shape[1], 3, 3), std=0.01), # bc01
                        b = Constant((512,), val=0.01),
                        printinfo=self.verbose
                        #output_shape = (b, 512, 28, 28)
                        )
        

        pool_2x2 = Pool(input=conv_3x3, 
                        #input_shape=conv_3x3.output_shape, # (b, 512, 28, 28)
                        poolsize=2, 
                        poolstride=2, 
                        poolpad=0,
                        mode = 'max',
                        printinfo=self.verbose
                        #output_shape = (b, 512, 14, 14)
                        )

        conv_3x3 = Conv(input=pool_2x2,
                        #input_shape=pool_2x2.output_shape, # (b, 512, 14, 14)
                        convstride=1,
                        padsize=1,
                        W = Normal((512, pool_2x2.output_shape[1], 3, 3), std=0.005), # bc01
                        b = Constant((512,)),
                        printinfo=self.verbose
                        #output_shape = (b, 512, 14, 14)
                        )
        conv_3x3 = Conv(input=conv_3x3,
                        #input_shape=conv_3x3.output_shape, # (b, 512, 14, 14)
                        convstride=1,
                        padsize=1,
                        W = Normal((512, conv_3x3.output_shape[1], 3, 3), std=0.005), # bc01
                        b = Constant((512,)),
                        printinfo=self.verbose
                        #output_shape = (b, 512, 14, 14)
                        )
        conv_3x3 = Conv(input=conv_3x3,
                        #input_shape=conv_3x3.output_shape, # (b, 512, 14, 14)
                        convstride=1,
                        padsize=1,
                        W = Normal((512, conv_3x3.output_shape[1], 3, 3), std=0.005), # bc01
                        b = Constant((512,)),
                        printinfo=self.verbose
                        #output_shape = (b, 512, 14, 14)
                        )
 
        pool_2x2 = Pool(input=conv_3x3, 
                        #input_shape=conv_3x3.output_shape, # (b, 512, 14, 14)
                        poolsize=2, 
                        poolstride=2, 
                        poolpad=0,
                        mode = 'max',
                        printinfo=self.verbose
                        #output_shape = (b, 512, 7, 7)
                        )

        flatten = Flatten(input = pool_2x2, #5
                        #input_shape = pool_2x2.output_shape, # (b, 512, 7, 7)
                        axis = 2, # expand dimensions after the first dimension
                        printinfo=self.verbose
                        #output_shape = (b, 25088)
                        )
        fc_4096 = FC(input= flatten, 
                        n_out=4096,
                        W = Normal((flatten.output_shape[1], 4096), std=0.001),
                        b = Constant((4096,),val=0.01),
                        printinfo=self.verbose
                        #input_shape = flatten.output_shape # (b, 25088)
                        )
        dropout= Dropout(input=fc_4096,
                        n_out=fc_4096.output_shape[1], 
                        prob_drop=0.5,
                        printinfo=self.verbose
                        #input_shape = fc_4096.output_shape # (b, 4096)
                        )
        fc_4096 = FC(input= dropout,  
                        n_out=4096,
                        W = Normal((dropout.output_shape[1], 4096), std=0.005),
                        b = Constant((4096,),val=0.01),
                        printinfo=self.verbose
                        #input_shape = dropout.output_shape # (b, 4096)
                        )
        dropout= Dropout(input=fc_4096, 
                        n_out=fc_4096.output_shape[1], 
                        prob_drop=0.5,
                        printinfo=self.verbose
                        #input_shape = fc_4096.output_shape # (b, 4096)
                        )
        softmax = Softmax(input=dropout,  
                        n_out=self.n_softmax_out,
                        W = Normal((dropout.output_shape[1], self.n_softmax_out), std=0.005),
                        b = Constant((self.n_softmax_out,),val=0),
                        printinfo=self.verbose
                        #input_shape = dropout.output_shape # (b, 4096)
                        )
        
        self.output_layer = softmax
        
        self.output = self.output_layer.output
        
        self.layers = get_layers(lastlayer = self.output_layer)
        
        self.layers = [layer for layer in self.layers \
            if layer.name not in ['LRN\t','Pool\t','Flatten\t','Dropout'+ str(0.5)]]
        
        self.params,self.weight_types = get_params(self.layers)
        
        # training related
        self.base_lr = np.float32(self.config['learning_rate'])
        self.shared_lr = theano.shared(self.base_lr)
        self.step_idx = 0
        self.mu = self.config['momentum'] # def: 0.9 # momentum
        self.eta = self.config['weight_decay'] #0.0002 # weight decay
        
        self.shared_x = theano.shared(np.zeros((
                                                3,
                                                self.input_width, 
                                                self.input_height,
                                                self.config['file_batch_size']
                                                ), 
                                                dtype=theano.config.floatX),  
                                                borrow=True)
                                              
        self.shared_y = theano.shared(np.zeros((self.config['file_batch_size'],), 
                                          dtype=int),   borrow=True)
        
        
        
        self.grads = T.grad(self.cost,self.params)
        
        subb_ind = T.iscalar('subb')  # sub batch index
        #print self.shared_x[:,:,:,subb_ind*self.batch_size:(subb_ind+1)*self.batch_size].shape.eval()
        self.subb_ind = subb_ind
        self.shared_x_slice = self.shared_x[:,:,:,subb_ind*self.batch_size:(subb_ind+1)*self.batch_size]
        self.shared_y_slice = self.shared_y[subb_ind*self.batch_size:(subb_ind+1)*self.batch_size]

Example 108

View license
def main_test_layers(model='relu'):
    X_train, y_train, X_val, y_val, X_test, y_test = \
                                    tl.files.load_mnist_dataset(shape=(-1,784))

    X_train = np.asarray(X_train, dtype=np.float32)
    y_train = np.asarray(y_train, dtype=np.int32)
    X_val = np.asarray(X_val, dtype=np.float32)
    y_val = np.asarray(y_val, dtype=np.int32)
    X_test = np.asarray(X_test, dtype=np.float32)
    y_test = np.asarray(y_test, dtype=np.int32)

    print('X_train.shape', X_train.shape)
    print('y_train.shape', y_train.shape)
    print('X_val.shape', X_val.shape)
    print('y_val.shape', y_val.shape)
    print('X_test.shape', X_test.shape)
    print('y_test.shape', y_test.shape)
    print('X %s   y %s' % (X_test.dtype, y_test.dtype))

    sess = tf.InteractiveSession()

    # placeholder
    x = tf.placeholder(tf.float32, shape=[None, 784], name='x')
    y_ = tf.placeholder(tf.int32, shape=[None, ], name='y_')

    # Note: the softmax is implemented internally in tl.cost.cross_entropy(y, y_)
    # to speed up computation, so we use identity in the last layer.
    # see tf.nn.sparse_softmax_cross_entropy_with_logits()
    if model == 'relu':
        network = tl.layers.InputLayer(x, name='input_layer')
        network = tl.layers.DropoutLayer(network, keep=0.8, name='drop1')
        network = tl.layers.DenseLayer(network, n_units=800,
                                        act = tf.nn.relu, name='relu1')
        network = tl.layers.DropoutLayer(network, keep=0.5, name='drop2')
        network = tl.layers.DenseLayer(network, n_units=800,
                                        act = tf.nn.relu, name='relu2')
        network = tl.layers.DropoutLayer(network, keep=0.5, name='drop3')
        network = tl.layers.DenseLayer(network, n_units=10,
                                        act = tf.identity,
                                        name='output_layer')
    elif model == 'dropconnect':
        network = tl.layers.InputLayer(x, name='input_layer')
        network = tl.layers.DropconnectDenseLayer(network, keep = 0.8,
                                                n_units=800, act = tf.nn.relu,
                                                name='dropconnect_relu1')
        network = tl.layers.DropconnectDenseLayer(network, keep = 0.5,
                                                n_units=800, act = tf.nn.relu,
                                                name='dropconnect_relu2')
        network = tl.layers.DropconnectDenseLayer(network, keep = 0.5,
                                                n_units=10,
                                                act = tf.identity,
                                                name='output_layer')

    # To print all attributes of a Layer.
    # attrs = vars(network)
    # print(', '.join("%s: %s\n" % item for item in attrs.items()))
    #
    # print(network.all_drop)     # {'drop1': 0.8, 'drop2': 0.5, 'drop3': 0.5}
    # print(drop1, drop2, drop3)  # Tensor("Placeholder_2:0", dtype=float32) Tensor("Placeholder_3:0", dtype=float32) Tensor("Placeholder_4:0", dtype=float32)
    # exit()

    y = network.outputs
    y_op = tf.argmax(tf.nn.softmax(y), 1)
    cost = tf.reduce_mean(tf.nn.sparse_softmax_cross_entropy_with_logits(y, y_))
    # Alternatively, you can use TensorLayer's function to compute cost:
    # cost = tl.cost.cross_entropy(y, y_)

    # You can add more penalty to the cost function as follow.
    # cost = cost + tl.cost.maxnorm_regularizer(1.0)(network.all_params[0]) + tl.cost.maxnorm_regularizer(1.0)(network.all_params[2])
    # cost = cost + tl.cost.lo_regularizer(0.0001)(network.all_params[0]) + tl.cost.lo_regularizer(0.0001)(network.all_params[2])
    # cost = cost + tl.cost.maxnorm_o_regularizer(0.001)(network.all_params[0]) + tl.cost.maxnorm_o_regularizer(0.001)(network.all_params[2])

    params = network.all_params
    # train
    n_epoch = 1
    batch_size = 128
    learning_rate = 0.0001
    print_freq = 10
    train_op = tf.train.AdamOptimizer(learning_rate, beta1=0.9, beta2=0.999,
                                epsilon=1e-08, use_locking=False).minimize(cost)

    sess.run(tf.initialize_all_variables()) # initialize all variables

    network.print_params()
    network.print_layers()

    print('   learning_rate: %f' % learning_rate)
    print('   batch_size: %d' % batch_size)

    for epoch in range(n_epoch):
        start_time = time.time()
        for X_train_a, y_train_a in tl.iterate.minibatches(X_train, y_train,
                                                    batch_size, shuffle=True):
            feed_dict = {x: X_train_a, y_: y_train_a}
            feed_dict.update( network.all_drop )    # enable dropout or dropconnect layers
            sess.run(train_op, feed_dict=feed_dict)

            # The optional feed_dict argument allows the caller to override the value of tensors in the graph. Each key in feed_dict can be one of the following types:
            # If the key is a Tensor, the value may be a Python scalar, string, list, or numpy ndarray that can be converted to the same dtype as that tensor. Additionally, if the key is a placeholder, the shape of the value will be checked for compatibility with the placeholder.
            # If the key is a SparseTensor, the value should be a SparseTensorValue.

        if epoch + 1 == 1 or (epoch + 1) % print_freq == 0:
            print("Epoch %d of %d took %fs" % (epoch + 1, n_epoch, time.time() - start_time))
            dp_dict = tl.utils.dict_to_one( network.all_drop ) # disable all dropout/dropconnect/denoising layers
            feed_dict = {x: X_train, y_: y_train}
            feed_dict.update(dp_dict)
            print("   train loss: %f" % sess.run(cost, feed_dict=feed_dict))
            dp_dict = tl.utils.dict_to_one( network.all_drop )
            feed_dict = {x: X_val, y_: y_val}
            feed_dict.update(dp_dict)
            print("   val loss: %f" % sess.run(cost, feed_dict=feed_dict))
            print("   val acc: %f" % np.mean(y_val ==
                                    sess.run(y_op, feed_dict=feed_dict)))
            try:
                # You can visualize the weight of 1st hidden layer as follow.
                tl.visualize.W(network.all_params[0].eval(), second=10,
                                        saveable=True, shape=[28, 28],
                                        name='w1_'+str(epoch+1), fig_idx=2012)
                # You can also save the weight of 1st hidden layer to .npz file.
                # tl.files.save_npz([network.all_params[0]] , name='w1'+str(epoch+1)+'.npz')
            except:
                raise Exception("You should change visualize_W(), if you want \
                            to save the feature images for different dataset")

    print('Evaluation')
    dp_dict = tl.utils.dict_to_one( network.all_drop )
    feed_dict = {x: X_test, y_: y_test}
    feed_dict.update(dp_dict)
    print("   test loss: %f" % sess.run(cost, feed_dict=feed_dict))
    print("   test acc: %f" % np.mean(y_test == sess.run(y_op,
                                                        feed_dict=feed_dict)))

    # Add ops to save and restore all the variables, including variables for training.
    # ref: https://www.tensorflow.org/versions/r0.8/how_tos/variables/index.html
    saver = tf.train.Saver()
    save_path = saver.save(sess, "model.ckpt")
    print("Model saved in file: %s" % save_path)


    # You can also save the parameters into .npz file.
    tl.files.save_npz(network.all_params , name='model.npz')
    # You can only save one parameter as follow.
    # tl.files.save_npz([network.all_params[0]] , name='model.npz')
    # Then, restore the parameters as follow.
    # load_params = tl.files.load_npz(path='', name='model.npz')
    # tl.files.assign_params(sess, load_params, network)

    # In the end, close TensorFlow session.
    sess.close()

Example 109

Project: fos-legacy
Source File: tracks.py
View license
    def __init__(self,fname,ang_table=None,colormap=None, line_width=3., shrink=None,subset=None,data_ext=None):

        self.position = (0,0,0)

        self.fname = fname
        
        self.manycolors = True
        
        self.bbox = None

        self.list_index = None

        self.affine = None

        self.data = None

        self.list_index = None

        self.rot_angle = 0

        self.colormap = None
                
        self.min = None
         
        self.max = None

        self.mean = None

        self.material_color = False

        self.fadeout = False

        self.fadein = False

        self.fadeout_speed = 0.

        self.fadein_speed = 0.

        self.min_length = 20.

        self.data_ext = data_ext

        self.angle = 0.

        self.angular_speed = .5

        self.line_width = line_width

        self.opacity = 1.

        self.near_pick = None

        self.far_pick = None

        self.near_pick_prev = None

        self.far_pick_prev = None

        self.picked_track = None

        self.pick_color = [1,1,0]

        self.brain_color = [1,1,1]

        self.yellow_indices = None

        self.dummy_data = False

        if subset != None:

            self.data_subset = subset #[0,20000]#None

        else:

            self.data_subset = None

        self.orbit_demo = False          

        self.orbit_anglez = 0.

        self.orbit_anglez_rate = 0.
        

        self.orbit_anglex = 0.

        self.orbit_anglex_rate = 0.


        self.orbit_angley = 0.

        self.orbit_angley_rate = 0.


        self.angle_table = ang_table

        
        self.angle_table_index = 0



        

        self.shrink = shrink

        self.picking_example = False

        if self.data_ext!=None:

            self.data=self.data_ext

        else:

            import dipy.io.trackvis as tv

            lines,hdr = tv.read(self.fname)

            ras = tv.aff_from_hdr(hdr)

            self.affine=ras

            tracks = [l[0] for l in lines]

            if self.yellow_indices != None :

                tracks = [t for t in tracks if tm.length(t) > 20]

            print 'tracks loaded'

            #self.data = [100*np.array([[0,0,0],[1,0,0],[2,0,0]]).astype(np.float32) ,100*np.array([[0,1,0],[0,2,0],[0,3,0]]).astype(np.float32)]#tracks[:20000]

            if self.dummy_data:

                self.data = [100*np.array([[0,0,0],[1,0,0],[2,0,0]]).astype(np.float32) ,100*np.array([[0,1,0],[0,2,0],[0,3,0]]).astype(np.float32)]

            if self.data_subset!=None:

                self.data = tracks[self.data_subset[0]:self.data_subset[1]]

            else:

                self.data = tracks



            data_stats = np.concatenate(tracks)

            self.min=np.min(data_stats,axis=0)

            self.max=np.max(data_stats,axis=0)

            self.mean=np.mean(data_stats,axis=0)

            if self.shrink != None:

                self.data = [ self.shrink*t  for t in self.data]

            del data_stats

            del lines

Example 110

Project: theano_alexnet
Source File: alex_net.py
View license
def compile_models(model, config, flag_top_5=False):

    x = model.x
    y = model.y
    rand = model.rand
    weight_types = model.weight_types

    cost = model.cost
    params = model.params
    errors = model.errors
    errors_top_5 = model.errors_top_5
    batch_size = model.batch_size

    mu = config['momentum']
    eta = config['weight_decay']

    # create a list of gradients for all model parameters
    grads = T.grad(cost, params)
    updates = []

    learning_rate = theano.shared(np.float32(config['learning_rate']))
    lr = T.scalar('lr')  # symbolic learning rate

    if config['use_data_layer']:
        raw_size = 256
    else:
        raw_size = 227

    shared_x = theano.shared(np.zeros((3, raw_size, raw_size,
                                       batch_size),
                                      dtype=theano.config.floatX),
                             borrow=True)
    shared_y = theano.shared(np.zeros((batch_size,), dtype=int),
                             borrow=True)

    rand_arr = theano.shared(np.zeros(3, dtype=theano.config.floatX),
                             borrow=True)

    vels = [theano.shared(param_i.get_value() * 0.)
            for param_i in params]

    if config['use_momentum']:

        assert len(weight_types) == len(params)

        for param_i, grad_i, vel_i, weight_type in \
                zip(params, grads, vels, weight_types):

            if weight_type == 'W':
                real_grad = grad_i + eta * param_i
                real_lr = lr
            elif weight_type == 'b':
                real_grad = grad_i
                real_lr = 2. * lr
            else:
                raise TypeError("Weight Type Error")

            if config['use_nesterov_momentum']:
                vel_i_next = mu ** 2 * vel_i - (1 + mu) * real_lr * real_grad
            else:
                vel_i_next = mu * vel_i - real_lr * real_grad

            updates.append((vel_i, vel_i_next))
            updates.append((param_i, param_i + vel_i_next))

    else:
        for param_i, grad_i, weight_type in zip(params, grads, weight_types):
            if weight_type == 'W':
                updates.append((param_i,
                                param_i - lr * grad_i - eta * lr * param_i))
            elif weight_type == 'b':
                updates.append((param_i, param_i - 2 * lr * grad_i))
            else:
                raise TypeError("Weight Type Error")

    # Define Theano Functions

    train_model = theano.function([], cost, updates=updates,
                                  givens=[(x, shared_x), (y, shared_y),
                                          (lr, learning_rate),
                                          (rand, rand_arr)])

    validate_outputs = [cost, errors]
    if flag_top_5:
        validate_outputs.append(errors_top_5)

    validate_model = theano.function([], validate_outputs,
                                     givens=[(x, shared_x), (y, shared_y),
                                             (rand, rand_arr)])

    train_error = theano.function(
        [], errors, givens=[(x, shared_x), (y, shared_y), (rand, rand_arr)])

    return (train_model, validate_model, train_error,
            learning_rate, shared_x, shared_y, rand_arr, vels)

Example 111

Project: fos-legacy
Source File: tracks.py
View license
    def __init__(self,fname,ang_table=None,colormap=None, line_width=3., shrink=None,subset=None,tracks=None,text=None,text_color=np.array([1,0,0])):

        self.position = (0,0,0)

        self.fname = fname
        
        self.manycolors = True
        
        self.bbox = None

        self.list_index = None

        self.affine = None

        self.data = None

        self.list_index = None

        self.rot_angle = 0

        self.colormap = None

        self.text = text
        
        self.min = None
         
        self.max = None

        self.mean = None

        self.material_color = False

        self.fadeout = False

        self.fadein = False

        self.fadeout_speed = 0.

        self.fadein_speed = 0.

        self.min_length = 20.

        self.angle = 0.

        self.angular_speed = .5

        self.line_width = line_width

        self.opacity = 1.

        self.near_pick = None

        self.far_pick = None

        self.near_pick_prev = None

        self.far_pick_prev = None

        self.picked_track = None

        self.pick_color = [1,1,0]

        self.brain_color = [1,1,1]

        self.yellow_indices = None

        self.dummy_data = False

        self.tracks = tracks

        if subset != None:

            self.data_subset = subset #[0,20000]#None

        else:

            self.data_subset = None

        self.orbit_demo = False          

        self.orbit_anglez = 0.

        self.orbit_anglez_rate = 10.
        

        self.orbit_anglex = 0.

        self.orbit_anglex_rate = 2.


        self.angle_table = ang_table

        
        self.angle_table_index = 0



        

        self.shrink = shrink

        self.picking_example = False


        self.partial_colors = False
        

        import dipy.io.trackvis as tv


        if self.tracks == None:

            lines,hdr = tv.read(self.fname)

            ras = tv.aff_from_hdr(hdr)

            self.affine=ras

            tracks = [l[0] for l in lines]

            del lines
        


        else:

            tracks = self.tracks


        if self.yellow_indices != None :

            tracks = [t for t in tracks if tm.length(t) > 20]

        print '%d tracks loaded' % len(tracks)

        #self.data = [100*np.array([[0,0,0],[1,0,0],[2,0,0]]).astype(np.float32) ,100*np.array([[0,1,0],[0,2,0],[0,3,0]]).astype(np.float32)]#tracks[:20000]

        if self.dummy_data:

            self.data = [100*np.array([[0,0,0],[1,0,0],[2,0,0]]).astype(np.float32) ,100*np.array([[0,1,0],[0,2,0],[0,3,0]]).astype(np.float32)]

        if self.data_subset!=None:

            self.data = tracks[self.data_subset[0]:self.data_subset[1]]

        else:

            self.data = tracks


        

        if self.shrink != None:

            self.data = [ self.shrink*t  for t in self.data]
            

            
        data_stats = np.concatenate(tracks)

        self.min=np.min(data_stats,axis=0)
         
        self.max=np.max(data_stats,axis=0)

        self.mean=np.mean(data_stats,axis=0)

        del data_stats

Example 112

Project: fos-legacy
Source File: tracks.py
View license
    def __init__(self,fname,colormap=None, line_width=1., shrink=None, thinning = 0,
                 angle_table = None, manycolors = False, brain_color=[1,1,1]):

        self.position = (0,0,0)

        self.fname = fname
        
        self.manycolors = manycolors

        #self.color = monocolor
        
        self.bbox = None

        self.list_index = None

        self.affine = None

        self.data = None

        self.list_index = None

        self.rot_angle = 0

        self.colormap = None
                
        self.min = None
         
        self.max = None

        self.mean = None

        self.material_color = False

        self.fadeout = False

        self.fadein = False

        self.fadeout_speed = 0.

        self.fadein_speed = 0.

        self.min_length = 20.

        self.angle = 0.

        self.angular_speed = .5

        self.line_width = line_width

        self.opacity = 1.

        self.opacity_rate = 0.

        self.near_pick = None

        self.far_pick = None

        self.near_pick_prev = None

        self.far_pick_prev = None

        self.picked_track = None

        self.pick_color = [1,1,0]

        #self.brain_color = [1,1,1] # white
        #self.brain_color = [.941,.862,.510] # buff
        self.brain_color = brain_color
        
        self.yellow_indices = None

        self.dummy_data = False

        self.data_subset = [0,20000]#None


        self.orbit_demo = False          

        self.orbit_anglez =0.

        self.orbit_anglez_rate = 10.
        

        self.orbit_anglex = 0.

        self.orbit_anglex_rate = 2.


        self.angle_table = angle_table

        if angle_table != None:
            print 'Tracks angle_table shape %s' % str(self.angle_table.shape)

        self.angle_table_index = 0

        #print 'angle_table_index %d' % self.angle_table_index

        self.shrink = shrink

        self.picking_example = False

        import dipy.io.trackvis as tv

        lines,hdr = tv.read(self.fname)

        ras = tv.aff_from_hdr(hdr
)
        self.affine=ras

        tracks = [l[0] for l in lines]

        print 'tracks %d loaded' % len(tracks)

        self.thinning = thinning

        if self.yellow_indices != None :

            tracks = [t for t in tracks if tm.length(t) > 20]

        if self.thinning != 0:

            tracks = [tracks[k] for k in range(0,len(tracks),self.thinning)]

        print '%d tracks active' % len(tracks)

        #self.data = [100*np.array([[0,0,0],[1,0,0],[2,0,0]]).astype(np.float32) ,100*np.array([[0,1,0],[0,2,0],[0,3,0]]).astype(np.float32)]#tracks[:20000]

        if self.dummy_data:

            self.data = [100*np.array([[0,0,0],[1,0,0],[2,0,0]]).astype(np.float32) ,100*np.array([[0,1,0],[0,2,0],[0,3,0]]).astype(np.float32)]

        if self.data_subset!=None:

            self.data = tracks[self.data_subset[0]:self.data_subset[1]]

        else:

            self.data = tracks

        if self.shrink != None:

            self.data = [ self.shrink*t  for t in self.data]
            

        data_stats = np.concatenate(tracks)

        self.min=np.min(data_stats,axis=0)
         
        self.max=np.max(data_stats,axis=0)

        self.mean=np.mean(data_stats,axis=0)

        del data_stats
        
        del lines

Example 113

Project: MemN2N-babi-python
Source File: util.py
View license
def build_model(general_config):
    """
    Build model

    NOTE: (for default config)
    1) Model's architecture (embedding B)
        LookupTable -> ElemMult -> Sum -> [ Duplicate -> { Parallel -> Memory -> Identity } -> AddTable ] -> LinearNB -> Softmax

    2) Memory's architecture
        a) Query module (embedding A)
            Parallel -> { LookupTable + ElemMult + Sum } -> Identity -> MatVecProd -> Softmax

        b) Output module (embedding C)
            Parallel -> { LookupTable + ElemMult + Sum } -> Identity -> MatVecProd
    """
    train_config = general_config.train_config
    dictionary   = general_config.dictionary
    use_bow      = general_config.use_bow
    nhops        = general_config.nhops
    add_proj     = general_config.add_proj
    share_type   = general_config.share_type
    enable_time  = general_config.enable_time
    add_nonlin   = general_config.add_nonlin

    in_dim    = train_config["in_dim"]
    out_dim   = train_config["out_dim"]
    max_words = train_config["max_words"]
    voc_sz    = train_config["voc_sz"]

    if not use_bow:
        train_config["weight"] = np.ones((in_dim, max_words), np.float32)
        for i in range(in_dim):
            for j in range(max_words):
                train_config["weight"][i][j] = (i + 1 - (in_dim + 1) / 2) * \
                                               (j + 1 - (max_words + 1) / 2)
        train_config["weight"] = \
            1 + 4 * train_config["weight"] / (in_dim * max_words)

    memory = {}
    model = Sequential()
    model.add(LookupTable(voc_sz, in_dim))
    if not use_bow:
        if enable_time:
            model.add(ElemMult(train_config["weight"][:, :-1]))
        else:
            model.add(ElemMult(train_config["weight"]))

    model.add(Sum(dim=1))

    proj = {}
    for i in range(nhops):
        if use_bow:
            memory[i] = MemoryBoW(train_config)
        else:
            memory[i] = MemoryL(train_config)

        # Override nil_word which is initialized in "self.nil_word = train_config["voc_sz"]"
        memory[i].nil_word = dictionary['nil']
        model.add(Duplicate())
        p = Parallel()
        p.add(memory[i])

        if add_proj:
            proj[i] = LinearNB(in_dim, in_dim)
            p.add(proj[i])
        else:
            p.add(Identity())

        model.add(p)
        model.add(AddTable())
        if add_nonlin:
            model.add(ReLU())

    model.add(LinearNB(out_dim, voc_sz, True))
    model.add(Softmax())

    # Share weights
    if share_type == 1:
        # Type 1: adjacent weight tying
        memory[0].emb_query.share(model.modules[0])
        for i in range(1, nhops):
            memory[i].emb_query.share(memory[i - 1].emb_out)

        model.modules[-2].share(memory[len(memory) - 1].emb_out)

    elif share_type == 2:
        # Type 2: layer-wise weight tying
        for i in range(1, nhops):
            memory[i].emb_query.share(memory[0].emb_query)
            memory[i].emb_out.share(memory[0].emb_out)

    if add_proj:
        for i in range(1, nhops):
            proj[i].share(proj[0])

    # Cost
    loss = CrossEntropyLoss()
    loss.size_average = False
    loss.do_softmax_bprop = True
    model.modules[-1].skip_bprop = True

    return memory, model, loss

Example 114

Project: dilation
Source File: test.py
View license
def test_image(options):
    options.feat_dir = join(options.feat_dir, options.feat_layer_name)
    if not exists(options.feat_dir):
        os.makedirs(options.feat_dir)

    label_margin = 186

    if options.up:
        zoom = 1
    else:
        zoom = 8

    if options.gpu >= 0:
        caffe.set_mode_gpu()
        caffe.set_device(options.gpu)
        print('Using GPU ', options.gpu)
    else:
        caffe.set_mode_cpu()
        print('Using CPU')

    mean_pixel = np.array(options.mean, dtype=np.float32)
    net = caffe.Net(options.deploy_net, options.weights, caffe.TEST)

    image_paths = [line.strip() for line in open(options.image_list, 'r')]
    image_names = [split(p)[1] for p in image_paths]
    input_dims = list(net.blobs['data'].shape)

    assert input_dims[0] == 1
    batch_size, num_channels, input_height, input_width = input_dims
    print('Input size:', input_dims)
    caffe_in = np.zeros(input_dims, dtype=np.float32)

    output_height = input_height - 2 * label_margin
    output_width = input_width - 2 * label_margin

    result_list = []
    feat_list = []

    for i in range(len(image_names)):
        print('Predicting', image_names[i])
        image = cv2.imread(image_paths[i]).astype(np.float32) - mean_pixel
        image_size = image.shape
        print('Image size:', image_size)
        image = cv2.copyMakeBorder(image, label_margin, label_margin,
                                   label_margin, label_margin,
                                   cv2.BORDER_REFLECT_101)
        num_tiles_h = image_size[0] // output_height + \
                      (1 if image_size[0] % output_height else 0)
        num_tiles_w = image_size[1] // output_width + \
                      (1 if image_size[1] % output_width else 0)
        prediction = []
        feat = []
        for h in range(num_tiles_h):
            col_prediction = []
            col_feat = []
            for w in range(num_tiles_w):
                offset = [output_height * h,
                          output_width * w]
                tile = image[offset[0]:offset[0] + input_height,
                             offset[1]:offset[1] + input_width, :]
                margin = [0, input_height - tile.shape[0],
                          0, input_width - tile.shape[1]]
                tile = cv2.copyMakeBorder(tile, margin[0], margin[1],
                                          margin[2], margin[3],
                                          cv2.BORDER_REFLECT_101)
                caffe_in[0] = tile.transpose([2, 0, 1])
                blobs = []
                if options.bin:
                    blobs = [options.feat_layer_name]
                out = net.forward_all(blobs=blobs, **{net.inputs[0]: caffe_in})
                prob = out['prob'][0]
                if options.bin:
                    col_feat.append(out[options.feat_layer_name][0])
                col_prediction.append(prob)
            col_prediction = np.concatenate(col_prediction, axis=2)
            if options.bin:
                col_feat = np.concatenate(col_feat, axis=2)
                feat.append(col_feat)
            prediction.append(col_prediction)
        prob = np.concatenate(prediction, axis=1)
        if options.bin:
            feat = np.concatenate(feat, axis=1)

        if zoom > 1:
            zoom_prob = util.interp_map(
                prob, zoom, image_size[1], image_size[0])
        else:
            zoom_prob = prob[:, :image_size[0], :image_size[1]]
        prediction = np.argmax(zoom_prob.transpose([1, 2, 0]), axis=2)
        if options.bin:
            out_path = join(options.feat_dir,
                            splitext(image_names[i])[0] + '.bin')
            print('Writing', out_path)
            write_array(out_path, feat.astype(np.float32))
            feat_list.append(out_path)
        out_path = join(options.result_dir,
                        splitext(image_names[i])[0] + '.png')
        print('Writing', out_path)
        cv2.imwrite(out_path, prediction)
        result_list.append(out_path)

    print('================================')
    print('All results are generated.')
    print('================================')

    result_list_path = join(options.result_dir, 'results.txt')
    print('Writing', result_list_path)
    with open(result_list_path, 'w') as fp:
        fp.write('\n'.join(result_list))
    if options.bin:
        feat_list_path = join(options.feat_dir, 'feats.txt')
        print('Writing', feat_list_path)
        with open(feat_list_path, 'w') as fp:
            fp.write('\n'.join(feat_list))

Example 115

Project: tagger
Source File: model.py
View license
    def build(self,
              dropout,
              char_dim,
              char_lstm_dim,
              char_bidirect,
              word_dim,
              word_lstm_dim,
              word_bidirect,
              lr_method,
              pre_emb,
              crf,
              cap_dim,
              training=True,
              **kwargs
              ):
        """
        Build the network.
        """
        # Training parameters
        n_words = len(self.id_to_word)
        n_chars = len(self.id_to_char)
        n_tags = len(self.id_to_tag)

        # Number of capitalization features
        if cap_dim:
            n_cap = 4

        # Network variables
        is_train = T.iscalar('is_train')
        word_ids = T.ivector(name='word_ids')
        char_for_ids = T.imatrix(name='char_for_ids')
        char_rev_ids = T.imatrix(name='char_rev_ids')
        char_pos_ids = T.ivector(name='char_pos_ids')
        tag_ids = T.ivector(name='tag_ids')
        if cap_dim:
            cap_ids = T.ivector(name='cap_ids')

        # Sentence length
        s_len = (word_ids if word_dim else char_pos_ids).shape[0]

        # Final input (all word features)
        input_dim = 0
        inputs = []

        #
        # Word inputs
        #
        if word_dim:
            input_dim += word_dim
            word_layer = EmbeddingLayer(n_words, word_dim, name='word_layer')
            word_input = word_layer.link(word_ids)
            inputs.append(word_input)
            # Initialize with pretrained embeddings
            if pre_emb and training:
                new_weights = word_layer.embeddings.get_value()
                print 'Loading pretrained embeddings from %s...' % pre_emb
                pretrained = {}
                emb_invalid = 0
                for i, line in enumerate(codecs.open(pre_emb, 'r', 'utf-8')):
                    line = line.rstrip().split()
                    if len(line) == word_dim + 1:
                        pretrained[line[0]] = np.array(
                            [float(x) for x in line[1:]]
                        ).astype(np.float32)
                    else:
                        emb_invalid += 1
                if emb_invalid > 0:
                    print 'WARNING: %i invalid lines' % emb_invalid
                c_found = 0
                c_lower = 0
                c_zeros = 0
                # Lookup table initialization
                for i in xrange(n_words):
                    word = self.id_to_word[i]
                    if word in pretrained:
                        new_weights[i] = pretrained[word]
                        c_found += 1
                    elif word.lower() in pretrained:
                        new_weights[i] = pretrained[word.lower()]
                        c_lower += 1
                    elif re.sub('\d', '0', word.lower()) in pretrained:
                        new_weights[i] = pretrained[
                            re.sub('\d', '0', word.lower())
                        ]
                        c_zeros += 1
                word_layer.embeddings.set_value(new_weights)
                print 'Loaded %i pretrained embeddings.' % len(pretrained)
                print ('%i / %i (%.4f%%) words have been initialized with '
                       'pretrained embeddings.') % (
                            c_found + c_lower + c_zeros, n_words,
                            100. * (c_found + c_lower + c_zeros) / n_words
                      )
                print ('%i found directly, %i after lowercasing, '
                       '%i after lowercasing + zero.') % (
                          c_found, c_lower, c_zeros
                      )

        #
        # Chars inputs
        #
        if char_dim:
            input_dim += char_lstm_dim
            char_layer = EmbeddingLayer(n_chars, char_dim, name='char_layer')

            char_lstm_for = LSTM(char_dim, char_lstm_dim, with_batch=True,
                                 name='char_lstm_for')
            char_lstm_rev = LSTM(char_dim, char_lstm_dim, with_batch=True,
                                 name='char_lstm_rev')

            char_lstm_for.link(char_layer.link(char_for_ids))
            char_lstm_rev.link(char_layer.link(char_rev_ids))

            char_for_output = char_lstm_for.h.dimshuffle((1, 0, 2))[
                T.arange(s_len), char_pos_ids
            ]
            char_rev_output = char_lstm_rev.h.dimshuffle((1, 0, 2))[
                T.arange(s_len), char_pos_ids
            ]

            inputs.append(char_for_output)
            if char_bidirect:
                inputs.append(char_rev_output)
                input_dim += char_lstm_dim

        #
        # Capitalization feature
        #
        if cap_dim:
            input_dim += cap_dim
            cap_layer = EmbeddingLayer(n_cap, cap_dim, name='cap_layer')
            inputs.append(cap_layer.link(cap_ids))

        # Prepare final input
        if len(inputs) != 1:
            inputs = T.concatenate(inputs, axis=1)

        #
        # Dropout on final input
        #
        if dropout:
            dropout_layer = DropoutLayer(p=dropout)
            input_train = dropout_layer.link(inputs)
            input_test = (1 - dropout) * inputs
            inputs = T.switch(T.neq(is_train, 0), input_train, input_test)

        # LSTM for words
        word_lstm_for = LSTM(input_dim, word_lstm_dim, with_batch=False,
                             name='word_lstm_for')
        word_lstm_rev = LSTM(input_dim, word_lstm_dim, with_batch=False,
                             name='word_lstm_rev')
        word_lstm_for.link(inputs)
        word_lstm_rev.link(inputs[::-1, :])
        word_for_output = word_lstm_for.h
        word_rev_output = word_lstm_rev.h[::-1, :]
        if word_bidirect:
            final_output = T.concatenate(
                [word_for_output, word_rev_output],
                axis=1
            )
            tanh_layer = HiddenLayer(2 * word_lstm_dim, word_lstm_dim,
                                     name='tanh_layer', activation='tanh')
            final_output = tanh_layer.link(final_output)
        else:
            final_output = word_for_output

        # Sentence to Named Entity tags - Score
        final_layer = HiddenLayer(word_lstm_dim, n_tags, name='final_layer',
                                  activation=(None if crf else 'softmax'))
        tags_scores = final_layer.link(final_output)

        # No CRF
        if not crf:
            cost = T.nnet.categorical_crossentropy(tags_scores, tag_ids).mean()
        # CRF
        else:
            transitions = shared((n_tags + 2, n_tags + 2), 'transitions')

            small = -1000
            b_s = np.array([[small] * n_tags + [0, small]]).astype(np.float32)
            e_s = np.array([[small] * n_tags + [small, 0]]).astype(np.float32)
            observations = T.concatenate(
                [tags_scores, small * T.ones((s_len, 2))],
                axis=1
            )
            observations = T.concatenate(
                [b_s, observations, e_s],
                axis=0
            )

            # Score from tags
            real_path_score = tags_scores[T.arange(s_len), tag_ids].sum()

            # Score from transitions
            b_id = theano.shared(value=np.array([n_tags], dtype=np.int32))
            e_id = theano.shared(value=np.array([n_tags + 1], dtype=np.int32))
            padded_tags_ids = T.concatenate([b_id, tag_ids, e_id], axis=0)
            real_path_score += transitions[
                padded_tags_ids[T.arange(s_len + 1)],
                padded_tags_ids[T.arange(s_len + 1) + 1]
            ].sum()

            all_paths_scores = forward(observations, transitions)
            cost = - (real_path_score - all_paths_scores)

        # Network parameters
        params = []
        if word_dim:
            self.add_component(word_layer)
            params.extend(word_layer.params)
        if char_dim:
            self.add_component(char_layer)
            self.add_component(char_lstm_for)
            params.extend(char_layer.params)
            params.extend(char_lstm_for.params)
            if char_bidirect:
                self.add_component(char_lstm_rev)
                params.extend(char_lstm_rev.params)
        self.add_component(word_lstm_for)
        params.extend(word_lstm_for.params)
        if word_bidirect:
            self.add_component(word_lstm_rev)
            params.extend(word_lstm_rev.params)
        if cap_dim:
            self.add_component(cap_layer)
            params.extend(cap_layer.params)
        self.add_component(final_layer)
        params.extend(final_layer.params)
        if crf:
            self.add_component(transitions)
            params.append(transitions)
        if word_bidirect:
            self.add_component(tanh_layer)
            params.extend(tanh_layer.params)

        # Prepare train and eval inputs
        eval_inputs = []
        if word_dim:
            eval_inputs.append(word_ids)
        if char_dim:
            eval_inputs.append(char_for_ids)
            if char_bidirect:
                eval_inputs.append(char_rev_ids)
            eval_inputs.append(char_pos_ids)
        if cap_dim:
            eval_inputs.append(cap_ids)
        train_inputs = eval_inputs + [tag_ids]

        # Parse optimization method parameters
        if "-" in lr_method:
            lr_method_name = lr_method[:lr_method.find('-')]
            lr_method_parameters = {}
            for x in lr_method[lr_method.find('-') + 1:].split('-'):
                split = x.split('_')
                assert len(split) == 2
                lr_method_parameters[split[0]] = float(split[1])
        else:
            lr_method_name = lr_method
            lr_method_parameters = {}

        # Compile training function
        print 'Compiling...'
        if training:
            updates = Optimization(clip=5.0).get_updates(lr_method_name, cost, params, **lr_method_parameters)
            f_train = theano.function(
                inputs=train_inputs,
                outputs=cost,
                updates=updates,
                givens=({is_train: np.cast['int32'](1)} if dropout else {})
            )
        else:
            f_train = None

        # Compile evaluation function
        if not crf:
            f_eval = theano.function(
                inputs=eval_inputs,
                outputs=tags_scores,
                givens=({is_train: np.cast['int32'](0)} if dropout else {})
            )
        else:
            f_eval = theano.function(
                inputs=eval_inputs,
                outputs=forward(observations, transitions, viterbi=True,
                                return_alpha=False, return_best_sequence=True),
                givens=({is_train: np.cast['int32'](0)} if dropout else {})
            )

        return f_train, f_eval

Example 116

Project: scikit-image
Source File: thresholding.py
View license
def threshold_isodata(image, nbins=256, return_all=False):
    """Return threshold value(s) based on ISODATA method.

    Histogram-based threshold, known as Ridler-Calvard method or inter-means.
    Threshold values returned satisfy the following equality:

    `threshold = (image[image <= threshold].mean() +`
                 `image[image > threshold].mean()) / 2.0`

    That is, returned thresholds are intensities that separate the image into
    two groups of pixels, where the threshold intensity is midway between the
    mean intensities of these groups.

    For integer images, the above equality holds to within one; for floating-
    point images, the equality holds to within the histogram bin-width.

    Parameters
    ----------
    image : (N, M) ndarray
        Input image.
    nbins : int, optional
        Number of bins used to calculate histogram. This value is ignored for
        integer arrays.
    return_all: bool, optional
        If False (default), return only the lowest threshold that satisfies
        the above equality. If True, return all valid thresholds.

    Returns
    -------
    threshold : float or int or array
        Threshold value(s).

    References
    ----------
    .. [1] Ridler, TW & Calvard, S (1978), "Picture thresholding using an
           iterative selection method"
           IEEE Transactions on Systems, Man and Cybernetics 8: 630-632,
           DOI:10.1109/TSMC.1978.4310039
    .. [2] Sezgin M. and Sankur B. (2004) "Survey over Image Thresholding
           Techniques and Quantitative Performance Evaluation" Journal of
           Electronic Imaging, 13(1): 146-165,
           http://www.busim.ee.boun.edu.tr/~sankur/SankurFolder/Threshold_survey.pdf
           DOI:10.1117/1.1631315
    .. [3] ImageJ AutoThresholder code,
           http://fiji.sc/wiki/index.php/Auto_Threshold

    Examples
    --------
    >>> from skimage.data import coins
    >>> image = coins()
    >>> thresh = threshold_isodata(image)
    >>> binary = image > thresh
    """
    hist, bin_centers = histogram(image.ravel(), nbins)

    # image only contains one unique value
    if len(bin_centers) == 1:
        if return_all:
            return bin_centers
        else:
            return bin_centers[0]

    hist = hist.astype(np.float32)

    # csuml and csumh contain the count of pixels in that bin or lower, and
    # in all bins strictly higher than that bin, respectively
    csuml = np.cumsum(hist)
    csumh = np.cumsum(hist[::-1])[::-1] - hist

    # intensity_sum contains the total pixel intensity from each bin
    intensity_sum = hist * bin_centers

    # l and h contain average value of all pixels in that bin or lower, and
    # in all bins strictly higher than that bin, respectively.
    # Note that since exp.histogram does not include empty bins at the low or
    # high end of the range, csuml and csumh are strictly > 0, except in the
    # last bin of csumh, which is zero by construction.
    # So no worries about division by zero in the following lines, except
    # for the last bin, but we can ignore that because no valid threshold
    # can be in the top bin. So we just patch up csumh[-1] to not cause 0/0
    # errors.
    csumh[-1] = 1
    l = np.cumsum(intensity_sum) / csuml
    h = (np.cumsum(intensity_sum[::-1])[::-1] - intensity_sum) / csumh

    # isodata finds threshold values that meet the criterion t = (l + m)/2
    # where l is the mean of all pixels <= t and h is the mean of all pixels
    # > t, as calculated above. So we are looking for places where
    # (l + m) / 2 equals the intensity value for which those l and m figures
    # were calculated -- which is, of course, the histogram bin centers.
    # We only require this equality to be within the precision of the bin
    # width, of course.
    all_mean = (l + h) / 2.0
    bin_width = bin_centers[1] - bin_centers[0]

    # Look only at thresholds that are below the actual all_mean value,
    # for consistency with the threshold being included in the lower pixel
    # group. Otherwise can get thresholds that are not actually fixed-points
    # of the isodata algorithm. For float images, this matters less, since
    # there really can't be any guarantees anymore anyway.
    distances = all_mean - bin_centers
    thresholds = bin_centers[(distances >= 0) & (distances < bin_width)]

    if return_all:
        return thresholds
    else:
        return thresholds[0]

Example 117

View license
def main_test_stacked_denoise_AE(model='relu'):
    X_train, y_train, X_val, y_val, X_test, y_test = \
                                tl.files.load_mnist_dataset(shape=(-1,784))

    X_train = np.asarray(X_train, dtype=np.float32)
    y_train = np.asarray(y_train, dtype=np.int64)
    X_val = np.asarray(X_val, dtype=np.float32)
    y_val = np.asarray(y_val, dtype=np.int64)
    X_test = np.asarray(X_test, dtype=np.float32)
    y_test = np.asarray(y_test, dtype=np.int64)

    print('X_train.shape', X_train.shape)
    print('y_train.shape', y_train.shape)
    print('X_val.shape', X_val.shape)
    print('y_val.shape', y_val.shape)
    print('X_test.shape', X_test.shape)
    print('y_test.shape', y_test.shape)
    print('X %s   y %s' % (X_test.dtype, y_test.dtype))

    sess = tf.InteractiveSession()

    x = tf.placeholder(tf.float32, shape=[None, 784], name='x')
    y_ = tf.placeholder(tf.int64, shape=[None, ], name='y_')

    if model == 'relu':
        act = tf.nn.relu
        act_recon = tf.nn.softplus
    elif model == 'sigmoid':
        act = tf.nn.sigmoid
        act_recon = act

    # Define network
    print("\nBuild Network")
    network = tl.layers.InputLayer(x, name='input_layer')
    # denoise layer for AE
    network = tl.layers.DropoutLayer(network, keep=0.5, name='denoising1')
    # 1st layer
    network = tl.layers.DropoutLayer(network, keep=0.8, name='drop1')
    network = tl.layers.DenseLayer(network, n_units=800, act = act, name=model+'1')
    x_recon1 = network.outputs
    recon_layer1 = tl.layers.ReconLayer(network, x_recon=x, n_units=784,
                                        act = act_recon, name='recon_layer1')
    # 2nd layer
    network = tl.layers.DropoutLayer(network, keep=0.5, name='drop2')
    network = tl.layers.DenseLayer(network, n_units=800, act = act, name=model+'2')
    recon_layer2 = tl.layers.ReconLayer(network, x_recon=x_recon1, n_units=800,
                                        act = act_recon, name='recon_layer2')
    # 3rd layer
    network = tl.layers.DropoutLayer(network, keep=0.5, name='drop3')
    network = tl.layers.DenseLayer(network, n_units=10,
                                            act = tf.identity,
                                            name='output_layer')

    # Define fine-tune process
    y = network.outputs
    y_op = tf.argmax(tf.nn.softmax(y), 1)
    ce = tf.reduce_mean(tf.nn.sparse_softmax_cross_entropy_with_logits(y, y_))
    cost = ce

    n_epoch = 200
    batch_size = 128
    learning_rate = 0.0001
    print_freq = 10

    train_params = network.all_params

        # train_op = tf.train.GradientDescentOptimizer(0.5).minimize(cost)
    train_op = tf.train.AdamOptimizer(learning_rate , beta1=0.9, beta2=0.999,
        epsilon=1e-08, use_locking=False).minimize(cost, var_list=train_params)

    # Initialize all variables including weights, biases and the variables in train_op
    sess.run(tf.initialize_all_variables())

    # Pre-train
    print("\nAll Network Params before pre-train")
    network.print_params()
    print("\nPre-train Layer 1")
    recon_layer1.pretrain(sess, x=x, X_train=X_train, X_val=X_val,
                            denoise_name='denoising1', n_epoch=200,
                            batch_size=128, print_freq=10, save=True,
                            save_name='w1pre_')
    print("\nPre-train Layer 2")
    recon_layer2.pretrain(sess, x=x, X_train=X_train, X_val=X_val,
                            denoise_name='denoising1', n_epoch=200,
                            batch_size=128, print_freq=10, save=False)
    print("\nAll Network Params after pre-train")
    network.print_params()

    # Fine-tune
    print("\nFine-tune Network")
    correct_prediction = tf.equal(tf.argmax(y, 1), y_)
    acc = tf.reduce_mean(tf.cast(correct_prediction, tf.float32))

    print('   learning_rate: %f' % learning_rate)
    print('   batch_size: %d' % batch_size)

    for epoch in range(n_epoch):
        start_time = time.time()
        for X_train_a, y_train_a in tl.iterate.minibatches(
                                    X_train, y_train, batch_size, shuffle=True):
            feed_dict = {x: X_train_a, y_: y_train_a}
            feed_dict.update( network.all_drop )     # enable noise layers
            feed_dict[set_keep['denoising1']] = 1    # disable denoising layer
            sess.run(train_op, feed_dict=feed_dict)

        if epoch + 1 == 1 or (epoch + 1) % print_freq == 0:
            print("Epoch %d of %d took %fs" % (epoch + 1, n_epoch, time.time() - start_time))
            train_loss, train_acc, n_batch = 0, 0, 0
            for X_train_a, y_train_a in tl.iterate.minibatches(
                                    X_train, y_train, batch_size, shuffle=True):
                dp_dict = tl.utils.dict_to_one( network.all_drop )    # disable noise layers
                feed_dict = {x: X_train_a, y_: y_train_a}
                feed_dict.update(dp_dict)
                err, ac = sess.run([cost, acc], feed_dict=feed_dict)
                train_loss += err
                train_acc += ac
                n_batch += 1
            print("   train loss: %f" % (train_loss/ n_batch))
            print("   train acc: %f" % (train_acc/ n_batch))
            val_loss, val_acc, n_batch = 0, 0, 0
            for X_val_a, y_val_a in tl.iterate.minibatches(
                                        X_val, y_val, batch_size, shuffle=True):
                dp_dict = tl.utils.dict_to_one( network.all_drop )    # disable noise layers
                feed_dict = {x: X_val_a, y_: y_val_a}
                feed_dict.update(dp_dict)
                err, ac = sess.run([cost, acc], feed_dict=feed_dict)
                val_loss += err
                val_acc += ac
                n_batch += 1
            print("   val loss: %f" % (val_loss/ n_batch))
            print("   val acc: %f" % (val_acc/ n_batch))
            try:
                # visualize the 1st hidden layer during fine-tune
                tl.visualize.W(network.all_params[0].eval(), second=10,
                            saveable=True, shape=[28, 28],
                            name='w1_'+str(epoch+1), fig_idx=2012)
            except:
                raise Exception("# You should change visualize_W(), if you want to save the feature images for different dataset")

    print('Evaluation')
    test_loss, test_acc, n_batch = 0, 0, 0
    for X_test_a, y_test_a in tl.iterate.minibatches(
                                X_test, y_test, batch_size, shuffle=True):
        dp_dict = tl.utils.dict_to_one( network.all_drop )    # disable noise layers
        feed_dict = {x: X_test_a, y_: y_test_a}
        feed_dict.update(dp_dict)
        err, ac = sess.run([cost, acc], feed_dict=feed_dict)
        test_loss += err
        test_acc += ac
        n_batch += 1
    print("   test loss: %f" % (test_loss/n_batch))
    print("   test acc: %f" % (test_acc/n_batch))
        # print("   test acc: %f" % np.mean(y_test == sess.run(y_op, feed_dict=feed_dict)))

    # Add ops to save and restore all the variables.
    # ref: https://www.tensorflow.org/versions/r0.8/how_tos/variables/index.html
    saver = tf.train.Saver()
    # you may want to save the model
    save_path = saver.save(sess, "model.ckpt")
    print("Model saved in file: %s" % save_path)
    sess.close()

Example 118

Project: vispy
Source File: generation.py
View license
def create_box(width=1, height=1, depth=1, width_segments=1, height_segments=1,
               depth_segments=1, planes=None):
    """ Generate vertices & indices for a filled and outlined box.

    Parameters
    ----------
    width : float
        Box width.
    height : float
        Box height.
    depth : float
        Box depth.
    width_segments : int
        Box segments count along the width.
    height_segments : float
        Box segments count along the height.
    depth_segments : float
        Box segments count along the depth.
    planes: array_like
        Any combination of ``{'-x', '+x', '-y', '+y', '-z', '+z'}``
        Included planes in the box construction.

    Returns
    -------
    vertices : array
        Array of vertices suitable for use as a VertexBuffer.
    faces : array
        Indices to use to produce a filled box.
    outline : array
        Indices to use to produce an outline of the box.
    """

    planes = (('+x', '-x', '+y', '-y', '+z', '-z')
              if planes is None else
              [d.lower() for d in planes])

    w_s, h_s, d_s = width_segments, height_segments, depth_segments

    planes_m = []
    if '-z' in planes:
        planes_m.append(create_plane(width, depth, w_s, d_s, '-z'))
        planes_m[-1][0]['position'][..., 2] -= height / 2
    if '+z' in planes:
        planes_m.append(create_plane(width, depth, w_s, d_s, '+z'))
        planes_m[-1][0]['position'][..., 2] += height / 2

    if '-y' in planes:
        planes_m.append(create_plane(height, width, h_s, w_s, '-y'))
        planes_m[-1][0]['position'][..., 1] -= depth / 2
    if '+y' in planes:
        planes_m.append(create_plane(height, width, h_s, w_s, '+y'))
        planes_m[-1][0]['position'][..., 1] += depth / 2

    if '-x' in planes:
        planes_m.append(create_plane(depth, height, d_s, h_s, '-x'))
        planes_m[-1][0]['position'][..., 0] -= width / 2
    if '+x' in planes:
        planes_m.append(create_plane(depth, height, d_s, h_s, '+x'))
        planes_m[-1][0]['position'][..., 0] += width / 2

    positions = np.zeros((0, 3), dtype=np.float32)
    texcoords = np.zeros((0, 2), dtype=np.float32)
    normals = np.zeros((0, 3), dtype=np.float32)

    faces = np.zeros((0, 3), dtype=np.uint32)
    outline = np.zeros((0, 2), dtype=np.uint32)

    offset = 0
    for vertices_p, faces_p, outline_p in planes_m:
        positions = np.vstack((positions, vertices_p['position']))
        texcoords = np.vstack((texcoords, vertices_p['texcoord']))
        normals = np.vstack((normals, vertices_p['normal']))

        faces = np.vstack((faces, faces_p + offset))
        outline = np.vstack((outline, outline_p + offset))
        offset += vertices_p['position'].shape[0]

    vertices = np.zeros(positions.shape[0],
                        [('position', np.float32, 3),
                         ('texcoord', np.float32, 2),
                         ('normal', np.float32, 3),
                         ('color', np.float32, 4)])

    colors = np.ravel(positions)
    colors = np.hstack((np.reshape(np.interp(colors,
                                             (np.min(colors),
                                              np.max(colors)),
                                             (0, 1)),
                                   positions.shape),
                        np.ones((positions.shape[0], 1))))

    vertices['position'] = positions
    vertices['texcoord'] = texcoords
    vertices['normal'] = normals
    vertices['color'] = colors

    return vertices, faces, outline

Example 119

Project: glumpy
Source File: primitives.py
View license
def tube(top_radius=1.0, base_radius=1.0, height=2.0, slices=32, caps=(True,True)):
    """
    Z-axis aligned tube centered at origin.

    Parameters
    ----------
    top_radius : float
        The radius at the top of the cylinder
    base_radius : float
        The radius at the base of the cylinder
    height : float
        The height of the cylinder
    slices : float
        The number of subdivisions around the Z axis.
    caps : bool, bool
        Whether to add caps at top and bottom
    """

    vtype = [('position', np.float32, 3),
             ('texcoord', np.float32, 2),
             ('normal',   np.float32, 3)]
    itype = np.uint32

    if top_radius <= 0.0:
        caps[0] = False
    if base_radius <= 0.0:
        caps[1] = False

    theta  = np.linspace(0,2*np.pi, slices+1, endpoint=True)
    cos_theta = np.cos(theta)
    sin_theta = np.sin(theta)

    n_side = slices+1
    n_cap  = slices+1
    n = 2*n_side + caps[0]*n_cap + caps[1]*n_cap

    i0 = 0           # Start index of base side
    i1 = i0 + n_side # Start index of top side

    vertices = np.zeros(n, dtype=vtype)
    sides    = vertices[:2*n_side]
    topside  = sides[:+n_side]
    baseside = sides[-n_side:]
    if caps[0] and caps[1]:
        C = vertices[-2*n_cap:]
        topcap   = C[:+n_cap]
        basecap  = C[-n_cap:]
        i2 = i1 + n_side # Start index of top cap
        i3 = i2 + n_cap  # Start index of base cap
    elif caps[0]:
        topcap = vertices[-n_cap:]
        basecap = None
        i2 = i1 + n_side # Start index of top cap
        i3 = None
    elif caps[1]:
        topcap = None
        basecap = vertices[-n_cap:]
        i2 = None
        i3 = i1 + n_side # Start index of top cap

    r = (top_radius-base_radius)
    hypothenus = np.sqrt(r*r + height*height)

    topside["position"][:,0] = cos_theta * top_radius
    topside["position"][:,1] = sin_theta * top_radius
    topside["position"][:,2] = +height/2.0
    topside["texcoord"][:,0] = np.linspace(0,1,slices+1,endpoint=True)
    topside["texcoord"][:,1] = 0
    topside["normal"][:,0] = cos_theta * height/hypothenus
    topside["normal"][:,1] = sin_theta * height/hypothenus
    topside["normal"][:,2] = -(top_radius-base_radius)/hypothenus


    baseside["position"][:,0] = cos_theta * base_radius
    baseside["position"][:,1] = sin_theta * base_radius
    baseside["position"][:,2] = -height/2.0
    baseside["texcoord"][:,0] = np.linspace(0,1,slices+1,endpoint=True)
    baseside["texcoord"][:,1] = 1
    baseside["normal"][:,0] = cos_theta * height/hypothenus
    baseside["normal"][:,1] = sin_theta * height/hypothenus
    baseside["normal"][:,2] = -(top_radius-base_radius)/hypothenus

    if caps[0]:
        topcap["position"][1:] = topside["position"][:-1]
        topcap["position"][0]  = 0, 0, +height/2.0
        topcap["texcoord"][1:,0] = (1+cos_theta[:-1])/2
        topcap["texcoord"][1:,1] = (1+sin_theta[:-1])/2
        topcap["texcoord"][0]  = 0.5, 0.5
        topcap["normal"]  = 0,0,+1
    if caps[1]:
        basecap["position"][1:] = baseside["position"][:-1]
        basecap["position"][0]  = 0, 0, -height/2.0
        basecap["texcoord"][1:,0] = (1+cos_theta[:-1])/2
        basecap["texcoord"][1:,1] = (1+sin_theta[:-1])/2
        basecap["texcoord"][0]  = 0.5, 0.5
        basecap["normal"]  = 0,0,-1

    indices = []
    # side triangles
    for i in range(slices):
        indices.extend([i0+i, i0+i+1, i1+i  ])
        indices.extend([i1+i, i1+i+1, i0+i+1])

    # caps triangles
    for i in range(slices+1):
        if caps[0]:
            indices.extend([i2,i2+1+i%slices,i2+1+(i+1)%slices])
        if caps[1]:
            indices.extend([i3,i3+1+i%slices,i3+1+(i+1)%slices])

    indices = np.array(indices, dtype=itype)
    return vertices.view(gloo.VertexBuffer), indices.view(gloo.IndexBuffer)

Example 120

Project: glumpy
Source File: agg_path_collection.py
View license
    def __init__(self, user_dtype=None, transform=None, viewport=None,
                 vertex=None, fragment=None, **kwargs):
        """
        Initialize the collection.

        Parameters
        ----------

        user_dtype: list
            The base dtype can be completed (appended) by the used_dtype. It
            only make sense if user also provide vertex and/or fragment shaders

        viewport: glumpy.Transforms
            The viewport to use to rende the collection

        transform: glumpy.Tranforms
            The default vertex shader apply the supplied transform to the
            vertices positions before computing the actual vertices positions
            for path thickness. Note that it is necessary to add the
            glumpy.transforms.Viewport transform at the end of the supplied transform.

        vertex: string
            Vertex shader code

        fragment: string
            Fragment  shader code

        caps : string
            'local', 'shared' or 'global'

        join : string
            'local', 'shared' or 'global'

        color : string
            'local', 'shared' or 'global'

        miter_limit : string
            'local', 'shared' or 'global'

        linewidth : string
            'local', 'shared' or 'global'

        antialias : string
            'local', 'shared' or 'global'
        """

        base_dtype = [ ('p0',         (np.float32, 3), '!local', (0,0,0)),
                       ('p1',         (np.float32, 3), '!local', (0,0,0)),
                       ('p2',         (np.float32, 3), '!local', (0,0,0)),
                       ('p3',         (np.float32, 3), '!local', (0,0,0)),
                       ('uv',         (np.float32, 2), '!local', (0,0)),

                       ('caps',       (np.float32, 2), 'global', (0,0)),
                       ('join',       (np.float32, 1), 'global', 0),
                       ('color',      (np.float32, 4), 'global', (0,0,0,1)),
                       ('miter_limit',(np.float32, 1), 'global', 4),
                       ('linewidth',  (np.float32, 1), 'global', 1),
                       ('antialias',  (np.float32, 1), 'global', 1) ]

        dtype = base_dtype
        if user_dtype:
            dtype.extend(user_dtype)

        if vertex is None:
            vertex = library.get('collections/agg-path.vert')
        if fragment is None:
            fragment = library.get('collections/agg-path.frag')

        Collection.__init__(self, dtype=dtype, itype=np.uint32, mode=gl.GL_TRIANGLES,
                            vertex=vertex, fragment=fragment, **kwargs)

        # Set hooks if necessary
        program = self._programs[0]
        if "transform" in program.hooks:
            if transform is not None:
                program["transform"] = transform
            else:
                program["transform"] = Position() + Viewport()

        if "viewport" in program.hooks:
            if viewport is not None:
                program["viewport"] = viewport
            else:
                program["viewport"] = Viewport()

Example 121

Project: hedge
Source File: box-in-box.py
View license
def main():
    from hedge.backends import guess_run_context
    rcon = guess_run_context(["cuda"])

    if rcon.is_head_rank:
        mesh = make_boxmesh()
        #from hedge.mesh import make_rect_mesh
        #mesh = make_rect_mesh(
        #       boundary_tagger=lambda fvi, el, fn, all_v: ["inflow"])
        mesh_data = rcon.distribute_mesh(mesh)
    else:
        mesh_data = rcon.receive_mesh()

    for order in [3]:
        from pytools import add_python_path_relative_to_script
        add_python_path_relative_to_script("..")

        from gas_dynamics_initials import UniformMachFlow
        box = UniformMachFlow(angle_of_attack=0)

        from hedge.models.gas_dynamics import GasDynamicsOperator
        op = GasDynamicsOperator(dimensions=3,
                gamma=box.gamma, mu=box.mu,
                prandtl=box.prandtl, spec_gas_const=box.spec_gas_const,
                bc_inflow=box, bc_outflow=box, bc_noslip=box,
                inflow_tag="inflow", outflow_tag="outflow", noslip_tag="noslip")

        discr = rcon.make_discretization(mesh_data, order=order,
                        debug=[
                            #"cuda_no_plan",
                            #"cuda_dump_kernels",
                            #"dump_dataflow_graph",
                            #"dump_optemplate_stages",
                            #"dump_dataflow_graph",
                            #"print_op_code",
                            "cuda_no_plan_el_local",
                            ],
                        default_scalar_type=numpy.float32,
                        tune_for=op.op_template())

        from hedge.visualization import SiloVisualizer, VtkVisualizer  # noqa
        #vis = VtkVisualizer(discr, rcon, "shearflow-%d" % order)
        vis = SiloVisualizer(discr, rcon)

        fields = box.volume_interpolant(0, discr)

        navierstokes_ex = op.bind(discr)

        max_eigval = [0]

        def rhs(t, q):
            ode_rhs, speed = navierstokes_ex(t, q)
            max_eigval[0] = speed
            return ode_rhs

        rhs(0, fields)

        if rcon.is_head_rank:
            print "---------------------------------------------"
            print "order %d" % order
            print "---------------------------------------------"
            print "#elements=", len(mesh.elements)

        from hedge.timestep import RK4TimeStepper
        stepper = RK4TimeStepper()

        # diagnostics setup ---------------------------------------------------
        from pytools.log import LogManager, add_general_quantities, \
                add_simulation_quantities, add_run_info

        logmgr = LogManager("navierstokes-%d.dat" % order, "w", rcon.communicator)
        add_run_info(logmgr)
        add_general_quantities(logmgr)
        add_simulation_quantities(logmgr)
        discr.add_instrumentation(logmgr)
        stepper.add_instrumentation(logmgr)

        logmgr.add_watches(["step.max", "t_sim.max", "t_step.max"])

        from pytools.log import LogQuantity

        class ChangeSinceLastStep(LogQuantity):
            """Records the change of a variable between a time step and the previous
               one"""

            def __init__(self, name="change"):
                LogQuantity.__init__(self, name, "1", "Change since last time step")

                self.old_fields = 0

            def __call__(self):
                result = discr.norm(fields - self.old_fields)
                self.old_fields = fields
                return result

        logmgr.add_quantity(ChangeSinceLastStep())

        # timestep loop -------------------------------------------------------
        try:
            from hedge.timestep import times_and_steps
            step_it = times_and_steps(
                    final_time=200,
                    #max_steps=500,
                    logmgr=logmgr,
                    max_dt_getter=lambda t: op.estimate_timestep(discr,
                        stepper=stepper, t=t, max_eigenvalue=max_eigval[0]))

            for step, t, dt in step_it:
                if step % 200 == 0:
                #if False:
                    visf = vis.make_file("box-%d-%06d" % (order, step))

                    #rhs_fields = rhs(t, fields)

                    vis.add_data(visf,
                            [
                                ("rho", discr.convert_volume(
                                    op.rho(fields), kind="numpy")),
                                ("e", discr.convert_volume(
                                    op.e(fields), kind="numpy")),
                                ("rho_u", discr.convert_volume(
                                    op.rho_u(fields), kind="numpy")),
                                ("u", discr.convert_volume(
                                    op.u(fields), kind="numpy")),

                                # ("rhs_rho", discr.convert_volume(
                                #     op.rho(rhs_fields), kind="numpy")),
                                # ("rhs_e", discr.convert_volume(
                                #     op.e(rhs_fields), kind="numpy")),
                                # ("rhs_rho_u", discr.convert_volume(
                                #     op.rho_u(rhs_fields), kind="numpy")),
                                ],
                            expressions=[
                                ("p", "(0.4)*(e- 0.5*(rho_u*u))"),
                                ],
                            time=t, step=step
                            )
                    visf.close()

                fields = stepper(fields, t, dt, rhs)

        finally:
            vis.close()
            logmgr.save()
            discr.close()

Example 122

Project: vispy
Source File: surface_plot.py
View license
    def set_data(self, x=None, y=None, z=None, colors=None):
        """Update the data in this surface plot.

        Parameters
        ----------
        x : ndarray | None
            1D array of values specifying the x positions of vertices in the
            grid. If None, values will be assumed to be integers.
        y : ndarray | None
            1D array of values specifying the x positions of vertices in the
            grid. If None, values will be assumed to be integers.
        z : ndarray
            2D array of height values for each grid vertex.
        colors : ndarray
            (width, height, 4) array of vertex colors.
        """
        if x is not None:
            if self._x is None or len(x) != len(self._x):
                self.__vertices = None
            self._x = x

        if y is not None:
            if self._y is None or len(y) != len(self._y):
                self.__vertices = None
            self._y = y

        if z is not None:
            if self._x is not None and z.shape[0] != len(self._x):
                raise TypeError('Z values must have shape (len(x), len(y))')
            if self._y is not None and z.shape[1] != len(self._y):
                raise TypeError('Z values must have shape (len(x), len(y))')
            self._z = z
            if (self.__vertices is not None and
                    self._z.shape != self.__vertices.shape[:2]):
                self.__vertices = None

        if self._z is None:
            return

        update_mesh = False
        new_vertices = False

        # Generate vertex and face array
        if self.__vertices is None:
            new_vertices = True
            self.__vertices = np.empty((self._z.shape[0], self._z.shape[1], 3),
                                       dtype=np.float32)
            self.generate_faces()
            self.__meshdata.set_faces(self.__faces)
            update_mesh = True

        # Copy x, y, z data into vertex array
        if new_vertices or x is not None:
            if x is None:
                if self._x is None:
                    x = np.arange(self._z.shape[0])
                else:
                    x = self._x
            self.__vertices[:, :, 0] = x.reshape(len(x), 1)
            update_mesh = True

        if new_vertices or y is not None:
            if y is None:
                if self._y is None:
                    y = np.arange(self._z.shape[1])
                else:
                    y = self._y
            self.__vertices[:, :, 1] = y.reshape(1, len(y))
            update_mesh = True

        if new_vertices or z is not None:
            self.__vertices[..., 2] = self._z
            update_mesh = True

        if colors is not None:
            self.__meshdata.set_vertex_colors(colors)
            update_mesh = True

        # Update MeshData
        if update_mesh:
            self.__meshdata.set_vertices(
                self.__vertices.reshape(self.__vertices.shape[0] *
                                        self.__vertices.shape[1], 3))
            MeshVisual.set_data(self, meshdata=self.__meshdata)

Example 123

Project: DeepDreamVideo
Source File: 2_dreaming_time.py
View license
def main(input, output, image_type, gpu, model_path, model_name, preview, octaves, octave_scale, iterations, jitter, zoom, stepsize, blend, layers, guide_image, start_frame, end_frame, verbose):
    make_sure_path_exists(input)
    make_sure_path_exists(output)

     # let max nr of frames
    nrframes =len([name for name in os.listdir(input) if os.path.isfile(os.path.join(input, name))])
    if nrframes == 0:
        print("no frames to process found")
        sys.exit(0)

    if preview is None: preview = 0
    if octaves is None: octaves = 4
    if octave_scale is None: octave_scale = 1.5
    if iterations is None: iterations = 5
    if jitter is None: jitter = 32
    if zoom is None: zoom = 1
    if stepsize is None: stepsize = 1.5
    if blend is None: blend = 0.5 #can be nr (constant), random, or loop
    if verbose is None: verbose = 1
    if layers is None: layers = 'customloop' #['inception_4c/output']
    if start_frame is None:
    	frame_i = 1
    else:
        frame_i = int(start_frame)
    if not end_frame is None:
    	nrframes = int(end_frame)+1
    else:
        nrframes = nrframes+1

    #Load DNN
    net_fn   = model_path + 'deploy.prototxt'
    param_fn = model_path + model_name #'bvlc_googlenet.caffemodel'

    if gpu is None:
        print("SHITTTTTTTTTTTTTT You're running CPU man =D")
    else:
        caffe.set_mode_gpu()
        caffe.set_device(int(args.gpu))
        print("GPU mode [device id: %s]" % args.gpu)
        print("using GPU, but you'd still better make a cup of coffee")

    # Patching model to be able to compute gradients.
    # Note that you can also manually add "force_backward: true" line to "deploy.prototxt".
    model = caffe.io.caffe_pb2.NetParameter()
    text_format.Merge(open(net_fn).read(), model)
    model.force_backward = True
    open('tmp.prototxt', 'w').write(str(model))

    net = caffe.Classifier('tmp.prototxt', param_fn,
                           mean = np.float32([104.0, 116.0, 122.0]), # ImageNet mean, training set dependent
                           channel_swap = (2,1,0)) # the reference model has channels in BGR order instead of RGB

    if verbose == 3:
        from IPython.display import clear_output, Image, display
        print("display turned on")
    frame = np.float32(PIL.Image.open(input + '/%08d.%s' % (frame_i, image_type) ))
    if preview is not 0:
        frame = np.float32(resizePicture(input + '/%08d.%s' % (frame_i, image_type), preview))

    now = time.time()
    totaltime = 0
    
    if blend == 'loop':
        blend_forward = True
        blend_at = 0.4
        blend_step = 0.1

    for i in xrange(frame_i, nrframes):
        print('Processing frame #{}').format(frame_i)

        #Choosing Layer
        if layers == 'customloop': #loop over layers as set in layersloop array
            endparam = layersloop[frame_i % len(layersloop)]
        else: #loop through layers one at a time until this specific layer
            endparam = layers[frame_i % len(layers)]

        #Choosing between normal dreaming, and guided dreaming
        if guide_image is None:
            frame = deepdream(net, frame, image_type=image_type, verbose=verbose, iter_n = iterations, step_size = stepsize, octave_n = octaves, octave_scale = octave_scale, jitter=jitter, end = endparam)
        else:
            guide = np.float32(PIL.Image.open(guide_image))
            print('Setting up Guide with selected image')
            guide_features = prepare_guide(net,PIL.Image.open(guide_image), end=endparam)

            frame = deepdream_guided(net, frame, image_type=image_type, verbose=verbose, iter_n = iterations, step_size = stepsize, octave_n = octaves, octave_scale = octave_scale, jitter=jitter, end = endparam, objective_fn=objective_guide, guide_features=guide_features,)

        saveframe = output + "/%08d.%s" % (frame_i, image_type)

        later = time.time()
        difference = int(later - now)
        totaltime += difference
        avgtime = (totaltime / i)
        # Stats (stolen + adapted from Samim: https://github.com/samim23/DeepDreamAnim/blob/master/dreamer.py)
        print '***************************************'
        print 'Saving Image As: ' + saveframe
        print 'Frame ' + str(i) + ' of ' + str(nrframes-1)
        print 'Frame Time: ' + str(difference) + 's'
        timeleft = avgtime * ((nrframes-1) - frame_i)        
        m, s = divmod(timeleft, 60)
        h, m = divmod(m, 60)
        print 'Estimated Total Time Remaining: ' + str(timeleft) + 's (' + "%d:%02d:%02d" % (h, m, s) + ')'
        print '***************************************'

        PIL.Image.fromarray(np.uint8(frame)).save(saveframe)
        newframe = input + "/%08d.%s" % (frame_i,image_type)

        if blend == 0:
            newimg = PIL.Image.open(newframe)
            if preview is not 0:
                newimg = resizePicture(newframe,preview)
            frame = newimg
        else:
       
            if blend == 'random':
            	blendval=randint(5,10)/10.
            elif blend == 'loop':
                if blend_at > 1 - blend_step: blend_forward = False
                elif blend_at <= 0.5: blend_forward = True
                if blend_forward: blend_at += blend_step
                else: blend_at -= blend_step
                blendval = blend_at
            else: blendval = float(blend)
            frame = morphPicture(saveframe,newframe,blendval,preview)

        frame = np.float32(frame)

        now = time.time()
        frame_i += 1

Example 124

Project: kaggle-Rain
Source File: NNregression_v2.py
View license
def do_regression(num_epochs=60, # No. of epochs to train
                  init_file=None,  # Saved parameters to initialise training
                  epoch_size=680780,  # Whole dataset size
                  valid_size=34848, # Size of validation holdout set
                  train_batch_multiple=10637,  # No. of minibatches per batch
                  valid_batch_multiple=1089,  # No. of minibatches per batch
                  train_minibatch_size=64,
                  valid_minibatch_size=32,
                  eval_multiple=50,  # No. of minibatches to ave. in report
                  save_model=True,
                  input_width=19,
                  rng_seed=100005,
                  cross_val=0,  # Cross-validation subset label
                  dataver=2,  # Label for different runs/architectures/etc
                  rate_init=1.0,
                  rate_decay=0.999983):

    ###################################################
    ################# 0. User inputs ##################
    ###################################################
    for i in range(1,len(sys.argv)):
        if sys.argv[i].startswith('-'):
            option = sys.argv[i][1:]
            if option == 'i': init_file = sys.argv[i+1]
            elif option[0:2] == 'v=' : dataver = int(option[2:])
            elif option[0:3] == 'cv=' : cross_val = int(option[3:])
            elif option[0:3] == 'rs=' : rng_seed = int(option[3:])
            elif option[0:3] == 'ri=' : rate_init = np.float32(option[3:])
            elif option[0:3] == 'rd=' : rate_decay = np.float32(option[3:])
                                
    print("Running with dataver %s" % (dataver))
    print("Running with cross_val %s" % (cross_val))
    
    
    ###################################################
    ############# 1. Housekeeping values ##############
    ###################################################
    # Batch size is possibly not equal to epoch size due to memory limits
    train_batch_size = train_batch_multiple*train_minibatch_size 
    assert epoch_size >= train_batch_size
    
    # Number of times we expect the training/validation generator to be called
    max_train_gen_calls = (num_epochs*epoch_size)//train_batch_size 

    # Number of evaluations (total minibatches / eval_multiple)
    num_eval = max_train_gen_calls*train_batch_multiple / eval_multiple
    
    
    ###################################################
    ###### 2. Define model and theano variables #######
    ###################################################
    if rng_seed is not None:
        print("Setting RandomState with seed=%i" % (rng_seed))
        rng = np.random.RandomState(rng_seed)
        set_rng(rng)
    
    print("Defining variables...")
    index = T.lscalar() # Minibatch index
    x = T.tensor3('x') # Inputs 
    y = T.fvector('y') # Target
    
    print("Defining model...")
    network_0 = build_1Dregression_v2(
                        input_var=x,
                        input_width=input_width,
                        h_num_units=[120,120,120],
                        h_grad_clip=1.0,
                        output_width=1
                        )
                        
    if init_file is not None:
        print("Loading initial model parametrs...")
        init_model = np.load(init_file)
        init_params = init_model[init_model.files[0]]           
        LL.set_all_param_values([network_0], init_params)
        
    
    ###################################################                                
    ################ 3. Import data ###################
    ###################################################
    # Loading data generation model parameters
    print("Defining shared variables...")
    train_set_y = theano.shared(np.zeros(1, dtype=theano.config.floatX),
                                borrow=True) 
    train_set_x = theano.shared(np.zeros((1,1,1), dtype=theano.config.floatX),
                                borrow=True)
    
    valid_set_y = theano.shared(np.zeros(1, dtype=theano.config.floatX),
                                borrow=True)
    valid_set_x = theano.shared(np.zeros((1,1,1), dtype=theano.config.floatX),
                                borrow=True)
    
    # Validation data (pick a single augmented instance, rand0 here)
    print("Creating validation data...")    
    chunk_valid_data = np.load(
        "./valid/data_valid_augmented_cv%s_t%s_rand0.npy" 
        % (cross_val, input_width)
        ).astype(theano.config.floatX)
    chunk_valid_answers = np.load(
        "./valid/data_valid_expected_cv%s.npy" 
        % (cross_val)
        ).astype(theano.config.floatX)     
    
    print "chunk_valid_answers.shape", chunk_valid_answers.shape
    print("Assigning validation data...")
    valid_set_y.set_value(chunk_valid_answers[:])
    valid_set_x.set_value(chunk_valid_data.transpose(0,2,1))
    
    # Create output directory
    if not os.path.exists("output_cv%s_v%s" % (cross_val, dataver)):
        os.makedirs("output_cv%s_v%s" % (cross_val, dataver))
    
    
    ###################################################                                
    ########### 4. Create Loss expressions ############
    ###################################################
    print("Defining loss expressions...")
    prediction_0 = LL.get_output(network_0) 
    train_loss = aggregate(T.abs_(prediction_0 - y.dimshuffle(0,'x')))
    
    valid_prediction_0 = LL.get_output(network_0, deterministic=True)
    valid_loss = aggregate(T.abs_(valid_prediction_0 - y.dimshuffle(0,'x')))
    
    
    ###################################################                                
    ############ 5. Define update method  #############
    ###################################################
    print("Defining update choices...")
    params = LL.get_all_params(network_0, trainable=True)
    learn_rate = T.scalar('learn_rate', dtype=theano.config.floatX)
    
    updates = lasagne.updates.adadelta(train_loss, params,
                                       learning_rate=learn_rate)
    
    
    ###################################################                                
    ######### 6. Define train/valid functions #########
    ###################################################    
    print("Defining theano functions...")
    train_model = theano.function(
        [index, learn_rate],
        train_loss,
        updates=updates,
        givens={
            x: train_set_x[(index*train_minibatch_size):
                            ((index+1)*train_minibatch_size)],
            y: train_set_y[(index*train_minibatch_size):
                            ((index+1)*train_minibatch_size)]  
        }
    )
    
    validate_model = theano.function(
        [index],
        valid_loss,
        givens={
            x: valid_set_x[index*valid_minibatch_size:
                            (index+1)*valid_minibatch_size],
            y: valid_set_y[index*valid_minibatch_size:
                            (index+1)*valid_minibatch_size]
        }
    )
    
    
    ###################################################                                
    ################ 7. Begin training ################
    ###################################################  
    print("Begin training...")
    sys.stdout.flush()
    
    cum_iterations = 0
    this_train_loss = 0.0
    this_valid_loss = 0.0
    best_valid_loss = np.inf
    best_iter = 0
    
    train_eval_scores = np.empty(num_eval)
    valid_eval_scores = np.empty(num_eval)
    eval_index = 0
    aug_index = 0
    
    for batch in xrange(max_train_gen_calls):
        start_time = time.time()        
        chunk_train_data = np.load(
            "./train/data_train_augmented_cv%s_t%s_rand%s.npy" %
            (cross_val, input_width, aug_index)
            ).astype(theano.config.floatX)
        chunk_train_answers = np.load(
            "./train/data_train_expected_cv%s.npy" % 
            (cross_val)
            ).astype(theano.config.floatX)     
            
        train_set_y.set_value(chunk_train_answers[:])
        train_set_x.set_value(chunk_train_data.transpose(0, 2, 1))
        
        # Iterate over minibatches in each batch
        for mini_index in xrange(train_batch_multiple):
            this_rate = np.float32(rate_init*(rate_decay**cum_iterations))
            this_train_loss += train_model(mini_index, this_rate)
            cum_iterations += 1
            
            # Report loss 
            if (cum_iterations % eval_multiple == 0):
                this_train_loss = this_train_loss / eval_multiple
                this_valid_loss = np.mean([validate_model(i) for
                                    i in xrange(valid_batch_multiple)])
                train_eval_scores[eval_index] = this_train_loss
                valid_eval_scores[eval_index] = this_valid_loss
                
                # Save report every five evaluations
                if ((eval_index+1) % 5 == 0):
                    np.savetxt(
                        "output_cv%s_v%s/training_scores.txt" %
                        (cross_val, dataver),
                         train_eval_scores, fmt="%.5f"
                         )
                    np.savetxt(
                        "output_cv%s_v%s/validation_scores.txt" %
                        (cross_val, dataver),
                         valid_eval_scores, fmt="%.5f"
                         )
                    np.savetxt(
                        "output_cv%s_v%s/last_learn_rate.txt" %
                        (cross_val, dataver),
                        [np.array(this_rate)], fmt="%.5f"
                        )
                
                # Save model if best validation score
                if (this_valid_loss < best_valid_loss):  
                    best_valid_loss = this_valid_loss
                    best_iter = cum_iterations-1
                    
                    if save_model:
                        np.savez("output_cv%s_v%s/model.npz" % 
                                 (cross_val, dataver),
                                 LL.get_all_param_values(network_0))
                    
                # Reset evaluation reports
                eval_index += 1
                this_train_loss = 0.0
                this_valid_loss = 0.0
                
        aug_index += 1
            
        end_time = time.time()
        print("Computing time for batch %d: %f" % (batch, end_time-start_time))
        
    print("Best validation loss %f after %d epochs" %
          (best_valid_loss, (best_iter*train_minibatch_size//epoch_size)))
    
    del train_set_x, train_set_y, valid_set_x, valid_set_y
    gc.collect()
    
    return None

Example 125

Project: DeepSurv
Source File: deep_surv.py
View license
    def __init__(self, n_in,
    learning_rate, hidden_layers_sizes = None,
    lr_decay = 0.0, momentum = 0.9,
    L2_reg = 0.0, L1_reg = 0.0,
    activation = lasagne.nonlinearities.rectify,
    dropout = None,
    batch_norm = False,
    standardize = False,
    ):
        """
        This class implements and trains a DeepSurv model.

        Parameters:
            n_in: number of input nodes.
            learning_rate: learning rate for training.
            lr_decay: coefficient for Power learning rate decay.
            L2_reg: coefficient for L2 weight decay regularization. Used to help
                prevent the model from overfitting.
            L1_reg: coefficient for L1 weight decay regularization
            momentum: coefficient for momentum. Can be 0 or None to disable.
            hidden_layer_sizes: a list of integers to determine the size of
                each hidden layer.
            activation: a lasagne activation class.
                Default: lasagne.nonlinearities.rectify
            batch_norm: True or False. Include batch normalization layers.
            dropout: if not None or 0, the percentage of dropout to include
                after each hidden layer. Default: None
            standardize: True or False. Include standardization layer after
                input layer.
        """

        self.X = T.fmatrix('x')  # patients covariates
        self.E = T.ivector('e') # the observations vector

        # Default Standardization Values: mean = 0, std = 1
        self.offset = theano.shared(numpy.zeros(shape = n_in, dtype=numpy.float32))
        self.scale = theano.shared(numpy.ones(shape = n_in, dtype=numpy.float32))

        network = lasagne.layers.InputLayer(shape=(None,n_in),
            input_var = self.X)

        if standardize:
            network = lasagne.layers.standardize(network,self.offset,
                                                self.scale,
                                                shared_axes = 0)
        self.standardize = standardize

        # Construct Neural Network
        for n_layer in (hidden_layers_sizes or []):
            if activation == lasagne.nonlinearities.rectify:
                W_init = lasagne.init.GlorotUniform()
            else:
                # TODO: implement other initializations
                W_init = lasagne.init.GlorotUniform()


            network = lasagne.layers.DenseLayer(
                network, num_units = n_layer,
                nonlinearity = activation,
                W = W_init
            )

            if batch_norm:
                network = lasagne.layers.batch_norm(network)

            if not dropout is None:
                network = lasagne.layers.DropoutLayer(network, p = dropout)

        # Combine Linear to output Log Hazard Ratio - same as Faraggi
        network = lasagne.layers.DenseLayer(
            network, num_units = 1,
            nonlinearity = lasagne.nonlinearities.linear,
            W = lasagne.init.GlorotUniform()
        )

        self.network = network
        self.params = lasagne.layers.get_all_params(self.network,
                                                    trainable = True)
        self.hidden_layers = lasagne.layers.get_all_layers(self.network)[1:]

        # Relevant Functions
        self.partial_hazard = T.exp(self.risk(deterministic = True)) # e^h(x)

        # Set Hyper-parameters:
        self.n_in = n_in
        self.learning_rate = learning_rate
        self.lr_decay = lr_decay
        self.L2_reg = L2_reg
        self.L1_reg = L1_reg
        self.momentum = momentum

Example 126

Project: mmdgm
Source File: GPUVAE_Z_X.py
View license
    def gen_xz(self, x, z, n_batch):
        
        x, z = ndict.ordereddicts((x, z))
        
        A = np.ones((1, n_batch)).astype(np.float32)
        
        for i in z: z[i] = z[i].astype(np.float32)
        for i in x: x[i] = x[i].astype(np.float32)
        
        _z = {}

        # If x['x'] was given but not z['z']: generate z ~ q(z|x)
        if x.has_key('x') and not z.has_key('z'):
            '''
            print x['x'].shape
            print x['mean_prior'].shape
            print A.shape
            '''
            q_mean, q_logvar = self.dist_qz['z'](*([x['x'], x['mean_prior']] + [A]))
            q_hidden = self.dist_qz['hidden'](*([x['x'], x['mean_prior']] + [A]))

            _z['mean'] = q_mean
            _z['logvar'] = q_logvar
            _z['hidden'] = q_hidden
            
            # Require epsilon
            if not z.has_key('eps'):
                eps = self.gen_eps(n_batch)['eps']
            
            z['z'] = q_mean + np.exp(0.5 * q_logvar) * eps
            
        elif not z.has_key('z'):
            if self.type_pz in ['gaussian','gaussianmarg']:
                z['z'] = np.random.standard_normal(size=(self.n_z, n_batch)).astype(np.float32)
            elif self.type_pz == 'laplace':
                z['z'] = np.random.laplace(size=(self.n_z, n_batch)).astype(np.float32)
            elif self.type_pz == 'studentt':
                z['z'] = np.random.standard_t(np.dot(np.exp(self.w['logv'].get_value()), A)).astype(np.float32)
            elif self.type_pz == 'mog':
                i = np.random.randint(self.n_mixture)
                loc = np.dot(self.w['mog_mean'+str(i)].get_value(), A)
                scale = np.dot(np.exp(.5*self.w['mog_logvar'+str(i)].get_value()), A)
                z['z'] = np.random.normal(loc=loc, scale=scale).astype(np.float32)
            else:
                raise Exception('Unknown type_pz')
        # Generate from p(x|z)
        
        if self.type_px == 'bernoulli':
            #print 'xz p'
            p = self.dist_px['x'](*([z['z']] + [A]))
            _z['x'] = p
            if not x.has_key('x'):
                #print 'xz ber'
                x['x'] = np.random.binomial(n=1,p=p)
        elif self.type_px == 'bounded01' or self.type_px == 'gaussian':
            x_mean, x_logvar = self.dist_px['x'](*([z['z']] + [A]))
            _z['x'] = x_mean
            if not x.has_key('x'):
                x['x'] = np.random.normal(x_mean, np.exp(x_logvar/2))
                if self.type_px == 'bounded01':
                    x['x'] = np.maximum(np.zeros(x['x'].shape), x['x'])
                    x['x'] = np.minimum(np.ones(x['x'].shape), x['x'])
        elif self.type_px == 'exponential':
            x_mean = self.dist_px['x'](*([z['z']] + [A]))
            _z['x'] = x_mean
            if not x.has_key('x'):
                x['x'] = np.random.exponential(x_mean)
        elif self.type_px == 'mixture':
            x_mean, x_logvar = self.dist_px['x'](*([z['z']] + [A]))
            _z['x'] = x_mean
            if not x.has_key('x'):
                normal_list = np.asarray([1,6,10,14,18])
                exponential_list = np.asarray([0,3,5,9,13,17,21,22,23,24,25,26,27])
                uniform_list = np.asarray([2,4,7,11,15,19])
                threemodal_list = np.asarray([8,12,16,20])
                x = np.zeros((x_mean.shape))
                x[exponential_list,:] = np.random.exponential(x_mean[exponential_list,:])
                x[normal_list,:] = np.random.normal(x_mean[normal_list, :], np.exp(x_logvar[normal_list, :]/2))
                x[uniform_list,:] = np.random.sample(x[uniform_list,:].shape) * 3.5 - 1.75
                #x[threemodal_list,:] = 
                x['x'] = x
        
        else: raise Exception("")
        
        return x, z, _z

Example 127

Project: apogee
Source File: isomodel.py
View license
    def __init__(self,
                 loggmin=None,loggmax=None,
                 imfmodel='lognormalChabrier2001',
                 Z=None,
                 expsfh=False,band='Ks',
                 dontgather=False,
                 basti=False,
                 parsec=True,
                 minage=None,
                 maxage=10.,
                 jkmin=None,
                 stage=None,
                 eta=None):
        """
        NAME:
           __init__
        PURPOSE:
           initialize isocmodel
        INPUT:
           Z= metallicity (if not set, use flat prior in Z over all Z; can be list)
           loggmin= if set, cut logg at this minimum
           loggmax= if set, cut logg at this maximum, if 'rc', then this is the function of teff and z appropriate for the APOGEE RC sample
           imfmodel= (default: 'lognormalChabrier2001') IMF model to use (see isodist.imf code for options)
           band= band to use for M_X (JHK)
           expsfh= if True, use an exponentially-declining star-formation history (can be set to a decay time-scale in Gyr)
           dontgather= if True, don't gather surrounding Zs
           basti= if True, use Basti isochrones (if False, use PARSEC)
           parsec= if True (=default), use PARSEC isochrones, if False, use Padova
           stage= if True, only use this evolutionary stage
           minage= (None) minimum log10 of age/yr
           maxage= (10.) maximum log10 of age/yr
           jkmin= (None) if set, only consider J-Ks greater than this
           eta= (None) mass-loss efficiency parameter
        OUTPUT:
           object
        HISTORY:
           2012-11-07 - Written - Bovy (IAS)
        """
        self._band= band
        self._loggmin= loggmin
        self._loggmax= loggmax
        if isinstance(expsfh,(int,float,numpy.float32,numpy.float64)):
            self._expsfh= expsfh
        elif expsfh:
            self._expsfh= 8.
        else:
            self._expsfh= False
        self._Z= Z
        self._imfmodel= imfmodel
        self._basti= basti
        self._eta= eta
        if isinstance(loggmax,str) and loggmax.lower() == 'rc':
            from apogee.samples.rc import loggteffcut
        #Read isochrones
        if basti:
            zs= numpy.array([0.0001,0.0003,0.0006,0.001,0.002,0.004,0.008,
                             0.01,0.0198,0.03,0.04])
        elif parsec:
            zs= numpy.arange(0.0005,0.06005,0.0005)
        else:
            zs= numpy.arange(0.0005,0.03005,0.0005)
        if Z is None:
            Zs= zs
        elif isinstance(Z,float):
            if basti or dontgather:
                Zs= [Z]
            elif Z < 0.001 or Z > 0.0295:
                Zs= [Z] 
            elif Z < 0.0015 or Z > 0.029:
                Zs= [Z-0.0005,Z,Z+0.0005] #build up statistics
            elif Z < 0.01:
                Zs= [Z-0.001,Z-0.0005,Z,Z+0.0005,Z+0.001] #build up statistics
            else:
                Zs= [Z-0.0005,Z,Z+0.0005] #build up statistics
        if basti:
            p= isodist.BastiIsochrone(Z=Zs,eta=eta)
        else:
            p= isodist.PadovaIsochrone(Z=Zs,parsec=parsec,eta=eta)
        if basti:
            #Force BaSTI to have equal age sampling
            lages= list(numpy.log10(numpy.arange(0.1,1.,0.1))+9.)
            lages.extend(list(numpy.log10(numpy.arange(1.0,10.5,0.5))+9.))
            lages= numpy.array(lages)
        #Get relevant data
        sample= []
        weights= []
        massweights= []
        loggs= []
        teffs= []
        pmasses= []
        plages= []
        pjks= []
        for logage in p.logages():
            if logage > maxage: continue
            if not minage is None and logage < minage: continue
            if basti and numpy.sum((logage == lages)) == 0: continue
            for zz in range(len(Zs)):
                thisiso= p(logage,Zs[zz],asrecarray=True,stage=stage)
                if len(thisiso.M_ini) == 0: continue
                #Calculate int_IMF for this IMF model
                if not imfmodel == 'lognormalChabrier2001': #That would be the default
                    if imfmodel == 'exponentialChabrier2001':
                        int_IMF= isodist.imf.exponentialChabrier2001(thisiso.M_ini,int=True)
                    elif imfmodel == 'kroupa2003':
                        int_IMF= isodist.imf.kroupa2003(thisiso.M_ini,int=True)
                    elif imfmodel == 'chabrier2003':
                        int_IMF= isodist.imf.chabrier2003(thisiso.M_ini,int=True)
                    else:
                        raise IOError("imfmodel option not understood (non-existing model)")
                elif basti:
                    int_IMF= isodist.imf.lognormalChabrier2001(thisiso.M_ini,int=True)
                else:
                    int_IMF= thisiso.int_IMF
                dN= (numpy.roll(int_IMF,-1)-int_IMF)/(int_IMF[-1]-int_IMF[0])/10**(logage-7.)
                dmass= thisiso.M_ini*(numpy.roll(int_IMF,-1)-int_IMF)/numpy.sum((thisiso.M_ini*(numpy.roll(int_IMF,-1)-int_IMF))[:-1])/10**(logage-7.)
                for ii in range(1,len(int_IMF)-1):
                    if basti:
                        JK= 0.996*(thisiso.J[ii]-thisiso.K[ii])+0.00923
                    else:
                        JK= thisiso.J[ii]-thisiso.Ks[ii]
                    if not jkmin is None and JK < jkmin: continue
                    if band.lower() == 'h':
                        if basti:
                            raise NotImplementedError("'H' not implemented for BaSTI yet")
                            J= JK+thisiso.K[ii]-0.044
                            H= J-(0.980*(thisiso.J[ii]-thisiso.H[ii])-0.045)
                        else:
                            H= thisiso.H[ii]
                    elif band.lower() == 'j':
                        if basti:
                            raise NotImplementedError("'J' not implemented for BaSTI yet")
                            J= JK+thisiso.K[ii]-0.044
                        else:
                            H= thisiso.J[ii]
                    elif band.lower() == 'k' or band.lower() == 'ks':
                        if basti:
                            H= thisiso.K[ii]-0.046
                        else:
                            H= thisiso.Ks[ii]
                    if JK < 0.3 \
                            or (isinstance(loggmax,str) and loggmax == 'rc' and (thisiso['logg'][ii] > loggteffcut(10.**thisiso['logTe'][ii],Zs[zz],upper=True))) \
                            or (not isinstance(loggmax,str) and not loggmax is None and thisiso['logg'][ii] > loggmax) \
                            or (not loggmin is None and thisiso['logg'][ii] < loggmin):
                        continue
                    if dN[ii] > 0.: 
                        sample.append([JK,H])
                        loggs.append([thisiso.logg[ii]])
                        teffs.append([10.**thisiso.logTe[ii]])
                        pmasses.append(thisiso.M_ini[ii])
                        plages.append(logage)
                        pjks.append(JK)
                        if basti: #BaSTI is sampled uniformly in age, not logage, but has a finer sampling below 1 Gyr
                            if logage < 9.:
                                if self._expsfh:
                                    weights.append(dN[ii]/5.*numpy.exp((10.**(logage-7.))/self._expsfh/100.)) #e.g., Binney (2010)
                                    massweights.append(dmass[ii]/5.*numpy.exp((10.**(logage-7.))/self._expsfh/100.)) #e.g., Binney (2010)
                                else:
                                    weights.append(dN[ii]/5.)
                                    massweights.append(dmass[ii]/5.)
                            else:
                                if self._expsfh:
                                    weights.append(dN[ii]*numpy.exp((10.**(logage-7.))/self._expsfh/100.)) #e.g., Binney (2010)
                                    massweights.append(dmass[ii]*numpy.exp((10.**(logage-7.))/self._expsfh/100.)) #e.g., Binney (2010)
                                else:
                                    weights.append(dN[ii])
                                    massweights.append(dmass[ii])
                        else:
                            if self._expsfh:
                                weights.append(dN[ii]*10**(logage-7.)*numpy.exp((10.**(logage-7.))/self._expsfh/100.)) #e.g., Binney (2010)
                                massweights.append(dmass[ii]*10**(logage-7.)*numpy.exp((10.**(logage-7.))/self._expsfh/100.)) #e.g., Binney (2010)
                            else:
                                weights.append(dN[ii]*10**(logage-7.))
                                massweights.append(dmass[ii]*10**(logage-7.))
                    else: 
                        continue #no use in continuing here   
        #Form array
        sample= numpy.array(sample)
        loggs= numpy.array(loggs)
        teffs= numpy.array(teffs)
        pmasses= numpy.array(pmasses)
        plages= numpy.array(plages)-9.
        pjks= numpy.array(pjks)
        weights= numpy.array(weights)
        massweights= numpy.array(massweights)
        #Cut out low weights
        if False:
            indx= (weights > 10.**-5.*numpy.sum(weights))
        else:
            indx= numpy.ones(len(weights),dtype='bool')
        self._sample= sample[indx,:]
        self._weights= weights[indx]
        self._massweights= massweights[indx]
        self._loggs= loggs[indx]
        self._teffs= teffs[indx]
        self._masses= pmasses[indx]
        self._lages= plages[indx]
        self._jks= pjks[indx]
        #Setup KDE
        self._kde= dens_kde.densKDE(self._sample,w=self._weights,
                                    h=2.*self._sample.shape[0]**(-1./5.),#h='scott',
                                    kernel='biweight',
                                    variable=True,variablenitt=3,
                                    variableexp=0.5)
        return None

Example 128

Project: scipy
Source File: pilutil.py
View license
def toimage(arr, high=255, low=0, cmin=None, cmax=None, pal=None,
            mode=None, channel_axis=None):
    """Takes a numpy array and returns a PIL image.

    The mode of the PIL image depends on the array shape and the `pal` and
    `mode` keywords.

    For 2-D arrays, if `pal` is a valid (N,3) byte-array giving the RGB values
    (from 0 to 255) then ``mode='P'``, otherwise ``mode='L'``, unless mode
    is given as 'F' or 'I' in which case a float and/or integer array is made.

    Notes
    -----
    For 3-D arrays, the `channel_axis` argument tells which dimension of the
    array holds the channel data.

    For 3-D arrays if one of the dimensions is 3, the mode is 'RGB'
    by default or 'YCbCr' if selected.

    The numpy array must be either 2 dimensional or 3 dimensional.

    """
    data = asarray(arr)
    if iscomplexobj(data):
        raise ValueError("Cannot convert a complex-valued array.")
    shape = list(data.shape)
    valid = len(shape) == 2 or ((len(shape) == 3) and
                                ((3 in shape) or (4 in shape)))
    if not valid:
        raise ValueError("'arr' does not have a suitable array shape for "
                         "any mode.")
    if len(shape) == 2:
        shape = (shape[1], shape[0])  # columns show up first
        if mode == 'F':
            data32 = data.astype(numpy.float32)
            image = Image.frombytes(mode, shape, data32.tostring())
            return image
        if mode in [None, 'L', 'P']:
            bytedata = bytescale(data, high=high, low=low,
                                 cmin=cmin, cmax=cmax)
            image = Image.frombytes('L', shape, bytedata.tostring())
            if pal is not None:
                image.putpalette(asarray(pal, dtype=uint8).tostring())
                # Becomes a mode='P' automagically.
            elif mode == 'P':  # default gray-scale
                pal = (arange(0, 256, 1, dtype=uint8)[:, newaxis] *
                       ones((3,), dtype=uint8)[newaxis, :])
                image.putpalette(asarray(pal, dtype=uint8).tostring())
            return image
        if mode == '1':  # high input gives threshold for 1
            bytedata = (data > high)
            image = Image.frombytes('1', shape, bytedata.tostring())
            return image
        if cmin is None:
            cmin = amin(ravel(data))
        if cmax is None:
            cmax = amax(ravel(data))
        data = (data*1.0 - cmin)*(high - low)/(cmax - cmin) + low
        if mode == 'I':
            data32 = data.astype(numpy.uint32)
            image = Image.frombytes(mode, shape, data32.tostring())
        else:
            raise ValueError(_errstr)
        return image

    # if here then 3-d array with a 3 or a 4 in the shape length.
    # Check for 3 in datacube shape --- 'RGB' or 'YCbCr'
    if channel_axis is None:
        if (3 in shape):
            ca = numpy.flatnonzero(asarray(shape) == 3)[0]
        else:
            ca = numpy.flatnonzero(asarray(shape) == 4)
            if len(ca):
                ca = ca[0]
            else:
                raise ValueError("Could not find channel dimension.")
    else:
        ca = channel_axis

    numch = shape[ca]
    if numch not in [3, 4]:
        raise ValueError("Channel axis dimension is not valid.")

    bytedata = bytescale(data, high=high, low=low, cmin=cmin, cmax=cmax)
    if ca == 2:
        strdata = bytedata.tostring()
        shape = (shape[1], shape[0])
    elif ca == 1:
        strdata = transpose(bytedata, (0, 2, 1)).tostring()
        shape = (shape[2], shape[0])
    elif ca == 0:
        strdata = transpose(bytedata, (1, 2, 0)).tostring()
        shape = (shape[2], shape[1])
    if mode is None:
        if numch == 3:
            mode = 'RGB'
        else:
            mode = 'RGBA'

    if mode not in ['RGB', 'RGBA', 'YCbCr', 'CMYK']:
        raise ValueError(_errstr)

    if mode in ['RGB', 'YCbCr']:
        if numch != 3:
            raise ValueError("Invalid array shape for mode.")
    if mode in ['RGBA', 'CMYK']:
        if numch != 4:
            raise ValueError("Invalid array shape for mode.")

    # Here we know data and mode is correct
    image = Image.frombytes(mode, shape, strdata)
    return image

Example 129

Project: kwiklib
Source File: kwik.py
View license
def add_spikes(fd, channel_group_id=None,
                time_samples=None, time_fractional=0,
                recording=0, cluster=0, cluster_original=0,
                features_masks=None, features=None, masks=None,
                waveforms_raw=None, waveforms_filtered=None,
                fill_empty=True):
    """fd is returned by `open_files`: it is a dict {type: tb_file_handle}."""
    if channel_group_id is None:
        channel_group_id = '0'
    kwik = fd.get('kwik', None)
    kwx = fd.get('kwx', None)
    # The KWIK needs to be there.
    assert kwik is not None
    # The channel group id containing the new cluster group must be specified.
    assert channel_group_id is not None

    spikes = kwik.root.channel_groups.__getattr__(channel_group_id).spikes
        
    time_samples = ensure_vector(time_samples)
    nspikes = len(time_samples)
    
    ds_features_masks = kwx.root.channel_groups.__getattr__(channel_group_id).features_masks
    ds_waveforms_raw = kwx.root.channel_groups.__getattr__(channel_group_id).waveforms_raw
    ds_waveforms_filtered = kwx.root.channel_groups.__getattr__(channel_group_id).waveforms_filtered
        
    nfeatures = ds_features_masks.shape[1]
    
    if features_masks is None:
        # Deal with masks only if they're supported in the file.
        if ds_features_masks.ndim == 2:
            features_masks = features
        # Otherwise, deal with masks.
        # Fill features masks only if features OR masks is specified.
        elif not(features is None and masks is None):
            # Unless features and masks are both None,
            # make sure features AND masks are filled.
            if features is None:
                features = np.zeros((nspikes, nfeatures), dtype=np.float32)
            if masks is None:
                masks = np.zeros((features.shape[0], nfeatures), dtype=np.float32)
                
            # Ensure both have 2 dimensions.
            if features.ndim == 1:
                features = np.expand_dims(features, axis=0)
            if masks.ndim == 1:
                masks = np.expand_dims(masks, axis=0)
                
            # Tile the masks if needed: same mask value on each channel.
            if masks.shape[1] < features.shape[1]:
                nfeatures_per_channel = features.shape[1] // masks.shape[1]
                masks = np.repeat(masks, nfeatures_per_channel, axis=1)
                
            # Concatenate features and masks.
            features_masks = np.dstack((features, masks))
    time_fractional = ensure_vector(time_fractional, size=nspikes)
    recording = ensure_vector(recording, size=nspikes)
    cluster = ensure_vector(cluster, size=nspikes)
    cluster_original = ensure_vector(cluster_original, size=nspikes)
    
    if fill_empty and features_masks is None:
        features_masks = empty_row(ds_features_masks, nrows=nspikes)
    if features_masks is not None and features_masks.ndim < ds_features_masks.ndim:
        features_masks = np.expand_dims(features_masks, axis=0)
    
    if fill_empty and waveforms_raw is None:
        waveforms_raw = empty_row(ds_waveforms_raw, nrows=nspikes)
    if waveforms_raw is not None and waveforms_raw.ndim < 3:
        waveforms_raw = np.expand_dims(waveforms_raw, axis=0)
        
    if fill_empty and waveforms_filtered is None:
        waveforms_filtered = empty_row(ds_waveforms_filtered, nrows=nspikes)
    if waveforms_filtered is not None and waveforms_filtered.ndim < 3:
        waveforms_filtered = np.expand_dims(waveforms_filtered, axis=0)
        
    # Make sure we add the correct number of rows to every object.
    assert len(time_samples) == nspikes
    assert len(time_fractional) == nspikes
    assert len(recording) == nspikes
    assert len(cluster) == nspikes
    assert len(cluster_original) == nspikes
    
    if features_masks is not None:
        assert features_masks.shape[0] == nspikes
    
    if waveforms_raw is not None:
        assert waveforms_raw.shape[0] == nspikes
    if waveforms_filtered is not None:
        assert waveforms_filtered.shape[0] == nspikes
        
    spikes.time_samples.append(time_samples)
    spikes.time_fractional.append(time_fractional)
    spikes.recording.append(recording)
    spikes.clusters.main.append(cluster)
    spikes.clusters.original.append(cluster_original)
    
    if features_masks is not None:
        ds_features_masks.append(features_masks)
    
    if waveforms_raw is not None:
        ds_waveforms_raw.append(waveforms_raw.astype(np.int16))
    if waveforms_filtered is not None:
        ds_waveforms_filtered.append(waveforms_filtered.astype(np.int16))
        

Example 130

Project: tensorlayer
Source File: tutorial_cifar10.py
View license
def main_test_cnn_naive():
    """Without any distorting, whitening and cropping for training data.
    This method work well for MNIST, but not CIFAR-10.
    """
    model_file_name = "model_cifar10_naive.ckpt"
    resume = False # load model, resume from previous checkpoint?

    X_train, y_train, X_test, y_test = tl.files.load_cifar10_dataset(
                                        shape=(-1, 32, 32, 3), plotable=False)

    X_train = np.asarray(X_train, dtype=np.float32)
    y_train = np.asarray(y_train, dtype=np.int32)
    X_test = np.asarray(X_test, dtype=np.float32)
    y_test = np.asarray(y_test, dtype=np.int32)

    print('X_train.shape', X_train.shape)   # (50000, 32, 32, 3)
    print('y_train.shape', y_train.shape)   # (50000,)
    print('X_test.shape', X_test.shape)     # (10000, 32, 32, 3)
    print('y_test.shape', y_test.shape)     # (10000,)
    print('X %s   y %s' % (X_test.dtype, y_test.dtype))

    sess = tf.InteractiveSession()

    batch_size = 128

    x = tf.placeholder(tf.float32, shape=[batch_size, 32, 32, 3])
                                # [batch_size, height, width, channels]
    y_ = tf.placeholder(tf.int32, shape=[batch_size,])

    network = tl.layers.InputLayer(x, name='input_layer')
    network = tl.layers.Conv2dLayer(network,
                        act = tf.nn.relu,
                        shape = [5, 5, 3, 64],  # 64 features for each 5x5x3 patch
                        strides=[1, 1, 1, 1],
                        padding='SAME',
                        name ='cnn_layer1')     # output: (?, 32, 32, 64)
    network = tl.layers.PoolLayer(network,
                        ksize=[1, 3, 3, 1],
                        strides=[1, 2, 2, 1],
                        padding='SAME',
                        pool = tf.nn.max_pool,
                        name ='pool_layer1')   # output: (?, 16, 16, 64)
    # local response normalization, you can also try batch normalization.
    # References: ImageNet Classification with Deep Convolutional Neural Networks
    #   it increases the accuracy but consume more time.
    #     https://www.tensorflow.org/versions/master/api_docs/python/nn.html#local_response_normalization
    network.outputs = tf.nn.lrn(network.outputs, 4, bias=1.0, alpha=0.001 / 9.0,
                                                        beta=0.75, name='norm1')
    network = tl.layers.Conv2dLayer(network,
                        act = tf.nn.relu,
                        shape = [5, 5, 64, 64], # 64 features for each 5x5 patch
                        strides=[1, 1, 1, 1],
                        padding='SAME',
                        name ='cnn_layer2')     # output: (?, 16, 16, 64)
    # Another local response normalization.
    network.outputs = tf.nn.lrn(network.outputs, 4, bias=1.0, alpha=0.001 / 9.0,
                                                        beta=0.75, name='norm2')
    network = tl.layers.PoolLayer(network,
                        ksize=[1, 3, 3, 1],
                        strides=[1, 2, 2, 1],
                        padding='SAME',
                        pool = tf.nn.max_pool,
                        name ='pool_layer2')   # output: (?, 8, 8, 64)

    network = tl.layers.FlattenLayer(network, name='flatten_layer')
                                                            # output: (?, 4096)
    network = tl.layers.DenseLayer(network, n_units=384,
                            act = tf.nn.relu, name='relu1') # output: (?, 384)
    network = tl.layers.DenseLayer(network, n_units=192,
                            act = tf.nn.relu, name='relu2') # output: (?, 192)
    network = tl.layers.DenseLayer(network, n_units=10,
                            act = tf.identity,
                            name='output_layer')            # output: (?, 10)

    y = network.outputs

    ce = tf.reduce_mean(tf.nn.sparse_softmax_cross_entropy_with_logits(y, y_))
    cost = ce

    correct_prediction = tf.equal(tf.cast(tf.argmax(y, 1), tf.int32), y_)
    acc = tf.reduce_mean(tf.cast(correct_prediction, tf.float32))

    # train
    n_epoch = 10000
    learning_rate = 0.0001
    print_freq = 1

    train_params = network.all_params
    train_op = tf.train.AdamOptimizer(learning_rate, beta1=0.9, beta2=0.999,
        epsilon=1e-08, use_locking=False).minimize(cost, var_list=train_params)

    sess.run(tf.initialize_all_variables())
    if resume:
        print("Load existing model " + "!"*10)
        saver = tf.train.Saver()
        saver.restore(sess, model_file_name)

    network.print_params()
    network.print_layers()

    print('   learning_rate: %f' % learning_rate)
    print('   batch_size: %d' % batch_size)

    for epoch in range(n_epoch):
        start_time = time.time()
        for X_train_a, y_train_a in tl.iterate.minibatches(
                                X_train, y_train, batch_size, shuffle=True):
            feed_dict = {x: X_train_a, y_: y_train_a}
            feed_dict.update( network.all_drop )        # enable all dropout/dropconnect/denoising layers
            sess.run(train_op, feed_dict=feed_dict)

        if epoch + 1 == 1 or (epoch + 1) % print_freq == 0:
            print("Epoch %d of %d took %fs" % (epoch + 1, n_epoch, time.time() - start_time))
            train_loss, train_acc, n_batch = 0, 0, 0
            for X_train_a, y_train_a in tl.iterate.minibatches(
                                    X_train, y_train, batch_size, shuffle=True):
                dp_dict = tl.utils.dict_to_one( network.all_drop )    # disable all dropout/dropconnect/denoising layers
                feed_dict = {x: X_train_a, y_: y_train_a}
                feed_dict.update(dp_dict)
                err, ac = sess.run([cost, acc], feed_dict=feed_dict)
                assert not np.isnan(err), 'Model diverged with cost = NaN'
                train_loss += err; train_acc += ac; n_batch += 1
            print("   train loss: %f" % (train_loss/ n_batch))
            print("   train acc: %f" % (train_acc/ n_batch))
            test_loss, test_acc, n_batch = 0, 0, 0
            for X_test_a, y_test_a in tl.iterate.minibatches(
                                    X_test, y_test, batch_size, shuffle=True):
                dp_dict = tl.utils.dict_to_one( network.all_drop )    # disable all dropout/dropconnect/denoising layers
                feed_dict = {x: X_test_a, y_: y_test_a}
                feed_dict.update(dp_dict)
                err, ac = sess.run([cost, acc], feed_dict=feed_dict)
                test_loss += err; test_acc += ac; n_batch += 1
            print("   test loss: %f" % (test_loss/ n_batch))
            print("   test acc: %f" % (test_acc/ n_batch))
            try:
                pass
                # tl.visualize.CNN2d(network.all_params[0].eval(), second=10, saveable=True, name='cnn1_'+str(epoch+1), fig_idx=2012)
            except:
                raise Exception("# You should change visualize.CNN(), \
                if you want to save the feature images for different dataset")

        if (epoch + 1) % 1 == 0:
            print("Save model " + "!"*10);
            saver = tf.train.Saver()
            save_path = saver.save(sess, model_file_name)

Example 131

Project: opendrift
Source File: opendrift3D.py
View license
    def vertical_mixing(self):
        """Mix particles vertically according to eddy diffusivity and buoyancy

            Buoyancy is expressed as terminal velocity, which is the
            steady-state vertical velocity due to positive or negative
            buoyant behaviour. It is usually a function of particle density,
            diameter, and shape.

            Vertical particle displacemend du to turbulent mixing is
            calculated using the "binned random walk scheme" (Thygessen and
            Aadlandsvik, 2007).
            The formulation of this scheme is copied from LADIM (IMR).
        """

        if self.config['processes']['turbulentmixing'] is False:
            logging.debug('Turbulent mixing deactivated.')
            return

        from opendrift.models import eddydiffusivity

        dz = self.config['turbulentmixing']['verticalresolution']
        dz = np.float32(dz)  # Convert to avoid error for older numpy
        dt_mix = self.config['turbulentmixing']['timestep']

        # minimum height/maximum depth for each particle
        Zmin = -1.*self.environment.sea_floor_depth_below_sea_level

        # place particle in center of bin
        surface = self.elements.z == 0
        self.elements.z[~surface] = np.round(self.elements.z[~surface]/dz)*dz

        #avoid that elements are below bottom
        bottom = np.where(self.elements.z < Zmin)
        self.elements.z[bottom] = np.round(Zmin/dz)*dz + dz/2.

        # Eventual model specific preparions
        self.prepare_vertical_mixing()

        # get profile of eddy diffusivity
        # get vertical eddy diffusivity from environment or specific model
        if (self.config['turbulentmixing']['diffusivitymodel'] ==
                'environment'):
            if 'ocean_vertical_diffusivity' in self.environment_profiles:
                Kprofiles = self.environment_profiles[
                    'ocean_vertical_diffusivity']
                logging.debug('use diffusivity from ocean model')
            else:
                # NB: using constant diffusivity, and value from first
                # element only - this should be checked/improved!
                Kprofiles = \
                    self.environment.ocean_vertical_diffusivity[0] * \
                    np.ones((len(self.environment_profiles['z']),
                             self.num_elements_active()))
                logging.debug('use constant diffusivity')
        else:
            logging.debug('use functional expression for diffusivity')
            Kprofiles = getattr(
                eddydiffusivity,
                self.config['turbulentmixing']['diffusivitymodel'])(self)

        # get profiles of salinity and temperature
        # (to save interpolation time in the inner loop)
        if 'TSprofiles' in self.config['turbulentmixing']:
            if self.config['turbulentmixing']['TSprofiles'] is True:
                Sprofiles = self.environment_profiles['sea_water_salinity']
                Tprofiles = \
                    self.environment_profiles['sea_water_temperature']

        # prepare vertical interpolation coordinates
        z_i = range(Kprofiles.shape[0])
        z_index = interp1d(-self.environment_profiles['z'],
                           z_i, bounds_error=False)

        # internal loop for fast time step of vertical mixing model
        # binned random walk needs faster time step compared
        # to horizontal advection
        ntimes_mix = int(self.time_step.total_seconds()/dt_mix)
        logging.debug('Vertical mixing module:')
        logging.debug('turbulent diffusion with binned random walk scheme')
        logging.debug('using ' + str(ntimes_mix) + ' fast time steps of dt=' +
                      str(dt_mix) + 's')
        for i in range(0, ntimes_mix):
            #remember which particles belong to the exact surface
            surface = self.elements.z == 0

            # update terminal velocity according to environmental variables
            if 'TSprofiles' in self.config['turbulentmixing']:
                if self.config['turbulentmixing']['TSprofiles'] is True:
                    self.update_terminal_velocity(Tprofiles=Tprofiles,
                                                  Sprofiles=Sprofiles,
                                                  z_index=z_index)
                else:
                    self.update_terminal_velocity()
            else:
                # this is faster, but ignores density gradients in
                # water column for the inner loop
                self.update_terminal_velocity()

            w = self.elements.terminal_velocity

            # diffusivity K at depth z
            zi = z_index(-self.elements.z)
            upper = np.maximum(np.floor(zi).astype(np.int), 0)
            lower = np.minimum(upper+1, Kprofiles.shape[0]-1)
            weight_upper = 1 - (zi - upper)
            K1 = Kprofiles[upper, range(Kprofiles.shape[1])] * \
                weight_upper + \
                Kprofiles[lower, range(Kprofiles.shape[1])] * \
                (1-weight_upper)

            # K at depth z-dz ; gradient of K is required for correct
            # solution with random walk scheme
            zi = z_index(-(self.elements.z-dz))
            upper = np.maximum(np.floor(zi).astype(np.int), 0)
            lower = np.minimum(upper+1, Kprofiles.shape[0]-1)
            weight_upper = 1 - (zi - upper)
            K2 = Kprofiles[upper, range(Kprofiles.shape[1])] * \
                weight_upper + \
                Kprofiles[lower, range(Kprofiles.shape[1])] * \
                (1-weight_upper)

            # calculate rise/sink probability dependent on K and w
            p = dt_mix * (2.0*K1 + dz*w)/(2.0*dz*dz)  # probability to rise
            q = dt_mix * (2.0*K2 - dz*w)/(2.0*dz*dz)  # probability to sink

            # check if probabilities are reasonable or wrong; which can happen if K is very high (K>0.1)
            wrong = p+q > 1.00002
            if wrong.sum() > 0:
                logging.info('WARNING! '+str(wrong.sum())+' elements have p+q>1; you might need a smaller mixing time step')
                # fixing p and q by scaling them to assure p+q<1:
                norm = p+q
                p[wrong] = p[wrong]/norm[wrong] 
                q[wrong] = q[wrong]/norm[wrong]

            # use probabilities to mix some particles up or down
            RandKick = np.random.random(self.num_elements_active())           
            up = np.where(RandKick < p)
            down = np.where(RandKick > 1.0 - q)           
            self.elements.z[up] = self.elements.z[up] + dz # move to layer above
            self.elements.z[down] = self.elements.z[down] - dz # move to layer underneath

            # put the particles that belong to the surface slick (if present) back to the surface
            self.elements.z[surface] = 0.

            #avoid that elements are below bottom
            bottom = np.where(self.elements.z < Zmin)
            self.elements.z[bottom] = np.round(Zmin/dz)*dz + dz/2.

            # Call surface interaction:
            # reflection at surface or formation of slick and wave mixing if implemented for this class
            self.surface_interaction(dt_mix)

Example 132

Project: visvis
Source File: baseTexture.py
View license
    def SetMap(self, *args):
        """ SetMap(*args)
        
        Set the colormap data. This method accepts several arguments:
        
        A list/tuple of tuples where each tuple represents a RGB or RGBA color.
        
        A dict with keys 'red', 'green', 'blue', 'alpha' (or only the first
        letter). Each dict should contain a list of 2-element tuples that
        specify index and color value. Indices should be between 0 and 1.
        
        A numpy array specifying the RGB or RGBA tuples.
        
        """
        
        # one argument given?
        if len(args)==1:
            args = args[0]
        
        # store
        self._current = args
        
        # init
        data = None
        
        # parse input
        
        if isinstance(args, dict):
            # DICT
            
            # Allow several color names
            for key in list(args.keys()):
                if key.lower() in ['r', 'red']:
                    args['r'] = args[key]
                elif key.lower() in ['g', 'green']:
                    args['g'] = args[key]
                if key.lower() in ['b', 'blue']:
                    args['b'] = args[key]
                if key.lower() in ['a', 'alpha']:
                    args['a'] = args[key]
            # Init data, alpha 1
            data2 = np.zeros((256,4),np.float32)
            data2[:,3] = 1.0
            # For each channel ...
            for i in range(4):
                channel = 'rgba'[i]
                if not channel in args:
                    continue
                # Get value list and check
                values = args[channel]
                if not hasattr(values,'__len__'):
                    raise ValueError('Invalid colormap.')
                # Init interpolation
                data = np.zeros((len(values),), dtype=np.float32)
                x = np.linspace(0.0, 1.0, 256)
                xp = np.zeros((len(values),), dtype=np.float32)
                # Insert values
                count = -1
                for el in values:
                    count += 1                    
                    if not hasattr(el,'__len__') or len(el) != 2:
                        raise ValueError('Colormap dict entries must have 2 elements.')
                    xp[count] = el[0]
                    data[count] = el[1]
                # Interpolate
                data2[:,i] = np.interp(x, xp, data)
            # Set
            data = data2
        
        elif isinstance(args, (tuple, list)):
            # LIST
            
            data = np.zeros((len(args),4), dtype=np.float32)
            data[:,3] = 1.0 # init alpha to be all ones
            count = -1
            for el in args:
                count += 1
                if not hasattr(el,'__len__') or len(el) not in [3,4]:
                    raise ValueError('Colormap entries must have 3 or 4 elements.')
                elif len(el)==3:
                    data[count,:] = el[0], el[1], el[2], 1.0
                elif len(el)==4:
                    data[count,:] = el[0], el[1], el[2], el[3]
        
        elif isinstance(args, np.ndarray):
            # ARRAY
            
            if args.ndim != 2 or args.shape[1] not in [3,4]:
                raise ValueError('Colormap entries must have 3 or 4 elements.')
            elif args.shape[1]==3:
                data = np.ones((args.shape[0], 4), dtype=np.float32)
                data[:, 0:3] = args
            elif args.shape[1]==4:
                data = args
            else:
                raise ValueError("Invalid argument to set colormap.")
        
        # Apply interpolation (if required)
        if data is not None:   
            if data.shape[0] == 256 and data.dtype == np.float32:
                data2 = data
            else:
                # interpolate first            
                x = np.linspace(0.0, 1.0, 256)
                xp = np.linspace(0.0, 1.0, data.shape[0])            
                data2 = np.zeros((256,4),np.float32)
                for i in range(4):
                    data2[:,i] = np.interp(x, xp, data[:,i])
            # store texture
            #self._data = data2
            self.SetData(data2)

Example 133

Project: instaseis
Source File: util.py
View license
def _validate_and_write_waveforms(st, callback, starttime, endtime, scale,
                                  source, receiver, db, label, format):
    if not label:
        label = ""
    else:
        label += "_"

    for tr in st:
        # Half the filesize but definitely sufficiently accurate.
        tr.data = np.require(tr.data, dtype=np.float32)

    if scale != 1.0:
        for tr in st:
            tr.data *= scale

    # Sanity checks. Raise internal server errors in case something fails.
    # This should not happen and should have been caught before.
    if endtime > st[0].stats.endtime:
        msg = ("Endtime larger than the extracted endtime: endtime=%s, "
               "largest db endtime=%s" % (
                _format_utc_datetime(endtime),
                _format_utc_datetime(st[0].stats.endtime)))
        callback((tornado.web.HTTPError(500, log_message=msg, reason=msg),
                  None))
        return
    if starttime < st[0].stats.starttime - 3600.0:
        msg = ("Starttime more than one hour before the starttime of the "
               "seismograms.")
        callback((tornado.web.HTTPError(500, log_message=msg, reason=msg),
                  None))
        return

    if isinstance(source, FiniteSource):
        mu = None
    else:
        mu = st[0].stats.instaseis.mu

    # Trim, potentially pad with zeroes.
    st.trim(starttime, endtime, pad=True, fill_value=0.0, nearest_sample=False)

    # Checked in another function and just a sanity check.
    assert format in ("miniseed", "saczip")

    if format == "miniseed":
        with io.BytesIO() as fh:
            st.write(fh, format="mseed")
            fh.seek(0, 0)
            binary_data = fh.read()
        callback((binary_data, mu))
    # Write a number of SAC files into an archive.
    elif format == "saczip":
        byte_strings = []
        for tr in st:
            # Write SAC headers.
            tr.stats.sac = obspy.core.AttribDict()
            # Write WGS84 coordinates to the SAC files.
            tr.stats.sac.stla = geocentric_to_elliptic_latitude(
                receiver.latitude)
            tr.stats.sac.stlo = receiver.longitude
            tr.stats.sac.stdp = receiver.depth_in_m
            tr.stats.sac.stel = 0.0
            if isinstance(source, FiniteSource):
                tr.stats.sac.evla = geocentric_to_elliptic_latitude(
                    source.hypocenter_latitude)
                tr.stats.sac.evlo = source.hypocenter_longitude
                tr.stats.sac.evdp = source.hypocenter_depth_in_m
                # Force source has no magnitude.
                if not isinstance(source, ForceSource):
                    tr.stats.sac.mag = source.moment_magnitude
                src_lat = source.hypocenter_latitude
                src_lng = source.hypocenter_longitude
            else:
                tr.stats.sac.evla = geocentric_to_elliptic_latitude(
                    source.latitude)
                tr.stats.sac.evlo = source.longitude
                tr.stats.sac.evdp = source.depth_in_m
                # Force source has no magnitude.
                if not isinstance(source, ForceSource):
                    tr.stats.sac.mag = source.moment_magnitude
                src_lat = source.latitude
                src_lng = source.longitude
            # Thats what SPECFEM uses for a moment magnitude....
            tr.stats.sac.imagtyp = 55
            # The event origin time relative to the reference which I'll
            # just assume to be the starttime here?
            tr.stats.sac.o = source.origin_time - starttime

            # Sac coordinates are elliptical thus it only makes sense to
            # have elliptical distances.
            dist_in_m, az, baz = gps2dist_azimuth(
                lat1=tr.stats.sac.evla,
                lon1=tr.stats.sac.evlo,
                lat2=tr.stats.sac.stla,
                lon2=tr.stats.sac.stlo)

            tr.stats.sac.dist = dist_in_m / 1000.0
            tr.stats.sac.az = az
            tr.stats.sac.baz = baz

            # XXX: Is this correct? Maybe better use some function in
            # geographiclib?
            tr.stats.sac.gcarc = locations2degrees(
                lat1=src_lat,
                long1=src_lng,
                lat2=receiver.latitude,
                long2=receiver.longitude)

            # Set two more headers. See #45.
            tr.stats.sac.lpspol = 1
            tr.stats.sac.lcalda = 0

            # Some provenance.
            tr.stats.sac.kuser0 = "InstSeis"
            tr.stats.sac.kuser1 = db.info.velocity_model[:8]
            tr.stats.sac.user0 = scale
            # Prefix version numbers to identify them at a glance.
            tr.stats.sac.kt7 = "A" + db.info.axisem_version[:7]
            tr.stats.sac.kt8 = "I" + __version__[:7]

            # Times have to be set by hand.
            t, _ = utcdatetime_to_sac_nztimes(tr.stats.starttime)
            for key, value in t.items():
                tr.stats.sac[key] = value

            with io.BytesIO() as temp:
                tr.write(temp, format="sac")
                temp.seek(0, 0)
                filename = "%s%s.sac" % (label, tr.id)
                byte_strings.append((filename, temp.read()))
        callback((byte_strings, mu))

Example 134

Project: tensorlayer
Source File: tutorial_mnist.py
View license
def main_test_layers(model='relu'):
    X_train, y_train, X_val, y_val, X_test, y_test = \
                                    tl.files.load_mnist_dataset(shape=(-1,784))

    X_train = np.asarray(X_train, dtype=np.float32)
    y_train = np.asarray(y_train, dtype=np.int32)
    X_val = np.asarray(X_val, dtype=np.float32)
    y_val = np.asarray(y_val, dtype=np.int32)
    X_test = np.asarray(X_test, dtype=np.float32)
    y_test = np.asarray(y_test, dtype=np.int32)

    print('X_train.shape', X_train.shape)
    print('y_train.shape', y_train.shape)
    print('X_val.shape', X_val.shape)
    print('y_val.shape', y_val.shape)
    print('X_test.shape', X_test.shape)
    print('y_test.shape', y_test.shape)
    print('X %s   y %s' % (X_test.dtype, y_test.dtype))

    sess = tf.InteractiveSession()

    # placeholder
    x = tf.placeholder(tf.float32, shape=[None, 784], name='x')
    y_ = tf.placeholder(tf.int32, shape=[None, ], name='y_')

    # Note: the softmax is implemented internally in tl.cost.cross_entropy(y, y_)
    # to speed up computation, so we use identity in the last layer.
    # see tf.nn.sparse_softmax_cross_entropy_with_logits()
    if model == 'relu':
        network = tl.layers.InputLayer(x, name='input_layer')
        network = tl.layers.DropoutLayer(network, keep=0.8, name='drop1')
        network = tl.layers.DenseLayer(network, n_units=800,
                                        act = tf.nn.relu, name='relu1')
        network = tl.layers.DropoutLayer(network, keep=0.5, name='drop2')
        network = tl.layers.DenseLayer(network, n_units=800,
                                        act = tf.nn.relu, name='relu2')
        network = tl.layers.DropoutLayer(network, keep=0.5, name='drop3')
        network = tl.layers.DenseLayer(network, n_units=10,
                                        act = tf.identity,
                                        name='output_layer')
    elif model == 'dropconnect':
        network = tl.layers.InputLayer(x, name='input_layer')
        network = tl.layers.DropconnectDenseLayer(network, keep = 0.8,
                                                n_units=800, act = tf.nn.relu,
                                                name='dropconnect_relu1')
        network = tl.layers.DropconnectDenseLayer(network, keep = 0.5,
                                                n_units=800, act = tf.nn.relu,
                                                name='dropconnect_relu2')
        network = tl.layers.DropconnectDenseLayer(network, keep = 0.5,
                                                n_units=10,
                                                act = tf.identity,
                                                name='output_layer')

    # To print all attributes of a Layer.
    # attrs = vars(network)
    # print(', '.join("%s: %s\n" % item for item in attrs.items()))
    #
    # print(network.all_drop)     # {'drop1': 0.8, 'drop2': 0.5, 'drop3': 0.5}
    # print(drop1, drop2, drop3)  # Tensor("Placeholder_2:0", dtype=float32) Tensor("Placeholder_3:0", dtype=float32) Tensor("Placeholder_4:0", dtype=float32)
    # exit()

    y = network.outputs
    y_op = tf.argmax(tf.nn.softmax(y), 1)
    cost = tf.reduce_mean(tf.nn.sparse_softmax_cross_entropy_with_logits(y, y_))
    # Alternatively, you can use TensorLayer's function to compute cost:
    # cost = tl.cost.cross_entropy(y, y_)

    # You can add more penalty to the cost function as follow.
    # cost = cost + tl.cost.maxnorm_regularizer(1.0)(network.all_params[0]) + tl.cost.maxnorm_regularizer(1.0)(network.all_params[2])
    # cost = cost + tl.cost.lo_regularizer(0.0001)(network.all_params[0]) + tl.cost.lo_regularizer(0.0001)(network.all_params[2])
    # cost = cost + tl.cost.maxnorm_o_regularizer(0.001)(network.all_params[0]) + tl.cost.maxnorm_o_regularizer(0.001)(network.all_params[2])

    params = network.all_params
    # train
    n_epoch = 100
    batch_size = 128
    learning_rate = 0.0001
    print_freq = 5
    train_op = tf.train.AdamOptimizer(learning_rate, beta1=0.9, beta2=0.999,
                                epsilon=1e-08, use_locking=False).minimize(cost)

    sess.run(tf.initialize_all_variables()) # initialize all variables

    network.print_params()
    network.print_layers()

    print('   learning_rate: %f' % learning_rate)
    print('   batch_size: %d' % batch_size)

    for epoch in range(n_epoch):
        start_time = time.time()
        for X_train_a, y_train_a in tl.iterate.minibatches(X_train, y_train,
                                                    batch_size, shuffle=True):
            feed_dict = {x: X_train_a, y_: y_train_a}
            feed_dict.update( network.all_drop )    # enable dropout or dropconnect layers
            sess.run(train_op, feed_dict=feed_dict)

            # The optional feed_dict argument allows the caller to override the value of tensors in the graph. Each key in feed_dict can be one of the following types:
            # If the key is a Tensor, the value may be a Python scalar, string, list, or numpy ndarray that can be converted to the same dtype as that tensor. Additionally, if the key is a placeholder, the shape of the value will be checked for compatibility with the placeholder.
            # If the key is a SparseTensor, the value should be a SparseTensorValue.

        if epoch + 1 == 1 or (epoch + 1) % print_freq == 0:
            print("Epoch %d of %d took %fs" % (epoch + 1, n_epoch, time.time() - start_time))
            dp_dict = tl.utils.dict_to_one( network.all_drop ) # disable all dropout/dropconnect/denoising layers
            feed_dict = {x: X_train, y_: y_train}
            feed_dict.update(dp_dict)
            print("   train loss: %f" % sess.run(cost, feed_dict=feed_dict))
            dp_dict = tl.utils.dict_to_one( network.all_drop )
            feed_dict = {x: X_val, y_: y_val}
            feed_dict.update(dp_dict)
            print("   val loss: %f" % sess.run(cost, feed_dict=feed_dict))
            print("   val acc: %f" % np.mean(y_val ==
                                    sess.run(y_op, feed_dict=feed_dict)))
            try:
                # You can visualize the weight of 1st hidden layer as follow.
                tl.visualize.W(network.all_params[0].eval(), second=10,
                                        saveable=True, shape=[28, 28],
                                        name='w1_'+str(epoch+1), fig_idx=2012)
                # You can also save the weight of 1st hidden layer to .npz file.
                # tl.files.save_npz([network.all_params[0]] , name='w1'+str(epoch+1)+'.npz')
            except:
                raise Exception("You should change visualize_W(), if you want \
                            to save the feature images for different dataset")

    print('Evaluation')
    dp_dict = tl.utils.dict_to_one( network.all_drop )
    feed_dict = {x: X_test, y_: y_test}
    feed_dict.update(dp_dict)
    print("   test loss: %f" % sess.run(cost, feed_dict=feed_dict))
    print("   test acc: %f" % np.mean(y_test == sess.run(y_op,
                                                        feed_dict=feed_dict)))

    # Add ops to save and restore all the variables, including variables for training.
    # ref: https://www.tensorflow.org/versions/r0.8/how_tos/variables/index.html
    saver = tf.train.Saver()
    save_path = saver.save(sess, "model.ckpt")
    print("Model saved in file: %s" % save_path)


    # You can also save the parameters into .npz file.
    tl.files.save_npz(network.all_params , name='model.npz')
    # You can only save one parameter as follow.
    # tl.files.save_npz([network.all_params[0]] , name='model.npz')
    # Then, restore the parameters as follow.
    # load_params = tl.files.load_npz(path='', name='model.npz')
    # tl.files.assign_params(sess, load_params, network)

    # In the end, close TensorFlow session.
    sess.close()

Example 135

Project: librosa
Source File: spectrum.py
View license
@cache(level=30)
def istft(stft_matrix, hop_length=None, win_length=None, window='hann',
          center=True, dtype=np.float32):
    """
    Inverse short-time Fourier transform (ISTFT).

    Converts a complex-valued spectrogram `stft_matrix` to time-series `y`
    by minimizing the mean squared error between `stft_matrix` and STFT of
    `y` as described in [1]_.

    In general, window function, hop length and other parameters should be same
    as in stft, which mostly leads to perfect reconstruction of a signal from
    unmodified `stft_matrix`.

    .. [1] D. W. Griffin and J. S. Lim,
        "Signal estimation from modified short-time Fourier transform,"
        IEEE Trans. ASSP, vol.32, no.2, pp.236–243, Apr. 1984.

    Parameters
    ----------
    stft_matrix : np.ndarray [shape=(1 + n_fft/2, t)]
        STFT matrix from `stft`

    hop_length  : int > 0 [scalar]
        Number of frames between STFT columns.
        If unspecified, defaults to `win_length / 4`.

    win_length  : int <= n_fft = 2 * (stft_matrix.shape[0] - 1)
        When reconstructing the time series, each frame is windowed
        and each sample is normalized by the sum of squared window
        according to the `window` function (see below).

        If unspecified, defaults to `n_fft`.

    window      : string, tuple, number, function, np.ndarray [shape=(n_fft,)]
        - a window specification (string, tuple, or number);
          see `scipy.signal.get_window`
        - a window function, such as `scipy.signal.hanning`
        - a user-specified window vector of length `n_fft`

        .. see also:: `filters.get_window`

    center      : boolean
        - If `True`, `D` is assumed to have centered frames.
        - If `False`, `D` is assumed to have left-aligned frames.

    dtype       : numeric type
        Real numeric type for `y`.  Default is 32-bit float.

    Returns
    -------
    y : np.ndarray [shape=(n,)]
        time domain signal reconstructed from `stft_matrix`

    See Also
    --------
    stft : Short-time Fourier Transform

    Notes
    -----
    This function caches at level 30.

    Examples
    --------
    >>> y, sr = librosa.load(librosa.util.example_audio_file())
    >>> D = librosa.stft(y)
    >>> y_hat = librosa.istft(D)
    >>> y_hat
    array([ -4.812e-06,  -4.267e-06, ...,   6.271e-06,   2.827e-07], dtype=float32)

    Exactly preserving length of the input signal requires explicit padding.
    Otherwise, a partial frame at the end of `y` will not be represented.

    >>> n = len(y)
    >>> n_fft = 2048
    >>> y_pad = librosa.util.fix_length(y, n + n_fft // 2)
    >>> D = librosa.stft(y_pad, n_fft=n_fft)
    >>> y_out = librosa.util.fix_length(librosa.istft(D), n)
    >>> np.max(np.abs(y - y_out))
    1.4901161e-07
    """

    n_fft = 2 * (stft_matrix.shape[0] - 1)

    # By default, use the entire frame
    if win_length is None:
        win_length = n_fft

    # Set the default hop, if it's not already specified
    if hop_length is None:
        hop_length = int(win_length // 4)

    ifft_window = get_window(window, win_length, fftbins=True)

    # Pad out to match n_fft
    ifft_window = util.pad_center(ifft_window, n_fft)

    n_frames = stft_matrix.shape[1]
    expected_signal_len = n_fft + hop_length * (n_frames - 1)
    y = np.zeros(expected_signal_len, dtype=dtype)
    ifft_window_sum = np.zeros(expected_signal_len, dtype=dtype)
    ifft_window_square = ifft_window * ifft_window

    for i in range(n_frames):
        sample = i * hop_length
        spec = stft_matrix[:, i].flatten()
        spec = np.concatenate((spec.conj(), spec[-2:0:-1]), 0)
        ytmp = ifft_window * fft.ifft(spec).real

        y[sample:(sample + n_fft)] = y[sample:(sample + n_fft)] + ytmp
        ifft_window_sum[sample:(sample + n_fft)] += ifft_window_square

    # Normalize by sum of squared window
    approx_nonzero_indices = ifft_window_sum > util.tiny(ifft_window_sum)
    y[approx_nonzero_indices] /= ifft_window_sum[approx_nonzero_indices]

    if center:
        y = y[int(n_fft // 2):-int(n_fft // 2)]

    return y

Example 136

View license
def main_test_cnn_layer():
    """Reimplementation of the TensorFlow official MNIST CNN tutorials:
        # https://www.tensorflow.org/versions/r0.8/tutorials/mnist/pros/index.html
        # https://github.com/tensorflow/tensorflow/blob/master/tensorflow/models/image/mnist/convolutional.py
        More TensorFlow official CNN tutorials can be found here:
        # tutorial_cifar10.py
        # https://www.tensorflow.org/versions/master/tutorials/deep_cnn/index.html
    """
    X_train, y_train, X_val, y_val, X_test, y_test = \
                    tl.files.load_mnist_dataset(shape=(-1, 28, 28, 1))

    X_train = np.asarray(X_train, dtype=np.float32)
    y_train = np.asarray(y_train, dtype=np.int64)
    X_val = np.asarray(X_val, dtype=np.float32)
    y_val = np.asarray(y_val, dtype=np.int64)
    X_test = np.asarray(X_test, dtype=np.float32)
    y_test = np.asarray(y_test, dtype=np.int64)

    print('X_train.shape', X_train.shape)
    print('y_train.shape', y_train.shape)
    print('X_val.shape', X_val.shape)
    print('y_val.shape', y_val.shape)
    print('X_test.shape', X_test.shape)
    print('y_test.shape', y_test.shape)
    print('X %s   y %s' % (X_test.dtype, y_test.dtype))

    sess = tf.InteractiveSession()

    # Define the batchsize at the begin, you can give the batchsize in x and y_
    # rather than 'None', this can allow TensorFlow to apply some optimizations
    # – especially for convolutional layers.
    batch_size = 128

    x = tf.placeholder(tf.float32, shape=[batch_size, 28, 28, 1])   # [batch_size, height, width, channels]
    y_ = tf.placeholder(tf.int64, shape=[batch_size,])

    network = tl.layers.InputLayer(x, name='input_layer')
    network = tl.layers.Conv2dLayer(network,
                        act = tf.nn.relu,
                        shape = [5, 5, 1, 32],  # 32 features for each 5x5 patch
                        strides=[1, 1, 1, 1],
                        padding='SAME',
                        name ='cnn_layer1')     # output: (?, 28, 28, 32)
    network = tl.layers.PoolLayer(network,
                        ksize=[1, 2, 2, 1],
                        strides=[1, 2, 2, 1],
                        padding='SAME',
                        pool = tf.nn.max_pool,
                        name ='pool_layer1',)   # output: (?, 14, 14, 32)
    network = tl.layers.Conv2dLayer(network,
                        act = tf.nn.relu,
                        shape = [5, 5, 32, 64], # 64 features for each 5x5 patch
                        strides=[1, 1, 1, 1],
                        padding='SAME',
                        name ='cnn_layer2')     # output: (?, 14, 14, 64)
    network = tl.layers.PoolLayer(network,
                        ksize=[1, 2, 2, 1],
                        strides=[1, 2, 2, 1],
                        padding='SAME',
                        pool = tf.nn.max_pool,
                        name ='pool_layer2',)   # output: (?, 7, 7, 64)
    network = tl.layers.FlattenLayer(network, name='flatten_layer')   # output: (?, 3136)
    network = tl.layers.DropoutLayer(network, keep=0.5, name='drop1') # output: (?, 3136)
    network = tl.layers.DenseLayer(network, n_units=256,
                                    act = tf.nn.relu, name='relu1')   # output: (?, 256)
    network = tl.layers.DropoutLayer(network, keep=0.5, name='drop2') # output: (?, 256)
    network = tl.layers.DenseLayer(network, n_units=10,
                                    act = tf.identity,
                                    name='output_layer')    # output: (?, 10)

    y = network.outputs

    ce = tf.reduce_mean(tf.nn.sparse_softmax_cross_entropy_with_logits(y, y_))
    cost = ce

    correct_prediction = tf.equal(tf.argmax(y, 1), y_)
    acc = tf.reduce_mean(tf.cast(correct_prediction, tf.float32))

    # train
    n_epoch = 200
    learning_rate = 0.0001
    print_freq = 10

    train_params = network.all_params
    train_op = tf.train.AdamOptimizer(learning_rate, beta1=0.9, beta2=0.999,
        epsilon=1e-08, use_locking=False).minimize(cost, var_list=train_params)

    sess.run(tf.initialize_all_variables())
    network.print_params()
    network.print_layers()

    print('   learning_rate: %f' % learning_rate)
    print('   batch_size: %d' % batch_size)

    for epoch in range(n_epoch):
        start_time = time.time()
        for X_train_a, y_train_a in tl.iterate.minibatches(
                                    X_train, y_train, batch_size, shuffle=True):
            feed_dict = {x: X_train_a, y_: y_train_a}
            feed_dict.update( network.all_drop )        # enable noise layers
            sess.run(train_op, feed_dict=feed_dict)

        if epoch + 1 == 1 or (epoch + 1) % print_freq == 0:
            print("Epoch %d of %d took %fs" % (epoch + 1, n_epoch, time.time() - start_time))
            train_loss, train_acc, n_batch = 0, 0, 0
            for X_train_a, y_train_a in tl.iterate.minibatches(
                                    X_train, y_train, batch_size, shuffle=True):
                dp_dict = tl.utils.dict_to_one( network.all_drop )    # disable noise layers
                feed_dict = {x: X_train_a, y_: y_train_a}
                feed_dict.update(dp_dict)
                err, ac = sess.run([cost, acc], feed_dict=feed_dict)
                train_loss += err; train_acc += ac; n_batch += 1
            print("   train loss: %f" % (train_loss/ n_batch))
            print("   train acc: %f" % (train_acc/ n_batch))
            val_loss, val_acc, n_batch = 0, 0, 0
            for X_val_a, y_val_a in tl.iterate.minibatches(
                                        X_val, y_val, batch_size, shuffle=True):
                dp_dict = tl.utils.dict_to_one( network.all_drop )    # disable noise layers
                feed_dict = {x: X_val_a, y_: y_val_a}
                feed_dict.update(dp_dict)
                err, ac = sess.run([cost, acc], feed_dict=feed_dict)
                val_loss += err; val_acc += ac; n_batch += 1
            print("   val loss: %f" % (val_loss/ n_batch))
            print("   val acc: %f" % (val_acc/ n_batch))
            try:
                tl.visualize.CNN2d(network.all_params[0].eval(),
                                    second=10, saveable=True,
                                    name='cnn1_'+str(epoch+1), fig_idx=2012)
            except:
                raise Exception("# You should change visualize.CNN(), if you want to save the feature images for different dataset")

    print('Evaluation')
    test_loss, test_acc, n_batch = 0, 0, 0
    for X_test_a, y_test_a in tl.iterate.minibatches(
                                X_test, y_test, batch_size, shuffle=True):
        dp_dict = tl.utils.dict_to_one( network.all_drop )    # disable noise layers
        feed_dict = {x: X_test_a, y_: y_test_a}
        feed_dict.update(dp_dict)
        err, ac = sess.run([cost, acc], feed_dict=feed_dict)
        test_loss += err; test_acc += ac; n_batch += 1
    print("   test loss: %f" % (test_loss/n_batch))
    print("   test acc: %f" % (test_acc/n_batch))

Example 137

Project: VIP
Source File: derotation.py
View license
def frame_rotate(array, angle, imlib='opencv', interpolation='bicubic', cxy=None):
    """ Rotates a frame.
    
    Parameters
    ----------
    array : array_like 
        Input frame, 2d array.
    angle : float
        Rotation angle.
    imlib : {'opencv', 'skimage'}, str optional
        Library used for image transformations. Opencv is faster than ndimage or
        skimage.
    interpolation : {'bicubic', 'bilinear', 'nearneig'}, optional
        'nneighbor' stands for nearest-neighbor interpolation,
        'bilinear' stands for bilinear interpolation,
        'bicubic' for interpolation over 4x4 pixel neighborhood.
        The 'bicubic' is the default. The 'nearneig' is the fastest method and
        the 'bicubic' the slowest of the three. The 'nearneig' is the poorer
        option for interpolation of noisy astronomical images.
    cxy : float, optional
        Coordinates X,Y  of the point with respect to which the rotation will be 
        performed. By default the rotation is done with respect to the center 
        of the frame; central pixel if frame has odd size.
        
    Returns
    -------
    array_out : array_like
        Resulting frame.
        
    """
    if not array.ndim == 2:
        raise TypeError('Input array is not a frame or 2d array.')
    array = np.float32(array)
    y, x = array.shape
    
    if not cxy:  
        cy, cx = frame_center(array)
    else:
        cx, cy = cxy

    if imlib not in ['skimage', 'opencv']:
        raise ValueError('Imlib not recognized, try opencv or ndimage')

    if imlib=='skimage' or no_opencv:
        if interpolation == 'bilinear':
            order = 1
        elif interpolation == 'bicubic':
            order = 3
        elif interpolation == 'nearneig':
            order = 0
        else:
            raise TypeError('Interpolation method not recognized.')

        min_val = np.min(array)
        im_temp = array - min_val
        max_val = np.max(im_temp)
        im_temp /= max_val

        array_out = rotate(im_temp, angle, order=order, center=cxy, cval=np.nan)

        array_out *= max_val
        array_out += min_val
        array_out = np.nan_to_num(array_out)

    else:
        if interpolation == 'bilinear':
            intp = cv2.INTER_LINEAR
        elif interpolation == 'bicubic':
            intp= cv2.INTER_CUBIC
        elif interpolation == 'nearneig':
            intp = cv2.INTER_NEAREST
        else:
            raise TypeError('Interpolation method not recognized.')

        M = cv2.getRotationMatrix2D((cx,cy), angle, 1)
        array_out = cv2.warpAffine(array.astype(np.float32), M, (x, y), flags=intp)
             
    return array_out

Example 138

Project: visvis
Source File: screenshot.py
View license
def screenshot(filename, ob=None, sf=2, bg=None, format=None, tension=-0.25):
    """ screenshot(filename, ob=None sf=2, bg=None, format=None)
    
    Make a screenshot and store it to a file, using cubic interpolation
    to increase the resolution (and quality) of the image.
    
    Parameters
    ----------
    filename : string
        The name of the file to store the screenshot to. If filename is None, 
        the interpolated image is returned as a numpy array.
    ob : Axes, AxesContainer, or Figure
        The object to take the screenshot of. The AxesContainer can be
        obtained using vv.gca().parent. It can be usefull to take a 
        screeshot of an axes including thickmarks and labels.
    sf : integer
        The scale factor. The image is increased in size with this factor,
        using a high quality interpolation method. A factor of 2 or 3
        is recommended; the image quality does not improve with higher
        factors. If using a sf larger than 1, the image is best saved in
        the jpg format.
    bg : 3-element tuple or char
        The color of the background. If bg is given, ob.bgcolor is set to
        bg before the frame is captured.
    format : string
        The format for the screenshot to be saved in. If not given, the
        format is deduced from the filename.
    
    Notes
    -----
    Uses vv.getframe(ob) to obtain the image in the figure or axes. 
    That image is interpolated with the given scale factor (sf) using 
    bicubic interpolation. Then  vv.imwrite(filename, ..) is used to 
    store the resulting image to a file.
    
    Rationale
    ---------
    We'd prefer storing screenshots of plots as vector (eps) images, but 
    the nature of OpenGl prevents this. By applying high quality 
    interpolation (using a cardinal spline), the resolution can be increased, 
    thereby significantly improving the visibility/smoothness for lines 
    and fonts. Use this to produce publication quality snapshots of your
    plots.
    
    """
    
    # The tension controls the
    # responsivenes of the filter. The more negative, the more overshoot,
    # but the more it is capable to make for example font glyphs smooth.
    # If tension is 0, the interpolator is a Catmull-Rom spline.
    
    # Scale must be integer
    s = int(sf)
    
    # Object given?
    if ob is None:
        ob = vv.gcf()
    
    # Get figure
    fig = ob
    if not hasattr(ob, 'DrawNow'):
        fig = ob.GetFigure()
    
    # Get object to set background of
    bgob = ob
    
    # Set background
    if bg and fig:
        bgOld = bgob.bgcolor
        bgob.bgcolor = bg
        fig.DrawNow()  
    
    # Obtain image      
    im1 = vv.getframe(ob)
    shape1 = im1.shape
    
    # Return background
    if bg and fig:
        bgob.bgcolor = bgOld
        fig.Draw()
    
    # Pad original image, so we have no trouble at the edges
    shape2 = shape1[0]+2, shape1[1]+2, 3
    im2 = np.zeros(shape2, dtype=np.float32) # Also make float
    im2[1:-1,1:-1,:] = im1
    im2[0,:,:] = im2[1,:,:]
    im2[-1,:,:] = im2[-2,:,:]
    im2[:,0,:] = im2[:,1,:]
    im2[:,-1,:] = im2[:,-2,:]
    
    # Create empty new image. It is sized by the scaleFactor, 
    # but the last row is not. 
    shape3 = (shape1[0]-1)*s+1, (shape1[1]-1)*s+1, 3    
    im3 = np.zeros(shape3, dtype=np.float32)
    
    # Fill in values!
    for dy in range(s+1):
        for dx in range(s+1):
            
            # Get interpolation fraction and coefs
            ty = float(dy)/s
            tx = float(dx)/s
            cy = getCardinalSplineCoefs(ty, tension)
            cx = getCardinalSplineCoefs(tx, tension)
            
            # Create tmp image to which we add the contributions
            # Note that this image is 1 pixel smaller in each dimension.
            # The last pixel is filled because dy and dx iterate INCLUDING s.
            shapeTmp = shape1[0]-1, shape1[1]-1, 3
            imTmp = np.zeros(shapeTmp, dtype=np.float32)
            
            # Collect all 16 neighbours and weight them apropriately
            for iy in range(4):
                for ix in range(4):
                    
                    # Get weight
                    w = cy[iy]*cx[ix]
                    if w==0:
                        continue
                    
                    # Get slice. Note that we start at 0,1,2,3 rather than
                    # -1,0,1,2, because we padded the image.
                    D = {0:-3,1:-2,2:-1,3:None}
                    slicey = slice(iy, D[iy])
                    slicex = slice(ix, D[ix])
                    
                    # Get contribution and add to temp image
                    imTmp += w * im2[slicey, slicex, :]
            
            # Store contributions            
            D = [-1 for tmp in range(s)]; D.append(None)
            slicey = slice(dy,D[dy],s)
            slicex = slice(dx,D[dx],s)
            im3[slicey, slicex, :] = imTmp
    
    
    # Correct for overshoot
    im3[im3>1]=1
    im3[im3<0]=0
    
    # Store image to file
    if filename is not None:
        vv.imwrite(filename, im3, format)
    else:
        return im3

Example 139

Project: medpy
Source File: texture.py
View license
def coarseness(image, voxelspacing = None, mask = slice(None)):
    r"""
    Takes a simple or multi-spectral image and returns the coarseness of the texture.
    
    Step1  At each pixel, compute six averages for the windows of size 2**k x 2**k,
            k=0,1,...,5, around the pixel. 
    Step2  At each pixel, compute absolute differences E between the pairs of non 
            overlapping averages in every directions.
    step3  At each pixel, find the value of k that maximises the difference Ek in either 
            direction and set the best size Sbest=2**k
    step4  Compute the coarseness feature Fcrs by averaging Sbest over the entire image.

    Parameters
    ----------
    image : array_like or list/tuple of array_like 
        A single image or a list/tuple of images (for multi-spectral case).
    voxelspacing : sequence of floats
        The side-length of each voxel.
    mask : array_like
        A binary mask for the image or a slice object
        
    Returns
    -------
    coarseness : float
        The size of coarseness of the given texture. It is basically the size of
        repeating elements in the image. 
        
    See Also
    --------
    
    
    """
    # Step1:  At each pixel (x,y), compute six averages for the windows
    # of size 2**k x 2**k, k=0,1,...,5, around the pixel.

    image = numpy.asarray(image, dtype=numpy.float32)
   
  
    # set default mask or apply given mask
    if not type(mask) is slice:
        if not type(mask[0] is slice):
            mask = numpy.array(mask, copy=False, dtype = numpy.bool)
    image = image[mask]
    
    # set default voxel spacing if not suppliec
    if None == voxelspacing:
        voxelspacing = tuple([1.] * image.ndim)
    
    if len(voxelspacing) != image.ndim:
        print "Voxel spacing and image dimensions do not fit."
        return None
    # set padding for image border control
    padSize = tuple((numpy.rint((2**5.0) * voxelspacing[jj]),0) for jj in xrange(image.ndim))        
    Apad = numpy.pad(image,pad_width=padSize, mode='reflect')

    # Allocate memory
    E = numpy.empty((6,image.ndim)+image.shape)

    # prepare some slicer 
    rawSlicer           = [slice(None)] * image.ndim
    slicerForImageInPad = [slice(padSize[d][0],None)for d in xrange(image.ndim)]

    for k in xrange(6):

        size_vs = tuple(numpy.rint((2**k) * voxelspacing[jj]) for jj in xrange(image.ndim))
        A = uniform_filter(Apad, size = size_vs, mode = 'mirror')

        # Step2: At each pixel, compute absolute differences E(x,y) between 
        # the pairs of non overlapping averages in the horizontal and vertical directions.
        for d in xrange(image.ndim):
            borders = numpy.rint((2**k) * voxelspacing[d])
            
            slicerPad_k_d   = slicerForImageInPad[:]
            slicerPad_k_d[d]= slice((padSize[d][0]-borders if borders < padSize[d][0] else 0),None)
            A_k_d           = A[slicerPad_k_d]

            AslicerL        = rawSlicer[:]
            AslicerL[d]     = slice(0, -borders)
            
            AslicerR        = rawSlicer[:]
            AslicerR[d]     = slice(borders, None)

            E[k,d,...] = numpy.abs(A_k_d[AslicerL] - A_k_d[AslicerR])

    # step3: At each pixel, find the value of k that maximises the difference Ek(x,y)
    # in either direction and set the best size Sbest(x,y)=2**k
    
    k_max = E.max(1).argmax(0)
    dim = E.argmax(1)
    dim_vox_space = numpy.asarray([voxelspacing[dim[k_max.flat[i]].flat[i]] for i in xrange(k_max.size)]).reshape(k_max.shape) 
    S = (2**k_max) * dim_vox_space

    # step4: Compute the coarseness feature Fcrs by averaging Sbest(x,y) over the entire image.
    return S.mean()

Example 140

Project: tensorlayer
Source File: tutorial_mnist.py
View license
def main_test_cnn_layer():
    """Reimplementation of the TensorFlow official MNIST CNN tutorials:
        # https://www.tensorflow.org/versions/r0.8/tutorials/mnist/pros/index.html
        # https://github.com/tensorflow/tensorflow/blob/master/tensorflow/models/image/mnist/convolutional.py
        More TensorFlow official CNN tutorials can be found here:
        # tutorial_cifar10.py
        # https://www.tensorflow.org/versions/master/tutorials/deep_cnn/index.html
    """
    X_train, y_train, X_val, y_val, X_test, y_test = \
                    tl.files.load_mnist_dataset(shape=(-1, 28, 28, 1))

    X_train = np.asarray(X_train, dtype=np.float32)
    y_train = np.asarray(y_train, dtype=np.int64)
    X_val = np.asarray(X_val, dtype=np.float32)
    y_val = np.asarray(y_val, dtype=np.int64)
    X_test = np.asarray(X_test, dtype=np.float32)
    y_test = np.asarray(y_test, dtype=np.int64)

    print('X_train.shape', X_train.shape)
    print('y_train.shape', y_train.shape)
    print('X_val.shape', X_val.shape)
    print('y_val.shape', y_val.shape)
    print('X_test.shape', X_test.shape)
    print('y_test.shape', y_test.shape)
    print('X %s   y %s' % (X_test.dtype, y_test.dtype))

    sess = tf.InteractiveSession()

    # Define the batchsize at the begin, you can give the batchsize in x and y_
    # rather than 'None', this can allow TensorFlow to apply some optimizations
    # – especially for convolutional layers.
    batch_size = 128

    x = tf.placeholder(tf.float32, shape=[batch_size, 28, 28, 1])   # [batch_size, height, width, channels]
    y_ = tf.placeholder(tf.int64, shape=[batch_size,])

    network = tl.layers.InputLayer(x, name='input_layer')
    network = tl.layers.Conv2dLayer(network,
                        act = tf.nn.relu,
                        shape = [5, 5, 1, 32],  # 32 features for each 5x5 patch
                        strides=[1, 1, 1, 1],
                        padding='SAME',
                        name ='cnn_layer1')     # output: (?, 28, 28, 32)
    network = tl.layers.PoolLayer(network,
                        ksize=[1, 2, 2, 1],
                        strides=[1, 2, 2, 1],
                        padding='SAME',
                        pool = tf.nn.max_pool,
                        name ='pool_layer1',)   # output: (?, 14, 14, 32)
    network = tl.layers.Conv2dLayer(network,
                        act = tf.nn.relu,
                        shape = [5, 5, 32, 64], # 64 features for each 5x5 patch
                        strides=[1, 1, 1, 1],
                        padding='SAME',
                        name ='cnn_layer2')     # output: (?, 14, 14, 64)
    network = tl.layers.PoolLayer(network,
                        ksize=[1, 2, 2, 1],
                        strides=[1, 2, 2, 1],
                        padding='SAME',
                        pool = tf.nn.max_pool,
                        name ='pool_layer2',)   # output: (?, 7, 7, 64)
    network = tl.layers.FlattenLayer(network, name='flatten_layer')   # output: (?, 3136)
    network = tl.layers.DropoutLayer(network, keep=0.5, name='drop1') # output: (?, 3136)
    network = tl.layers.DenseLayer(network, n_units=256,
                                    act = tf.nn.relu, name='relu1')   # output: (?, 256)
    network = tl.layers.DropoutLayer(network, keep=0.5, name='drop2') # output: (?, 256)
    network = tl.layers.DenseLayer(network, n_units=10,
                                    act = tf.identity,
                                    name='output_layer')    # output: (?, 10)

    y = network.outputs

    ce = tf.reduce_mean(tf.nn.sparse_softmax_cross_entropy_with_logits(y, y_))
    cost = ce

    correct_prediction = tf.equal(tf.argmax(y, 1), y_)
    acc = tf.reduce_mean(tf.cast(correct_prediction, tf.float32))

    # train
    n_epoch = 200
    learning_rate = 0.0001
    print_freq = 10

    train_params = network.all_params
    train_op = tf.train.AdamOptimizer(learning_rate, beta1=0.9, beta2=0.999,
        epsilon=1e-08, use_locking=False).minimize(cost, var_list=train_params)

    sess.run(tf.initialize_all_variables())
    network.print_params()
    network.print_layers()

    print('   learning_rate: %f' % learning_rate)
    print('   batch_size: %d' % batch_size)

    for epoch in range(n_epoch):
        start_time = time.time()
        for X_train_a, y_train_a in tl.iterate.minibatches(
                                    X_train, y_train, batch_size, shuffle=True):
            feed_dict = {x: X_train_a, y_: y_train_a}
            feed_dict.update( network.all_drop )        # enable noise layers
            sess.run(train_op, feed_dict=feed_dict)

        if epoch + 1 == 1 or (epoch + 1) % print_freq == 0:
            print("Epoch %d of %d took %fs" % (epoch + 1, n_epoch, time.time() - start_time))
            train_loss, train_acc, n_batch = 0, 0, 0
            for X_train_a, y_train_a in tl.iterate.minibatches(
                                    X_train, y_train, batch_size, shuffle=True):
                dp_dict = tl.utils.dict_to_one( network.all_drop )    # disable noise layers
                feed_dict = {x: X_train_a, y_: y_train_a}
                feed_dict.update(dp_dict)
                err, ac = sess.run([cost, acc], feed_dict=feed_dict)
                train_loss += err; train_acc += ac; n_batch += 1
            print("   train loss: %f" % (train_loss/ n_batch))
            print("   train acc: %f" % (train_acc/ n_batch))
            val_loss, val_acc, n_batch = 0, 0, 0
            for X_val_a, y_val_a in tl.iterate.minibatches(
                                        X_val, y_val, batch_size, shuffle=True):
                dp_dict = tl.utils.dict_to_one( network.all_drop )    # disable noise layers
                feed_dict = {x: X_val_a, y_: y_val_a}
                feed_dict.update(dp_dict)
                err, ac = sess.run([cost, acc], feed_dict=feed_dict)
                val_loss += err; val_acc += ac; n_batch += 1
            print("   val loss: %f" % (val_loss/ n_batch))
            print("   val acc: %f" % (val_acc/ n_batch))
            try:
                tl.visualize.CNN2d(network.all_params[0].eval(),
                                    second=10, saveable=True,
                                    name='cnn1_'+str(epoch+1), fig_idx=2012)
            except:
                raise Exception("# You should change visualize.CNN(), if you want to save the feature images for different dataset")

    print('Evaluation')
    test_loss, test_acc, n_batch = 0, 0, 0
    for X_test_a, y_test_a in tl.iterate.minibatches(
                                X_test, y_test, batch_size, shuffle=True):
        dp_dict = tl.utils.dict_to_one( network.all_drop )    # disable noise layers
        feed_dict = {x: X_test_a, y_: y_test_a}
        feed_dict.update(dp_dict)
        err, ac = sess.run([cost, acc], feed_dict=feed_dict)
        test_loss += err; test_acc += ac; n_batch += 1
    print("   test loss: %f" % (test_loss/n_batch))
    print("   test acc: %f" % (test_acc/n_batch))

Example 141

Project: hwrt
Source File: segmentation.py
View license
def get_segmentation(recording,
                     single_clf,
                     single_stroke_clf,
                     stroke_segmented_classifier):
    """
    Get a list of segmentations of recording with the probability of the
    segmentation being correct.

    Parameters
    ----------
    recording : A list of lists
        Each sublist represents a stroke
    single_clf : object
        A classifier for single symbols
    single_stroke_clf : object
        A classifier which decides if a single stroke is a complete symbol
    stroke_segmented_classifier : object
        Classifier which decides if two strokes belong to one symbol or not

    Returns
    -------
    list of tuples :
        Segmentations together with their probabilities. Each probability
        has to be positive and the sum may not be bigger than 1.0.

    Examples
    --------
    >>> stroke1 = [{'x': 0, 'y': 0, 'time': 0}, {'x': 12, 'y': 12, 'time': 1}]
    >>> stroke2 = [{'x': 0, 'y': 10, 'time': 2}, {'x': 12, 'y': 0, 'time': 3}]
    >>> stroke3 = [{'x': 14, 'y': 0, 'time': 5}, {'x': 14, 'y': 12, 'time': 6}]
    >>> #get_segmentation([stroke1, stroke2, stroke3], single_clf)
    [
      ([[0, 1], [2]], 0.8),
      ([[0], [1,2]], 0.1),
      ([[0,2], [1]], 0.05)
    ]
    """
    mst_wood = get_mst_wood(recording, single_clf)
    return [(normalize_segmentation([mst['strokes'] for mst in mst_wood]),
             1.0)]

    # HandwrittenData(json.dumps(recording)).show()
    # return [([[i for i in range(len(recording))]], 1.0)]
    # #mst_wood = break_mst(mst, recording)  # TODO

    # for i in range(0, 2**len(points)):
    #     segmentation = get_segmentation_from_mst(mst, i)
    # TODO

    X_symbol = [get_median_stroke_distance(recording)]

    # Pre-segment to 8 strokes
    # TODO: Take first 4 strokes and add strokes within their bounding box
    # TODO: What if that is more then 8 strokes?
    # -> Geometry
    #    Build tree structure. A stroke `c` is the child of another stroke `p`,
    #    if the bounding box of `c` is within the bounding box of `p`.
    #       Problem: B <-> 13
    g_top_segmentations = [([], 1.0)]  # g_top_segmentations

    # range(int(math.ceil(float(len(recording))/8))):
    for chunk_part in mst_wood:
        # chunk = recording[8*chunk_part:8*(chunk_part+1)]
        chunk = [recording[stroke] for stroke in chunk_part['strokes']]

        # Segment after pre-segmentation
        prob = [[1.0 for _ in chunk] for _ in chunk]
        for strokeid1, strokeid2 in itertools.product(range(len(chunk)),
                                                      range(len(chunk))):
            if strokeid1 == strokeid2:
                continue
            X = get_stroke_features(chunk, strokeid1, strokeid2)
            X += X_symbol
            X = numpy.array([X], dtype=numpy.float32)
            prob[strokeid1][strokeid2] = stroke_segmented_classifier(X)

        # Top segmentations
        ts = list(partitions.get_top_segmentations(prob, 500))
        for i, segmentation in enumerate(ts):
            symbols = apply_segmentation(chunk, segmentation)
            min_top2 = partitions.TopFinder(1, find_min=True)
            for i, symbol in enumerate(symbols):
                predictions = single_clf.predict(symbol)
                min_top2.push("value-%i" % i,
                              predictions[0]['probability'] +
                              predictions[1]['probability'])
            ts[i][1] *= list(min_top2)[0][1]
        # for i, segmentation in enumerate(ts):
        #     ts[i][0] = update_segmentation_data(ts[i][0], 8*chunk_part)
        g_top_segmentations = merge_segmentations(g_top_segmentations,
                                                  ts,
                                                  chunk_part['strokes'])
    return [(normalize_segmentation(seg), probability)
            for seg, probability in g_top_segmentations]

Example 142

Project: spinn
Source File: stack.py
View license
    def make_backprop_scan(self, error_signal,
                           extra_cost_inputs=None,
                           compute_embedding_gradients=True):
        """
        Compile the backpropagation graph for a ThinStack instance.

        This method analyzes the feedforward graph of a ThinStack instance and
        analytically calculates its gradient with respect to the stack's
        parameters of interest.

        This method must be called by the client after constructing ThinStack.
        After creating the gradient graph, the outputs of the graph (gradient
        batches) are installed into an instance member self.gradients. This
        object is a dictionary mapping from stack shared variables (parameters)
        to their symbolic gradients.

        Args:
            error_signal: The external gradient d(cost)/d(stack top). A Theano
                batch of size `batch_size * model_dim`.
            extra_cost_inputs: Other symbolic variables which may be involved
                in the symbolic cost expression / the error signal that should
                be passed as extra inputs to a backpropagation scan.
            compute_embedding_gradients: Calculate gradients for each embedding
                vector.
        """

        # The meat of the implementation of this method is in the inner
        # function step_b. What comes before is mainly preparatory code; what
        # comes after is the scan invocation that unrolls the step_b function.

        assert hasattr(self, "stack_2_ptrs"), \
            ("self._make_scan (forward pass) must be defined before "
             "self.make_backprop_scan is called")

        # We need to add extra updates to the `_zero_updates` member, so we
        # must be called before `_zero_updates` is read.
        assert self._zero is None, \
            ("Can only install backprop on a fresh ThinStack. Don't call "
             "ThinStack.zero before setting up backprop.")

        if (compute_embedding_gradients
            and self._embedding_projection_network not in [None, util.IdentityLayer]):
            raise ValueError(
                "Do not support backprop for both an embedding projection "
                "layer and individual embeddings.")

        if self.use_input_batch_norm:
            raise ValueError(
                "Thin-stack backprop not supported with input batch-norm. Jon "
                "worked on BN gradients for 3 days without success, and then "
                "dropped it.")

        if extra_cost_inputs is None:
            extra_cost_inputs = []

        # Analytically calculate backpropagation graphs for an abstract single
        # timestep. The list 'wrt' contains all variables for which we should
        # collect gradients (this list of course depends on the nature of the
        # forward-prop activation function). The latter 3 returns are actually
        # graph "thunks" which return a concrete backpropagation graph when
        # provided with actual symbolic inputs / error signals.
        #
        # These graphs will be repeatedly invoked and chained together in the
        # code further below.
        wrt, f_proj_delta, f_shift_delta, f_reduce_delta = \
                self._make_backward_graphs()
        wrt_shapes = [wrt_i.get_value().shape for wrt_i in wrt]

        # Build shared variables for accumulating wrt deltas.
        wrt_vars = [theano.shared(np.zeros(wrt_shape, dtype=np.float32),
                                  name=self._prefix + "bwd/wrt/%s" % wrt_i)
                    for wrt_i, wrt_shape in zip(wrt, wrt_shapes)]
        # All of these need to be zeroed out in between batches.
        self._zero_updates += wrt_vars

        # Also accumulate embedding gradients separately
        if compute_embedding_gradients:
            dE = theano.shared(np.zeros(self.embeddings.get_value().shape,
                                        dtype=np.float32),
                               name=self._prefix + "bwd/wrt/embeddings")
            self._zero_updates.append(dE)
        else:
            # Make dE a dummy variable.
            dE = T.zeros((1,))

        # Useful batch zero-constants.
        zero_stack = T.zeros((self.batch_size, self.model_dim))
        zero_extra_inps = [T.zeros((self.batch_size, extra_shape[-1]))
                           for extra_shape in self.recurrence.extra_outputs]
        # Zero Jacobian matrices for masked reductions during backprop. May not
        # be used.
        zero_jac_wrts = [T.zeros((self.batch_size,) + wrt_shape)
                         for wrt_shape in wrt_shapes]

        DUMMY = util.zeros_nobroadcast((1,))

        def lookup(t_f, stack_fwd, stack_2_ptrs_t, buffer_cur_t,
                  stack_bwd_t, extra_bwd):
            """Retrieve all relevant bwd inputs/outputs at time `t`."""

            grad_cursor = t_f * self.batch_size + self._stack_shift
            main_grad = cuda.AdvancedSubtensor1Floats("B_maingrad")(
                stack_bwd_t, grad_cursor)
            extra_grads = tuple([
                cuda.AdvancedSubtensor1Floats("B_extragrad_%i" % i)(
                    extra_bwd_i, grad_cursor)
                for i, extra_bwd_i in enumerate(extra_bwd)])

            # Find the timesteps of the two elements involved in the potential
            # reduce at this timestep.
            t_c1 = (t_f - 1.0) * self.batch_size + self._stack_shift
            t_c2 = stack_2_ptrs_t

            # Find the two elements involved in the potential reduce.
            c1 = cuda.AdvancedSubtensor1Floats("B_stack1")(stack_fwd, t_c1)
            c2 = cuda.AdvancedSubtensor1Floats("B_stack2")(stack_fwd, t_c2)

            buffer_top_t = cuda.AdvancedSubtensor1Floats("B_buffer_top")(
                self.buffer_t, buffer_cur_t + self._buffer_shift)

            # Retrieve extra inputs from auxiliary stack(s).
            extra_inps_t = tuple([
                cuda.AdvancedSubtensor1Floats("B_extra_inp_%i" % i)(
                    extra_inp_i, t_c1)
                for extra_inp_i in self.final_aux_stacks])

            inputs = (c1, c2, buffer_top_t) + extra_inps_t
            grads = (main_grad,) + extra_grads
            return t_c1, t_c2, inputs, grads

        def step_b(# sequences
                   t_f, transitions_t_f, stack_2_ptrs_t, buffer_cur_t,
                   # accumulators
                   dE,
                   # rest (incl. outputs_info, non_sequences)
                   *rest):

            # Separate the accum arguments from the non-sequence arguments.
            n_wrt = len(wrt_shapes)
            n_extra_bwd = len(self.recurrence.extra_outputs)
            wrt_deltas = rest[:n_wrt]
            stack_bwd_t = rest[n_wrt]
            extra_bwd = rest[n_wrt + 1:n_wrt + 1 + n_extra_bwd]
            id_buffer, stack_final = \
                rest[n_wrt + 1 + n_extra_bwd:n_wrt + 1 + n_extra_bwd + 2]

            # At first iteration, drop the external error signal into the main
            # backward stack.
            stack_bwd_next = ifelse(T.eq(t_f, self.seq_length),
                                    T.set_subtensor(stack_bwd_t[-self.batch_size:], error_signal),
                                    stack_bwd_t)


            # Retrieve all relevant inputs/outputs at this timestep.
            t_c1, t_c2, inputs, grads = \
                lookup(t_f, stack_final, stack_2_ptrs_t, buffer_cur_t,
                       stack_bwd_next, extra_bwd)
            main_grad = grads[0]

            # Calculate deltas for this timestep.
            r_delta_inp, r_delta_wrt = f_reduce_delta(inputs, grads)
            # NB: main_grad is not passed to shift function.
            s_delta_inp, s_delta_wrt = f_shift_delta(inputs, grads[1:])

            # Check that delta function outputs match (at least in number).
            assert len(r_delta_inp) == len(s_delta_inp), \
                "%i %i" % (len(r_delta_inp), len(s_delta_inp))
            assert len(r_delta_wrt) == len(s_delta_wrt), \
                "%i %i" % (len(r_delta_wrt), len(s_delta_wrt))
            assert len(r_delta_inp) == 3 + len(self.aux_stacks), \
                "%i %i" % (len(r_delta_inp), 3 + len(self.aux_stacks))
            assert len(r_delta_wrt) == len(wrt)

            # Retrieve embedding indices on buffer at this timestep.
            # (Necessary for sending embedding gradients.)
            buffer_ids_t = cuda.AdvancedSubtensor1Floats("B_buffer_ids")(
                    id_buffer, buffer_cur_t + self._buffer_shift)

            # Prepare masks for op-wise gradient accumulation.
            # TODO: Record actual transitions (e.g. for model 1S and higher)
            # and repeat those here
            mask = transitions_t_f
            masks = [mask, mask.dimshuffle(0, "x"),
                     mask.dimshuffle(0, "x", "x")]

            # Accumulate gradients from the embedding projection network into
            # the stack op gradients.
            if f_proj_delta is not None:
                # Look up raw buffer top for this timestep -- i.e., buffer top
                # *before* the op at this timestep was performed. This was the
                # input to the projection network at this timestep.
                proj_input = cuda.AdvancedSubtensor1Floats("B_raw_buffer_top")(
                    self._raw_buffer_t, buffer_cur_t + self._buffer_shift)

                proj_inputs = (proj_input,)
                if self.use_input_dropout:
                    embedding_dropout_mask = cuda.AdvancedSubtensor1Floats("B_buffer_dropout")(
                        self._embedding_dropout_masks, buffer_cur_t + self._buffer_shift)
                    proj_inputs = (proj_input, embedding_dropout_mask)

                # Compute separate graphs based on gradient from above.
                # NB: We discard the delta_inp return here. The delta_inp
                # should actually be passed back to the raw embedding
                # parameters, but we don't have any reason to support this in
                # practice. (Either we backprop to embeddings or project them
                # and learn the projection -- not both.)
                if r_delta_inp[2] is not None:
                    _, m_proj_delta_wrt = f_proj_delta(proj_inputs,
                                                       (r_delta_inp[2],))
                    r_delta_wrt = util.merge_update_lists(r_delta_wrt, m_proj_delta_wrt)

                # If we shifted (moved the buffer top onto the stack), the
                # gradient from above is a combination of the accumulated stack
                # gradient (main_grad) and any buffer top deltas from the shift
                # function (e.g. tracking LSTM gradient).
                embedding_grad = main_grad
                if s_delta_inp[2] is not None:
                    embedding_grad += s_delta_inp[2]
                _, p_proj_delta_wrt = f_proj_delta(proj_inputs,
                                                   (embedding_grad,))
                s_delta_wrt = util.merge_update_lists(s_delta_wrt, p_proj_delta_wrt)

            ###########
            # STEP 1. #
            ###########
            # Accumulate deltas onto the graph inputs, switching over
            # shift/reduce decisions.
            stacks = (stack_bwd_next, stack_bwd_next,
                      (compute_embedding_gradients and dE) or None)
            cursors = (t_c1, t_c2,
                       (compute_embedding_gradients and buffer_ids_t) or None)
            # Handle potential aux bwd stacks.
            stacks += extra_bwd
            cursors += ((t_c1,)) * len(extra_bwd)
            new_stacks = {}
            for stack, cursor, r_delta, s_delta in zip(stacks, cursors, r_delta_inp, s_delta_inp):
                if stack is None or cursor is None:
                    continue
                elif r_delta is None and s_delta is None:
                    # Disconnected gradient.
                    continue

                base = new_stacks.get(stack, stack)
                mask_i = masks[(r_delta or s_delta).ndim - 1]
                if r_delta is None:
                    delta = (1. - mask_i) * s_delta
                elif s_delta is None:
                    delta = mask_i * r_delta
                else:
                    delta = mask_i * r_delta + (1. - mask_i) * s_delta

                # Run subtensor update on associated structure using the
                # current cursor.
                new_stack = cuda.AdvancedIncSubtensor1Floats(inplace=True)(
                    base, delta, cursor)
                new_stacks[stack] = new_stack

            ###########
            # STEP 2. #
            ###########
            # Accumulate deltas with respect to the variables in the wrt store,
            # again switching over shift/reduce decision.
            new_wrt_deltas = {}
            wrt_data = enumerate(zip(wrt, zero_jac_wrts, wrt_deltas,
                                     r_delta_wrt, s_delta_wrt))
            for i, (wrt_var, wrt_zero, accum_delta, r_delta, s_delta) in wrt_data:
                if r_delta is None and s_delta is None:
                    # Disconnected gradient.
                    continue

                # Check that tensors returned by delta functions match shape
                # expectations.
                assert r_delta is None or accum_delta.ndim == r_delta.ndim - 1, \
                    "%s %i %i" % (wrt_var.name, accum_delta.ndim, r_delta.ndim)
                assert s_delta is None or accum_delta.ndim == s_delta.ndim - 1, \
                    "%s %i %i" % (wrt_var.name, accum_delta.ndim, s_delta.ndim)

                mask_i = masks[(r_delta or s_delta).ndim - 1]
                if r_delta is None:
                    delta = T.switch(mask_i, wrt_zero, s_delta)
                elif s_delta is None:
                    delta = T.switch(mask_i, r_delta, wrt_zero)
                else:
                    delta = T.switch(mask_i, r_delta, s_delta)
                # TODO: Is this at all efficient? (Bring back GPURowSwitch?)
                delta = delta.sum(axis=0)
                # TODO: we want this to be inplace
                new_wrt_deltas[accum_delta] = accum_delta + delta

            # On shift ops, backprop the stack_bwd error onto the embedding
            # projection network / embedding parameters.
            # TODO make sparse?
            if compute_embedding_gradients:
                new_stacks[dE] = cuda.AdvancedIncSubtensor1Floats(inplace=True)(
                    new_stacks.get(dE, dE), (1. - masks[1]) * main_grad, buffer_ids_t)

            updates = dict(new_wrt_deltas.items() + new_stacks.items())
            updates = util.prepare_updates_dict(updates)

            return updates

        # TODO: These should come from forward pass -- not fixed -- in model
        # 1S, etc.
        transitions_f = T.cast(self.transitions.dimshuffle(1, 0),
                               dtype=theano.config.floatX)

        ts_f = T.cast(T.arange(1, self.seq_length + 1), dtype=theano.config.floatX)

        # Representation of buffer using embedding indices rather than values
        id_buffer = T.cast(self.X.flatten(), theano.config.floatX)
        # Build sequence of buffer pointers, where buf_ptrs[i] indicates the
        # buffer pointer values *before* computation at timestep *i* proceeds.
        # (This means we need to slice off the last actual buf_ptr output and
        # prepend a dummy.)
        buf_ptrs = T.concatenate([T.zeros((1, self.batch_size,)),
                                  self.buf_ptrs[:-1]], axis=0)

        sequences = [ts_f, transitions_f, self.stack_2_ptrs, buf_ptrs]
        outputs_info = []

        # Shared variables: Accumulated wrt deltas and bwd stacks.
        non_sequences = [dE] + wrt_vars
        non_sequences += [self.stack_bwd] + self.aux_bwd_stacks
        # More auxiliary data
        non_sequences += [id_buffer, self.final_stack]

        # More helpers (not referenced directly in code, but we need to include
        # them as non-sequences to satisfy scan strict mode)
        aux_data = [self.stack, self.buffer_t] + self.aux_stacks + self.final_aux_stacks
        aux_data += [self.X, self.transitions, self._raw_buffer_t]
        if self.use_input_dropout:
            aux_data.append(self._embedding_dropout_masks)
        aux_data += self._vs.vars.values() + extra_cost_inputs
        non_sequences += list(set(aux_data))

        bscan_ret, self.bscan_updates = theano.scan(
                step_b, sequences, outputs_info, non_sequences,
                go_backwards=True,
                n_steps=self.seq_length,
#                strict=True,
                name=self._prefix + "stack_bwd")

        self.gradients = {wrt_i: self.bscan_updates.get(wrt_var)
                          for wrt_i, wrt_var in zip(wrt, wrt_vars)}
        if compute_embedding_gradients:
            self.embedding_gradients = self.bscan_updates[dE]

Example 143

Project: aplpy
Source File: rgb.py
View license
def make_rgb_cube(files, output, north=False, system=None, equinox=None):
    '''
    Make an RGB data cube from a list of three FITS images.

    This method can read in three FITS files with different
    projections/sizes/resolutions and uses Montage to reproject
    them all to the same projection.

    Two files are produced by this function. The first is a three-dimensional
    FITS cube with a filename give by `output`, where the third dimension
    contains the different channels. The second is a two-dimensional FITS
    image with a filename given by `output` with a `_2d` suffix. This file
    contains the mean of the different channels, and is required as input to
    FITSFigure if show_rgb is subsequently used to show a color image
    generated from the FITS cube (to provide the correct WCS information to
    FITSFigure).

    Parameters
    ----------

    files : tuple or list
       A list of the filenames of three FITS filename to reproject.
       The order is red, green, blue.

    output : str
       The filename of the output RGB FITS cube.

    north : bool, optional
       By default, the FITS header generated by Montage represents the
       best fit to the images, often resulting in a slight rotation. If
       you want north to be straight up in your final mosaic, you should
       use this option.

    system : str, optional
       Specifies the system for the header (default is EQUJ).
       Possible values are: EQUJ EQUB ECLJ ECLB GAL SGAL

    equinox : str, optional
       If a coordinate system is specified, the equinox can also be given
       in the form YYYY. Default is J2000.
    '''

    # Check whether the Python montage module is installed. The Python module
    # checks itself whether the Montage command-line tools are available, and
    # if they are not then importing the Python module will fail.
    try:
        import montage_wrapper as montage
    except ImportError:
        raise Exception("Both the Montage command-line tools and the"
                        " montage-wrapper Python module are required"
                        " for this function")

    # Check that input files exist
    for f in files:
        if not os.path.exists(f):
            raise Exception("File does not exist : " + f)

    # Create work directory
    work_dir = tempfile.mkdtemp()

    raw_dir = '%s/raw' % work_dir
    final_dir = '%s/final' % work_dir

    images_raw_tbl = '%s/images_raw.tbl' % work_dir
    header_hdr = '%s/header.hdr' % work_dir

    # Create raw and final directory in work directory
    os.mkdir(raw_dir)
    os.mkdir(final_dir)

    # Create symbolic links to input files
    for i, f in enumerate(files):
        os.symlink(os.path.abspath(f), '%s/image_%i.fits' % (raw_dir, i))

    # List files and create optimal header
    montage.mImgtbl(raw_dir, images_raw_tbl, corners=True)
    montage.mMakeHdr(images_raw_tbl, header_hdr, north_aligned=north, system=system, equinox=equinox)

    # Read header in with astropy.io.fits
    header = fits.Header.fromtextfile(header_hdr)

    # Find image dimensions
    nx = int(header['NAXIS1'])
    ny = int(header['NAXIS2'])

    # Generate empty datacube
    image_cube = np.zeros((len(files), ny, nx), dtype=np.float32)

    # Loop through files
    for i in range(len(files)):

        # Reproject channel to optimal header
        montage.reproject('%s/image_%i.fits' % (raw_dir, i),
                          '%s/image_%i.fits' % (final_dir, i),
                          header=header_hdr, exact_size=True, bitpix=-32)

        # Read in and add to datacube
        image_cube[i, :, :] = fits.getdata('%s/image_%i.fits' % (final_dir, i))

    # Write out final cube
    fits.writeto(output, image_cube, header, clobber=True)

    # Write out collapsed version of cube
    fits.writeto(output.replace('.fits', '_2d.fits'), \
                   np.mean(image_cube, axis=0), header, clobber=True)

    # Remove work directory
    shutil.rmtree(work_dir)

Example 144

Project: pyvision
Source File: img.py
View license
    def __init__(self,data,bw_annotate=False):
        '''
        Create an image from a file or a PIL Image, OpenCV Image, or numpy array.
         
        @param data: this can be a numpy array, PIL image, or opencv image.
        @param bw_annotate: generate a black and white image to make color annotations show up better
        @return: an Image object instance
        '''

        self.filename = None
        self.pil = None
        self.matrix2d = None
        self.matrix3d = None
        self.opencv = None
        self.opencv2 = None
        self.opencv2bw = None
        self.annotated = None
        self.bw_annotate = bw_annotate
        
        # Convert floating point ipl images to numpy arrays
        if isinstance(data,cv.iplimage) and data.nChannels == 3 and data.depth == 32:
            w,h = cv.GetSize(data)
            data = np.frombuffer(data.tostring(),dtype=np.float32)
            data.shape = (h,w,3)
            data = data.transpose((2,1,0))
            data = data[::-1,:,:]
            
        # Convert floating point ipl images to numpy arrays
        if isinstance(data,cv.iplimage) and data.nChannels == 1 and data.depth == 32:
            w,h = cv.GetSize(data)
            data = np.frombuffer(data.tostring(),dtype=np.float32)
            data.shape = (h,w)
            data = data.transpose((2,1,0))
            data = data[::-1,:,:]
            
        # Numpy format
        if isinstance(data,numpy.ndarray) and len(data.shape) == 2 and data.dtype != np.uint8:
            self.type=TYPE_MATRIX_2D
            self.matrix2d = data
            
            self.width,self.height = self.matrix2d.shape
            self.channels = 1
            
            if self.matrix2d.dtype == numpy.float32:
                self.depth=32
            elif self.matrix2d.dtype == numpy.float64:
                self.depth=64
            else:
                raise TypeError("Unsuppoted format for ndarray images: %s"%self.matrix2d.dtype)
            
        # OpenCV2 gray scale format
        elif isinstance(data,numpy.ndarray) and len(data.shape) == 2 and data.dtype == np.uint8:
            self.type=TYPE_OPENCV2BW
            self.opencv2bw = data
            
            self.height,self.width = self.opencv2bw.shape
            self.channels = 1
            self.depth=8
            
        # Numpy color format    
        elif isinstance(data,numpy.ndarray) and len(data.shape) == 3 and data.shape[0]==3 and data.dtype != np.uint8:
            self.type=TYPE_MATRIX_RGB
            self.matrix3d = data
            self.channels=3
            self.width = self.matrix3d.shape[1]
            self.height = self.matrix3d.shape[2]
            # set the types
            if self.matrix3d.dtype == numpy.float32:
                self.depth=32
            elif self.matrix3d.dtype == numpy.float64:
                self.depth=64
            else:
                raise TypeError("Unsuppoted format for ndarray images: %s"%self.matrix2d.dtype)
            
        # OpenCV2 color format    
        elif isinstance(data,numpy.ndarray) and len(data.shape) == 3 and data.shape[2]==3 and data.dtype == np.uint8:
            self.type=TYPE_OPENCV2
            self.opencv2 = data
            self.channels=3
            self.width = self.opencv2.shape[1]
            self.height = self.opencv2.shape[0]
            self.depth=8
            
        # Load as a pil image    
        elif isinstance(data,PIL.Image.Image) or type(data) == str:
            if type(data) == str:
                # Assume this is a filename
                # TODO: Removing the filename causes errors in other unittest.
                #       Those errors should be corrected.
                self.filename = data
                data = PIL.Image.open(data)
                    
            self.type=TYPE_PIL
            self.pil = data
            self.width,self.height = self.pil.size
                        
            if self.pil.mode == 'L':
                self.channels = 1
            elif self.pil.mode == 'RGB':
                self.channels = 3
            #elif self.pil.mode == 'RGBA':
                # 
            #    self.pil = self.pil.convert('RGB')
            #    self.channels = 3
            else:
                self.pil.convert('RGB')
                self.channels = 3
                #   raise TypeError("Unsuppoted format for PIL images: %s"%self.pil.mode)
            
            self.depth = 8
        
        # opencv format             
        elif isinstance(data,cv.iplimage):
            self.type=TYPE_OPENCV
            self.opencv=data 
            
            self.width = data.width
            self.height = data.height
                        
            assert data.nChannels in (1,3)
            self.channels = data.nChannels 
            
            assert data.depth == 8
            self.depth = data.depth   

        # unknown type
        else:
            raise TypeError("Could not create from type: %s %s"%(data,type(data)))
        
        self.size = (self.width,self.height)
        self.data = data

Example 145

Project: visvis
Source File: text_freetype.py
View license
    def Compile(self, textObject):
        
        # Get text string with escaped text converted to Unicode
        tt = self.ConvertEscapedText(textObject.text)
        
        # Init arrays
        vertices = np.zeros((len(tt)*4,3), dtype=np.float32)
        texcoords = np.zeros((len(tt)*4,2), dtype=np.float32)
        
        # Prepare
        textObject._shift = None
        pen = [0,0]
        prev = None
        
        # Calculate font size
        # note: I can also imagine doing it the other way around; i.e. textSize 
        # becomes as FreeType sees it, and we scale the fonts in the 
        # prerendered text renderer. We'd have to change all the 
        # uses fontSize though. So this is simply easier.
        fontSize = textObject.fontSize * 1.4 * TEX_SCALE
        fig = textObject.GetFigure()
        if fig:
            fontSize *= fig._relativeFontSize
        
        # Store tex_scale, integer fontsize and residu
        textObject._tex_scale = TEX_SCALE
        textObject._actualFontSize = int(round(fontSize)) 
        
        fonts = [(self.GetFont(textObject.fontName, textObject._actualFontSize), 0, False, False)]
        font, voffset, bold, italic = fonts[-1]
        escaped = False
        for i,charcode in enumerate(tt):
            if not escaped:
                if charcode == '_':
                    font = self.GetFont(textObject.fontName, font.size*0.7, bold, italic)
                    voffset -= 0.1*font.size
                    prev = None  # Disable kerning
                    continue
                elif charcode == '^':
                    font = self.GetFont(textObject.fontName, font.size*0.7, bold, italic)
                    voffset += 0.5*font.size
                    prev = None
                    continue
                elif charcode == '\x06':
                    italic = True
                    font = self.GetFont(textObject.fontName, font.size, bold, italic)
                    continue
                elif charcode == '\x07':
                    bold = True
                    font = self.GetFont(textObject.fontName, font.size, bold, italic)
                    continue
                elif charcode == '{':
                    fonts.append((font, voffset, bold, italic))
                    continue
                elif charcode == '}':
                    if len(fonts) > 1:
                        fonts.pop()
                        font, voffset, bold, italic = fonts[-1]
                        continue
                elif charcode == '\\':
                    if i < len(tt)-1 and tt[i+1] in '_^{}':
                        escaped = True
                        continue
            glyph = font[charcode]
            if glyph is None:
                continue # Character not available
            kerning = glyph.get_kerning(prev)
            x0 = pen[0] + glyph.offset[0] + kerning
            y0 = pen[1] + glyph.offset[1] + voffset
            x1 = x0 + glyph.size[0]
            y1 = y0 - glyph.size[1]
            u0 = glyph.texcoords[0]
            v0 = glyph.texcoords[1]
            u1 = glyph.texcoords[2]
            v1 = glyph.texcoords[3]

            index     = i*4
            _vertices  = [[x0,y0,1],[x0,y1,1],[x1,y1,1], [x1,y0,1]]
            _texcoords = [[u0,v0],[u0,v1],[u1,v1], [u1,v0]]

            vertices[i*4:i*4+4] = _vertices
            texcoords[i*4:i*4+4] = _texcoords
            pen[0] = pen[0]+glyph.advance[0]/64.0 + kerning
            pen[1] = pen[1]+glyph.advance[1]/64.0
            prev = charcode
            font, voffset, bold, italic = fonts[-1]
            escaped = False
        
        # Flip and shift vertices
        vertices /= TEX_SCALE
        vertices *= (1, -1, 1)
        vertices += (0, font.ascender/TEX_SCALE, 0)
        
        # Store width
        if False: #glyph is not None: # Not used
            textObject.width = pen[0]-glyph.advance[0]/64.0+glyph.size[0] if tt else 0
        
        # Update dynamic texture
        self.atlas.upload()
        
        # Store data. 
        textObject._SetCompiledData(vertices, texcoords)

Example 146

Project: tf_resnet_cifar
Source File: main_train.py
View license
def train_and_val():
  with tf.Graph().as_default():
    # train/test phase indicator
    phase_train = tf.placeholder(tf.bool, name='phase_train')

    # learning rate is manually set
    learning_rate = tf.placeholder(tf.float32, name='learning_rate')
    tf.scalar_summary('learning_rate', learning_rate)

    # global step
    global_step = tf.Variable(0, trainable=False, name='global_step')

    # train/test inputs
    train_image_batch, train_label_batch = m.make_train_batch(
      FLAGS.train_tf_path, FLAGS.train_batch_size)
    val_image_batch, val_label_batch = m.make_validation_batch(
      FLAGS.val_tf_path, FLAGS.val_batch_size)
    image_batch, label_batch = control_flow_ops.cond(phase_train,
      lambda: (train_image_batch, train_label_batch),
      lambda: (val_image_batch, val_label_batch))

    # model outputs
    logits = m.residual_net(image_batch, FLAGS.residual_net_n, 10, phase_train)

    # total loss
    loss = m.loss(logits, label_batch)
    accuracy = m.accuracy(logits, label_batch)
    tf.scalar_summary('train_loss', loss)
    tf.scalar_summary('train_accuracy', accuracy)

    # train one step
    train_op = m.train_op(loss, global_step, learning_rate)

    # saver
    saver = tf.train.Saver(tf.all_variables())

    # start session
    sess = tf.Session(config=tf.ConfigProto(log_device_placement=False))

    # summary writer
    summary_op = tf.merge_all_summaries()
    summary_writer = tf.train.SummaryWriter(
      FLAGS.log_dir, graph=sess.graph)

    # initialize parameters or load from a checkpoint
    if FLAGS.load_dir != '':
      # load from checkpoint
      checkpoint = tf.train.get_checkpoint_state(FLAGS.load_dir)
      model_checkpoint_path = checkpoint.model_checkpoint_path
      if checkpoint and model_checkpoint_path:
        saver.restore(sess, model_checkpoint_path)
        print('Model restored from %s' % model_checkpoint_path)
      else:
        raise 'Load directory provided by no checkpoint found'
    else:
      init_op = tf.initialize_all_variables()
      print('Initializing...')
      sess.run(init_op, {phase_train.name: True})

    print('Start training...')
    # train loop
    tf.train.start_queue_runners(sess=sess)
    curr_lr = 0.0
    for step in xrange(FLAGS.max_steps):
      # # set learning rate manually
      # if step <= 5000:
      #   _lr = 1e-2
      # elif step <= 32000:
      #   _lr = 1e-1
      # elif step <= 48000:
      #   _lr = 1e-2
      # else:
      #   _lr = 1e-3
      # set learning rate manually
      if step <= 48000:
        _lr = 1e-2
      else:
        _lr = 1e-3
      if curr_lr != _lr:
        curr_lr = _lr
        print('Learning rate set to %f' % curr_lr)

      # train
      fetches = [train_op, loss]
      if step > 0 and step % FLAGS.summary_interval == 0:
        fetches += [accuracy, summary_op]
      sess_outputs = sess.run(
        fetches, {phase_train.name: True, learning_rate.name: curr_lr})

      # summary
      if step > 0 and step % FLAGS.summary_interval == 0:
        train_loss_value, train_acc_value, summary_str = sess_outputs[1:]
        print('[%s] Iteration %d, train loss = %f, train accuracy = %f' %
            (datetime.now(), step, train_loss_value, train_acc_value))
        summary_writer.add_summary(summary_str, step)

      # validation
      if step > 0 and step % FLAGS.val_interval == 0:
        print('Evaluating...')
        n_val_samples = 10000
        val_batch_size = FLAGS.val_batch_size
        n_val_batch = int(n_val_samples / val_batch_size)
        val_logits = np.zeros((n_val_samples, 10), dtype=np.float32)
        val_labels = np.zeros((n_val_samples), dtype=np.int64)
        val_losses = []
        for i in xrange(n_val_batch):
          fetches = [logits, label_batch, loss]
          session_outputs = sess.run(
            fetches, {phase_train.name: False})
          val_logits[i*val_batch_size:(i+1)*val_batch_size, :] = session_outputs[0]
          val_labels[i*val_batch_size:(i+1)*val_batch_size] = session_outputs[1]
          val_losses.append(session_outputs[2])
        pred_labels = np.argmax(val_logits, axis=1)
        val_accuracy = np.count_nonzero(
          pred_labels == val_labels) / n_val_samples
        val_loss = float(np.mean(np.asarray(val_losses)))
        print('Test accuracy = %f' % val_accuracy)
        val_summary = tf.Summary()
        val_summary.value.add(tag='val_accuracy', simple_value=val_accuracy)
        val_summary.value.add(tag='val_loss', simple_value=val_loss)
        summary_writer.add_summary(val_summary, step)

      # save variables
      if step > 0 and step % FLAGS.save_interval == 0:
        checkpoint_path = os.path.join(FLAGS.log_dir, 'checkpoint')
        saver.save(sess, checkpoint_path, global_step=step)
        print('Checkpoint saved at %s' % checkpoint_path)

Example 147

View license
    def setup_train(self):

        # dimensions: (batch, time, notes, input_data) with input_data as in architecture
        self.input_mat = T.btensor4()
        # dimensions: (batch, time, notes, onOrArtic) with 0:on, 1:artic
        self.output_mat = T.btensor4()
        
        self.epsilon = np.spacing(np.float32(1.0))

        def step_time(in_data, *other):
            other = list(other)
            split = -len(self.t_layer_sizes) if self.dropout else len(other)
            hiddens = other[:split]
            masks = [None] + other[split:] if self.dropout else []
            new_states = self.time_model.forward(in_data, prev_hiddens=hiddens, dropout=masks)
            return new_states
        
        def step_note(in_data, *other):
            other = list(other)
            split = -len(self.p_layer_sizes) if self.dropout else len(other)
            hiddens = other[:split]
            masks = [None] + other[split:] if self.dropout else []
            new_states = self.pitch_model.forward(in_data, prev_hiddens=hiddens, dropout=masks)
            return new_states
        
        # We generate an output for each input, so it doesn't make sense to use the last output as an input.
        # Note that we assume the sentinel start value is already present
        # TEMP CHANGE: NO SENTINEL
        input_slice = self.input_mat[:,0:-1]
        n_batch, n_time, n_note, n_ipn = input_slice.shape
        
        # time_inputs is a matrix (time, batch/note, input_per_note)
        time_inputs = input_slice.transpose((1,0,2,3)).reshape((n_time,n_batch*n_note,n_ipn))
        num_time_parallel = time_inputs.shape[1]
        
        # apply dropout
        if self.dropout > 0:
            time_masks = theano_lstm.MultiDropout( [(num_time_parallel, shape) for shape in self.t_layer_sizes], self.dropout)
        else:
            time_masks = []

        time_outputs_info = [initial_state_with_taps(layer, num_time_parallel) for layer in self.time_model.layers]
        time_result, _ = theano.scan(fn=step_time, sequences=[time_inputs], non_sequences=time_masks, outputs_info=time_outputs_info)
        
        self.time_thoughts = time_result
        
        # Now time_result is a list of matrix [layer](time, batch/note, hidden_states) for each layer but we only care about 
        # the hidden state of the last layer.
        # Transpose to be (note, batch/time, hidden_states)
        last_layer = get_last_layer(time_result)
        n_hidden = last_layer.shape[2]
        time_final = get_last_layer(time_result).reshape((n_time,n_batch,n_note,n_hidden)).transpose((2,1,0,3)).reshape((n_note,n_batch*n_time,n_hidden))
        
        # note_choices_inputs represents the last chosen note. Starts with [0,0], doesn't include last note.
        # In (note, batch/time, 2) format
        # Shape of start is thus (1, N, 2), concatenated with all but last element of output_mat transformed to (x, N, 2)
        start_note_values = T.alloc(np.array(0,dtype=np.int8), 1, time_final.shape[1], 2 )
        correct_choices = self.output_mat[:,1:,0:-1,:].transpose((2,0,1,3)).reshape((n_note-1,n_batch*n_time,2))
        note_choices_inputs = T.concatenate([start_note_values, correct_choices], axis=0)
        
        # Together, this and the output from the last LSTM goes to the new LSTM, but rotated, so that the batches in
        # one direction are the steps in the other, and vice versa.
        note_inputs = T.concatenate( [time_final, note_choices_inputs], axis=2 )
        num_timebatch = note_inputs.shape[1]
        
        # apply dropout
        if self.dropout > 0:
            pitch_masks = theano_lstm.MultiDropout( [(num_timebatch, shape) for shape in self.p_layer_sizes], self.dropout)
        else:
            pitch_masks = []

        note_outputs_info = [initial_state_with_taps(layer, num_timebatch) for layer in self.pitch_model.layers]
        note_result, _ = theano.scan(fn=step_note, sequences=[note_inputs], non_sequences=pitch_masks, outputs_info=note_outputs_info)
        
        self.note_thoughts = note_result
        
        # Now note_result is a list of matrix [layer](note, batch/time, onOrArticProb) for each layer but we only care about 
        # the hidden state of the last layer.
        # Transpose to be (batch, time, note, onOrArticProb)
        note_final = get_last_layer(note_result).reshape((n_note,n_batch,n_time,2)).transpose(1,2,0,3)
        
        # The cost of the entire procedure is the negative log likelihood of the events all happening.
        # For the purposes of training, if the ouputted probability is P, then the likelihood of seeing a 1 is P, and
        # the likelihood of seeing 0 is (1-P). So the likelihood is (1-P)(1-x) + Px = 2Px - P - x + 1
        # Since they are all binary decisions, and are all probabilities given all previous decisions, we can just
        # multiply the likelihoods, or, since we are logging them, add the logs.
        
        # Note that we mask out the articulations for those notes that aren't played, because it doesn't matter
        # whether or not those are articulated.
        # The padright is there because self.output_mat[:,:,:,0] -> 3D matrix with (b,x,y), but we need 3d tensor with 
        # (b,x,y,1) instead
        active_notes = T.shape_padright(self.output_mat[:,1:,:,0])
        mask = T.concatenate([T.ones_like(active_notes),active_notes], axis=3)
        
        loglikelihoods = mask * T.log( 2*note_final*self.output_mat[:,1:] - note_final - self.output_mat[:,1:] + 1 + self.epsilon )
        self.cost = T.neg(T.sum(loglikelihoods))
        
        updates, _, _, _, _ = create_optimization_updates(self.cost, self.params, method="adadelta")
        self.update_fun = theano.function(
            inputs=[self.input_mat, self.output_mat],
            outputs=self.cost,
            updates=updates,
            allow_input_downcast=True)

        self.update_thought_fun = theano.function(
            inputs=[self.input_mat, self.output_mat],
            outputs= ensure_list(self.time_thoughts) + ensure_list(self.note_thoughts) + [self.cost],
            allow_input_downcast=True)

Example 148

Project: chainer
Source File: n_step_lstm.py
View license
    def forward(self, inputs):
        (hx, cx), inputs = _split(inputs, 2)
        ws, inputs = _split(inputs, self.n_layers * 8)
        bs, inputs = _split(inputs, self.n_layers * 8)
        x_list = inputs

        hx = cuda.cupy.ascontiguousarray(hx)
        cx = cuda.cupy.ascontiguousarray(cx)

        x_desc = cudnn.create_tensor_nd_descriptor(x_list[0][..., None])

        length = len(x_list)
        n_units = hx.shape[2]

        xs = cuda.cupy.concatenate(x_list, axis=0)
        ys = cuda.cupy.empty((len(xs), n_units), dtype=xs.dtype)

        handle = cudnn.get_handle()
        self.handle = handle

        rnn_desc = cudnn.create_rnn_descriptor(
            n_units, self.n_layers, self.states.desc,
            libcudnn.CUDNN_LINEAR_INPUT, libcudnn.CUDNN_UNIDIRECTIONAL,
            libcudnn.CUDNN_LSTM, libcudnn.CUDNN_DATA_FLOAT)
        self.rnn_desc = rnn_desc

        c_x_descs = _make_tensor_descriptor_array(x_list)
        hx_desc = cudnn.create_tensor_nd_descriptor(hx)
        cx_desc = cudnn.create_tensor_nd_descriptor(cx)

        weights_size = libcudnn.getRNNParamsSize(
            handle, rnn_desc.value, x_desc.value, libcudnn.CUDNN_DATA_FLOAT)
        w = cuda.cupy.empty((weights_size // 4, 1, 1), dtype=numpy.float32)
        w_desc = cudnn.create_filter_descriptor(w)

        for layer in six.moves.range(self.n_layers):
            for lin_layer_id in six.moves.range(8):
                mat = cudnn.get_rnn_lin_layer_matrix_params(
                    handle, rnn_desc, layer, x_desc, w_desc, w,
                    lin_layer_id)
                m = mat.reshape(mat.size)
                m[...] = ws[layer * 8 + lin_layer_id].ravel()
                bias = cudnn.get_rnn_lin_layer_bias_params(
                    handle, rnn_desc, layer, x_desc, w_desc, w,
                    lin_layer_id)
                b = bias.reshape(bias.size)
                b[...] = bs[layer * 8 + lin_layer_id]
        self.w = w
        self.w_desc = w_desc

        sections = numpy.cumsum([len(x) for x in x_list[:-1]])
        y_list = cuda.cupy.split(ys, sections)

        c_y_descs = _make_tensor_descriptor_array(y_list)
        hy = cuda.cupy.empty_like(hx)
        cy = cuda.cupy.empty_like(cx)
        hy_desc = cudnn.create_tensor_nd_descriptor(hy)
        cy_desc = cudnn.create_tensor_nd_descriptor(cy)

        work_size = libcudnn.getRNNWorkspaceSize(
            handle, rnn_desc.value, length, c_x_descs.data)
        workspace = cuda.cupy.empty((work_size,), dtype='b')
        self.workspace = workspace

        if not self.train:
            libcudnn.RNNForwardInference(
                handle, rnn_desc.value, length,
                c_x_descs.data, xs.data.ptr, hx_desc.value, hx.data.ptr,
                cx_desc.value, cx.data.ptr, w_desc.value, w.data.ptr,
                c_y_descs.data, ys.data.ptr, hy_desc.value, hy.data.ptr,
                cy_desc.value, cy.data.ptr, workspace.data.ptr, work_size)

        else:
            reserve_size = libcudnn.getRNNTrainingReserveSize(
                handle, rnn_desc.value, length, c_x_descs.data)
            self.reserve_space = cuda.cupy.empty((reserve_size,), dtype='b')
            libcudnn.RNNForwardTraining(
                handle, rnn_desc.value, length,
                c_x_descs.data, xs.data.ptr, hx_desc.value, hx.data.ptr,
                cx_desc.value, cx.data.ptr, w_desc.value, w.data.ptr,
                c_y_descs.data, ys.data.ptr, hy_desc.value, hy.data.ptr,
                cy_desc.value, cy.data.ptr,
                workspace.data.ptr, work_size,
                self.reserve_space.data.ptr, reserve_size)

        self.c_y_descs = c_y_descs
        self.ys = ys
        self.c_x_descs = c_x_descs

        return tuple([hy, cy] + y_list)

Example 149

Project: lowrank-gru
Source File: RNN_memory.py
View license
def main(datafile, n_iter, n_batch, n_hidden, time_steps, learning_rate, savefile, scale_penalty, use_scale,
         model, n_hidden_lstm, loss_function, cost_every_t, n_gru_lr_proj, initial_b_u, load_iter):


    # --- Manage data --------------------
    #f = file('/u/shahamar/complex_RNN/trainingRNNs-master/memory_data_300.pkl', 'rb')
    #f = file('./memory_data_20.pkl', 'rb')
    #f = file('./memory_data_500.pkl', 'rb')
    f = file(datafile, 'rb')
    dict = cPickle.load(f)
    f.close()


    train_x = dict['train_x'] 
    train_y = dict['train_y']
    test_x = dict['test_x'] 
    test_y = dict['test_y'] 


    #import pdb; pdb.set_trace()
    #cPickle.dump(dict, file('/u/shahamar/complex_RNN/trainingRNNs-master/memory_data.pkl', 'wb'))

    #import pdb; pdb.set_trace()

    n_train = train_x.shape[1]
    n_test = test_x.shape[1]
    n_input = train_x.shape[2]
    n_output = train_y.shape[2]

    num_batches = int(n_train / n_batch)
    time_steps = train_x.shape[0]
    
   #######################################################################

    gradient_clipping = np.float32(1)

    if (model == 'LSTM'):   
        inputs, parameters, costs = LSTM(n_input, n_hidden_lstm, n_output,
                                         out_every_t=cost_every_t, loss_function=loss_function)
        gradients = T.grad(costs[0], parameters)
        gradients = [T.clip(g, -gradient_clipping, gradient_clipping) for g in gradients]
    
    elif (model == 'GRU'): 
        inputs, parameters, costs = GRU(n_input, n_hidden_lstm, n_output, 
                                         out_every_t=cost_every_t, loss_function=loss_function, initial_b_u=initial_b_u)
        gradients = T.grad(costs[0], parameters)
        gradients = [T.clip(g, -gradient_clipping, gradient_clipping) for g in gradients]
   
    elif (model == 'GRU_LR'):
        inputs, parameters, costs = GRU_LR(n_input, n_hidden_lstm, n_output, n_gru_lr_proj,
                                         out_every_t=cost_every_t, loss_function=loss_function, initial_b_u=initial_b_u)
        gradients = T.grad(costs[0], parameters)
        #gradients = [T.clip(g, -gradient_clipping, gradient_clipping) for g in gradients]
        gradients = clip_gradients_norm(gradients, gradient_clipping, parameters, fix_nan = True)
    
    elif (model == 'GRU_LRDiag'):
        inputs, parameters, costs = GRU_LRDiag(n_input, n_hidden_lstm, n_output, n_gru_lr_proj,
                                         out_every_t=cost_every_t, loss_function=loss_function, initial_b_u=initial_b_u)
        gradients = T.grad(costs[0], parameters)
        #gradients = [T.clip(g, -gradient_clipping, gradient_clipping) for g in gradients]
        gradients = clip_gradients_norm(gradients, gradient_clipping, parameters, fix_nan = True)

    elif (model == 'complex_RNN'):
        inputs, parameters, costs = complex_RNN(n_input, n_hidden, n_output, scale_penalty,
                                                out_every_t=cost_every_t, loss_function=loss_function)
        if use_scale is False:
            parameters.pop()
        gradients = T.grad(costs[0], parameters)

    elif (model == 'complex_RNN_LSTM'):
        inputs, parameters, costs = complex_RNN_LSTM(n_input, n_hidden, n_hidden_lstm, n_output, scale_penalty,
                                                     out_every_t=cost_every_t, loss_function=loss_function)

    elif (model == 'IRNN'):
        inputs, parameters, costs = IRNN(n_input, n_hidden, n_output,
                                         out_every_t=cost_every_t, loss_function=loss_function)
        gradients = T.grad(costs[0], parameters)
        gradients = [T.clip(g, -gradient_clipping, gradient_clipping) for g in gradients]

    elif (model == 'RNN'):
        inputs, parameters, costs = tanhRNN(n_input, n_hidden, n_output,
                                            out_every_t=cost_every_t, loss_function=loss_function)
        gradients = T.grad(costs[0], parameters)
        gradients = [T.clip(g, -gradient_clipping, gradient_clipping) for g in gradients]

    else:
        print >> sys.stderr, "Unsuported model:", model
        return
 

   




    s_train_x = theano.shared(train_x)
    s_train_y = theano.shared(train_y)

    s_test_x = theano.shared(test_x)
    s_test_y = theano.shared(test_y)


    # --- Compile theano functions --------------------------------------------------

    index = T.iscalar('i')

    #updates, rmsprop = rms_prop(learning_rate, parameters, gradients)
    updates, rmsprop = rms_prop(learning_rate, parameters, gradients, epsilon=1e-8)

    givens = {inputs[0] : s_train_x[:, n_batch * index : n_batch * (index + 1), :],
              inputs[1] : s_train_y[:, n_batch * index : n_batch * (index + 1), :]}

    givens_test = {inputs[0] : s_test_x,
                   inputs[1] : s_test_y}
    
   
    
    #train = theano.function([index], costs[0], givens=givens, updates=updates)
    train = theano.function([index], [costs[0], costs[2]], givens=givens, updates=updates)
    #test = theano.function([], costs[1], givens=givens_test)
    test = theano.function([], [costs[1], costs[2]], givens=givens_test)

    # --- Training Loop ---------------------------------------------------------------

    # f1 = file('/data/lisatmp3/shahamar/memory/RNN_100.pkl', 'rb')
    # data1 = cPickle.load(f1)
    # f1.close()
    # train_loss = data1['train_loss']
    # test_loss = data1['test_loss']
    # best_params = data1['best_params']
    # best_test_loss = data1['best_test_loss']

    # for i in xrange(len(parameters)):
    #     parameters[i].set_value(data1['parameters'][i])

    # for i in xrange(len(parameters)):
    #     rmsprop[i].set_value(data1['rmsprop'][i])

    corrected_accuracy = lambda a: 4.0 / 5.0 * (time_steps * (a - 1.0) + 20.0 * a - 10.0)

    train_loss = []
    test_loss = []
    best_params = [p.get_value() for p in parameters]
    best_test_loss = 1e6

    if load_iter > 0:
	loaded_data = cPickle.load(file(savefile, 'rb'))
        loaded_parameters = loaded_data['parameters']
	for j, p in enumerate(parameters):
        	p.set_value(loaded_parameters[j])
        loaded_rmsprop = loaded_data['rmsprop']
        for j, r in enumerate(rmsprop):
		r.set_value(loaded_rmsprop[j])
        train_loss = loaded_data['train_loss']
	test_loss = loaded_data['test_loss']
        best_params = loaded_data['best_params']
        best_test_loss =  loaded_data['best_test_loss']
        model = loaded_data['model']
        time_steps = loaded_data['time_steps']

    for i in xrange(load_iter, n_iter):
#        start_time = timeit.default_timer()
     #   pdb.set_trace()

        if (n_iter % int(num_batches) == 0):
            #import pdb; pdb.set_trace()
            inds = np.random.permutation(int(n_train))
            data_x = s_train_x.get_value()
            s_train_x.set_value(data_x[:,inds,:])
            data_y = s_train_y.get_value()
            s_train_y.set_value(data_y[:,inds,:])


        mse, accuracy = train(i % int(num_batches))
        train_loss.append(mse)
        print >> sys.stderr, "Iteration:", i
        #print >> sys.stderr, "mse:", mse
        print >> sys.stderr, "ce:", mse, "accuracy:", accuracy, "corrected accuracy:", corrected_accuracy(accuracy), "/ 8"
        print >> sys.stderr, ''

        #if (i % 50==0):
        if (i % 500==0):
            mse, accuracy = test()
            print >> sys.stderr, ''
            print >> sys.stderr, "TEST"
            #print >> sys.stderr, "mse:", mse
            print >> sys.stderr, "ce:", mse, "accuracy:", accuracy, "corrected accuracy:", corrected_accuracy(accuracy), " / 8"
            print >> sys.stderr, ''
            test_loss.append(mse)

            if mse < best_test_loss:
                best_params = [p.get_value() for p in parameters]
                best_test_loss = mse
                print >> sys.stderr, "BEST!"

            save_vals = {'parameters': [p.get_value() for p in parameters],
                         'rmsprop': [r.get_value() for r in rmsprop],
                         'train_loss': train_loss,
                         'test_loss': test_loss,
                         'best_params': best_params,
                         'best_test_loss': best_test_loss,
                         'model': model,
                         'time_steps': time_steps}

            cPickle.dump(save_vals,
                         file(savefile, 'wb'),
                         cPickle.HIGHEST_PROTOCOL)

Example 150

Project: scikit-learn
Source File: olivetti_faces.py
View license
def fetch_olivetti_faces(data_home=None, shuffle=False, random_state=0,
                         download_if_missing=True):
    """Loader for the Olivetti faces data-set from AT&T.

    Read more in the :ref:`User Guide <olivetti_faces>`.

    Parameters
    ----------
    data_home : optional, default: None
        Specify another download and cache folder for the datasets. By default
        all scikit learn data is stored in '~/scikit_learn_data' subfolders.

    shuffle : boolean, optional
        If True the order of the dataset is shuffled to avoid having
        images of the same person grouped.

    download_if_missing : optional, True by default
        If False, raise a IOError if the data is not locally available
        instead of trying to download the data from the source site.

    random_state : optional, integer or RandomState object
        The seed or the random number generator used to shuffle the
        data.

    Returns
    -------
    An object with the following attributes:

    data : numpy array of shape (400, 4096)
        Each row corresponds to a ravelled face image of original size 64 x 64 pixels.

    images : numpy array of shape (400, 64, 64)
        Each row is a face image corresponding to one of the 40 subjects of the dataset.

    target : numpy array of shape (400, )
        Labels associated to each face image. Those labels are ranging from
        0-39 and correspond to the Subject IDs.

    DESCR : string
        Description of the modified Olivetti Faces Dataset.

    Notes
    ------

    This dataset consists of 10 pictures each of 40 individuals. The original
    database was available from (now defunct)

        http://www.cl.cam.ac.uk/research/dtg/attarchive/facedatabase.html

    The version retrieved here comes in MATLAB format from the personal
    web page of Sam Roweis:

        http://www.cs.nyu.edu/~roweis/

    """
    data_home = get_data_home(data_home=data_home)
    if not exists(data_home):
        makedirs(data_home)
    filepath = _pkl_filepath(data_home, TARGET_FILENAME)
    if not exists(filepath):
        print('downloading Olivetti faces from %s to %s'
              % (DATA_URL, data_home))
        fhandle = urlopen(DATA_URL)
        buf = BytesIO(fhandle.read())
        mfile = loadmat(buf)
        faces = mfile['faces'].T.copy()
        joblib.dump(faces, filepath, compress=6)
        del mfile
    else:
        faces = joblib.load(filepath)
    # We want floating point data, but float32 is enough (there is only
    # one byte of precision in the original uint8s anyway)
    faces = np.float32(faces)
    faces = faces - faces.min()
    faces /= faces.max()
    faces = faces.reshape((400, 64, 64)).transpose(0, 2, 1)
    # 10 images per class, 400 images total, each class is contiguous.
    target = np.array([i // 10 for i in range(400)])
    if shuffle:
        random_state = check_random_state(random_state)
        order = random_state.permutation(len(faces))
        faces = faces[order]
        target = target[order]
    return Bunch(data=faces.reshape(len(faces), -1),
                 images=faces,
                 target=target,
                 DESCR=MODULE_DOCS)