def MSE(sampleSize):
    totalSE = 0.0
    for i in range(sampleSize):
        x = random() * 6.0
        y = random() * 6.0
        error = targetFunction(x,y) - f(x,y)
        totalSE = totalSE + error * error
    print 'The estimated MSE: ', (totalSE / sampleSize)

def MSE(sampleSize):
    totalSE = 0.0
    for i in range(sampleSize):
        in1 = random() * 6.0
        in2 = random() * 6.0
        error = targetFunction(in1,in2) - f(in1,in2)
        totalSE = totalSE + error * error
    print('The estimated MSE: ', (totalSE / sampleSize))

def test_com_jacobian(dq_norm=1e-3, q=None):
    if q is None:
        q = hrp.dof_llim + random(56) * (hrp.dof_ulim - hrp.dof_llim)
    dq = random(56) * dq_norm
    com = hrp.compute_com(q)
    J_com = hrp.compute_com_jacobian(q)
    expected = com + dot(J_com, dq)
    actual = hrp.compute_com(q + dq)
    assert norm(actual - expected) < 2 * dq_norm ** 2
    return J_com

def check_on_random_instance():
    """

    Check our criterion with an LP solver on a random instance of the problem.
    Returns the random point, the criterion's outcome on this point (true or
    false) and whether an LP solution was found (positive or negative). If the
    criterion is correct, true == positive and false == negative.

    """

    K1, K2, K3, C1, C2 = 2 * pylab.random(5) - 1
    px, py = X / (X + Y), Y / (X + Y)

    D4max = .5 * min(2, 1 + C1, 1 + C2, 2 + C1 + C2)
    D4min = .5 * max(-1 + C1, C1 + C2, -1 + C2, 0)
    D4 = D4min + (D4max - D4min) * pylab.random()
    D1, D2, D3 = .5 * (1 + C1) - D4, -.5 * (C1 + C2) + D4, .5 * (1 + C2) - D4

    c = cvxopt.matrix(pylab.array([[1.]] * 8))   # score vector
    G = cvxopt.matrix(pylab.array([
        [+1, 0., 0., 0., 0., 0., 0., 0.],
        [-1, 0., 0., 0., 0., 0., 0., 0.],
        [0., +1, 0., 0., 0., 0., 0., 0.],
        [0., -1, 0., 0., 0., 0., 0., 0.],
        [0., 0., +1, 0., 0., 0., 0., 0.],
        [0., 0., -1, 0., 0., 0., 0., 0.],
        [0., 0., 0., +1, 0., 0., 0., 0.],
        [0., 0., 0., -1, 0., 0., 0., 0.],
        [0., 0., 0., 0., +1, 0., 0., 0.],
        [0., 0., 0., 0., -1, 0., 0., 0.],
        [0., 0., 0., 0., 0., +1, 0., 0.],
        [0., 0., 0., 0., 0., -1, 0., 0.],
        [0., 0., 0., 0., 0., 0., +1, 0.],
        [0., 0., 0., 0., 0., 0., -1, 0.],
        [0., 0., 0., 0., 0., 0., 0., +1],
        [0., 0., 0., 0., 0., 0., 0., -1]]))
    h = cvxopt.matrix(pylab.array([[1.]] * 16))  # h - G x >= 0
    A = cvxopt.matrix(pylab.array([
        [D1, D2, D3, D4, 0, 0, 0, 0],
        [0, 0, 0, 0, D1, D2, D3, D4],
        [-py * D1, +py * D2, +py * D3, -py * D4,
         +px * D1, +px * D2, -px * D3, -px * D4]]))
    b = cvxopt.matrix(pylab.array([K1, K2, K3]))
    sol = cvxopt.solvers.lp(c, G, h, A, b)

    K3min = -1 + py * abs(K1 - C1) + px * abs(K2 - C2)
    K3max = 1 - py * abs(K1 + C1) - px * abs(K2 + C2)

    is_true = K3min <= K3 <= K3max
    is_positive = sol['x'] is not None
    return is_true, is_positive, (K1, K2, K3, C1, C2)

def main():
	shifts = [
		[-1,  1], [0,  1], [1,  1],
		[-1,  0],          [1,  0],
		[-1, -1], [0, -1], [1, -1]
	]

	num_atoms = 100
	num_dims = 2 # dimensions
	coords = pl.random((num_atoms, num_dims))
	chosen = pl.random_integers(num_atoms) # from 1 to num_atoms
	chosen -= 1 # from 0 to num_atoms - 1

	for i in range(len(shifts)):
		coords = pl.vstack((coords, coords[:num_atoms] + shifts[i]))
	num_atoms *= 9 # after 8 shifts added

	max_distance = 0.9
	for i in range(num_atoms):
		if i != chosen:
			dx = coords[chosen, 0] - coords[i, 0]
			dy = coords[chosen, 1] - coords[i, 1]
			distance = pl.sqrt(dx*dx + dy*dy)
			if distance < max_distance:
				pl.plot([coords[i, 0]], [coords[i, 1]], "bo")
			else:
				pl.plot([coords[i, 0]], [coords[i, 1]], "ko")

	# plot last for visibility
	pl.plot([coords[chosen, 0]], [coords[chosen, 1]], "ro")
	pl.grid(True)
	pl.show()

 def get_reading (self):
     reading = []
     t = self.target
     for s in self.sensors:
         r = sqrt ((t[0] - s[0])**2 + (t[1] - s[1])**2) + (random () - 0.5)
         reading.append ([s[0], s[1], r])
     return reading

def main():
	num_atoms = 64
	num_dims = 2 # dimensions
	coords = pl.random((num_atoms, num_dims))

	axis_limits = [0.0, 1.0]
	points = plot_atoms(coords, num_atoms, axis_limits)

	update_limit = 16
	update_count = 0
	while True:
		chosen = pl.random_integers(num_atoms) - 1 # [0, num_atoms - 1]
		new_x, new_y = new_xy(coords, chosen, axis_limits)

		energy_old, energy_new =  energies(coords, chosen, num_atoms, new_x, new_y)

		if energy_new < energy_old:
			coords[chosen, 0] = new_x
			coords[chosen, 1] = new_y
			points[chosen].set_data([new_x, new_y])
			update_count += 1

		if not update_count < update_limit:
			pl.draw()
			update_count = 0

	"""

    def beeswarm(self, data, position, ratio=2.):
        r"""Naive plotting of the data points

        We assume gaussian distribution so we expect fewers dots
        far from the mean/median. We'd like those dots to be close to the
        axes. conversely, we expect lots of dots centered around the mean, in
        which case, we'd like them to be spread in the box. We uniformly
        distribute position using

        .. math::

            X = X + \dfrac{ U()-0.5 }{ratio} \times factor

        but the factor is based on an arctan function:

        .. math::

            factor = 1 - \arctan( \dfrac{X - \mu }{\pi/2})

        The farther the data is from the mean :math:`\mu`,
        the closest it is to the axes that goes through the box.

        """
        N = len(data)
        m = np.median(data)
        sd = np.std(data)
        # arctan function to have a tapering window
        factor = 1. - np.abs(np.arctan((data-m)/sd)/1.570796)  # pi/2

        newdata = position + (pylab.random(N) - 0.5)/float(ratio) * factor
        return newdata

    def __init__(self, Fs= 16000, TinSec= 10):
        '''
        Fs: 取樣頻率，預設值為 16000，
        TinSec: 保存語音長度，預設值為 10 sec
        '''
        print('RyAudio use %s'%pa.get_portaudio_version_text())
        
        self.Fs= Fs
        self.spBufferSize=           1024
        self.fftWindowSize= self.spBufferSize

        self.aP= pa.PyAudio()
        self.iS= pa.Stream(PA_manager= self.aP, input= True, rate= self.Fs, channels= 1, format= pa.paInt16)
        self.oS= pa.Stream(PA_manager= self.aP, output= True, rate= self.Fs, channels= 1, format= pa.paInt16)
        self.iTh= None
        self.oTh= None

        #self.sound=     None
        #self.soundTime= 0
        self.gettingSound= True
        self.playingSound= True

        self.t= 0
        self.b= None  # byte string
        self.x= None  # ndarray
        self.fft= None
        self.f0= 0#None
        self.en= 0#None
        self.fm= 0#None # frequency mean
        self.fv= 0#None # frequency var
        self.fs= 0#None # frequency std
        self.enP= 0#None # AllPass
        self.enPL= 0#None # LowPass
        self.enPH= 0#None # HighPass

        self.entropy= 0#None

        self.frameI=   0

        #self.frameN= self.spBufferSize/4  #1024/4 = 256
        self.TinSec= TinSec #10 # sec
        self.frameN= self.Fs*self.TinSec/self.spBufferSize #self.spBufferSize/4  #1024/4 = 256
        self.frameN= int(self.frameN)

        self.specgram= pl.random([self.frameN, self.spBufferSize/2])

        self.xBuf= pl.random([self.frameN, self.spBufferSize])

 def mc_counts2samples( self, counts ) :
     values=[]
     bin=-180
     for value in counts:
         for i in arange(value) :
             values.append(bin+random()*self.step)
         bin +=  self.step
     return array(values)

def sqr_cplx(z):
    """
    Fonction racine d'un complexe. Prend aléatoirement la racine
    positive ou négative
    """
    r, theta = cart2pol(z)
    r = pl.sqrt(r)
    theta = theta/2 + int(2*pl.random())*pl.pi
    return pol2cart(r, theta)

 def next(self):
     Tnext = ((self.Konstant * self.t1) * 2) - self.t0
     if len(self.values) % 100 > 70:
         self.values.append(pylab.random() * 2 - 1)
     else:
         self.values.append(Tnext)
     self.t0 = self.t1
     self.t1 = Tnext
     return self.values[-1]

 def markov_trajectory(self,distribution,state=None):
     if state == None :
         state = self.start_point(distribution)
     trj = []
     for i in range(self.frames):
         state = self.markov_step(state,distribution)
         angle = self.mod2pi((state+random())*self.step-180)
         trj.append(angle)
     return array(trj)

def random_euler_angles():
    r1,r2,r3 = pylab.random(3)
    q1 = pylab.sqrt(1.0-r1)*pylab.sin(2.0*pylab.pi*r2)
    q2 = pylab.sqrt(1.0-r1)*pylab.cos(2.0*pylab.pi*r2)
    q3 = pylab.sqrt(r1)*pylab.sin(2.0*pylab.pi*r3)
    q4 = pylab.sqrt(r1)*pylab.cos(2.0*pylab.pi*r3)
    phi = math.atan2(2.0*(q1*q2+q3*q4), 1.0-2.0*(q2**2+q3**2))
    theta = math.asin(2.0*(q1*q3-q4*q2))
    psi = math.atan2(2.0*(q1*q4+q2*q3), 1.0-2.0*(q3**2+q4**2))
    return [phi,theta,psi]

 def initialize(self):
     self.state = pylab.zeros([self.n, self.n])
     for x in xrange(self.n):
         for y in xrange(self.n):
             self.state[x, y] = 1 if pylab.random() < self.init_p else 0
     self.state[self.n/2, self.n/2] = 2
     self.next_state = pylab.zeros([self.n, self.n])
     for x in xrange(self.n):
         for y in xrange(self.n):
         	if self.state[x, y] == 1 or self.state[x, y] == 2:
         		self.total_num_trees += 1

def GenerateRandomTrajectory(ncurve, ndof, bound):
    def vector2string(v):
        ndof = len(v)
        s = str(ndof)
        for a in v:
            s += ' %f' % a
        return s

    p0a = vector2string(random(ndof) * 2 * bound - bound)
    p0b = vector2string(random(ndof) * 2 * bound - bound)
    p1a = vector2string(random(ndof) * 2 * bound - bound)
    p1b = vector2string(random(ndof) * 2 * bound - bound)
    s = '%d' % ncurve
    s += '\n1.0 ' + p0a + ' ' + p0b
    for k in range(ncurve - 1):
        a = random(ndof) * 2 * bound - bound
        b = random(ndof) * 2 * bound - bound
        c = 2 * b - a
        pa = vector2string(a)
        pb = vector2string(b)
        pc = vector2string(c)
        s += ' ' + pa + ' ' + pb + '\n1.0 ' + pb + ' ' + pc
    s += ' ' + p1a + ' ' + p1b
    Tv, p0v, p1v, p2v, p3v = string2p(s)
    return BezierToTrajectoryString(Tv, p0v, p1v, p2v, p3v)

    def __init__(self,n=100,state=0.5,J=1.,T=2.,):
        """Constructor..."""
        self.n = n
        self.spins = -1 * pl.ones((n,n))
        self.J = J
        self.T=T
        self.init = state
        self.config = "n=%d,init=%f,J=%f,T=%f" % (self.n,self.init,self.J,self.T)

        for i in xrange(n):
            for j in xrange(n):
                if pl.random() < (1 + state) / 2.:
                    self.spins[i,j] = 1

def packing(dataset, radius=4, nifti=False, randoffset=False):
	"""return a hexagonal close sphere packing grid for a PyMVPA fMRI dataset
    
    Keyword arguments:
    radius-- radius in voxels of the spheres to pack (default 4)
    nifti-- write out a seed voxel mask as a nifti
    randomoffset-- random jitter of the seed voxel grid

    """
    
	from pylab import find, random
	from numpy import ones, zeros, arange, sqrt, remainder
	from mvpa2.suite import fmri_dataset, Dataset
	import os

	if randoffset:
		ro = random(3)
	else:
		ro = zeros(3)
	
	minco = dataset.fa.voxel_indices.min(0)
	maxco = dataset.fa.voxel_indices.max(0)
	rect = ones(dataset.a.voxel_dim)
	
	fac = sqrt(6)*2*radius/3
	for iz,z in enumerate(arange(minco[2], maxco[2], fac)):
		for iy,y in enumerate(arange(minco[1], maxco[1], fac)):
			for x in arange(minco[0], maxco[0], 2*radius):
				hx = x + remainder(iy, 2)*radius + ro[0]*radius
				hy = y + remainder(iz, 2)*fac + ro[1]*radius
				hz = z + ro[2]*radius
				if hz <= maxco[2]:
					rect [hx, hy, hz] += 1

	maskedrect = dataset.mapper.forward1(rect)
	roiIndex = find((maskedrect == 2))
	print  'number of seed voxel: '+str(len(roiIndex))
	
	maskedrectds = Dataset([maskedrect])
	maskedrectds.a = dataset.a.copy()
	maskedrectds.fa = dataset.fa.copy()
    
	if nifti:
		from nibabel import Nifti1Image
		Nifti1Image(maskedrectds.O.squeeze(),
					dataset.a.imghdr.get_best_affine()
					).to_filename(os.path.join('sparse'+str(int(radius))+'.nii.gz'))

	return roiIndex, maskedrectds

    def step(self):
        """Run one step of the Metropolis algorithm"""
        site = self.getsite()
        
        #Now need to work out the energy difference between the lattices.
        #This is used to determine the probability that we'll keep the
        #new configuration.
        Ediff = -2 * self.Hij(site)
        Sdiff = -2 * self.spins[site]

        if self.probaccept(Ediff) > pl.random():
            self.spinflip(site)
            return (Ediff, Sdiff)
        else:
            return (0.,0.)

def generate_waveforms(
    N_harmonics=[8,16,32,64],
    path='/home/ritz/mix/audio samples instruments/basic-waveforms/',
    name='square_rnd_phase-%03dh-G2-(i).wav'):
    '''Generate a series of basic additive waveforms with a varying number of harmonics.'''
    f0 = 49.0
    w0 = pl.round(sr / f0)
    f0 = sr / w0
    pl.clf()
    for ii, N in enumerate(N_harmonics):
        H = 1.0 + 2 * pl.arange(N)
        waev = beep(pl.c_[f0 * H, pl.random(N) * tau, 1 / H], 8 * w0)
        pl.subplot(3,4,ii+1)
        pl.plot(waev)
        write(waev, path + (name % N))