def xover(chrom,N,p):
    """Single point crossover with probability N,precision p
    """
    N = round(chrom.shape[0]*N)
    index1 = scipy.arange(chrom.shape[0])
    index2 = scipy.unique(scipy.around(scipy.rand(chrom.shape[0],)*chrom.shape[0]))[0:chrom.shape[0]/2]
    sel1,sel2 = [],[]
    for i in range(len(index1)):
        if index1[i] not in index2:
            sel1.append(index1[i])
        else:
            sel2.append(index1[i])
    select1 = sel1[0:min([int(round(len(sel1)*N)),int(round(len(sel2)*N))])]
    select2 = sel2[0:min([int(round(len(sel1)*N)),int(round(len(sel2)*N))])]       
    
    # set xover points
    xoverpnt = scipy.around(scipy.rand(len(select1),)*(chrom.shape[1]-1))
    
    # perform xover
    nchrom = copy.deepcopy(chrom)
    for i in range(len(select1)):
        try:
            slice1 = chrom[select1[i],0:int(xoverpnt[i])]
            slice2 = chrom[select2[i],0:int(xoverpnt[i])]
            nchrom[select2[i],0:int(xoverpnt[i])] = slice1
            nchrom[select1[i],0:int(xoverpnt[i])] = slice2
        except:
            nchrom = nchrom
    
    return nchrom

 def __init__(self, basef, 
              translate=True, 
              rotate=False, 
              conditioning=None, 
              asymmetry=None,
              oscillate=False, 
              penalize=None,
              ):
     FunctionEnvironment.__init__(self, basef.xdim, basef.xopt)
     self.desiredValue = basef.desiredValue            
     self.toBeMinimized = basef.toBeMinimized
     
     if translate:            
         self.xopt = (rand(self.xdim) - 0.5) * 9.8
         
     self._diags = eye(self.xdim)            
     self._R = eye(self.xdim)            
     self._Q = eye(self.xdim)            
     
     if conditioning is not None:
         self._diags = generateDiags(conditioning, self.xdim)
     if rotate:
         self._R = orth(rand(basef.xdim, basef.xdim))        
         if conditioning:
             self._Q = orth(rand(basef.xdim, basef.xdim))
                     
     tmp = lambda x: dot(self._Q, dot(self._diags, dot(self._R, x-self.xopt)))
     if asymmetry is not None:
         tmp2 = tmp
         tmp = lambda x: asymmetrify(tmp2(x), asymmetry)
     if oscillate:
         tmp3 = tmp
         tmp = lambda x: oscillatify(tmp3(x))
     
     self.f = lambda x: basef.f(tmp(x))

    def test_blackbox(self):
        for A, b in self.cases:
            x = solve(A, b, verb=False, maxiter=A.shape[0])
            assert(norm(b - A*x)/norm(b - A*rand(b.shape[0],)) < 1e-4)

        # Special tests
        # (1) Make sure BSR format is preserved, and B is multiple vecs
        A, b = self.cases[-1]
        (x, ml) = solve(A, b, return_solver=True, verb=False)
        assert(ml.levels[0].B.shape[1] == 3)
        assert(ml.levels[0].A.format == 'bsr')

        # (2) Run with solver and make sure that solution is still good
        x = solve(A, b, existing_solver=ml, verb=False)
        assert(norm(b - A*x)/norm(b - A*rand(b.shape[0],)) < 1e-4)

        # (3) Convert to CSR, make sure B is a single vector
        (x, ml) = solve(A.tocsr(), b, return_solver=True, verb=False)
        assert(ml.levels[0].B.shape[1] == 1)
        assert(ml.levels[0].A.format == 'csr')

        # (4) Run with x0, maxiter and tol
        x = solve(A, b, existing_solver=ml, x0=zeros_like(b), tol=1e-8,
                  maxiter=300, verb=False)
        assert(norm(b - A*x)/norm(b - A*rand(b.shape[0],)) < 1e-7)

        # (5) Run nonsymmetric example, make sure BH isn't None
        A, b = self.cases[2]
        (x, ml) = solve(A, b, return_solver=True, verb=False,
                        maxiter=A.shape[0])
        assert(ml.levels[0].BH is not None)

 def run(self):
  # Para cada particula 
  for i in scipy.arange(self.ns):
   # atualiza melhor posicao da particula
   self.fit[i],aux = self.avalia_aptidao(self.pop[i])
   self.pop[i] = aux.copy()
   # atualiza melhor posicao da particula
   self.bfp_fitness[i],aux = self.avalia_aptidao(self.bfp[i])
   self.bfp[i] = aux.copy()
   if self.debug:
    print "self.fit[{0}] = {1} bfp_fitness = {2}".format(i,self.fit[i],self.bfp_fitness[i])
   if self.bfp_fitness[i] < self.fit[i]:
    self.bfp[i] = self.pop[i].copy()
    self.bfp_fitness[i] = self.fit[i].copy()
  
  # Atualiza melhor posicao global
  idx = self.bfp_fitness.argmax()
  curr_best_global_fitness = self.bfp_fitness[idx]
  curr_best_global = self.bfp[idx].copy()
  if curr_best_global_fitness > self.bfp_fitness.max():
    self.bfg = curr_best_global
 
  for i in scipy.arange(self.ns):
   # Atualiza velocidade
   self.v[i] = self.w*self.v[i] 
   self.v[i] = self.v[i] + self.c1*scipy.rand()*( self.bfp[i] - self.pop[i]) 
   self.v[i] = self.v[i] + self.c2*scipy.rand()*(self.bfg - self.pop[i])
   # Atualiza posicao
   self.pop[i] = self.pop[i] + self.v[i]
   
  # calcula fitness 
   self.fit[i],aux = self.avalia_aptidao(self.pop[i])
   self.pop[i] = aux.copy()

    def test_DAD(self):
        A = poisson((50, 50), format="csr")

        x = rand(A.shape[0])
        b = rand(A.shape[0])

        D = diag_sparse(1.0 / sqrt(10 ** (12 * rand(A.shape[0]) - 6))).tocsr()
        D_inv = diag_sparse(1.0 / D.data)

        DAD = D * A * D

        B = ones((A.shape[0], 1))

        # TODO force 2 level method and check that result is the same
        kwargs = {"max_coarse": 1, "max_levels": 2, "coarse_solver": "splu"}

        sa = smoothed_aggregation_solver(D * A * D, D_inv * B, **kwargs)

        residuals = []
        x_sol = sa.solve(b, x0=x, maxiter=10, tol=1e-12, residuals=residuals)

        avg_convergence_ratio = (residuals[-1] / residuals[0]) ** (1.0 / len(residuals))

        # print "Diagonal Scaling Test:   %1.3e,  %1.3e" %
        # (avg_convergence_ratio, 0.25)
        assert avg_convergence_ratio < 0.25

    def test_symmetry(self):
        # Test that a basic V-cycle yields a symmetric linear operator.  Common
        # reasons for failure are problems with using the same rho for the
        # pres/post-smoothers and using the same block_D_inv for
        # pre/post-smoothers.

        n = 500
        A = poisson((n,), format="csr")
        smoothers = [
            ("gauss_seidel", {"sweep": "symmetric"}),
            ("schwarz", {"sweep": "symmetric"}),
            ("block_gauss_seidel", {"sweep": "symmetric"}),
            "jacobi",
            "block_jacobi",
        ]
        rng = arange(1, n + 1, dtype="float").reshape(-1, 1)
        Bs = [ones((n, 1)), hstack((ones((n, 1)), rng))]

        for smoother in smoothers:
            for B in Bs:
                ml = smoothed_aggregation_solver(A, B, max_coarse=10, presmoother=smoother, postsmoother=smoother)
                P = ml.aspreconditioner()
                x = rand(n)
                y = rand(n)
                assert_approx_equal(dot(P * x, y), dot(x, P * y))

    def test_poisson(self):
        cases = []

        cases.append((500,))
        cases.append((250, 250))
        cases.append((25, 25, 25))

        for case in cases:
            A = poisson(case, format='csr')

            np.random.seed(0)  # make tests repeatable

            x = sp.rand(A.shape[0])
            b = A*sp.rand(A.shape[0])  # zeros_like(x)

            ml = ruge_stuben_solver(A, max_coarse=50)

            res = []
            x_sol = ml.solve(b, x0=x, maxiter=20, tol=1e-12,
                             residuals=res)
            del x_sol

            avg_convergence_ratio = (res[-1]/res[0])**(1.0/len(res))

            assert(avg_convergence_ratio < 0.20)

 def test_improve_candidates(self):
     ##
     # test improve_candidates for the Poisson problem and elasticity, where rho_scale is 
     # the amount that each successive improve_candidates option should improve convergence
     # over the previous improve_candidates option.
     improve_candidates_list = [None, [('block_gauss_seidel', {'iterations' : 4, 'sweep':'symmetric'})] ]
     # make tests repeatable
     numpy.random.seed(0) 
     
     cases = []
     A_elas,B_elas = linear_elasticity( (60,60), format='bsr')
     #                Matrix                              Candidates    rho_scale
     cases.append( (poisson( (61,61),  format='csr'), ones((61*61,1)), 0.9 ) )
     cases.append( (A_elas,                           B_elas,       0.9 ) )
     for (A,B,rho_scale) in cases:
         last_rho = -1.0
         x0 = rand(A.shape[0],1) 
         b = rand(A.shape[0],1)
         for improve_candidates in improve_candidates_list:
             ml = smoothed_aggregation_solver(A, B, max_coarse=10, improve_candidates=improve_candidates)
             residuals=[]
             x_sol = ml.solve(b,x0=x0,maxiter=20,tol=1e-10, residuals=residuals)
             rho = (residuals[-1]/residuals[0])**(1.0/len(residuals))
             if last_rho == -1.0:
                 last_rho = rho
             else:
                 # each successive improve_candidates option should be an improvement on the previous
                 # print "\nimprove_candidates Test: %1.3e, %1.3e, %d\n"%(rho,rho_scale*last_rho,A.shape[0])
                 assert(rho < rho_scale*last_rho)
                 last_rho = rho

def driver():
    """ 
    Driver function for testing Laguerre polynomials
    """

    from scipy import rand, randn
    from numpy import ceil

    all_tests = ValidationContainer()

    """ Laguerre case """
    N = int(ceil(50*rand()))
    all_tests.extend(quadrature_test(N,alpha=0.))
    all_tests.extend(evaluation_test(N,alpha=0.))

    """ Random generalized case 1 """
    N = int(ceil(50*rand()))
    alpha = 10*rand()
    all_tests.extend(quadrature_test(N, alpha=alpha))
    all_tests.extend(evaluation_test(N, alpha=alpha))

    """ Random generalized case 2 """
    N = int(ceil(50*rand()))
    alpha = 10*rand()
    all_tests.extend(quadrature_test(N, alpha=alpha))
    all_tests.extend(evaluation_test(N, alpha=alpha))

    return all_tests

	def __init__(self, parent=None, width = 10, height = 12, dpi = 100, sharex = None, sharey = None):
		self.fig = Figure(figsize = (width, height), dpi=dpi, facecolor = '#FFFFFF')

		self.ax = Axes3D(self.fig)
#		n = 100
#		for c, zl, zh in [('r', -50, -25), ('b', -30, -5)]:
#		    xs = randrange(n, 23, 32)
#		    ys = randrange(n, 0, 100)
#		    zs = randrange(n, zl, zh)
		self.ax.scatter3D(S.rand(200), S.rand(200), S.rand(200))#, c = c,  alpha = 0.8)

		self.ax.set_xlabel('X Label')
		self.ax.set_ylabel('Y Label')
		self.ax.set_zlabel('Z Label')

#		self.ax = self.fig.add_subplot(111, sharex = sharex, sharey = sharey)
#		self.fig.subplots_adjust(left=0.1, bottom=0.1, right=0.9, top=0.9)
		self.xtitle="x-Axis"
		self.ytitle="y-Axis"
		self.PlotTitle = "Plot"
		self.grid_status = True
		self.xaxis_style = 'linear'
		self.yaxis_style = 'linear'
		self.format_labels()
		self.ax.hold(True)
		FigureCanvas.__init__(self, self.fig)
		#self.fc = FigureCanvas(self.fig)
		FigureCanvas.setSizePolicy(self,
			QtGui.QSizePolicy.Expanding,
			QtGui.QSizePolicy.Expanding)
		FigureCanvas.updateGeometry(self)

    def test_poisson(self):
        cases = []
        
        # perturbed Laplacian
        A = poisson( (50,50), format='csr' )
        Ai = A.copy(); Ai.data = Ai.data + 1e-5j*rand(Ai.nnz)
        cases.append((Ai, 0.25))
        
        # imaginary Laplacian
        Ai = 1.0j*A
        cases.append((Ai, 0.25))
        
        ## JBS:  Not sure if this is a valid test case
        ## imaginary shift 
        #Ai = A + 1.1j*scipy.sparse.eye(A.shape[0], A.shape[1])
        #cases.append((Ai,0.8))

        for A,rratio in cases:
            [asa,work] = adaptive_sa_solver(A, num_candidates = 1, symmetry='symmetric')
            #sa  = smoothed_aggregation_solver(A, B = ones((A.shape[0],1)) )
    
            b = zeros((A.shape[0],))
            x0 = rand(A.shape[0],) + 1.0j*rand(A.shape[0],)
    
            residuals0 = []
            residuals1 = []

            sol0 = asa.solve(b, x0=x0, maxiter=20, tol=1e-10, residuals=residuals0)
            #sol1 =  sa.solve(b, x0=x0, maxiter=20, tol=1e-10, residuals=residuals1)
           
            conv_asa = (residuals0[-1]/residuals0[0])**(1.0/len(residuals0))
            #conv_sa  = (residuals1[-1]/residuals1[0])**(1.0/len(residuals1))
            
            assert( conv_asa < rratio )

 def run(self):
 
  for i in scipy.arange(self.ns):
   # Atualiza velocidade
   
   self.v[i] = self.w*self.v[i] 
   self.v[i] = self.v[i] + self.c1*scipy.rand()*( self.bfp[i] - self.pop[i]) 
   self.v[i] = self.v[i] + self.c2*scipy.rand()*(self.bfg - self.pop[i])
   for j in range(Dim):
    if self.v[i][j] >= 52.6:
     self.v[i][j] = 52.6
    elif self.v[i][j] <= -52.6:
     self.v[i][j] = -52.6
   # Atualiza posicao
   self.pop[i] = self.pop[i] + self.v[i]
   
   self.fit[i],self.pop[i] = self.avalia_aptidao(self.pop[i])
   
   # Atualiza melhor posicao da particula
   if self.fit[i] < self.bfp_fitness[i]:
    self.bfp[i] = self.pop[i]
    self.bfp_fitness[i] = self.fit[i]
   if self.debug:
    print "self.fit[{0}] = {1} bfp_fitness = {2}".format(i,self.fit[i],self.bfp_fitness[i])
  # Atualiza melhor posicao global
   if  self.bfp_fitness[i] < self.bfg_fitness:
    self.bfg_fitness = self.bfp_fitness[i].copy()
    self.bfg = self.bfp[i].copy()

 def __new__(self, nbinaries=1e6):
     arr = sp.ones(nbinaries, dtype=[(name, 'f8') for name in ['period', 'mass_ratio', 'eccentricity', 'phase', 'theta', 'inclination']])
     arr['eccentricity'] = 0.
     arr['phase'] = sp.rand(nbinaries)
     arr['theta'] = sp.rand(nbinaries) * 2 * sp.pi
     arr['inclination'] = sp.arccos(sp.rand(nbinaries) * 2. - 1.)
     return arr.view(OrbitalParameters)

 def setUp(self):
     SP.random.seed(1)
     nr = 3
     nc = 5
     n_dim1 = 8
     n_dim2 = 12
     # truncation of soft kronecker
     self.n_trunk = 10
     Xr = SP.rand(nr, n_dim1)
     Xc = SP.rand(nc, n_dim2)
     Cr = dlimix.CCovSqexpARD(n_dim1)
     Cr.setX(Xr)
     Cc = dlimix.CCovLinearARD(n_dim2)
     Cc.setX(Xc)
     self.C = dlimix.CKroneckerCF()
     self.C.setRowCovariance(Cr)
     self.C.setColCovariance(Cc)
     # set kronecker index
     self.kronecker_index = dlimix.CKroneckerCF.createKroneckerIndex(nc, nr)
     self.n = self.C.Kdim()
     self.n_dim = self.C.getNumberDimensions()
     self.name = "CKroneckerCF"
     self.n_params = self.C.getNumberParams()
     params = SP.exp(SP.randn(self.n_params))
     self.C.setParams(params)

 def runSequences(self, num_actions=1, num_features=1, num_states=1,
                  num_interactions=10000, gamma=None, _lambda=None, lr=None, r_states=None):
     if r_states is None:
         r_states = [rand(num_features) for _ in range(num_states)]
     else:
         num_features = len(r_states[0])
         num_states = len(r_states)
     state_seq = [choice(r_states) for  _ in range(num_interactions)]
     action_seq = [randint(0, num_actions - 1) for  _ in range(num_interactions)]
     rewards = [ones(num_interactions), rand(num_interactions), action_seq, [s[0] for s in state_seq]]
     datas = [list(zip(state_seq, action_seq, r)) for r in rewards]        
     res = []        
     for algo in self.algos:
         res.append((algo.__name__, []))
         for d in datas:
             l = algo(num_actions, num_features)
             if gamma is not None:       
                 l.rewardDiscount = gamma
             if _lambda is not None:
                 l._lambda = _lambda
             if lr is not None:
                 l.learningRate = lr                     
             self.trainWith(l, d)
             res[-1][-1].append([dot(l._theta, s) for s in r_states])
     return res

    def __init__(self, type='random', pars=parameters()):

        if type == 'random':
            ee = (rand(pars['Ne'], pars['Ne']) < pars['p_ee'])
            ei = (rand(pars['Ne'], pars['Ni']) < pars['p_ei'])
            ii = (rand(pars['Ni'], pars['Ni']) < pars['p_ii'])
            ie = (rand(pars['Ni'], pars['Ne']) < pars['p_ie'])
            self.A = vstack((hstack((ee, ei)), hstack((ie, ii))))
            self.A[range(pars['Ne'] + pars['Ni']), range(pars['Ne'] + pars['Ni'])] = 0  # remove selfloops

        elif type == 'none':
            self.A = zeros((pars['N'], pars['N']))  # no connectivity

        elif type == 'uni_torus':  # torus with uniform connectivity profile
            self.A = zeros((pars['N'], pars['N']))

            # construct matrix of pairwise distance
            distMat = zeros((pars['N'], pars['N']))
            for n1 in range(pars['N']):
                coord1 = linear2grid(n1, pars['N_col'])
                for n2 in arange(n1 + 1, pars['N']):
                    coord2 = linear2grid(n2, pars['N_col']) - coord1  # this sets neuron n1 to the origin
                    distMat[n1, n2] = toric_length(coord2, pars['N_row'], pars['N_col'])
            distMat = distMat + distMat.transpose()

            # construct adjajency matrix
            for n1 in range(pars['N']):
                neighbor_ids = nonzero(distMat[:, n1] < pars['sigma_con'])[0]
                random.shuffle(neighbor_ids)
                idx = neighbor_ids[0:min([pars['ncon'], len(neighbor_ids)])]
                self.A[idx, n1] = 1
        else:
            print "type " + type + " not yet implemented"

def initial_cond(coords, mass, dipole, temp, F):
    cm_coords = coords - tile(center_of_mass(coords, mass), (coords.shape[0], 1))

    print "computing inertia tensor and principal axes of inertia"

    mol_I, mol_Ix = eig(inertia_tensor(cm_coords, mass))
    mol_I.sort()

    print "principal moments of inertia are: ", mol_I

    # compute the ratio of the dipole energy to the
    # rotational energy

    print "x = (mu*F / kB*T_R) = ", norm(dipole) * F / kB_au / temp

    # random initial angular velocity vector
    # magnitude set so that 0.5 * I * w**2.0 = kT
    w_mag = sqrt(2.0 * kB_au * temp / mol_I.mean())
    w0 = 2.0 * rand(3) - 1.0
    w0 = w0 / norm(w0) * w_mag

    # random initial orientation / random unit quaternion
    q0 = 2.0 * rand(4) - 1.0
    q0 = q0 / norm(q0)

    return q0, w0

    def test_symmetry(self):
        # Test that a basic V-cycle yields a symmetric linear operator.  Common
        # reasons for failure are problems with using the same rho for the
        # pres/post-smoothers and using the same block_D_inv for
        # pre/post-smoothers.

        n = 500
        A = poisson((n,), format='csr')
        smoothers = [('gauss_seidel', {'sweep': 'symmetric'}),
                     ('schwarz', {'sweep': 'symmetric'}),
                     ('block_gauss_seidel', {'sweep': 'symmetric'}),
                     'jacobi', 'block_jacobi']
        Bs = [ones((n, 1)),
              hstack((ones((n, 1)),
                      arange(1, n + 1, dtype='float').reshape(-1, 1)))]

        for smoother in smoothers:
            for B in Bs:
                ml = rootnode_solver(A, B, max_coarse=10,
                                     presmoother=smoother,
                                     postsmoother=smoother)
                P = ml.aspreconditioner()
                x = rand(n,)
                y = rand(n,)
                assert_approx_equal(dot(P * x, y), dot(x, P * y))

    def test_basic(self):
        """check that method converges at a reasonable rate"""

        for A, B, c_factor, symmetry, smooth in self.cases:
            A = csr_matrix(A)

            ml = rootnode_solver(A, B, symmetry=symmetry, smooth=smooth,
                                 max_coarse=10)

            numpy.random.seed(0)  # make tests repeatable

            x = rand(A.shape[0]) + 1.0j * rand(A.shape[0])
            b = A * rand(A.shape[0])
            residuals = []

            x_sol = ml.solve(b, x0=x, maxiter=20, tol=1e-10,
                             residuals=residuals)

            avg_convergence_ratio =\
                (residuals[-1] / residuals[0]) ** (1.0 / len(residuals))

            # print "Complex Test:   %1.3e,  %1.3e,  %d,  %1.3e" % \
            #    (avg_convergence_ratio, c_factor, len(ml.levels),
            #     ml.operator_complexity())
            assert(avg_convergence_ratio < c_factor)

 def test_dimensions(self):
     # 1D
     img = scipy.rand(10)
     template = scipy.random.randint(0, 2, (3))
     result = ght(img, template)
     self.assertEqual(result.ndim, 1, "Computing ght with one-dimensional input data failed.")
     # 2D
     img = scipy.rand(10, 11)
     template = scipy.random.randint(0, 2, (3, 4))
     result = ght(img, template)
     self.assertEqual(result.ndim, 2, "Computing ght with two-dimensional input data failed.")
     # 3D
     img = scipy.rand(10, 11, 12)
     template = scipy.random.randint(0, 2, (3, 4, 5))
     result = ght(img, template)
     self.assertEqual(result.ndim, 3, "Computing ght with three-dimensional input data failed.")
     # 4D
     img = scipy.rand(10, 11, 12, 13)
     template = scipy.random.randint(0, 2, (3, 4, 5, 6))
     result = ght(img, template)
     self.assertEqual(result.ndim, 4, "Computing ght with four-dimensional input data failed.")
     # 5D
     img = scipy.rand(3, 4, 3, 4, 3)
     template = scipy.random.randint(0, 2, (2, 2, 2, 2, 2))
     result = ght(img, template)
     self.assertEqual(result.ndim, 5, "Computing ght with five-dimensional input data failed.")