def test_smolyak(self):

        import numpy

        f = lambda x: numpy.row_stack([
            x[0,:] * x[1,:]**0.5,
            x[1,:] * x[1,:] - x[0,:] * x[0,:]
        ])

        bounds = numpy.row_stack([
            [0.5,0.1],
            [2,3]
        ])

        from dolo.numeric.smolyak import SmolyakGrid
        sg = SmolyakGrid(bounds,3)

        values = f(sg.grid)
        sg.fit_values(values)
        theta_0 = sg.theta.copy()

        def fobj(theta):
            sg.theta = theta
            return sg(sg.grid)

        fobj(theta_0)

def mirrorpad(img, width):                            
    """Return the image resulting from padding width amount of pixels on each                                                  
    sides of the image img.  The padded values are mirror image with respect to                                                
    the borders of img.  Width can be an integer or a tuple of the form (north,                                                
    south, east, west).                                                                                                        
    """                                                                                                                        
    n, s, e, w = _get4widths(width)                                                                                            
    row = N.shape(img)[0]                                                                                                      
                                                                                                                                                               
    col = N.shape(img)[1]                                                                                                      
    
    if n != 0:                                                                                                                 
        north = img[:n,:]                                                                                                      
        img = N.row_stack((north[::-1,:], img))                                                                                
                                                                                       
    if s != 0:                                                                                                                 
        south = img[-s:,:]                                                                                                     
        img = N.row_stack((img, south[::-1,:]))                                                                                
        
    if e != 0:                                                                                                                 
        east = img[:,-e:]                                                                                                      
        img = N.column_stack((img, east[:,::-1]))                                                                              
        
    if w != 0:                                                                                                                 
        west = img[:,:w]                                                                                                       
        img = N.column_stack((west[:,::-1], img))                                                                              
        
        
        
        return img

    def __init__(self,d,l, dtype=None):

        self.d = d
        self.l = l

        [self.smolyak_points, self.smolyak_indices] = smolyak_grids(d,l)

        self.u_grid = array( self.smolyak_points.T, order='C')

        self.isup = max(max(self.smolyak_indices))
        self.n_points = len(self.smolyak_points)
        #self.grid = self.real_gri

        Ts = chebychev( self.smolyak_points.T, self.n_points - 1 )
        C = []
        for comb in self.smolyak_indices:
            p = reduce( mul, [Ts[comb[i],i,:] for i in range(self.d)] )
            C.append(p)
        C = numpy.row_stack( C )

        # C is such that :  X = theta * C
        self.__C__ = C
        self.__C_inv__ = numpy.linalg.inv(C)  # would be better to store preconditioned matrix

        self.bounds = numpy.row_stack([(0,1)]*d)

    def train(self, training_data_array): 
        data = training_data_array[0] #dict
            # Step 2: Forward propagation
        a1 = np.mat(data['y0']).T
            # 400*1 matrix
        z2 = np.dot(self.theta1, np.row_stack((1, a1)))
            # num_hidden_nodes*1 matrix
        a2 = self.sigmoid(z2)

        z3 = np.dot(self.theta2, np.row_stack((1, a2)))
            # 10*1 matrix
        a3 = self.sigmoid(z3)

            # Step 3: Back propagation
        y = [0] * 10 # y is a python list for easy initialization and is later turned into an np matrix (2 lines down).
        y[data['label']] = 1
            # 1*10 list
                      
        delta3 = a3 - np.mat(y).T
            # 10*1 matrix
        z2plus = np.row_stack((0, z2))
            # (num_hidden_nodes+1)*1 matrix
        delta2 = np.multiply(np.dot(self.theta2.T, delta3), self.sigmoid_prime(z2plus))
            # (num_hidden_nodes+1)*1 matrix
        delta2 = delta2[1:,0]
            # num_hidden_nodes*1 matrix

            # Step 4: Update weights
        self.theta1 -= self.LEARNING_RATE * np.dot(delta2, np.row_stack((1, a1)).T)
        self.theta2 -= self.LEARNING_RATE * np.dot(delta3, np.row_stack((1, a2)).T)

def find_goto_point_surface_1(grid,pt,display_list=None):
    ''' returns p_erratic,p_edge,surface_height
        p_erratic - point where the erratic's origin should go to. (in thok0 frame)
        p_edge - point on the edge closest to pt.
        p_erratic,p_edge are 3x1 matrices.
        surface_height - in thok0 coord frame.
    '''
    pt_thok = pt
    
    # everything is happening in the thok coord frame.
    close_pt,approach_vector = find_approach_direction(grid,pt_thok,display_list)
    if close_pt == None:
        return None,None,None

    # move perpendicular to approach direction.
#    lam = -(close_pt[0:2,0].T*approach_vector)[0,0]
#    lam = min(lam,-0.3) # atleast 0.3m from the edge.
    lam = -0.4
    goto_pt = close_pt[0:2,0] + lam*approach_vector # this is where I want the utm
                                                    # to be
    err_to_thok = tr.erraticTglobal(tr.globalTthok0(np.matrix([0,0,0]).T))
    goto_pt_erratic = -err_to_thok[0,0]*approach_vector + goto_pt # this is NOT
                # general. It uses info about the two frames. If frames move, bad
                # things can happen.

    if display_list != None:
        display_list.append(pu.CubeCloud(np.row_stack((goto_pt,close_pt[2,0])),color=(0,250,250), size=(0.012,0.012,0.012)))
        display_list.append(pu.Line(np.row_stack((goto_pt,close_pt[2,0])).A1,close_pt.A1,color=(255,20,0)))

    p_erratic = np.row_stack((goto_pt_erratic,np.matrix((close_pt[2,0]))))
    print 'p_erratic in thok0:', p_erratic.T

    return p_erratic,close_pt,close_pt[2,0]

 def predict_test(batch, model, params,  **kwparams):
   """ some code duplication here with forward pass, but I think we want the freedom in future """
   F = np.row_stack(x['image']['feat'] for x in batch)
   lda_enabled = params.get('lda',0)
   L = np.zeros((params.get('image_encoding_size',128),lda_enabled))
   if lda_enabled:
      L = np.row_stack(x['image']['topics'] for x in batch)
   We = model['We']
   if lda_enabled:
     Wlda = model['Wlda']
     lda = L.dot(Wlda)
   be = model['be']
   Xe = F.dot(We) + be # Xe becomes N x image_encoding_size
   #print('L shape', L.shape)
   #print('Wlda shape', Wlda.shape)
   generator_str = params['generator']
   Generator = decodeGenerator(generator_str)
   Ys = []
   guide_input = params.get('guide',None)
   for i,x in enumerate(batch):
     Xi = Xe[i,:]
     if guide_input=='cca':
       guide = get_guide(guide_input,F[i,:],kwparams.get('ccaweights'))
     else:
       guide = get_guide(guide_input,F[i,:],L=L[i,:])
     if (lda_enabled and guide_input=="image") or (lda_enabled and not guide_input):
       guide = lda[i,:]
       print 'guide = lda'
     gen_Y = Generator.predict(Xi, guide, model, model['Ws'], params, **kwparams)
     Ys.append(gen_Y)
   return Ys

    def add_margins(self, matrix, regions=2):
        nrow, ncol = np.shape(matrix)

        assert regions<=2
        for i in range(regions):
            nc = ncol/regions
            nr = nrow

            if i==0:
                m = matrix[:, 0:nc]
            else:
                m = matrix[:, nc:]

            left_col = np.ones((nr+2,1))
            bottom_row = np.ones((1, nc))
            top_row = np.reshape([[0,1] for _ in range(nc/2)], (1, nc))
            right_row = np.reshape([[0,1] for _ in range(nr/2+1)], (nr+2, 1))
            m = np.row_stack((top_row, m))
            m = np.row_stack((m, bottom_row))
            m = np.column_stack((left_col, m))
            m = np.column_stack((m, right_row))
            if i == 0:
                ans = m
            else:
                ans = np.column_stack((ans, m))
        return ans

def getstats_base(X, linds):
    sval = {}
    for l_ind in linds:
        print(l_ind, X['model_state']['layers'][l_ind]['name'])
        layer = X['model_state']['layers'][l_ind]
        w = layer['weights'][0]
        karray = stats.kurtosis(w)
        kall = stats.kurtosis(w.ravel())
        cf0 = np.corrcoef(w)
        cf0t = np.corrcoef(w.T)
        wmean = w.mean(1)
        w2mean = (w**2).mean(1)
        lname = X['model_state']['layers'][l_ind]['name']
        sval[lname] = {'karray': karray, 'kall': kall, 'corr0': cf0, 'corr0_t': cf0t,
                            'wmean': wmean, 'w2mean': w2mean}

        if 'filterSize' in X['model_state']['layers'][l_ind]:
            fs = X['model_state']['layers'][l_ind]['filterSize'][0]
            ws = w.shape
            w = w.reshape((ws[0] / (fs**2), fs, fs, ws[1]))
            mat = np.row_stack([np.row_stack([w[i, j, :, :] for i in range(w.shape[0])]).T for j in range(w.shape[1])] )
            cf = np.corrcoef(mat.T)
            cft = np.corrcoef(mat)
            mat2 = np.row_stack([np.row_stack([w[i, :, :, j] for i in range(w.shape[0])]).T for j in range(w.shape[3])] )
            cf2 = np.corrcoef(mat2.T)
            cf2t = np.corrcoef(mat2)
            sval[lname].update({'corr': cf, 'corr2': cf2, 'corr_t': cft, 'corr2_t': cf2t})

    return sval

def setup_data(db,h,w):
    corners = np.array([[0,0],[0,h-1],[w-1,0],[w-1,h-1]])
    im1_pts = np.array(db_unpack(db)[0])
    im1_pts = np.row_stack([im1_pts,corners])
    im2_pts = np.array(db_unpack(db)[1])
    im2_pts = np.row_stack([im2_pts,corners])
    return im1_pts,im2_pts

 def solveMartix(self):
     MNA1 = np.column_stack((self.C, self.D))
     MNA2 = np.column_stack((self.G, self.B))
     MNA = np.row_stack((MNA2, MNA1))
     RHS = np.row_stack((self.U, self.I))
     x = np.linalg.lstsq(MNA, RHS)
     self.x = x[0]

 def printMartix(self):
     MNA1 = np.column_stack((self.C, self.D))
     MNA2 = np.column_stack((self.G, self.B))
     MNA = np.row_stack((MNA2, MNA1))
     RHS = np.row_stack((self.U, self.I))
     print "MNA", MNA
     print "RHS", RHS

def plot_attention_loss(train_plot, test_plot) :
	
	testing_freq = (test_plot['iter'][1]-test_plot['iter'][0]) / (train_plot['iter'][1]-train_plot['iter'][0])
	loss_y = train_plot['iter']
	
	pl.figure(1)
	test_loss = test_plot['dir_loss']
	train_loss = train_plot['dir_loss']
	test_loss = numpy.row_stack(test_loss)
	test_loss = numpy.tile(test_loss, (1, testing_freq))
	test_loss = list(test_loss.flatten())
	test_loss += [test_loss[-1]] * max(0,len(train_loss) - len(test_loss))
	test_loss = test_loss[:len(train_loss)]
	pl.plot(loss_y, train_loss, 'k-', label='Train', linewidth=0.75, marker='x')
	pl.plot(loss_y, test_loss,  'r-', label='Test' , linewidth=0.75)
	pl.legend(loc='best')
	pl.xlabel('Iter')
	pl.title('Direction Loss')

	pl.figure(2)
	test_loss = test_plot['cls_loss']
	train_loss = train_plot['cls_loss']
	test_loss = numpy.row_stack(test_loss)
	test_loss = numpy.tile(test_loss, (1, testing_freq))
	test_loss = list(test_loss.flatten())
	test_loss += [test_loss[-1]] * max(0,len(train_loss) - len(test_loss))
	test_loss = test_loss[:len(train_loss)]
	pl.plot(loss_y, train_loss, 'k-', label='Train', linewidth=0.75, marker='x')
	pl.plot(loss_y, test_loss,  'r-', label='Test' , linewidth=0.75)
	pl.legend(loc='best')
	pl.xlabel('Iter')
	pl.title('Classification Loss')
	pl.show()

    def test(self, ratings):
        L = ratings.kv_dict.items()
        R = np.matrix([I[1] for I in L]).T 
        sqr_err = list()
       
        au = len(ratings.rows) - len(self.user_factor)
        self.user_factor = np.row_stack([self.user_factor, np.matrix(np.zeros([au, self.factor]))])
        self.user_bias = np.row_stack([self.user_bias, np.matrix(np.zeros([au, 1]))])
 
        ai = len(ratings.cols) - len(self.item_factor)
        self.item_factor = np.row_stack([self.item_factor, np.matrix(np.zeros([ai, self.factor]))])
        self.item_bias = np.row_stack([self.item_bias, np.matrix(np.zeros([ai, 1]))])
 
        for s in range(0, len(L), self.batch_size):
            mini_batch = L[s : s + self.batch_size]
            r = R[s : s + self.batch_size]
            
            uid = [I[0][0] for I in mini_batch]
            iid = [I[0][1] for I in mini_batch]
 
            base_line = self.mu + self.user_bias[uid] + self.item_bias[iid]
            r_pred = base_line + np.sum(np.multiply(self.user_factor[uid], self.item_factor[iid]), 1)

            err = r - r_pred
            sqr_err.append(float(err.T * err) / self.batch_size)
        
        rmse = math.sqrt(np.mean(np.matrix(sqr_err)))
        # sys.stderr.write('RMSE: %f\n' % (rmse))
        
        return rmse

	def save_Ian_E_H_Yen(self,folder):
		if not self.B_lower is None:
			print ('self.B_lower is not None, you should convert your problem with convertToOnesideInequalitySystem first')
			raise
		if not np.all(self.lowerbounds==0):
			print ('lower bound constraint on variables should at 0')
			raise
		
		import os		
		Aeq=self.Aequalities.tocoo()
		tmp=np.row_stack(([Aeq.shape[0],Aeq.shape[1],0.0],np.column_stack((Aeq.row+1,Aeq.col+1,Aeq.data))))
		np.savetxt(os.path.join(folder,'Aeq'),tmp,fmt='%d %d %f')
		np.savetxt(os.path.join(folder,'beq'),self.Bequalities,fmt='%f')
		np.savetxt(os.path.join(folder,'c'),self.costsvector,fmt='%f')
		nbvariables=self.costsvector.size
		upperbounded=np.nonzero(~np.isinf(self.upperbounds))[0]
		nbupperbounded=len(upperbounded)
		Aineq2=sparse.coo_matrix((np.ones(nbupperbounded),(np.arange(nbupperbounded),upperbounded)),(nbupperbounded,nbvariables))
		Aineq=scipy.sparse.vstack((self.Ainequalities,Aineq2)).tocoo()
		bupper=np.hstack((self.B_upper,self.upperbounds[upperbounded]))
		tmp=np.row_stack(([Aineq.shape[0],Aineq.shape[1],0.0],np.column_stack((Aineq.row+1,Aineq.col+1,Aineq.data))))
		np.savetxt(os.path.join(folder,'A'),tmp,fmt='%d %d %f')	
		np.savetxt(os.path.join(folder,'b'),bupper,fmt='%f')	
	
		with open(os.path.join(folder,'meta'), 'w') as f:
			f.write('nb	%d\n'%nbvariables)
			f.write('nf	%d\n'%0)
			f.write('mI	%d\n'%Aineq.shape[0])
			f.write('mE	%d\n'%Aeq.shape[0])

def plot_attention_loss(train_plot, test_plot):

    testing_freq = (test_plot["iter"][1] - test_plot["iter"][0]) / (train_plot["iter"][1] - train_plot["iter"][0])
    loss_y = train_plot["iter"]

    pl.figure(1)
    test_loss = test_plot["dir_loss"]
    train_loss = train_plot["dir_loss"]
    test_loss = numpy.row_stack(test_loss)
    test_loss = numpy.tile(test_loss, (1, testing_freq))
    test_loss = list(test_loss.flatten())
    test_loss += [test_loss[-1]] * max(0, len(train_loss) - len(test_loss))
    test_loss = test_loss[: len(train_loss)]
    pl.plot(loss_y, train_loss, "k-", label="Train", linewidth=0.75, marker="x")
    pl.plot(loss_y, test_loss, "r-", label="Test", linewidth=0.75)
    pl.legend(loc="best")
    pl.xlabel("Iter")
    pl.title("Direction Loss")

    pl.figure(2)
    test_loss = test_plot["cls_loss"]
    train_loss = train_plot["cls_loss"]
    test_loss = numpy.row_stack(test_loss)
    test_loss = numpy.tile(test_loss, (1, testing_freq))
    test_loss = list(test_loss.flatten())
    test_loss += [test_loss[-1]] * max(0, len(train_loss) - len(test_loss))
    test_loss = test_loss[: len(train_loss)]
    pl.plot(loss_y, train_loss, "k-", label="Train", linewidth=0.75, marker="x")
    pl.plot(loss_y, test_loss, "r-", label="Test", linewidth=0.75)
    pl.legend(loc="best")
    pl.xlabel("Iter")
    pl.title("Classification Loss")
    pl.show()

    def __call__(self, event):
        if event.button == 1:
            print 'You press button 1 to zomm, to mask use right button'
        if event.button >=2:
            if event.inaxes!=self.graph.axes: return
#            self.bt.append(event.button)
            self.xs.append(event.xdata)
            
            if len(self.xs)%2 == 0:  # only odd numbers 
                if event.button == 3:
                    self.wei.append(0)
                    print 'You press button 3 to mask, maksed regions is', self.xs[-2],self.xs[-1],self.wei[-1]
                    self.tmp=np.column_stack((self.xs[-2],self.xs[-1],self.wei[-1]))
                    if (self.masked[0][0])==0:
                         self.masked=np.row_stack((self.tmp))
                    else:
                        self.masked=np.row_stack((self.masked,self.tmp))
                    plot(self.xso[(self.xso>=self.xs[-2]) & (self.xso<=self.xs[-1])],self.yso[(self.xso>=self.xs[-2]) & (self.xso<=self.xs[-1])],color='red',lw='2')
                    self.graph.figure.canvas.draw()
                if event.button == 2:
                    self.wei.append(2)
                    print 'You press button 2 attributed weight 2 for this region', self.xs[-2],self.xs[-1],self.wei[-1]
                    self.tmp=np.column_stack((self.xs[-2],self.xs[-1],self.wei[-1]))
                    if (self.masked[0][0])==0:
                         self.masked=np.row_stack((self.tmp))
                    else:
                        self.masked=np.row_stack((self.masked,self.tmp))
                    plot(self.xso[(self.xso>=self.xs[-2]) & (self.xso<=self.xs[-1])],self.yso[(self.xso>=self.xs[-2]) & (self.xso<=self.xs[-1])],color='green',lw='2')                   
                    self.graph.figure.canvas.draw()

    def get_scan(self, avoid_duplicate=False, avg=1, remove_graze=True):
        ''' avoid_duplicate - prevent duplicate scans which will happen if get_scan is
            called at a rate faster than the scanning rate of the hokuyo.
            avoid_duplicate == True, get_scan will block till new scan is received.
            (~.2s for urg and 0.05s for utm)
        '''
        l = []
        l2 = []
        for i in range(avg):
            hscan = self.hokuyo.get_scan(avoid_duplicate)
            l.append(hscan.ranges)
            l2.append(hscan.intensities)

        ranges_mat = np.row_stack(l)
        ranges_mat[np.where(ranges_mat==0)] = -1000. # make zero pointvery negative
        ranges_avg = (ranges_mat.sum(0)/avg)
        if self.flip:
            ranges_avg = np.fliplr(ranges_avg)
        hscan.ranges = ranges_avg
    
        intensities_mat = np.row_stack(l2)
        if self.hokuyo_type == 'utm':
            hscan.intensities = (intensities_mat.sum(0)/avg)

        if remove_graze:
            if self.hokuyo_type=='utm':
                hp.remove_graze_effect_scan(hscan)
            else:
                print 'hokuyo_scan.Hokuyo.get_scan: hokuyo type is urg, but remove_graze is True'

        return hscan

def multi_scale_test(image, max_im_shrink):
    # shrink detecting and shrink only detect big face
    st = 0.5 if max_im_shrink >= 0.75 else 0.5 * max_im_shrink
    det_s = detect_face(image, st)
    if max_im_shrink > 0.75:
        det_s = np.row_stack((det_s,detect_face(image,0.75)))
    index = np.where(np.maximum(det_s[:, 2] - det_s[:, 0] + 1, det_s[:, 3] - det_s[:, 1] + 1) > 30)[0]
    det_s = det_s[index, :]
    # enlarge one times
    bt = min(2, max_im_shrink) if max_im_shrink > 1 else (st + max_im_shrink) / 2
    det_b = detect_face(image, bt)

    # enlarge small iamge x times for small face
    if max_im_shrink > 1.5:
        det_b = np.row_stack((det_b,detect_face(image,1.5)))
    if max_im_shrink > 2:
        bt *= 2
        while bt < max_im_shrink: # and bt <= 2:
            det_b = np.row_stack((det_b, detect_face(image, bt)))
            bt *= 2

        det_b = np.row_stack((det_b, detect_face(image, max_im_shrink)))

    # enlarge only detect small face
    if bt > 1:
        index = np.where(np.minimum(det_b[:, 2] - det_b[:, 0] + 1, det_b[:, 3] - det_b[:, 1] + 1) < 100)[0]
        det_b = det_b[index, :]
    else:
        index = np.where(np.maximum(det_b[:, 2] - det_b[:, 0] + 1, det_b[:, 3] - det_b[:, 1] + 1) > 30)[0]
        det_b = det_b[index, :]

    return det_s, det_b

def _build_series(series, dim_names, comment, delta_name, delta_unit):
    from glue.ligolw import array as ligolw_array
    Attributes = ligolw.sax.xmlreader.AttributesImpl
    elem = ligolw.LIGO_LW(
            Attributes({u"Name": unicode(series.__class__.__name__)}))
    if comment is not None:
        elem.appendChild(ligolw.Comment()).pcdata = comment
    elem.appendChild(ligolw.Time.from_gps(series.epoch, u"epoch"))
    elem.appendChild(ligolw_param.Param.from_pyvalue(u"f0", series.f0,
                                                     unit=u"s^-1"))
    delta = getattr(series, delta_name)
    if numpy.iscomplexobj(series.data.data):
        data = numpy.row_stack((numpy.arange(len(series.data.data)) * delta,
                             series.data.data.real, series.data.data.imag))
    else:
        data = numpy.row_stack((numpy.arange(len(series.data.data)) * delta,
                                series.data.data))
    a = ligolw_array.Array.build(series.name, data, dim_names=dim_names)
    a.Unit = str(series.sampleUnits)
    dim0 = a.getElementsByTagName(ligolw.Dim.tagName)[0]
    dim0.Unit = delta_unit
    dim0.Start = series.f0
    dim0.Scale = delta
    elem.appendChild(a)
    return elem

def softmaxCostAndGradient(y1, target, w2):
    """ Softmax cost function for word2vec models """
    ###################################################################
    # Implement the cost and gradients for one predicted word vector  #
    # and one target word vector as a building block for word2vec     #
    # models, assuming the softmax prediction function and cross      #
    # entropy loss.                                                   #
    # Inputs:                                                         #
    #   - predicted: numpy ndarray, predicted word vector (\hat{r} in #
    #           the written component)                                #
    #   - target: integer, the index of the target word               #
    #   - outputVectors: "output" vectors for all tokens              #
    # Outputs:                                                        #
    #   - cost: cross entropy cost for the softmax word prediction    #
    #   - gradPred: the gradient with respect to the predicted word   #
    #           vector                                                #
    #   - grad: the gradient with respect to all the other word       # 
    #           vectors                                               #
    # We will not provide starter code for this function, but feel    #
    # free to reference the code you previously wrote for this        #
    # assignment!                                                     #
    ###################################################################
    out = np.dot(w2,y1)
    y2 = softmaxactual(out)
    cost = -math.log(y2[target])
    y2[target] -= 1.0
    d2= y2
    grad2 = np.dot(np.row_stack(d2),np.column_stack(y1))
    d1 = np.dot(np.transpose(w2), np.row_stack(d2))
    return cost, d1, grad2