def velovect(u1,u2,d,minvel=1e-40,nvect=None,scalevar=None,scale=100,color='k',fig=None):
    '''Plots normalized velocity vectors'''


    if fig==None:
        ax=plt.gca()
    else:
        ax=fig.ax

    CC=d.getCenterPoints()
    n=np.sqrt(u1**2+u2**2)
    # remove zero velocity:
    m=n<minvel
    vr=np.ma.filled(np.ma.masked_array(u1/n,m),0.)
    vz=np.ma.filled(np.ma.masked_array(u2/n,m),0.)
    if scalevar != None:
        vr = vr*scalevar
        vz = vz*scalevar
    if nvect==None:
        Q=ax.quiver(CC[:,0],CC[:,1],vr,vz,pivot='middle',width=1e-3,minlength=0.,scale=scale,
                    headwidth=6)
    else:
        # regrid the data:
        tmp0=np.complex(0,nvect[0])
        tmp1=np.complex(0,nvect[1])
        grid_r, grid_z = np.mgrid[ax.get_xlim()[0]:ax.get_xlim()[1]:tmp0, ax.get_ylim()[0]:ax.get_ylim()[1]:tmp1]
        grid_vr = griddata(CC, vr, (grid_r, grid_z), method='nearest')
        grid_vz = griddata(CC, vz, (grid_r, grid_z), method='nearest')
        Q=ax.quiver(grid_r,grid_z,grid_vr,grid_vz,pivot='middle',width=2e-3,minlength=minvel,scale=scale,
                    headwidth=10,headlength=10,color=color,edgecolor=color,rasterized=True)

    plt.draw()
    return Q     

    def __init__(self, array, epsilon=-1, resample=-1):
        try:
            self.points = np.copy(array)
            if epsilon > 0:
                self.points = cv2.approxPolyDP(self.points, epsilon, True)
            if resample > 0:
                pass
            self.moments = cv2.moments(self.points, False)
            zeroth = self.moments["m00"]
            self.x = self.moments["m10"] / zeroth
            self.y = self.moments["m01"] / zeroth
            self.area = cv2.contourArea(self.points)
            self.perim = cv2.arcLength(self.points,True)
            self.hull = cv2.convexHull(self.points)
            self.area_hull = cv2.contourArea(self.hull)
            self.perim_hull = cv2.arcLength(self.hull,True)
            rect = cv2.minAreaRect(self.points)
            box = cv.BoxPoints(rect)
            a = np.complex(box[0][0], box[0][1])
            b = np.complex(box[1][0], box[1][1])
            c = np.complex(box[2][0], box[2][1])
            self.w, self.h = sorted([np.abs(a-b), np.abs(c-b)])
            self.is_valid = True
        except ZeroDivisionError:
            self.area = 0
            self.perim = 0
            self.hull = 0
            self.area_hull = 0
            self.perim_hull = 0

            self.w, self.h = 0,0
            self.is_valid = False

 def get_transfer_matrix_mine(self, index, energy):
     if energy == self.data[index + 1].height:
         km_ratio = np.sqrt(
             self.data[index].mass
             * np.complex(energy - self.data[index].height)
             / self.data[index + 1].mass
             / FAKE_ZERO
         )
     else:
         km_ratio = np.sqrt(
             self.data[index].mass
             * np.complex(energy - self.data[index].height)
             / self.data[index + 1].mass
             / np.complex(energy - self.data[index + 1].height)
         )
     kcp = 1 + km_ratio
     kcm = 1 - km_ratio
     d_this = self.get_delta(index, energy)
     d_next = self.get_delta(index + 1, energy)
     kd = self.get_wave_number(index, energy) * self.data[index].width
     # print "d_this =", d_this, "d_next =", d_next, "kd", kd, "1j * ( d_this + d_next - kd) =", 1j * ( d_this + d_next - kd)
     return np.matrix(
         [
             [kcp / 2 * np.exp(1j * (d_this - d_next - kd)), kcm / 2 * np.exp(1j * (d_this + d_next - kd))],
             [kcp / 2 * np.exp(-1j * (d_this + d_next - kd)), kcm / 2 * np.exp(-1j * (d_this - d_next - kd))],
         ]
     )

def uniform_seed_grid():

    #read bvals,gradients and data   
    fimg,fbvals, fbvecs = get_data('small_64D')    
    bvals=np.load(fbvals)
    gradients=np.load(fbvecs)
    img =ni.load(fimg)    
    data=img.get_data()
    
    x,y,z,g=data.shape   

    M=np.mgrid[.5:x-.5:np.complex(0,x),.5:y-.5:np.complex(0,y),.5:z-.5:np.complex(0,z)]
    M=M.reshape(3,x*y*z).T

    print(M.shape)
    print(M.dtype)

    for m in M: 
        print(m)
    gqs = GeneralizedQSampling(data,bvals,gradients)
    iT=iter(EuDX(gqs.QA,gqs.IN,seeds=M))    
    T=[]
    for t in iT:
        T.append(i)
    
    print('lenT',len(T))
    assert_equal(len(T), 1221)

def init(f, g):
    k = 0
    while k < len(f):
        f[k] = np.complex(k, k + 1)
        g[k] = np.complex(k, 2 * k + 1)
        k += 1
    return

    def pherr(self, i, phi):
        """add phase error to antenna i
        """

        for j in xrange(self.na):
            self.data[i,j] = self.data[i,j] * n.exp(n.complex(0.,n.sin(phi)))
            self.data[j,i] = self.data[j,i] * n.exp(n.complex(0.,n.sin(-phi)))

    def general_work(self, input_items, output_items):
        nread = self.nitems_read(0)  # number of items read on port 0
        ninput_items = len(input_items[0])
        noutput_items = len(output_items[0])
        nitems_to_consume = min(ninput_items, noutput_items)
        in0 = input_items[0]
        out0 = output_items[0]

        for ii in range(nitems_to_consume):
            x = in0[ii].real
            y = in0[ii].imag

            if x == 0 and y == 0:
                out0[ii] = numpy.complex(0)
            else:
                r = numpy.sqrt(numpy.square(x) + numpy.square(y))
                theta = numpy.arctan2(x, y) - self.angle

                a = r * numpy.cos(theta)
                b = r * numpy.sin(theta)

                out0[ii] = numpy.complex(a, b)

        self.consume(0, nitems_to_consume)
        return nitems_to_consume

def sdb2sri(Sdb, Sdeg):
	# convert DB/DEG to real/imag
	num_freqs = len(Sdb)
	Sri = np.zeros( (num_freqs, 2, 2), dtype=complex)

	for idx in range(len(Sdb)):
		db_mat = Sdb[idx]
		S11_db = db_mat[0][0]
		S12_db = db_mat[0][1]
		S21_db = db_mat[1][0]
		S22_db = db_mat[1][1]

		deg_mat = Sdeg[idx]
		S11_deg = deg_mat[0][0]
		S12_deg = deg_mat[0][1]
		S21_deg = deg_mat[1][0]
		S22_deg = deg_mat[1][1]

		S11 = 10**(S11_db/20) * np.complex( np.cos(S11_deg*np.pi/180), np.sin(S11_deg*np.pi/180) )
		S12 = 10**(S12_db/20) * np.complex( np.cos(S12_deg*np.pi/180), np.sin(S12_deg*np.pi/180) )
		S21 = 10**(S21_db/20) * np.complex( np.cos(S21_deg*np.pi/180), np.sin(S21_deg*np.pi/180) )
		S22 = 10**(S22_db/20) * np.complex( np.cos(S22_deg*np.pi/180), np.sin(S22_deg*np.pi/180) )

		Sri[idx][0][0] = S11
		Sri[idx][0][1] = S12
		Sri[idx][1][0] = S21
		Sri[idx][1][1] = S22

	return Sri

def getGroundState(psi0,V,a,b,dt,omega=1.0,delta=0.0,epsilon=0.048,phi=4.0/3.0,*args,**kwargs):
    dx=(b-a)/(psi0.shape[0])
    t0=time.clock()
    eigE,eigV,eigVdagger=getEigenHam2(a,b,psi0.shape[0],V,omega=omega,delta=delta,epsilon=epsilon,phi=phi,*args,**kwargs)
    t1=time.clock()
    print 'Got EigenHam in '+str(t1-t0)+' seconds!'
    psi0=psi0/np.sqrt(np.sum(psi0*np.conj(psi0)*dx))
    psi0eigB=np.einsum('ijk,ik->ij',eigVdagger,psi0)
    t0=time.clock()
    psi1eigB=splitStepPropagatorEigB2(psi0eigB,dt*np.complex(0.0,-1.0),a,b,eigE,eigV,eigVdagger)
    t1=time.clock()
    psi1eigB=psi1eigB/np.sqrt(np.sum(psi1eigB*np.conj(psi1eigB)*dx))
    t2=time.clock()
    print 'Completed one time step in '+str(t1-t0)+' seconds!'
    print 'Then normalized wavefunction in '+str(t2-t1)+' seconds!'
    diff=np.sum(np.abs(psi0eigB*np.conj(psi0eigB)*dx-psi1eigB*np.conj(psi1eigB)*dx))
    i=0
    while diff>0.1e-5:
        psi0eigB=psi1eigB
        psi1eigB=splitStepPropagatorEigB2(psi0eigB,dt*np.complex(0.0,-1.0),a,b,eigE,eigV,eigVdagger)
        psi1eigB=psi1eigB/np.sqrt(np.sum(psi1eigB*np.conj(psi1eigB)*dx))
        diffLast=diff
        diff=np.sum(np.abs(psi0eigB*np.conj(psi0eigB)*dx-psi1eigB*np.conj(psi1eigB)*dx))
        if diffLast<diff:
            print 'Not converging! Difference went up from %f to %f' %(diffLast,diff)
            break
        i+=1
    print i
    print diff
    psi1=np.einsum('ijk,ik->ij',eigV,psi1eigB)
    psi1=psi1.transpose()
    return psi1

def test_disbesldv():
    qxqyv = besselaesnew.disbesldv(2.0, 1.0, np.complex(-3.0, -1.0), np.complex(2.0, 2.0),
              [0.0, 2.0, 11.0], 1, 1, 3)
    assert_allclose(qxqyv[0], np.array([-0.17013114606375021, -0.18423853257632447, -0.17315784943727297]))
    assert_allclose(qxqyv[2], np.array([2.7440507429637e-002, 8.880686745447e-002, 3.426560831291e-002]))
    assert_allclose(qxqyv[1], np.array([-0.10412493484448178, -0.10844664064434061, -0.10447761803194042]))
    assert_allclose(qxqyv[3], np.array([0.10617613097471285, 0.11627387807684744, 0.10674211206906066]))

    def plot_error_map(self, point_error, ele_pos):
        """
        Creates plot of mean error calculated separately for every point of
        estimation space

        Parameters
        ----------
        point_error: numpy array
            Error of reconstruction calculated at every point of reconstruction
            space.
        ele_pos: numpy array
            Positions of electrodes.

        Returns
        -------
        mean_error: numpy array
            Accuracy mask.
        """
        ele_x, ele_y = ele_pos[:, 0], ele_pos[:, 1]
        x, y = np.mgrid[self.kcsd_xlims[0]:self.kcsd_xlims[1]:
                        np.complex(0, point_error.shape[1]),
                        self.kcsd_ylims[0]:self.kcsd_ylims[1]:
                        np.complex(0, point_error.shape[2])]
        mean_error = self.sigmoid_mean(point_error)
        plt.figure(figsize=(12, 7))
        ax1 = plt.subplot(111, aspect='equal')
        levels = np.linspace(0, 1., 25)
        im = ax1.contourf(x, y, mean_error, cmap='Greys')
        plt.colorbar(im)#im, fraction=0.046, pad=0.06)
        plt.scatter(ele_x, ele_y, 10, c='k')
        ax1.set_xlabel('Depth x [mm]')
        ax1.set_ylabel('Depth y [mm]')
        ax1.set_title('Sigmoidal mean point error')
        plt.show()
        return mean_error

def hermitian_subspace_basis(n):
    """
    returns a basis set for the real subspace of Hermitian n*n matrices
    """
    sqrt2over2 = numpy.sqrt(2.)/2.
    assert(n)
    basis = []
    for row in range(0,n):
        for col in range(row,n):
            if row == col:
                mat = numpy.zeros((n,n), dtype=complex)
                mat[row][col] = 1.
                basis.append(mat)
            else:
                mat = numpy.zeros((n,n), dtype=complex)
                mat[col][row] = sqrt2over2
                mat[row][col] = sqrt2over2
                basis.append(mat)                
                mat = numpy.zeros((n,n), dtype=complex)
                mat[row][col] = numpy.complex(0.,sqrt2over2)
                mat[col][row] = numpy.complex(0.,-sqrt2over2)
                basis.append(mat)

    # unit vectors ?
#    for M in basis:
#        for F in basis:
#           print hs.dot(F,M)
        # assert hs.dot(M,M) == 1.
              
    return basis            

def make_invariant(points):
	print '*' * 50
	print 'make_invariant called'
	print '*' * 50
	print

	N = len(points)

	points = map(lambda x : np.complex(x[0], x[1]), points)
	FD = np.fft.fft(points)
	# translational invariance
	translational_invar = FD[0]
	FD = map(lambda x : x - translational_invar, FD)
	# scalar invariance 
	val = ( ( FD[1] )*( FD[1].conjugate() ) ).real
	FD = map(lambda x : x / val, FD)
	# rotational invariance
	final_points = []
	for i in range(N):
		prev_ = FD[(i-1)%N]
		next_ = FD[(i+1)%N]
		angle = ( 2 * math.pi * i )/N
		cmplx1 = np.complex(math.cos(angle), math.sin(angle))
		cmplx2 = cmplx1.conjugate()
		u_i = prev_*cmplx2 + next_*cmplx1
		final_points.append(u_i)

	phi = math.atan( final_points[1].imag / final_points[1].real )
	cmplx3 = np.complex(math.cos(phi), math.sin(phi))

	final_points = map( lambda x : [ (x*cmplx3).real, (x*cmplx3).imag], final_points)	

	print final_points
	return final_points

 def __init__(self,s1pPath,delay=0.0):
     assert s1pPath.endswith('.s1p')
     self.title = re.split('/',s1pPath)[-1]
     self.frequencies = [] # in Hz
     self.s11Values = [] # as complex, linear values
     
     touchstoneFileHandle = open(s1pPath,'r',1)
     for line in touchstoneFileHandle:
         if line[0] == '!':
             pass # comment, ignore for the moment
         elif line[0] == '#':
             assert line.endswith('hz S db R 50\n')
         else: # should be a data line
             splittedLine = re.split('\s+',line)
             floatItems = []
             for item in splittedLine:
                 try:
                     floatItems.append(float(item))
                 except ValueError:
                     pass
             
             assert len(floatItems) == 3 # we only support S1P files for the moment
             self.frequencies.append(floatItems[0])
             
             s11angle = numpy.exp(numpy.complex(0,numpy.deg2rad(floatItems[2])))
             s11magnitude = 10.0**(floatItems[1]/20.0)
             
             self.s11Values.append(s11magnitude*s11angle)
     
     self.frequencies = numpy.array(self.frequencies)
     self.s11Values = numpy.array(self.s11Values)
     # compensate electrical delay
     self.s11Values = numpy.exp(numpy.complex(0,1)*2*2*numpy.pi*self.frequencies*delay) * self.s11Values

def load_data(filename,y1_col,y2_col,sformat='realimag',phase_conversion = 1,ampformat='lin',fdata_unit=1.,delimiter=None):
    '''
    sformat = 'realimag' or 'ampphase'
    ampformat = 'lin' or 'log'
    '''
    f = open(filename)
    lines = f.readlines()
    f.close()
    z_data = []
    f_data = []

    if sformat=='realimag':
        for line in lines:
            if ((line!="\n") and (line[0]!="#") and (line[0]!="!")) :
                lineinfo = line.split(delimiter)
                f_data.append(float(lineinfo[0])*fdata_unit)
                z_data.append(np.complex(float(lineinfo[y1_col]),float(lineinfo[y2_col])))
    elif sformat=='ampphase' and ampformat=='lin':
        for line in lines:
            if ((line!="\n") and (line[0]!="#") and (line[0]!="!") and (line[0]!="M") and (line[0]!="P")):
                lineinfo = line.split(delimiter)
                f_data.append(float(lineinfo[0])*fdata_unit)
                z_data.append(float(lineinfo[y1_col])*np.exp( np.complex(0.,phase_conversion*float(lineinfo[y2_col]))))
    elif sformat=='ampphase' and ampformat=='log':
        for line in lines:
            if ((line!="\n") and (line[0]!="#") and (line[0]!="!") and (line[0]!="M") and (line[0]!="P")):
                lineinfo = line.split(delimiter)
                f_data.append(float(lineinfo[0])*fdata_unit)
                linamp = 10**(float(lineinfo[y1_col])/20.)
                z_data.append(linamp*np.exp( np.complex(0.,phase_conversion*float(lineinfo[y2_col]))))
    else:
        print "ERROR"
    return np.array(f_data), np.array(z_data)

def IntegralLapLineDipoleDis(zin, z1, z2, del0, ra, order):
    cg = np.full(order + 2, np.complex(0.0, 0.0))

    z = (2.0 * zin - (z1 + z2)) / (z2 - z1)
    zplus1 = z + 1.0
    zmin1 = z - 1.0

    # Determine coefficients of powers of Delta for [ (Delta-Delta_0)/a ] ^ p
    for m in range(0, order + 1):
        cg[m] = RBINOM[order, m] * (-del0) ** (order - m) / ra ** order

    zterm1 = np.complex(0.0, 0.0)
    zterm2 = np.complex(0.0, 0.0)
    for n in range(1, order + 1):
        zterm1 = zterm1 + cg[n] * float(n) * z ** (n - 1)
    for n in range(0, order + 1):
        zterm2 = zterm2 + cg[n] * z ** n

    qmtot = np.complex(0.0, 0.0)
    for m in range(2, order + 1):
        qm = np.complex(0.0, 0.0)
        for n in range(1, int(m / 2) + 1):
            qm = qm + float(m - 2 * n + 1) * z ** (m - 2 * n) / float(2 * n - 1)
        qmtot = qmtot + 2.0 * cg[m] * qm

    wdis = (
        zterm1 * np.log(zmin1 / zplus1) + zterm2 * (1.0 / zmin1 - 1.0 / zplus1) + qmtot
    ) / (2.0 * np.pi)
    return wdis

    def check_numpy_scalar_argument_return_generic(self):
        f = PyCFunction('foo')
        f += Variable('a1', numpy.int_, 'in, out')
        f += Variable('a2', numpy.float_, 'in, out')
        f += Variable('a3', numpy.complex_, 'in, out')
        foo = f.build()
        args = 2, 1.2, 1+2j
        results = numpy.int_(2), numpy.float_(1.2), numpy.complex(1+2j)
        assert_equal(foo(*args),results)
        args = [2], [1.2], [1+2j]
        assert_equal(foo(*args),results)
        args = [2], [1.2], [1,2]
        assert_equal(foo(*args),results)

        f = PyCFunction('foo')
        f += Variable('a1', 'npy_int', 'in, out')
        f += Variable('a2', 'npy_float', 'in, out')
        f += Variable('a3', 'npy_complex', 'in, out')
        foo = f.build()
        args = 2, 1.2, 1+2j
        results = numpy.int_(2), numpy.float_(1.2), numpy.complex(1+2j)
        assert_equal(foo(*args),results)
        args = [2], [1.2], [1+2j]
        assert_equal(foo(*args),results)
        args = [2], [1.2], [1,2]
        assert_equal(foo(*args),results)

def getftzer(Jzer,ngrid=128,Rpix=100):
    '''
    ; Compute the Fourier Transform of Zernike mode

    ; ngrid = 128 ; grid half-size, pixels
    ; Rpix = 100 ; pupil radius in pixels

    :param Jzer:
    :return:
    '''

    x = np.arange(-ngrid,ngrid)
    y = np.arange(-ngrid,ngrid)
    theta = np.arctan2(x,y)

    n,m = zern_num(Jzer)
    f = np.roll(np.roll(dist(2*ngrid),
                        ngrid,
                        axis=0),
                ngrid,
                axis=1)/(2*ngrid)*Rpix
    f[ngrid][ngrid] = 1e-3

    ftmod = np.sqrt(n+1.0)*jv(n+1,2*np.pi*f)/(np.pi*f)

    if m == 0:
        zz = ftmod*np.complex(0,1.)**(n/2.)
    else:
        if (Jzer%2 == 0):
            fact=np.sqrt(2.)*np.cos(m*theta)
        else:
            fact=np.sqrt(2.)*np.sin(m*theta)
        zz = ftmod*fact*(-1)**((n-m/2.))*np.complex(0,1.)**m

    return zz

def interpolation2D(inputdata,
                    vsize=200, hsize=400, 
                    method='cubic'):
    """
    use different 2D interpolation methods from scipy
    
    visize and hsize are the output dim after interpolation
    
    methods: nearest, linear, cubic
    
    some bugs exist. This part does't work well.
    
    """
    inputdata = np.array(inputdata)
    datas = inputdata.shape

    grid_x, grid_y = np.mgrid[0:datas[0]:np.complex(0,vsize), 0:datas[1]:np.complex(0,hsize)]

    val = np.zeros([datas[0]*datas[1]])
    pos = np.zeros([datas[0]*datas[1], 2])

    for i in range(datas[0]):
        for j in range(datas[1]):
            pos[i*datas[1]+j,0] = i
            pos[i*datas[1]+j,1] = j
            val[i*datas[1]+j] = inputdata[i,j]

    grid_z = griddata(pos, val, (grid_x, grid_y), method=method)

    return grid_z

 def xdatan(self):
     if self.is_numerical:
         return self.xdata
     else:
         return (numpy.complex(self.center),
                 numpy.complex(self.xcoefficient),
                 numpy.int(self.ramification_index))