def apply(self, input_, application_call):
     if self._training_mode:
         mean, stdev = self._compute_training_statistics(input_)
     else:
         mean, stdev = self._prepare_population_statistics()
     # Useful for filtration of calls that were already made in
     # training mode when doing graph transformations.
     # Very important to cast to bool, as self._training_mode is
     # normally a list (to support nested context managers), which would
     # otherwise get passed by reference and be remotely mutated.
     application_call.metadata['training_mode'] = bool(self._training_mode)
     # Useful for retrieving a list of updates for population
     # statistics. Ditch the broadcastable first axis, though, to
     # make it the same dimensions as the population mean and stdev
     # shared variables.
     application_call.metadata['offset'] = mean[0]
     application_call.metadata['divisor'] = stdev[0]
     # Give these quantities roles in the graph.
     _add_role_and_annotate(mean, BATCH_NORM_OFFSET,
                            [self, application_call])
     _add_role_and_annotate(stdev, BATCH_NORM_DIVISOR,
                            [self, application_call])
     scale = _add_batch_axis(self.scale)
     shift = _add_batch_axis(self.shift)
     # Heavy lifting is done by the Theano utility function.
     normalized = bn.batch_normalization(input_, scale, shift, mean, stdev,
                                         mode=('low_mem'
                                               if self.conserve_memory
                                               else 'high_mem'))
     return normalized

 def conv_bn(inputs, gamma, beta, mean, std):
     return bn.batch_normalization(inputs,
                                   gamma.dimshuffle('x', 0, 'x', 'x'),
                                   beta.dimshuffle('x', 0, 'x', 'x'),
                                   mean.dimshuffle('x', 0, 'x', 'x'),
                                   std.dimshuffle('x', 0, 'x', 'x'),
                                   mode=mode)

def test_bn():
    def bn_ref(x, G, B, M, V):
        n = (x - M) / V
        return n * G + B

    numpy.random.seed(1234)
    X = 1 + numpy.random.random([10, 20]).astype("float32")
    B = 1 + numpy.random.random([20]).astype("float32")
    G = 1 + numpy.random.random([20]).astype("float32")
    M = 1 + numpy.random.random([20]).astype("float32")
    V = 1 + numpy.random.random([20]).astype("float32")

    x = theano.tensor.matrix("x")
    b = theano.tensor.vector("b")
    g = theano.tensor.vector("g")
    m = theano.tensor.vector("m")
    v = theano.tensor.vector("v")

    bn_ref_op = bn_ref(x, g, b, m, v)
    f_ref = theano.function([x, b, g, m, v], [bn_ref_op])
    res_ref = f_ref(X, G, B, M, V)
    for mode in ["low_mem", "high_mem"]:
        bn_op = batch_normalization(x, g, b, m, v, mode=mode)
        f = theano.function([x, b, g, m, v], [bn_op])
        res = f(X, G, B, M, V)
        utt.assert_allclose(res_ref, res)

        def bn(inputs, gamma, beta, mean, std):
            return batch_normalization(inputs, gamma, beta, mean, std, mode=mode)

        utt.verify_grad(bn, [X, G, B, M, V])

    bn_ref_op = bn_ref(x, g, b, x.mean(axis=0, keepdims=True), x.std(axis=0, keepdims=True))
    f_ref = theano.function([x, b, g], [bn_ref_op])
    res_ref = f_ref(X, G, B)
    for mode in ["low_mem", "high_mem"]:
        bn_op = batch_normalization(x, g, b, x.mean(axis=0, keepdims=True), x.std(axis=0, keepdims=True), mode=mode)
        f = theano.function([x, b, g], [bn_op])
        res = f(X, G, B)
        utt.assert_allclose(res_ref, res)

        def bn(inputs, gamma, beta, mean, std):
            return batch_normalization(inputs, gamma, beta, mean, std, mode=mode)

        utt.verify_grad(batch_normalization, [X, G, B, X.mean(axis=0)[numpy.newaxis], X.std(axis=0)[numpy.newaxis]])

def test_batch_normalization():

    def bn_ref(x, G, B, M, V):
        n = (x - M) / V
        return n * G + B

    numpy.random.seed(1234)
    X = 1 + numpy.random.random([10, 20]).astype('float32')
    B = 1 + numpy.random.random([20]).astype('float32')
    G = 1 + numpy.random.random([20]).astype('float32')
    M = 1 + numpy.random.random([20]).astype('float32')
    V = 1 + numpy.random.random([20]).astype('float32')

    x = theano.tensor.matrix('x')
    b = theano.tensor.vector('b')
    g = theano.tensor.vector('g')
    m = theano.tensor.vector('m')
    v = theano.tensor.vector('v')

    bn_ref_op = bn_ref(x, g, b, m, v)
    f_ref = theano.function([x, b, g, m, v], [bn_ref_op])
    res_ref = f_ref(X, G, B, M, V)
    for mode in ['low_mem', 'high_mem']:
        bn_op = bn.batch_normalization(x, g, b, m, v, mode=mode)
        f = theano.function([x, b, g, m, v], [bn_op])
        res = f(X, G, B, M, V)
        utt.assert_allclose(res_ref, res)

        def bn_f(inputs, gamma, beta, mean, std):
            return bn.batch_normalization(inputs, gamma, beta, mean, std, mode=mode)
        utt.verify_grad(bn_f, [X, G, B, M, V])

    bn_ref_op = bn_ref(x, g, b, x.mean(axis=0, keepdims=True), x.std(axis=0, keepdims=True))
    f_ref = theano.function([x, b, g], [bn_ref_op])
    res_ref = f_ref(X, G, B)
    for mode in ['low_mem', 'high_mem']:
        bn_op = bn.batch_normalization(x, g, b, x.mean(axis=0, keepdims=True), x.std(axis=0, keepdims=True), mode=mode)
        f = theano.function([x, b, g], [bn_op])
        res = f(X, G, B)
        utt.assert_allclose(res_ref, res)

        def bn_f(inputs, gamma, beta, mean, std):
            return bn.batch_normalization(inputs, gamma, beta, mean, std, mode=mode)
        utt.verify_grad(bn_f, [X, G, B,
                               X.mean(axis=0)[numpy.newaxis], X.std(axis=0)[numpy.newaxis]])

 def _inference(self, input_):
     output = bn.batch_normalization(input_,
         self.gamma.dimshuffle(*self.pattern),
         self.beta.dimshuffle(*self.pattern),
         self.pop_means.dimshuffle(*self.pattern),
         tensor.sqrt(self.pop_vars.dimshuffle(*self.pattern) +
                     self.epsilon),
         mode='low_mem')
     return output

 def conv_bn(inputs, gamma, beta, mean, std):
     return batch_normalization(
         inputs,
         gamma.dimshuffle("x", 0, "x", "x"),
         beta.dimshuffle("x", 0, "x", "x"),
         mean.dimshuffle("x", 0, "x", "x"),
         std.dimshuffle("x", 0, "x", "x"),
         mode=mode,
     )

 def get_hidden_values(self, input):
     lin_output = T.dot(input, self.W) + self.b
     bn_output = batch_normalization(
         inputs=lin_output,
         gamma=self.gamma_h,
         beta=self.beta_h,
         mean=lin_output.mean((0,), keepdims=True),
         std=lin_output.std((0,), keepdims=True),
         mode="low_mem",
     )
     return self.actv_fcn(bn_output)

 def get_reconstructed_input(self, hidden):
     lin_output = T.dot(hidden, self.W_prime) + self.b_prime
     bn_output = batch_normalization(
         inputs=lin_output,
         gamma=self.gamma_o,
         beta=self.beta_o,
         mean=lin_output.mean((0,), keepdims=True),
         std=lin_output.std((0,), keepdims=True),
         mode="low_mem",
     )
     return self.actv_fcn(bn_output)

 def _training(self, input_):
     self.batch_means = input_.mean(axis=self.axes, keepdims=False,
                                    dtype=floatX)
     self.batch_vars = input_.var(axis=self.axes, keepdims=False)
     output = bn.batch_normalization(input_,
         self.gamma.dimshuffle(*self.pattern),
         self.beta.dimshuffle(*self.pattern),
         self.batch_means.dimshuffle(*self.pattern),
         tensor.sqrt(self.batch_vars.dimshuffle(*self.pattern) +
                     self.epsilon),
         mode='low_mem')
     return output

def test_bn_feature_maps():
    def bn_ref(x, G, B, M, V):
        n = (x - M) / V
        return n * G + B

    numpy.random.seed(1234)
    X = 1 + numpy.random.random([10, 20, 4, 4]).astype("float32")
    B = 1 + numpy.random.random([20]).astype("float32")
    G = 1 + numpy.random.random([20]).astype("float32")
    M = 1 + numpy.random.random([20]).astype("float32")
    V = 1 + numpy.random.random([20]).astype("float32")

    x = theano.tensor.tensor4("x")
    b = theano.tensor.vector("b")
    g = theano.tensor.vector("g")
    m = theano.tensor.vector("m")
    v = theano.tensor.vector("v")

    bn_ref_op = bn_ref(
        x,
        g.dimshuffle("x", 0, "x", "x"),
        b.dimshuffle("x", 0, "x", "x"),
        m.dimshuffle("x", 0, "x", "x"),
        v.dimshuffle("x", 0, "x", "x"),
    )
    f_ref = theano.function([x, b, g, m, v], [bn_ref_op])
    res_ref = f_ref(X, G, B, M, V)

    for mode in ["low_mem", "high_mem"]:
        bn_op = batch_normalization(
            x,
            g.dimshuffle("x", 0, "x", "x"),
            b.dimshuffle("x", 0, "x", "x"),
            m.dimshuffle("x", 0, "x", "x"),
            v.dimshuffle("x", 0, "x", "x"),
            mode=mode,
        )
        f = theano.function([x, b, g, m, v], [bn_op])
        res = f(X, G, B, M, V)
        utt.assert_allclose(res_ref, res)

        def conv_bn(inputs, gamma, beta, mean, std):
            return batch_normalization(
                inputs,
                gamma.dimshuffle("x", 0, "x", "x"),
                beta.dimshuffle("x", 0, "x", "x"),
                mean.dimshuffle("x", 0, "x", "x"),
                std.dimshuffle("x", 0, "x", "x"),
                mode=mode,
            )

        utt.verify_grad(conv_bn, [X, G, B, M, V])

def bn_layer(x, a, b, normParam, params, phase):

    ''' Apply BN.    

    # phase = 0 : BN eval with m1v1, BN ups weighter average 
    # phase = 1 : BN eval with m2v2, no BN ups

    '''
        
    minAlpha = params.movingAvMin
    iterStep = params.movingAvStep                  
    # compute mean & variance    
    if params.model == 'convnet':
        mean1 = T.mean(x, axis = (0, 2, 3))
        var1 = T.var(x, axis = (0, 2, 3))
    else:
        mean1 = T.mean(x, axis = 0)
        var1 = T.var(x, axis = 0)

    # moving average as a proxi for validation model 
    alpha = (1.-phase)*T.maximum(minAlpha, 1./normParam['iter'])                     
    mean2 = (1.-alpha)*normParam['mean'] + alpha*mean1 
    var2 = (1.-alpha)*normParam['var'] + alpha*var1   

    mean = (1.-phase)*mean2 + phase*mean1 
    var = (1.-phase)*var1 + phase*var1
    std = T.sqrt(var+eps)

    # apply transformation: 
    if params.model == 'convnet':
        x = bn.batch_normalization(x, a.dimshuffle('x', 0, 'x', 'x'), b.dimshuffle('x', 0, 'x', 'x'), 
                                mean.dimshuffle('x', 0, 'x', 'x'), std.dimshuffle('x', 0, 'x', 'x'), mode='high_mem')
    else:    
        x = bn.batch_normalization(x, a, b, mean, std) 
    updateBN = [mean2, var2, mean1, var1, normParam['iter']+iterStep]  
    return x, updateBN

def test_bn_feature_maps():

    def bn_ref(x, G, B, M, V):
        n = (x - M) / V
        return n * G + B

    numpy.random.seed(1234)
    X = 1 + numpy.random.random([2, 3, 4, 4]).astype('float32')
    B = 1 + numpy.random.random([3]).astype('float32')
    G = 1 + numpy.random.random([3]).astype('float32')
    M = 1 + numpy.random.random([3]).astype('float32')
    V = 1 + numpy.random.random([3]).astype('float32')

    x = theano.tensor.tensor4('x')
    b = theano.tensor.vector('b')
    g = theano.tensor.vector('g')
    m = theano.tensor.vector('m')
    v = theano.tensor.vector('v')

    bn_ref_op = bn_ref(x,
                       g.dimshuffle('x', 0, 'x', 'x'),
                       b.dimshuffle('x', 0, 'x', 'x'),
                       m.dimshuffle('x', 0, 'x', 'x'),
                       v.dimshuffle('x', 0, 'x', 'x'))
    f_ref = theano.function([x, b, g, m, v], [bn_ref_op])
    res_ref = f_ref(X, G, B, M, V)

    for mode in ['low_mem', 'high_mem']:
        bn_op = bn.batch_normalization(x,
                                       g.dimshuffle('x', 0, 'x', 'x'),
                                       b.dimshuffle('x', 0, 'x', 'x'),
                                       m.dimshuffle('x', 0, 'x', 'x'),
                                       v.dimshuffle('x', 0, 'x', 'x'),
                                       mode=mode)
        f = theano.function([x, b, g, m, v], [bn_op])
        res = f(X, G, B, M, V)
        utt.assert_allclose(res_ref, res)

        def conv_bn(inputs, gamma, beta, mean, std):
            return bn.batch_normalization(inputs,
                                          gamma.dimshuffle('x', 0, 'x', 'x'),
                                          beta.dimshuffle('x', 0, 'x', 'x'),
                                          mean.dimshuffle('x', 0, 'x', 'x'),
                                          std.dimshuffle('x', 0, 'x', 'x'),
                                          mode=mode)
        utt.verify_grad(conv_bn, [X, G, B, M, V])

	def __init__(self, x, n_in, n_out, dropout_on,
				layer=0, act=T.nnet.sigmoid,
				w = None, b = None, dropout_rate=0.3):
		if w==None:
			w = theano.shared(
				value=w_init(n_in, n_out),
				name='w'+str(layer),
				borrow=True
			)

		if b==None:
			b = theano.shared(
				value=b_init(n_out),
				name='b'+str(layer),
				borrow=True
			)

		self.w = w
		self.b = b
		self.gamma = theano.shared(value = numpy.ones((n_out,), 
						dtype=theano.config.floatX), name='gamma')
		self.beta = theano.shared(value = numpy.zeros((n_out,), 
						dtype=theano.config.floatX), name='beta')

		rng = np.random.RandomState(42)
		srng = RandomStreams(rng.randint(10**9))
		mask = srng.binomial(n=1, p=1-dropout_rate, size=x.shape)
		cast_mark = T.cast(mask, theano.config.floatX)

		drop_input = T.switch(dropout_on, x*cast_mark,x*(1-dropout_rate))
		lin_output = T.dot(drop_input, self.w) + self.b

		bn_output = batch_normalization(inputs = lin_output,
			gamma = self.gamma, beta = self.beta, 
			mean = lin_output.mean((0,), keepdims=True),
			std = lin_output.std((0,), keepdims = True),
						mode='low_mem')

		self.output = (
			bn_output if act is None
			else act(bn_output)
		)

		self.params = [self.w, self.b]

def test_BNComposite():
    try:
        orig = theano.config.compute_test_value

        theano.config.compute_test_value = 'raise'

        def bn_ref(x, G, B, M, V):
            n = (x - M) / V
            return n * G + B

        numpy.random.seed(1234)
        X = 1 + numpy.random.random([10, 20]).astype('float32')
        B = 1 + numpy.random.random([20]).astype('float32')
        G = 1 + numpy.random.random([20]).astype('float32')
        M = 1 + numpy.random.random([20]).astype('float32')
        V = 1 + numpy.random.random([20]).astype('float32')

        x = theano.tensor.matrix('x')
        b = theano.tensor.vector('b')
        g = theano.tensor.vector('g')
        m = theano.tensor.vector('m')
        v = theano.tensor.vector('v')

        x.tag.test_value = numpy.random.rand(2, 2).astype(theano.config.floatX)
        b.tag.test_value = numpy.random.rand(2).astype(theano.config.floatX)
        g.tag.test_value = numpy.random.rand(2).astype(theano.config.floatX)
        m.tag.test_value = numpy.random.rand(2).astype(theano.config.floatX)
        v.tag.test_value = numpy.random.rand(2).astype(theano.config.floatX)

        bn_ref_op = bn_ref(x, g, b, m, v)
        f_ref = theano.function([x, b, g, m, v], [bn_ref_op])
        res_ref = f_ref(X, G, B, M, V)
        for mode in ['low_mem', 'high_mem']:
            bn_op = bn.batch_normalization(x, g, b, m, v, mode=mode)
            f = theano.function([x, b, g, m, v], [bn_op])
            res = f(X, G, B, M, V)
            utt.assert_allclose(res_ref, res)
    finally:
        theano.config.compute_test_value = orig

    def __init__(self, rng, input, n_in, n_out, stochastic=False, binary=True, W=None, b=None,
                 activation=T.nnet.relu):
        """
        Typical hidden layer of a MLP: units are fully-connected and have
        sigmoidal activation function. Weight matrix W is of shape (n_in,n_out)
        and the bias vector b is of shape (n_out,).
        NOTE : The nonlinearity used here is ReLU

        :type rng: numpy.random.RandomState
        :param rng: a random number generator used to initialize weights

        :type input: theano.tensor.dmatrix
        :param input: a symbolic tensor of shape (n_examples, n_in)

        :type n_in: int
        :param n_in: dimensionality of input

        :type n_out: int
        :param n_out: number of hidden units
        
        :type binary: boolean
        :param binary: indicates whether to implement Binary Connect binarization for weights
        
        :type stochastic: boolean
        :param stochastic: indicate whether to implement a stochatic or deterministic Binary Connect layer
        
        :type activation: theano.Op or function
        :param activation: Non linearity to be applied in the hidden
                           layer
        """
        self.input = input

        # `W` is initialized with `W_values` which is uniformely sampled
        # from sqrt(-6./(n_in+n_hidden)) and sqrt(6./(n_in+n_hidden))
        # for tanh activation function
        # the output of uniform if converted using asarray to dtype
        # theano.config.floatX so that the code is runable on GPU
        # Note : optimal initialization of weights is dependent on the
        #        activation function used (among other things).
        #        For example, results presented in [Xavier10] suggest that you
        #        should use 4 times larger initial weights for sigmoid
        #        compared to tanh
        #        We have no info for other function, so we use the same as
        #        tanh.
        if W is None:
            W_values = numpy.asarray(
                rng.uniform(
                    low=-numpy.sqrt(6. / (n_in + n_out)),
                    high=numpy.sqrt(6. / (n_in + n_out)),
                    size=(n_in, n_out)
                ),
                dtype=theano.config.floatX
            )
           
            if activation == theano.tensor.nnet.sigmoid:
                W_values *= 4

            W = theano.shared(value=W_values, name='W', borrow=True)

        if b is None:
            b_values = numpy.zeros((n_out,), dtype=theano.config.floatX)
            b = theano.shared(value=b_values, name='b', borrow=True)

       
        self.high = numpy.float32(numpy.sqrt(6. / (n_in + n_out)))
        self.W0 = numpy.float32(self.high/2)

        
        # binarize weights either deterministically or stochastically, if indicated
        def hard_sigma(w):
            p=T.clip((w+1)/2,0,1)
            return p
        
        if stochastic:
            p = hard_sigma(W/self.W0)
            p_mask = T.cast(numpy.random.binomial(n=1, p=p.eval(), size=(n_in, n_out)), theano.config.floatX)
            Wb = T.switch(p_mask,self.W0,-self.W0).eval()
        else:        
            Wb = T.switch(T.ge(W.get_value(),0),self.W0,-self.W0).eval()
       
        if binary:
            Wb = theano.shared(Wb, name='Wb', borrow=True)
            self.Wb = Wb
        else:
            self.Wb=W
            
        self.W = W
        self.b = b
        self.n_in=n_in
         
        self.gamma = theano.shared(value = numpy.ones((n_out,), dtype=theano.config.floatX), name='gamma')
        self.beta = theano.shared(value = numpy.zeros((n_out,), dtype=theano.config.floatX), name='beta')

        lin_output = T.dot(input, self.Wb) + self.b
        
        # batch normalization at output
        bn_output = batch_normalization(inputs = lin_output,
            gamma = self.gamma, beta = self.beta, mean = lin_output.mean((0,), keepdims=True),
            std = lin_output.std((0,), keepdims = True),
                        mode='low_mem')
#.........这里部分代码省略.........

 def bn_f(inputs, gamma, beta, mean, std):
     return bn.batch_normalization(inputs, gamma, beta, mean, std, mode=mode)

    def __init__(self, input, rng, n_in, n_out, stochastic=False, binary=True):
        # initialize with 0 the weights W as a matrix of shape (n_in, n_out)
        self.high = 2.#numpy.float32(numpy.sqrt(6. / (float(n_in) + float(n_out) )))
        self.W0 = numpy.float32(self.high/2.)        
        #self.W = theano.shared(value=numpy.zeros((n_in, n_out),
        #                         dtype=theano.config.floatX),
        #                         name='W', borrow=True)
        #self.high = numpy.float32(2.)
        #self.W0 = numpy.float32(self.high/2.)
        W_values = numpy.asarray(
            rng.uniform(
                low= -1.,#numpy.sqrt(6. / (n_in + n_out)),
                high= 1.,#numpy.sqrt(6. / (n_in + n_out)),
                size=(n_in, n_out)
            ),
            dtype=theano.config.floatX
        )
        #srng = RandomStreams(seed=420) 
        srng = theano.sandbox.rng_mrg.MRG_RandomStreams(420)
        self.W = theano.shared(value=W_values, name='W', borrow=True)

        # initialize the baises b as a vector of n_out 0s
        self.b = theano.shared(value=numpy.zeros((n_out,),
                                dtype=theano.config.floatX),
                                name='b', borrow=True)

        self.gamma = theano.shared(value = numpy.ones((n_out,), dtype=theano.config.floatX), name='gamma')
        self.beta = theano.shared(value = numpy.zeros((n_out,), dtype=theano.config.floatX), name='beta')
        
        def hard_sigma(w):
            p=T.clip((w+1)/2,0,1)
            return p

        if binary:
            Wb = hard_sigma(self.W/self.W0)
            if stochastic:
                #Wb = T.cast(numpy.random.binomial(n=1, p=T.ge(Wb), size=(n_in, n_out)),  theano.config.floatX)
                Wb = srng.binomial(n=1, p=Wb, size=(n_in, n_out) )

            else:
                Wb = T.round(Wb)

            # Leave below alone
            Wb = T.switch(Wb,self.W0, -self.W0)
            #Wb = T.cast(T.switch(Wb,self.W0, -self.W0), dtype=theano.config.floatX)
            #Wb = theano.shared(Wb.eval(), name='Wb', borrow=True)
            self.Wb = Wb

        else:
            self.Wb = self.W

        # parameters of the model
        self.linear = T.dot(input, self.Wb) + self.b

        bn_output = batch_normalization(inputs = self.linear,
                    gamma = self.gamma, beta = self.beta, mean = self.linear.mean((0,), keepdims=True),
                    std = T.ones_like(self.linear.var((0,), keepdims = True)), mode='high_mem')

        self.linear_output = bn_output
        self.y_pred = T.argmax(bn_output, axis=1)

        if binary:
            self.params = [self.W, self.Wb, self.gamma, self.beta, self.b]

        elif binary==False:
            self.params = [self.Wb, self.gamma, self.beta, self.b]

        self.len_params = len(self.params)
        self.n_in=n_in
        # keep track of model input
        self.input = input

    def __init__(self, rng, input, n_in, n_out, W=None, b=None,
                 activation=T.tanh):
        """
        Typical hidden layer of a MLP: units are fully-connected and have
        sigmoidal activation function. Weight matrix W is of shape (n_in,n_out)
        and the bias vector b is of shape (n_out,).
        NOTE : The nonlinearity used here is tanh

        Hidden unit activation is given by: tanh(dot(input,W) + b)

        :type rng: numpy.random.RandomState
        :param rng: a random number generator used to initialize weights

        :type input: theano.tensor.dmatrix
        :param input: a symbolic tensor of shape (n_examples, n_in)

        :type n_in: int
        :param n_in: dimensionality of input

        :type n_out: int
        :param n_out: number of hidden units

        :type activation: theano.Op or function
        :param activation: Non linearity to be applied in the hidden
                           layer
        """
        self.input = input

        # `W` is initialized with `W_values` which is uniformely sampled
        # from sqrt(-6./(n_in+n_hidden)) and sqrt(6./(n_in+n_hidden))
        # for tanh activation function
        # the output of uniform if converted using asarray to dtype
        # theano.config.floatX so that the code is runable on GPU
        # Note : optimal initialization of weights is dependent on the
        #        activation function used (among other things).
        #        For example, results presented in [Xavier10] suggest that you
        #        should use 4 times larger initial weights for sigmoid
        #        compared to tanh
        #        We have no info for other function, so we use the same as
        #        tanh.

        if W is None:
            W_values = numpy.asarray(
                rng.uniform(
                    low=-numpy.sqrt(6. / (n_in + n_out)),
                    high=numpy.sqrt(6. / (n_in + n_out)),
                    size=(n_in, n_out)
                ),
                dtype=theano.config.floatX
            )
            if activation == theano.tensor.nnet.sigmoid:
                W_values *= 4

            W = theano.shared(value=W_values, name='W', borrow=True)

        if b is None:
            b_values = numpy.zeros((n_out,), dtype=theano.config.floatX)
            b = theano.shared(value=b_values, name='b', borrow=True)

        self.gamma = theano.shared(value = numpy.ones((n_out,), dtype=theano.config.floatX), name='gamma')
        self.beta = theano.shared(value = numpy.zeros((n_out,), dtype=theano.config.floatX), name='beta')

        self.W = W
        self.b = b

        lin_output = T.dot(input, self.W) + self.b
        
        bn_output = batch_normalization(inputs = lin_output,
            gamma = self.gamma, beta = self.beta, mean = lin_output.mean((0,), keepdims=True),
            std = lin_output.std((0,), keepdims = True),
                        mode='low_mem')

        self.output = (T.clip(bn_output,0,20) if activation is 'relu' else activation(bn_output))        
        # self.output = (
        #     lin_output if activation is None
        #     else activation(lin_output)
        # )
        # parameters of the model
        self.params = [self.W, self.b, self.gamma, self.beta]