def write_dp_item(coe, dp, palette, width, input_, pp):
    assert(1 <= width <= 31)
    coe.write(dp, low(pp))
    coe.write(dp+1, eval('0b' + np.binary_repr(palette, 3)
                              + np.binary_repr(-width, 5)))
    coe.write(dp+2, high(pp))
    coe.write(dp+3, input_)

	def write_adc(self,addr,data):
		SCLK = 0x200
		CS = self.chip_select
		IDLE = SCLK
		SDA_SHIFT = 8
		self.snap.write_int('adc16_controller',IDLE,offset=0,blindwrite=True)
		for i in range(8):
			addr_bit = (addr>>(8-i-1))&1
			state = (addr_bit<<SDA_SHIFT) | CS
			self.snap.write_int('adc16_controller',state,offset=0,blindwrite=True)
			logging.debug("Printing address state written to adc16_controller, offset=0, clock low")
			logging.debug(np.binary_repr(state,width=32))
	#		print(np.binary_repr(state,width=32))
			state = (addr_bit<<SDA_SHIFT) | CS | SCLK
			self.snap.write_int('adc16_controller',state,offset=0,blindwrite=True)
			logging.debug("Printing address state written to adc16_controller, offset=0, clock high")
			logging.debug(np.binary_repr(state,width=32))
	#		print(np.binary_repr(state,width=32))
		for j in range(16):
			data_bit = (data>>(16-j-1))&1
			state = (data_bit<<SDA_SHIFT) | CS
			self.snap.write_int('adc16_controller',state,offset=0,blindwrite=True)
			logging.debug("Printing data state written to adc16_controller, offset=0, clock low")
			logging.debug(np.binary_repr(state,width=32))
	#		print(np.binary_repr(state,width=32))
			state =( data_bit<<SDA_SHIFT) | CS | SCLK	
			self.snap.write_int('adc16_controller',state,offset=0,blindwrite=True)		
			logging.debug("Printing data address state written to adc16_controller, offset=0, clock high")
			logging.debug(np.binary_repr(state,width=32))
	#		print(np.binary_repr(state,width=32))
		
		self.snap.write_int('adc16_controller',IDLE,offset=0,blindwrite=True)

def embed(cover,secret,pos,skip):
    file=open("in.txt","w")
    multiple=False
    coverMatrix=pgm_to_mat(cover)
    secretMatrix=pgm_to_mat(secret)
    stegoMatrix=np.zeros(np.shape(coverMatrix), dtype=np.complex_)
    np.copyto(stegoMatrix,coverMatrix)
    dummy=""
    if(skip<1):
        skip=1
        multiple=True
    for a in range(0,len(secretMatrix)):
        for b in range(0,len(secretMatrix)):
            dummy+=np.binary_repr(secretMatrix[a][b],width=8)
            #file.write(np.binary_repr(secretMatrix[a][b],width=8)+"\n")
    index=0
    for a in range(0,len(stegoMatrix)*len(stegoMatrix),skip):
        rown=int(a % len(stegoMatrix))
        coln=int(a / len(stegoMatrix))
        if(index>=len(dummy)):
            break
        stegoMatrix[coln][rown] = ( int(coverMatrix[coln][rown]) & ~(1 << hash(coln,rown,pos) )) | (int(dummy[index],2) << hash(coln,rown,pos))
        index += 1
        if(multiple):
            stegoMatrix[coln][rown] = (int(stegoMatrix[coln][rown]) & ~(1 << (3-hash(coln, rown, pos)))) | ( int(dummy[index], 2) << (3-hash(coln, rown, pos)))
            index += 1
        file.write(np.binary_repr(int(stegoMatrix[coln][rown]), 8) + "\n")
    return stegoMatrix

def encode_gray(number):
    
    s1 = numpy.binary_repr(number);
    s2 = "0" + numpy.binary_repr(number)[:-1];
    
    
    return "".join(["1" if s1[i] != s2[i] else "0" for i in range(len(s1))]);

def convolve(image,psf,doPSF=True,edgeCheck=False):
    """
    A reasonably fast convolution routine that supports re-entry with a
    pre-FFT'd PSF. Returns the convolved image and the FFT'd PSF.
    """
    datadim1 = image.shape[0]
    datadim2 = image.shape[1]
    if datadim1!=datadim2:
        ddim = max(datadim1,datadim2)
        s = numpy.binary_repr(ddim-1)
        s = s[:-1]+'0' # Guarantee that padding is used
    else:
        ddim = datadim1
        s = numpy.binary_repr(ddim-1)
    if s.find('0')>0:
        size = 2**len(s)
        if edgeCheck==True and size-ddim<8:
            size*=2
        boxd = numpy.zeros((size,size))
        r = size-datadim1
        r1 = r2 = r/2
        if r%2==1:
            r1 = r/2+1
        c = size-datadim2
        c1 = c2 = c/2
        if c%2==1:
            c1 = c/2+1
        boxdslice = (slice(r1,datadim1+r1),slice(c1,datadim2+c1))
        boxd[boxdslice] = image
    else:
        boxd = image

    if doPSF:
        # Pad the PSF to the image size
        boxp = boxd*0.
        if boxd.shape[0]==psf.shape[0]:
            boxp = psf.copy()
        else:
            r = boxp.shape[0]-psf.shape[0]
            r1 = r/2+1
            c = boxp.shape[1]-psf.shape[1]
            c1 = c/2+1
            boxpslice = (slice(r1,psf.shape[0]+r1),slice(c1,psf.shape[1]+c1))
            boxp[boxpslice] = psf.copy()
        # Store the transform of the image after the first iteration
        a = (numpy.fft.rfft2(boxp))
    else:
        a = psf
        # PSF transform and multiplication
    b = a*numpy.fft.rfft2(boxd)
    # Inverse transform, including phase-shift to put image back in center;
    #   this removes the requirement to do 2x zero-padding so makes things
    #   go a bit quicker.
    b = numpy.fft.fftshift(numpy.fft.irfft2(b)).real
    # If the image was padded, remove the padding
    if s.find('0')>0:
        b = b[boxdslice]

    return b,a

def displaystates(psi, N=5, pop=False):
  ''' print the population of each state '''
  for i in range(len(psi)):
    if np.around(abs(psi[i]), 5) > 0:
        if pop:
            print np.binary_repr(i).zfill(N), ": ", np.around(psi[i],3)
        else:
            print np.binary_repr(i).zfill(N), ": ", np.around(abs(psi[i]**2),3), np.around(psi[i], 5)

def read_binaries(tks1, tks2, encoding):
  dtks1 = decode(tks1, encoding)
  if not dtks1 or not validate_length(dtks1):
    error("The first samples have different sizes")
    return
  dtks2 = decode(tks2, encoding)
  if not dtks2 or not validate_length(dtks2):
    error("The second samples have different sizes")
    return
  btks1 = [ "".join([np.binary_repr(ord(c), width=8) for c in tk ]) for tk in dtks1 ]
  btks2 = [ "".join([np.binary_repr(ord(c), width=8) for c in tk ]) for tk in dtks2 ]
  return btks1, btks2, "01"

    def loadFIRcoeffs(self):
        N_freqs = len(map(float, unicode(self.textedit_DACfreqs.toPlainText()).split()))
        taps = 26
        
        for ch in range(N_freqs):
            # If the resonator's attenuation is >=99 then its FIR should be zeroed
            if self.zeroChannels[ch]:
                lpf = numpy.array([0.]*taps)*(2**11-1)
                print 'deleted ch ',ch
            else:
                lpf = numpy.array(self.fir)*(2**11-1)
                print ch
                #lpf = numpy.array([1.]+[0]*(taps-1))*(2**11-1)
        #    26 tap, 25 us matched fir
                #lpf = numpy.array([0.0875788844768 , 0.0840583257978 , 0.0810527406206 , 0.0779008825067 , 0.075106964962 , 0.0721712998256 , 0.0689723729398 , 0.066450095496 , 0.0638302570705 , 0.0613005685486 , 0.0589247737004 , 0.0565981917436 , 0.0544878914297 , 0.0524710948658 , 0.0503447054014 , 0.0483170854189 , 0.0463121066637 , 0.044504238059 , 0.0428469827102 , 0.0410615366471 , 0.0395570640218 , 0.0380071830756 , 0.0364836787854 , 0.034960959124 , 0.033456372241 , 0.0321854467182])*(2**11-1)
                #26 tap, 20 us matched fir
                #lpf = numpy.array([ 0.102806030245 , 0.097570344415 , 0.0928789946181 , 0.0885800360545 , 0.0841898850361 , 0.079995145104 , 0.0761649967857 , 0.0724892663141 , 0.0689470889358 , 0.0657584886557 , 0.0627766233242 , 0.0595952531565 , 0.0566356208278 , 0.053835736579 , 0.0510331408751 , 0.048623806127 , 0.0461240096904 , 0.0438134132285 , 0.0418265743203 , 0.0397546477453 , 0.0377809254888 , 0.0358044897245 , 0.0338686929847 , 0.0321034547839 , 0.0306255734188 , 0.0291036235859 ])*(2**11-1)
                #26 tap, 30 us matched fir
                #lpf = numpy.array([ 0.0781747107378 , 0.0757060398243 , 0.0732917718492 , 0.0708317694778 , 0.0686092845217 , 0.0665286923521 , 0.0643467681477 , 0.0621985982971 , 0.0600681642401 , 0.058054873199 , 0.0562486467178 , 0.0542955553149 , 0.0527148880657 , 0.05096365681 , 0.0491121116212 , 0.0474936094733 , 0.0458638771941 , 0.0443219286645 , 0.0429290438102 , 0.0415003391096 , 0.0401174498302 , 0.0386957715665 , 0.0374064708747 , 0.0362454802408 , 0.0350170176804 , 0.033873302383 ])*(2**11-1)
                #lpf = lpf[::-1]
                #    26 tap, lpf, 250 kHz,
                #lpf = numpy.array([-0 , 0.000166959420533 , 0.00173811663844 , 0.00420937801998 , 0.00333739357391 , -0.0056305703275 , -0.0212738104942 , -0.0318529375832 , -0.0193635986879 , 0.0285916612022 , 0.106763943766 , 0.18981814328 , 0.243495321192 , 0.243495321192 , 0.18981814328 , 0.106763943766 , 0.0285916612022 , -0.0193635986879 , -0.0318529375832 , -0.0212738104942 , -0.0056305703275 , 0.00333739357391 , 0.00420937801998 , 0.00173811663844 , 0.000166959420533 , -0])*(2**11-1)
                #    26 tap, lpf, 125 kHz.
                #lpf = numpy.array([0 , -0.000431898216436 , -0.00157886921107 , -0.00255492263971 , -0.00171727439076 , 0.00289724121972 , 0.0129123447233 , 0.0289345497995 , 0.0500906370566 , 0.0739622085341 , 0.0969821586979 , 0.115211955161 , 0.125291869266 , 0.125291869266 , 0.115211955161 , 0.0969821586979 , 0.0739622085341 , 0.0500906370566 , 0.0289345497995 , 0.0129123447233 , 0.00289724121972 , -0.00171727439076 , -0.00255492263971 , -0.00157886921107 , -0.000431898216436 , -0])*(2**11-1)
                #    Generic 40 tap matched filter for 25 us lifetime pulse
                #lpf = numpy.array([0.153725595011 , 0.141052390733 , 0.129753816201 , 0.119528429291 , 0.110045314901 , 0.101336838027 , 0.0933265803805 , 0.0862038188673 , 0.0794067694409 , 0.0729543134914 , 0.0674101836798 , 0.0618283869464 , 0.0567253144676 , 0.0519730940444 , 0.047978953698 , 0.043791412767 , 0.0404560656757 , 0.0372466775252 , 0.0345000956808 , 0.0319243455811 , 0.0293425115323 , 0.0268372778298 , 0.0245216835234 , 0.0226817116475 , 0.0208024488535 , 0.0189575043357 , 0.0174290665862 , 0.0158791788119 , 0.0144611054123 , 0.0132599563305 , 0.0121083419203 , 0.0109003580368 , 0.0100328742978 , 0.00939328253743 , 0.00842247241585 , 0.00789304712484 , 0.00725494259117 , 0.00664528407122 , 0.00606688645845 , 0.00552041438208])*(2**11-1)                
                #lpf = lpf[::-1]

            for n in range(taps/2):
                coeff0 = int(lpf[2*n])
                coeff1 = int(lpf[2*n+1])
                coeff0 = numpy.binary_repr(int(lpf[2*n]), 12)
                coeff1 = numpy.binary_repr(int(lpf[2*n+1]), 12)
                coeffs = int(coeff1+coeff0, 2)
                coeffs_bin = struct.pack('>l', coeffs)
                register_name = 'FIR_b' + str(2*n) + 'b' + str(2*n+1)
                self.roach.write(register_name, coeffs_bin)
                self.roach.write_int('FIR_load_coeff', (ch<<1) + (1<<0))
                self.roach.write_int('FIR_load_coeff', (ch<<1) + (0<<0))
        
        # Inactive channels will also be zeroed.
        lpf = numpy.array([0.]*taps)
        for ch in range(N_freqs, 256):
            for n in range(taps/2):
                #coeffs = struct.pack('>h', lpf[2*n]) + struct.pack('>h', lpf[2*n+1])
                coeffs = struct.pack('>h', lpf[2*n+1]) + struct.pack('>h', lpf[2*n])
                register_name = 'FIR_b' + str(2*n) + 'b' + str(2*n+1)
                self.roach.write(register_name, coeffs)
                self.roach.write_int('FIR_load_coeff', (ch<<1) + (1<<0))
                self.roach.write_int('FIR_load_coeff', (ch<<1) + (0<<0))
                
        print 'done loading fir.'
        self.status_text.setText('FIRs loaded')

 def send_32(self, inbits):
     print "inbits", inbits
     programming_bits = np.bitwise_and(inbits,65535) #bits 0 16 
     address_branch = (np.bitwise_and(inbits,8323072)>>16) #bits 17 to 22 
     print "send stuff"
     print "address_branch", np.binary_repr(address_branch)
     print "programming_bits", np.binary_repr(programming_bits)
     print "address_branch", (address_branch)
     print "programming_bits", (programming_bits<<7)
     final_address = (programming_bits<<7) + (address_branch) + 2**31
     print "final address", final_address
     biasusb_wrap.send_32(int(final_address))

 def get_key_slow(self, iarr, level=None):
     if level is None:
         level = self.level
     i1, i2, i3 = iarr
     rep1 = np.binary_repr(i1, width=self.level)
     rep2 = np.binary_repr(i2, width=self.level)
     rep3 = np.binary_repr(i3, width=self.level)
     inter = np.zeros(self.level*3, dtype='c')
     inter[self.dim_slices[0]] = rep1
     inter[self.dim_slices[1]] = rep2
     inter[self.dim_slices[2]] = rep3
     return int(inter.tostring(), 2)

    def ExecuteQP (self, qubits, funcao, memory, customMatrices = []): # Method for build the structures (Pages, VPPs and sizesList) of the Quantum Processes.
        if qubits > 5:
            sizeVPP = 4
            qtdPages = qubits/sizeVPP
            rest = qubits%sizeVPP
            if rest > 1:
                qtdPages += 1
        else:
            qtdPages = 1
            sizeVPP = qubits
            rest = 0
        
        Pages = []
        opIndex = 0
        sizesList = []
        for pageId in range (qtdPages):
            if rest == 1:
                qtdFunctions = sizeVPP + 1
                rest = 0
            elif rest > 1 and pageId == qtdPages - 1:
                qtdFunctions = rest
            else:
                qtdFunctions = sizeVPP
        
            Lvpp = []
            for VPPIndex in range (2**qtdFunctions): ## Creates each VPP of a Lvpp
                listOp = self.StringToList(funcao, qtdFunctions, opIndex)
                param1 = numpy.binary_repr(VPPIndex, qtdFunctions) ## First parameter of each function to fill the QPPs

                zero = numpy.complex(0)
                pos = 0
                list = []
                for tupleIndex in range (2**qtdFunctions): ## Creates each tuple of a VPP
                    param2 = numpy.binary_repr(tupleIndex,qtdFunctions)
                    temp = numpy.complex(1)
                    op = 0
            
                    while temp != zero and op < qtdFunctions:
                        temp = temp * self.getValue(listOp[op], int(param1[op:op+1:]), int(param2[op:op+1:]))
                        op += 1
                        
                    if temp != zero:
                        list.append([temp,pos])
                    pos += 1
                
                Lvpp.append(list)
            Pages.append(Lvpp)
            opIndex += qtdFunctions
            sizesList.append(2**(qubits-opIndex))

        self.ApplyValuesForQP(Pages,sizesList,memory,numpy.complex(1),0,0,numpy.binary_repr(0,qubits), qubits)

 def _send_32(self, in_bits, debug=False):
     programming_bits = np.bitwise_and(in_bits,65535) #bits 0 16 
     address_branch = (np.bitwise_and(in_bits,8323072)>>16) #bits 17 to 22 
     final_address = (programming_bits<<7) + (address_branch) + 2**31
     if debug:
         print "in_bits", in_bits
         print "send stuff"
         print "address_branch", np.binary_repr(address_branch)
         print "programming_bits", np.binary_repr(programming_bits)
         print "address_branch", (address_branch)
         print "programming_bits", (programming_bits<<7)
         print "final address", final_address
     self._client.send(str([0,final_address]))
     time.sleep(0.001)

 def get_instruction(self):
     """
     A generator to get the instruction binary code for
     respective assembly code
     """
     variable_address = 16
     for line in self._lines:
         if '@' in line:
             # for A instruction
             instruction = '0'
             # if @21 then its direct accessing
             if line[1:].isdigit():
                 instruction += numpy.binary_repr(int(line[1:]), 15)
             # Predefined variables
             elif line[1:] in Parser.SYMBOLS:
                 instruction += numpy.binary_repr(int(Parser.SYMBOLS[line[1:]]), 15)
             # user defined variables
             else:
                 Parser.SYMBOLS[line[1:]] = variable_address
                 instruction += numpy.binary_repr(variable_address, 15)
                 variable_address += 1
             yield instruction
         else:
             # for C instruction all the null cases are equal to '0'
             # hence initialized it with zero
              
             dest, rest, comp, jump = '0', '0', '0', '0'
             # to separate destination and rest of the code
             dest_rest = (['0', '0'] + list(line.split('=')))  
             dest_rest.reverse()
             rest, dest = dest_rest[:2]
             
             # C Instruction fixed starting values
             instruction = '111'
             
             # from the rest to get computation and jump
             comp_jump = (['0', '0'] + list(rest.split(';'))) 
             comp_jump.reverse()
             
             # if there is no JMP instruction and else
             if len(list(rest.split(';'))) == 1:
                 comp, jump = comp_jump[0:2]
             else:
                 jump, comp = comp_jump[0:2]
                 
             instruction += Parser.INSTRUCTIONS[comp] +\
                             Parser.DESTINATION[dest] +\
                             Parser.JUMP[jump]
             yield instruction

def prep(image,psf):
    datadim1 = image.shape[0]
    datadim2 = image.shape[1]
    if datadim1!=datadim2:
        ddim = max(datadim1,datadim2)
        s = numpy.binary_repr(ddim-1)
        s = s[:-1]+'0' # Guarantee that padding is used
    else:
        ddim = datadim1
        s = numpy.binary_repr(ddim-1)
    if s.find('0')>0:
        size = 2**len(s)
        boxd = numpy.zeros((size,size))
        r = size-datadim1
        r1 = r2 = r/2
        if r%2==1:
            r1 = r/2+1
        c = size-datadim2
        c1 = c2 = c/2
        if c%2==1:
            c1 = c/2+1
        boxdslice = (slice(r1,datadim1+r1),slice(c1,datadim2+c1))
        boxd[boxdslice] = image
    else:
        boxd = image

    boxp = boxd*0.
    if boxd.shape[0]==psf.shape[0]:
        boxp = psf.copy()
    else:
        r = boxp.shape[0]-psf.shape[0]
        r1 = r/2+1
        c = boxp.shape[1]-psf.shape[1]
        c1 = c/2+1
        boxpslice = (slice(r1,psf.shape[0]+r1),slice(c1,psf.shape[1]+c1))
        boxp[boxpslice] = psf.copy()

    from pyfft.cuda import Plan
    import pycuda.driver as cuda
    from pycuda.tools import make_default_context
    import pycuda.gpuarray as gpuarray
    cuda.init()
    context = make_default_context()
    stream = cuda.Stream()

    plan = Plan(boxp.shape,stream=stream)
    gdata = gpuarray.to_gpu(boxp.astype(numpy.complex64))
    plan.execute(gdata)
    return gdata,boxd.shape,boxdslice,plan,stream

    def col_names(self):
        mode = self.parameters_dict.StateReadout.readout_mode     
        names = np.array(range(self.output_size())[::-1])+1
         
        if mode == 'pmt':
            if self.output_size==1:
                dependents = [('', 'prob dark ', '')]
            else:
                dependents = [('', 'num dark {}'.format(x), '') for x in names ]
                
        if mode == 'pmt_states':
            if self.output_size==1:
                dependents = [('', 'prob dark ', '')]
            else:
                dependents = [('', ' {} dark ions'.format(x-1), '') for x in names ]
                
        if mode == 'pmt_parity':
            if self.output_size==1:
                dependents = [('', 'prob dark ', '')]
            else:
                dependents = [('', ' {} dark ions'.format(x-1), '') for x in names[1:] ]
                
            dependents.append(('', 'Parity', ''))        
                
        if mode == 'camera':
            dependents = [('', ' prob ion {}'.format(x), '') for x in range(self.output_size())]
            
        if mode == 'camera_states':
            num_of_ions=int(self.parameters_dict.IonsOnCamera.ion_number)
            names = range(2**num_of_ions)
            dependents=[]
            for name in names:
                temp= np.binary_repr(name,width=num_of_ions)
                temp = self.binary_to_state(temp)
                temp=('', 'Col {}'.format(temp), '')
                dependents.append(temp)
        
        if mode == 'camera_parity':
            num_of_ions=int(self.parameters_dict.IonsOnCamera.ion_number)
            names = range(2**num_of_ions)
            dependents=[]
            for name in names:
                temp= np.binary_repr(name,width=num_of_ions)
                temp = self.binary_to_state(temp)
                temp=('', 'Col {}'.format(temp), '')
                dependents.append(temp)
            dependents.append(('', 'Parity', ''))

        return  dependents

def intToArray(i, length=0):
    """Convert an unsigned integer to a binary array.
    Args:
        i: unsigned integer
        length: padding to length (default: 0)
    Returns:
        binary array
    """
    if length > 0:
        s = np.binary_repr(i, width=length)
    else:
        s = np.binary_repr(i)
    m = np.fromstring(s, 'u1') - ord('0')
    m = np.flipud(m)
    return m

def getSecret(stegoMatrix,pos,skip,x=None):
    file = open("out.txt", "w")
    index=0
    if(x==None):
        secret=int(math.sqrt(len(stegoMatrix)*len(stegoMatrix)/(8*skip)))
    else:
        secret=x
    dummy=""
    multiple=False
    if(skip<1):
        skip=1
        multiple=True
    secretMatrix=np.zeros((secret,secret),dtype=np.complex_)
    for a in range(0,len(stegoMatrix)):
        for b in range(0,len(stegoMatrix)):
            c=np.binary_repr(int(stegoMatrix[a][b]),8)
            if(index%skip==0):
                dummy+=c[7-hash(a,b,pos)]
                if(multiple):
                    dummy += c[4+hash(a, b, pos)]
            file.write(c+"\t"+"\n")
            index+=1

    sindex=0
    for a in range(0,min(len(dummy),secret*secret*8),8):
        secretMatrix[int(sindex/secret)][int(sindex%secret)]=int(dummy[a:a+8],2)
        sindex+=1
    return secretMatrix

def all_bit_strings(bits, dtype='uint8'):
    """
    Create a matrix of all binary strings of a given width as the rows.

    Parameters
    ----------
    bits : int
        The number of bits to count through.

    dtype : str or dtype object
        The dtype of the returned array.

    Returns
    -------
    bit_strings : ndarray, shape (2 ** bits, bits)
        The numbers from 0 to 2 ** bits - 1 as binary numbers, most
        significant bit first.

    Notes
    -----
    Obviously the memory requirements of this are exponential in the first
    argument, so use with caution.
    """
    return np.array([map(int, np.binary_repr(i, width=bits))
                     for i in xrange(0, 2 ** bits)], dtype=dtype)

  def _getAllowedShapes(self, shape):
    ''' Return set of allowed shapes that can be squeezed into given shape.

        Examples
        --------
        >>> PB = ParamBag() # fixing K,D doesn't matter
        >>> PB._getAllowedShapes(())
        set([()])
        >>> PB._getAllowedShapes((1))
        set([(), (1,)])
        >>> PB._getAllowedShapes((23))
        set([(23)])
        >>> PB._getAllowedShapes((3,1))
        set([(3), (3,1)])
        >>> PB._getAllowedShapes((1,1))
        set([(), (1,), (1,1)])
    '''
    allowedShapes = set()
    if len(shape) == 0:
      allowedShapes.add(tuple())
      return allowedShapes
    shapeVec = np.asarray(shape, dtype=np.int32)
    onesMask = shapeVec == 1
    keepMask = np.logical_not(onesMask)
    nOnes = sum(onesMask)
    for b in range(2**nOnes):
      bStr = np.binary_repr(b)
      bStr = '0'*(nOnes - len(bStr)) + bStr
      keepMask[onesMask] = np.asarray([int(x) > 0 for x in bStr])
      curShape = shapeVec[keepMask]
      allowedShapes.add(tuple(curShape))
    return allowedShapes

def usm(seqs_list):
    d = {}
    unique = list(set(''.join(seqs_list)))
    unique.sort()
    number_of_bits=int(math.ceil(np.log2(len(unique)))) 
    for number, char in enumerate(unique): d[char] = number
    mat=[]
    for i,seq in enumerate(seqs_list):
        list_b = [np.binary_repr(numb, width=number_of_bits) for numb in [d[char] for char in seqs_list[i]]]
        matrix_usmc = np.zeros([len(seq)+2,number_of_bits*2])
        #Forward Coordinates
        matrix_usmc[0,:number_of_bits] = matrix_usmc[len(seq)+1,number_of_bits:]= np.random.rand(number_of_bits)
        for j in range(1,len(seq)+1):
            binary = np.array([int(x) for x in str(list_b[j-1])])
            matrix_usmc[j,:number_of_bits] = matrix_usmc[j-1,:number_of_bits]+(0.5*(1-matrix_usmc[j-1,:number_of_bits]))*binary -(0.5*matrix_usmc[j-1,:number_of_bits])*(1-binary)
        #Backward Coordinates
        for j in reversed(range(1,len(seq)+1)):
            binary = np.array([int(x) for x in str(list_b[j-1])])
            matrix_usmc[j,number_of_bits:] = matrix_usmc[j+1,number_of_bits:]+(0.5*(1-matrix_usmc[j+1,number_of_bits:]))*binary - (0.5*matrix_usmc[j+1,number_of_bits:])*(1-binary)
        matrix_usmc = matrix_usmc[1:-1,:]
        mat.append(matrix_usmc)
    mat= np.array(mat)
    matrix = np.zeros([len(seqs_list), len(seqs_list)])
    for i, j in itertools.combinations(range(0,len(seqs_list)),2):
          matrix[i][j] = matrix[j][i] = calc_usm_d(mat[i][:,:number_of_bits],mat[i][:,number_of_bits:],mat[j][:,:number_of_bits],mat[j][:,number_of_bits:])
    return matrix