def load_shapes( self, narrow_shape, wide_shape ):
     
     self.ref_signals = []
     
     for i in range(10):
         for j in range(10):
             
             signal = scipy.array([])
             
             for k in range(5):
                 
                 if ( codes[i][k] == 'n' ):
                     signal = scipy.append( signal, narrow_shape )
                 else:
                     signal = scipy.append( signal, wide_shape )
                 
                 if ( codes[j][k] == 'n' ):
                     signal = scipy.append( signal, scipy.zeros(len(narrow_shape)) )
                 else:
                     signal = scipy.append( signal, scipy.zeros(len(wide_shape)) ) 
             
             signal /= norm(signal);
             
             self.ref_signals.append( signal )
             self.ref_values.append( value_map(i,j) )

def mix_prop_test():
    energyfile = "energies_0_K"
    vibfiles = ["0_thermo","1.54_thermo","1.3_thermo","2.3_thermo","1_thermo"]
    efile = open(energyfile,'r')
    efile.readline()
    comp = sp.array([])
    energy0 = sp.array([])
    for line in efile:
        comp = sp.append(comp,float(line.split()[0]))
        energy0 = sp.append(energy0,float(line.split()[1]))
    efile.close()
    phase = []
    for i in range(len(comp)):
        phase.append(Phase(comp[i],energy0[i],vibfiles[i]))
    vib_prop = []
    for i in range(len(phase)):
        vib_prop.append(phase[i].E_thermal)
    vib_mix_prop = binary_mixing_property(vib_prop,comp,2)

    # Text output of data
    for i in range(len(phase[0].T)):
        text = str(phase[0].T[i])
        for j in range(len(phase)):
            text += " "+str(vib_mix_prop[j,i])
        print text
    sys.exit()

    # Graphical output of data
    for i in range(len(phase)):
        plt.plot(phase[i].T,vib_mix_prop[i])
    plt.xlabel('T (K)')
    plt.ylabel('Delta E_vib (meV/cation)')
    plt.show()

def calculate_difference(frames_energy):
    """calculate difference of energies
    
    Implementation following paper "A Highly Robust Audio Fingerprinting System"
    
    :math:`F(n,m)=1` if :math:`E(n,m)-E(n,m+1)-(E(n-1,m)-E(n-1,m+1))>0`
    
    :math:`F(n,m)=0` if :math:`E(n,m)-E(n,m+1)-(E(n-1,m)-E(n-1,m+1))\leq 0`
    
    :param frames_energy: frames of energys
    :type frames_energy: scipy.array
    :return: fingerint
    """
    log.debug('Fingerprinting: calculate_difference')

    # fingerprint vector
    fingerprint = scipy.array([], dtype=int)
    
    # first frame is defined as previous frame
    prev_frame = frames_energy[0]
    print str(len(frames_energy))
    del frames_energy[0]
    
    for n, frame in enumerate(frames_energy):
        # every energy of frequency bands until length-1
        for m in range(len(frame)-1):
            # calculate difference with formula from paper            
            if (frame[m]-frame[m+1]-(prev_frame[m]-prev_frame[m+1]) > 0):
                fingerprint = scipy.append(fingerprint, 1)
            else:
                fingerprint = scipy.append(fingerprint, 0)
            
        prev_frame = frame
        
    return fingerprint

    def f(self,x,t):
        N = len(x)/2
        xdot = pl.array([])

        # modulus the x for periodicity.
        x[N:2*N]= x[N:2*N]%self.d
        # HERE ---->> 1Dify
        for i in range(N):
            temp = 0.0
            for j in range(N):
                if i == j:
                    continue
                #repulsive x interparticle force of j on i
                temp += self.qq*(x[N+i]-x[N+j])/(pl.sqrt((x[N+i]-x[N+j])**2)**3)

                # Add the forces from the two ancored charges
                # First one at x=0
                temp += self.qq*(x[N+i]-0.0)/(pl.sqrt((x[N+i]-0.0)**2)**3)
                # Second one at the other boundary i.e. self.d
                temp += self.qq*(x[N+i]-self.d)/(pl.sqrt((x[N+i]-self.d)**2)**3)

            # periodic force on particle i
            temp += self.As[i]*pl.sin(x[N+i])*pl.cos(t)
            temp -= self.beta*x[i]
            xdot = pl.append(xdot,temp)
        for i in range(N):
            xdot = pl.append(xdot,x[i])
        return xdot

def GP_train(x, y, cov_par, cov_func = None, cov_typ ='SE', \
             cov_fixed = None, prior = None, \
             MF = None, MF_par = None, MF_args = None, \
             MF_fixed = None):
    '''    
    Max likelihood optimization of GP hyper-parameters. Calls
    GP_negloglik. Takes care of merging / splitting the fixed /
    variable and cov / MF parameters
    '''
    if MF != None:
        merged_par = scipy.append(cov_par, MF_par)
        n_MF_par = len(MF_par)
        fixed = scipy.append(scipy.zeros(len(cov_par), 'bool'), \
                             scipy.zeros(n_MF_par, 'bool'))
        if (cov_fixed != None): fixed[0:-n_MF_par] = cov_fixed
        if (MF_fixed != None): fixed[-n_MF_par:] = MF_fixed
        if MF_args == None: MF_args = x[:]
    else:
        merged_par = cov_par[:]
        n_MF_par = 0
        fixed = scipy.zeros(len(cov_par), 'bool')
        if cov_fixed != None: fixed[:] = cov_fixed
    var_par_in = merged_par[fixed == False]
    fixed_par = merged_par[fixed == True]
    args = (x, y, cov_func, cov_typ, MF, n_MF_par, MF_args, fixed, \
            fixed_par, prior)
    var_par_out = \
        sop.fmin(GP_negloglik, var_par_in, args, disp = 0)
    par_out = scipy.copy(merged_par)
    par_out[fixed == False] = var_par_out
    par_out[fixed == True] = fixed_par
    if MF != None:
        return par_out[:-n_MF_par], par_out[-n_MF_par:]
    else:
        return par_out

 def onclick(self, event):
     print 'button=%d, x=%d, y=%d, xdata=%f, ydata=%f'%(event.button, event.x, event.y, event.xdata, event.ydata)
     x=SP.append(self.x, event.xdata)
     self.y=SP.append(self.y, event.ydata)
     self.x= x[:,SP.newaxis]
     self.on_show()
     self.status_text.setText("New data point: x=%f, y=%f"%(event.xdata, event.ydata))

    def update_step(self, input_signal=None, teaching_signal=None):
        """update the network with the given input and teach_output, input_signal and teaching_signal must be a column vector
        notice that input_signal is u(n+1) and output is output(n+1) 
        this step makes state(n) -> state(n+1)
        the x_history is a list of state's state_history , every item is a row vector like (100L,)"""

        if input_signal != None:
            assert input_signal.shape == (self.input_unit_amount, 1)
        if teaching_signal != None:
            assert teaching_signal.shape == (self.output_unit_amount, 1)

        if self.feedback_matrix != None and self.input_matrix != None:
            self.state = self.unit_type_ufunc(sp.dot(self.input_matrix, input_signal) + sp.dot(self.internal_matrix, self.state) + sp.dot(self.feedback_matrix, self.output))
            if teaching_signal == None:
                self.output = sp.dot(self.output_matrix, sp.append(input_signal.T,self.state.T).T)
            else:
                self.output = teaching_signal
        elif self.feedback_matrix != None:
            self.state = self.unit_type_ufunc(sp.dot(self.internal_matrix, self.state) + sp.dot(self.feedback_matrix, self.output))
            if teaching_signal == None:
                self.output = sp.dot(self.output_matrix, self.state)
            else:
                self.output = teaching_signal
        else:
            self.state = self.unit_type_ufunc(sp.dot(self.input_matrix, input_signal) + sp.dot(self.internal_matrix, self.state))
        if input_signal != None:
            self.state_history.append(sp.append(input_signal.T, self.state.T))
        else:
            self.state_history.append(self.state.reshape(-1))
        self.output_history.append(self.output)

 def create_adjacency_matrix(self, data=None):
     r"""
     Returns a weighted adjacency matrix, in CSR format based on the
     product of weight values sharing an interface.
     """
     #
     if data is None:
         data = self.data_vector
     #
     weights = data[self._cell_interfaces[:, 0]]
     weights = 2*weights * data[self._cell_interfaces[:, 1]]
     #
     # clearing any zero-weighted connections
     indices = weights > 0
     interfaces = self._cell_interfaces[indices]
     weights = weights[indices]
     #
     # getting cell connectivity info
     row = interfaces[:, 0]
     col = interfaces[:, 1]
     #
     # append row & col to each other, and weights to itself
     row = sp.append(row, interfaces[:, 1])
     col = sp.append(col, interfaces[:, 0])
     weights = sp.append(weights, weights)
     #
     # Generate sparse adjacency matrix in 'coo' format and convert to csr
     num_blks = self.nx*self.nz
     matrix = sprs.coo_matrix((weights, (row, col)), (num_blks, num_blks))
     #
     return matrix.tocsr()

    def load_shapes(self, narrow_shape, wide_shape):

        self.signal_group = []

        for i in range(10):
            for j in range(10):

                signal = scipy.array([])

                for k in range(5):

                    if codes[i][k] == "n":
                        signal = scipy.append(signal, narrow_shape)
                    else:
                        signal = scipy.append(signal, wide_shape)

                    if codes[j][k] == "n":
                        signal = scipy.append(signal, scipy.zeros(len(narrow_shape)))
                    else:
                        signal = scipy.append(signal, scipy.zeros(len(wide_shape)))

                signal /= pnorm(signal)

                self.signal_group.append(signal)
                self.value_group.append(value_map(i, j))

    def __init__(self, x, y, z, f, boundary = 'natural', dx=0, dy=0, dz=0, bounds_error=True, fill_value=scipy.nan):
        if dx != 0 or dy != 0 or dz != 0:
            raise NotImplementedError(
                "Trispline derivatives are not implemented, do not use tricubic "
                "interpolation if you need to compute magnetic fields!"
            )

        
        self._x = scipy.array(x,dtype=float)
        self._y = scipy.array(y,dtype=float)
        self._z = scipy.array(z,dtype=float)

        self._xlim = scipy.array((x.min(), x.max()))
        self._ylim = scipy.array((y.min(), y.max()))
        self._zlim = scipy.array((z.min(), z.max()))
        self.bounds_error = bounds_error
        self.fill_value = fill_value

        if f.shape != (self._x.size,self._y.size,self._z.size):
            raise ValueError("dimensions do not match f")
            
        if _tricub.ismonotonic(self._x) and _tricub.ismonotonic(self._y) and _tricub.ismonotonic(self._z):
            self._x = scipy.insert(self._x,0,2*self._x[0]-self._x[1])
            self._x = scipy.append(self._x,2*self._x[-1]-self._x[-2])
            self._y = scipy.insert(self._y,0,2*self._y[0]-self._y[1])
            self._y = scipy.append(self._y,2*self._y[-1]-self._y[-2])
            self._z = scipy.insert(self._z,0,2*self._z[0]-self._z[1])
            self._z = scipy.append(self._z,2*self._z[-1]-self._z[-2])

        
        self._f = scipy.zeros(scipy.array(f.shape)+(2,2,2))
        self._f[1:-1,1:-1,1:-1] = scipy.array(f) # place f in center, so that it is padded by unfilled values on all sides
        
        if boundary == 'clamped':
            # faces
            self._f[(0,-1),1:-1,1:-1] = f[(0,-1),:,:] 
            self._f[1:-1,(0,-1),1:-1] = f[:,(0,-1),:]
            self._f[1:-1,1:-1,(0,-1)] = f[:,:,(0,-1)]
            #verticies
            self._f[(0,0,-1,-1),(0,-1,0,-1),1:-1] = f[(0,0,-1,-1),(0,-1,0,-1),:] 
            self._f[(0,0,-1,-1),1:-1,(0,-1,0,-1)] = f[(0,0,-1,-1),:,(0,-1,0,-1)]
            self._f[1:-1,(0,0,-1,-1),(0,-1,0,-1)] = f[:,(0,0,-1,-1),(0,-1,0,-1)]
            #corners
            self._f[(0,0,0,0,-1,-1,-1,-1),(0,0,-1,-1,0,0,-1,-1),(0,-1,0,-1,0,-1,0,-1)] = f[(0,0,0,0,-1,-1,-1,-1),(0,0,-1,-1,0,0,-1,-1),(0,-1,0,-1,0,-1,0,-1)]
        elif boundary == 'natural':
            # faces
            self._f[(0,-1),1:-1,1:-1] = 2*f[(0,-1),:,:] - f[(1,-2),:,:]
            self._f[1:-1,(0,-1),1:-1] = 2*f[:,(0,-1),:] - f[:,(1,-2),:]
            self._f[1:-1,1:-1,(0,-1)] = 2*f[:,:,(0,-1)] - f[:,:,(1,-2)]
            #verticies
            self._f[(0,0,-1,-1),(0,-1,0,-1),1:-1] = 4*f[(0,0,-1,-1),(0,-1,0,-1),:] - f[(1,1,-2,-2),(0,-1,0,-1),:] - f[(0,0,-1,-1),(1,-2,1,-2),:] - f[(1,1,-2,-2),(1,-2,1,-2),:]
            self._f[(0,0,-1,-1),1:-1,(0,-1,0,-1)] = 4*f[(0,0,-1,-1),:,(0,-1,0,-1)] - f[(1,1,-2,-2),:,(0,-1,0,-1)] - f[(0,0,-1,-1),:,(1,-2,1,-2)] - f[(1,1,-2,-2),:,(1,-2,1,-2)]
            self._f[1:-1,(0,0,-1,-1),(0,-1,0,-1)] = 4*f[:,(0,0,-1,-1),(0,-1,0,-1)] - f[:,(1,1,-2,-2),(0,-1,0,-1)] - f[:,(0,0,-1,-1),(1,-2,1,-2)] - f[:,(1,1,-2,-2),(1,-2,1,-2)]
            #corners
            self._f[(0,0,0,0,-1,-1,-1,-1),(0,0,-1,-1,0,0,-1,-1),(0,-1,0,-1,0,-1,0,-1)] = 8*f[(0,0,0,0,-1,-1,-1,-1),(0,0,-1,-1,0,0,-1,-1),(0,-1,0,-1,0,-1,0,-1)] -f[(1,1,1,1,-2,-2,-2,-2),(0,0,-1,-1,0,0,-1,-1),(0,-1,0,-1,0,-1,0,-1)] -f[(0,0,0,0,-1,-1,-1,-1),(1,1,-2,-2,1,1,-2,-2),(0,-1,0,-1,0,-1,0,-1)] -f[(0,0,0,0,-1,-1,-1,-1),(0,0,-1,-1,0,0,-1,-1),(1,-2,1,-2,1,-2,1,-2)] -f[(1,1,1,1,-2,-2,-2,-2),(1,1,-2,-2,1,1,-2,-2),(0,-1,0,-1,0,-1,0,-1)] -f[(0,0,0,0,-1,-1,-1,-1),(1,1,-2,-2,1,1,-2,-2),(1,-2,1,-2,1,-2,1,-2)] -f[(1,1,1,1,-2,-2,-2,-2),(0,0,-1,-1,0,0,-1,-1),(1,-2,1,-2,1,-2,1,-2)] -f[(1,1,1,1,-2,-2,-2,-2),(1,1,-2,-2,1,1,-2,-2),(1,-2,1,-2,1,-2,1,-2)]

        self._regular = False
        if _tricub.isregular(self._x) and _tricub.isregular(self._y) and _tricub.isregular(self._z):
            self._regular = True

def read_mat_files(features_basename, labels_fname, camname_fname, actname_fname, partiname_fname):
    """docstring for read_mat_file"""
    
    print "reading features"
    tic = time.time()
    
    f = h5py.File(features_basename + '_part1.mat', 'r')
    ff = f["myData"]
    features1 = ff[:,0:N_1STCHUNK].T
    features2 = ff[:,N_1STCHUNK+1:N_2NDCHUNK].T
    features3 = ff[:,N_2NDCHUNK+1:].T
    features = sp.append(features1, features2,1)
    features = sp.append(features, features3,1)
#    features = sp.array(ff).T
    import ipdb; ipdb.set_trace()
    for nn in range(2,N_PARTS+1):
        f = h5py.File(features_basename + '_part' + str(nn)+ '.mat', 'r')
        import ipdb; ipdb.set_trace()
        ff = f["myData"]
        import ipdb; ipdb.set_trace()
        temp = sp.array(ff).T
        import ipdb; ipdb.set_trace()
        features = sp.append(features, temp,1)
        print nn
    
    print "time taken :", time.time() - tic, 'seconds'

    print "reading participant names"
    tic = time.time()
    partiNames = io.loadmat(partiname_fname)['myPartis']
    partiNames = sp.array([str(partiNames[i][0][0]) for i in xrange(partiNames.shape[0])])
    print "time taken :", time.time() - tic, 'seconds'

    print "reading labels"
    tic = time.time()
    labels = io.loadmat(labels_fname)['labels']
    labels = sp.array([str(labels[i][0][0]) for i in xrange(labels.shape[0])])
    print "time taken :", time.time() - tic, 'seconds'
    
    print "reading camera names"
    tic = time.time()
    camNames = io.loadmat(camname_fname)['myCams']
    camNames = sp.array([str(camNames[i][0][0]) for i in xrange(camNames.shape[0])])
    print "time taken :", time.time() - tic, 'seconds'

    
    print "reading action names"
    tic = time.time()
    actNames = io.loadmat(actname_fname)['myActs']
    actNames = sp.array([str(actNames[i][0][0]) for i in xrange(actNames.shape[0])])
    print "time taken :", time.time() - tic, 'seconds'
            


    # few sanity checks
    #assert(not features.isnan().any())
    #assert(not features.isinf().any())
    
    return features, labels, camNames, actNames, partiNames

 def matrix_rows(self):
     if self.node==None or self.Phi==None:
         raise NameError('DoF of point %(self.id)d has not been set.')
     
     cols = scipy.array([i*3 for i in self.nodes])
     data = scipy.array([phi for phi in self.Phi])
     
     cols = scipy.append(cols, 3*self.node)
     data = self.binding_weight*scipy.append(data, -1)
     
     return [[cols, data]], [0,1,2], None

def prob3():
	rate1,sig1 = wavfile.read('chopinw.wav')
	n = sig1.shape[0]
	rate2,sig2 = wavfile.read('balloon.wav')
	m = sig2.shape[0]
	sig1 = sp.append(sig1,sp.zeros((m,2)))
	sig2 = sp.append(sig2,sp.zeros((n,2)))
	f1 = sp.fft(sig1)
	f2 = sp.fft(sig2)
	out = sp.ifft((f1*f2))
	out = sp.real(out)
	scaled = sp.int16(out/sp.absolute(out).max() * 32767)
	wavfile.write('test.wav',rate1,scaled)

    def TOBEFIXED_build_k_coo_sub(self):
        import scipy
        import scipy.sparse as ss
        import alg3dpy.scipy_sparse as assparse
        #FIXME not considering pos to build the matrix!!!
        self.k_coo_sub = {}
        for sub in self.subcases.values():
            dim = 6*len( self.k_pos ) - \
                    len( self.index_to_delete[ sub.id ] )
            data = scipy.zeros(0, dtype='float64')
            row =  scipy.zeros(0, dtype='int64')
            col =  scipy.zeros(0, dtype='int64')
            for elem in self.elemdict.values():
                numg = len( elem.grids )
                for i in xrange( numg ):
                    gi = elem.grids[ i ]
                    offseti = gi.k_offset[ sub.id ]
                    consi = set( [] )
                    if sub.id in gi.cons.keys():
                        consi = gi.cons[ sub.id ]
                    for j in xrange( numg ):
                        gj = elem.grids[ j ]
                        offsetj = gj.k_offset[ sub.id ]
                        consj = set( [] )
                        if sub.id in gj.cons.keys():
                            consj = gj.cons[ sub.id ]
                        cons = consi | consj
                        index_to_delete = [ (c-1) for c in cons ]
                        k_grid = assparse.in_sparse(elem.k,i*6,i*6+5,j*6,j*6+5)

                        if len(index_to_delete) < 6:
                            k = k_grid
                            for d, r, c, in zip( k.data , k.row, k.col ):
                                #FIXME remove the search below
                                if not r in index_to_delete:
                                    sub_r = 0
                                    for k in index_to_delete:
                                        if r > k:
                                            sub_r += 1
                                    if not c in index_to_delete:
                                        sub_c = 0
                                        for m in index_to_delete:
                                            if c > m:
                                                sub_c += 1
                                        newr = r + gi.pos * 6 - offseti - sub_r
                                        newc = c + gj.pos * 6 - offsetj - sub_c
                                        data = scipy.append( data , d    )
                                        row  = scipy.append( row  , newr )
                                        col  = scipy.append( col  , newc )
            k_coo = ss.coo_matrix(( data, (row,col) ), shape=(dim,dim))
            self.k_coo_sub[sub.id] = k_coo

def center_orbit(sol,N,Dim):
    how_much = 20
    #sol[:,N:(2*N)]=sol[:,N:(2*N)]%(2.0*sp.pi)

    print('looking for orbit for with smallest velocity amplitude. This only works for 1D systems')

    # make an array of the mean velocity squared values. The array will be orderd as the solution is
    # ordered
    mean_vel_arr = sp.array([])
    # we also need an array of mean_positions to figure out where these things are
    mean_pos_arr = sp.array([])
    count = 0
    while count < N:
        mean_vel_arr = sp.append(mean_vel_arr,(sol[(-len(sol)/how_much):,count]**2).mean())
        print('mean vel appended: ' +str((sol[(-len(sol)/how_much):,count]**2).mean()))

        # so this is a little trickey... 
        # We also need the standard deveation becuase if the paticle is oscilating at the boundary
        # so it is right around 0 AND 2pi then its mean position is pi. We can use the standard
        # deviation to tell weather or not it is actualy at pi or at the boundary. The standard
        # deveation should certinaly be less than pi/2 unless there is only 1 particle.
        cur_mean_pos = (sol[(-len(sol)/how_much):,N+count]).mean()
        cur_mean_std = (sol[(-len(sol)/how_much):,N+count]).std()
        if (cur_mean_std > 3.0):
            print('foud particle oscilating at boundary')
            print('standard deviation: ' +str(cur_mean_std))
            cur_mean_pos = 0.0
        mean_pos_arr = sp.append(mean_pos_arr,cur_mean_pos)
        print('mean pos appended: ' +str(cur_mean_pos))

        count+=1

    print('mean pos arr: ' +str(mean_pos_arr))
    print('mean vel arr: ' +str(mean_vel_arr))

    # which particle is the one with the smallest v^2 trajectory? the_one will be the index of this
    # particle
    the_one = sp.argmin(mean_vel_arr)
    print('orbit with smallest velocity amplitued: '+str(the_one))
    print('mean vel of the_one: ' +str(mean_vel_arr[the_one]))
    print('mean pos of the_one: ' +str(mean_pos_arr[the_one]))

    # Now we need to shift everything to get it into the center. 
    # there are a few ways to do this. We are going to try this one but it might not be the best one

    shift = sp.pi-mean_pos_arr[the_one]
    # now shift everything by the right amount
    sol[:,N:] = (sol[:,N:]+shift)%(2.0*sp.pi)

    return sol

def createNormalizedDataSets():
    xw1 = norm(loc=0.3, scale=.15).rvs(20)
    yw1 = norm(loc=0.3, scale=.15).rvs(20)

    xw2 = norm(loc=0.7, scale=.15).rvs(20)
    yw2 = norm(loc=0.7, scale=.15).rvs(20)

    xw3 = norm(loc=0.2, scale=.15).rvs(20)
    yw3 = norm(loc=0.8, scale=.15).rvs(20)

    x = sp.append(sp.append(xw1, xw2), xw3)
    y = sp.append(sp.append(yw1, yw2), yw3)

    return x, y

 def shift(self, items):
     """Shift all buffers up or down a defined number of items on offset axis.
     Negative values indicate backward shift."""
     if items == 0:
         return
     self.offset += items
     for buffername, _  in self.bufferlist:
         buf = getattr(self, buffername)
         assert abs(items) <= len(buf), "Cannot shift further than length of buffer."
         fill = zeros((abs(items), len(buf[0])))
         if items < 0:
             buf[:] = append(buf[-items:], fill, 0)
         else:
             buf[:] = append(fill ,buf[0:-items] , 0)

 def square_wave(self,input):
     if type(input)==float:
         if input%(2.0*pl.pi)<= pl.pi:
             output = 1.0
         if input%(2.0*pl.pi)> pl.pi:
             output = -1.0
     if (type(input)==pl.ndarray)or(type(input)==list):
         output = pl.array([])
         for i,j in enumerate(input):
             if j%(2.0*pl.pi)<= pl.pi:
                 output = pl.append(output,1.0)
             if j%(2.0*pl.pi)> pl.pi:
                 output = pl.append(output,-1.0)
         print('square wave output: ' + str(output)) 
     return output

def nullSpaceBasis(A):
    """
    This funciton will find the basis of the null space of the matrix A.

    Inputs:
        A: The matrix you want the basis for
    Outputs:
        A numpy matrix containing the vectors as row vectors.

    Notes:
        If A is an empty matrix, an empty matrix is returned.

    """
    if A:
        U,s, Vh = la.svd(A)
        vecs = np.array([])
        toAppend = A.shape[1] -s.size
        s = sp.append(s,sp.zeros((1,toAppend)))
        for i in range(0,s.size):
            if s[i]==0:
                vecs = Vh[-toAppend:,:]
        if vecs.size ==0:
            vecs = sp.zeros((1,A.shape[1]))
        return sp.mat(vecs)
    else:
        return sp.zeros((0,0))

def coordinates_to_voxel_idx(coords_xyz, masker):
	# transform to homogeneous coordinates
	coords_h_xyz = sp.append(coords_xyz, ones([1,coords_xyz.shape[1]]),axis=0)
	
	# apply inverse affine transformation to get homogeneous coordinates in voxel space
	inv_transf = sp.linalg.inv(masker.volume.get_affine())
	coords_h_voxel_space = inv_transf.dot(coords_h_xyz)
	coords_h_voxel_space = sp.rint(coords_h_voxel_space).astype(int)
	
	# remove homogeneous dimension
	coords_voxel_space = coords_h_voxel_space[0:-1,:]
	
	# convert coordinates to idcs in a flattened voxel space
	flattened_idcs = sp.ravel_multi_index(coords_voxel_space, masker.dims)
	
	# check if there is any study data for the flattened idcs
	voxel_idcs = sp.zeros((1,len(flattened_idcs)),dtype=int64)
	for i in range(0,len(flattened_idcs)):
		idcs = find(masker.in_mask == flattened_idcs[i])
		if len(idcs > 0):
			voxel_idcs[0,i] = find(masker.in_mask == flattened_idcs[i])
		else:
			voxel_idcs[0,i] = nan
			
	return voxel_idcs