def Energy_condensate_full(Q, F1, x, y, H, mu, kappa, Ns) : 
    
    if Q==0 and F1 ==0 :
        return 1e14
    
    m = find_minimum(Q, F1, mu, kappa, Ns)
    if m[0] < H/2 :
        return 1e14
    result = 0
    for n1 in range(Ns) : 
        for n2 in range(Ns) : 
            M, dim = HamiltonianMatrix(n1, n2, Q, F1, 0, H, mu, kappa, Ns, 'T1') 
        
            B = np.identity(dim)
            B[dim/2:dim, dim/2:dim] = -np.identity(dim/2)
        
            eig = np.absolute(np.real(lin.eigvals(np.dot(B,M))))
        
            result += sum(eig)/2
            
            vec = [x[Ns * n1 + n2], np.conjugate(y[Ns * ((Ns - n1) % Ns) + (Ns - n2) % Ns])]
            
            result += np.dot(vec, np.dot(np.conj(vec).T, M))
            
    return result - 3 * Ns ** 2 * (np.abs(F1)**2 - np.abs(Q)**2)/2 - Ns**2 * mu*(1. + kappa) + Ns * H

    def test_linear_2p2s_with_dof_at_1_correlation_matrix(self):
        dof_index = 0
        output_index = 0

        model_instance = dict(models.model_data.model_structure)
        model_instance["parameters"] = numpy.array([2.0, 4.0])
        model_instance["inputs"] = numpy.array([[1.0, 10], [2.0, 20], [3.0, 30]])

        problem_instance = dict(models.model_data.problem_structure)
        problem_instance["output_indices"] = [output_index]
        problem_instance["inputs"] = model_instance["inputs"]
        problem_instance["parameters"] = [model_instance["parameters"][dof_index]]
        problem_instance["parameter_indices"] = [dof_index]

        dof = [1.0]
        # one standard deviation
        offset = -1
        measured = numpy.array([[1.0], [2.0], [3.0]]) + offset
        problem_instance["outputs"] = measured
        sens = [J_linear_2p2s([dof[0], model_instance["parameters"][1]], x) for x in model_instance["inputs"]]
        # identity covariance matrix of observation errors
        cov_obs_errs = numpy.identity(3)
        actual = metrics.confidence_measures.compute_covariance_matrix(sens, cov_obs_errs)
        expected = numpy.array([[1**2+2**2+3**2, 0], [0, 0]])
        [[self.assertEqual(act, exp) for act, exp in zip(acts, exps)] for acts, exps in zip(actual, expected)]
        
        # diagonal covariance matrix of observation errors
        multiplier = 2
        cov_obs_errs = numpy.identity(3) * multiplier
        actual = metrics.confidence_measures.compute_covariance_matrix(sens, cov_obs_errs)
        expected = numpy.array([[(1**2+2**2+3**2)/multiplier, 0], [0, 0]])
        [[self.assertEqual(act, exp) for act, exp in zip(acts, exps)] for acts, exps in zip(actual, expected)]

	def updateParameters(self, articlePicked, click,  userID):	
		self.counter +=1
		self.Wlong = vectorize(self.W)
		featureDimension = len(articlePicked.featureVector)
		T_X = vectorize(np.outer(articlePicked.featureVector, self.W.T[userID])) 
		self.A += np.outer(T_X, T_X)	
		self.b += click*T_X
		self.AInv = np.linalg.inv(self.A)
		self.UserTheta = matrixize(np.dot(self.AInv, self.b), len(articlePicked.featureVector)) 

		Xi_Matirx = np.zeros(shape = (featureDimension, self.userNum))
		Xi_Matirx.T[userID] = articlePicked.featureVector
		W_X = vectorize( np.dot(np.transpose(self.UserTheta), Xi_Matirx))
		self.batchGradient +=evaluateGradient(W_X, click, self.Wlong, self.lambda_, self.regu  )

		if self.counter%self.windowSize ==0:
			self.Wlong -= 1/(float(self.counter/self.windowSize)+1)*self.batchGradient
			self.W = matrixize(self.Wlong, self.userNum)
			self.W = normalize(self.W, axis=0, norm='l1')
			#print 'SVD', self.W
			self.batchGradient = np.zeros(self.userNum*self.userNum)
			# Use Ridge regression to fit W
		'''
		plt.pcolor(self.W_b)
		plt.colorbar
		plt.show()
		'''
		if self.W.T[userID].any() <0 or self.W.T[userID].any()>1:
			print self.W.T[userID]

		self.CoTheta = np.dot(self.UserTheta, self.W)
		self.BigW = np.kron(np.transpose(self.W), np.identity(n=len(articlePicked.featureVector)))
		self.CCA = np.dot(np.dot(self.BigW , self.AInv), np.transpose(self.BigW))
		self.BigTheta = np.kron(np.identity(n=self.userNum) , self.UserTheta)

def compose_matrix(scale=None, shear=None, angles=None, translate=None,perspective=None):
    """Return 4x4 transformation matrix from sequence of
    transformations.

    Code modified from the work of Christoph Gohlke link provided here
    http://www.lfd.uci.edu/~gohlke/code/transformations.py.html

    This is the inverse of the decompose_matrix function.

    Parameters
    -------------    
    scale : vector of 3 scaling factors
    shear : list of shear factors for x-y, x-z, y-z axes
    angles : list of Euler angles about static x, y, z axes
    translate : translation vector along x, y, z axes
    perspective : perspective partition of matrix

    Returns
    ---------
    matrix : 4x4 array
    

    Examples
    ----------
    >>> import math
    >>> import numpy as np
    >>> import dipy.core.geometry as gm
    >>> scale = np.random.random(3) - 0.5
    >>> shear = np.random.random(3) - 0.5
    >>> angles = (np.random.random(3) - 0.5) * (2*math.pi)
    >>> trans = np.random.random(3) - 0.5
    >>> persp = np.random.random(4) - 0.5
    >>> M0 = gm.compose_matrix(scale, shear, angles, trans, persp)
    """
    M = np.identity(4)
    if perspective is not None:
        P = np.identity(4)
        P[3, :] = perspective[:4]
        M = np.dot(M, P)
    if translate is not None:
        T = np.identity(4)
        T[:3, 3] = translate[:3]
        M = np.dot(M, T)
    if angles is not None:
        R = euler_matrix(angles[0], angles[1], angles[2], 'sxyz')
        M = np.dot(M, R)
    if shear is not None:
        Z = np.identity(4)
        Z[1, 2] = shear[2]
        Z[0, 2] = shear[1]
        Z[0, 1] = shear[0]
        M = np.dot(M, Z)
    if scale is not None:
        S = np.identity(4)
        S[0, 0] = scale[0]
        S[1, 1] = scale[1]
        S[2, 2] = scale[2]
        M = np.dot(M, S)
    M /= M[3, 3]
    return M

	def maximum_muscle_force(OptimalLengths, Angle, AngularVelocity, R, OptimalForces):
		"""
		OptimalLengths must be a 1X4 Matrix with optimal muscle lengths for 4 muscles.
		Angle should be the current angle. AngularVelocity should be the current 
		angular velocity. R is the 1X4 moment arm matrix. And OptimalForces should be
		a 4X1 matrix of forces produced by each muscle when at optimal lengths and 
		velocities.
		"""
		# Force-Length Considerations
		CurrentMuscleLengths = OptimalLengths.T - R.T*Angle
		NormalizedMuscleLengths = np.matrix([(CurrentMuscleLengths[i,0]/OptimalLengths.T[i,0])-1 for i in range(4)])
		# We subtract one from every Normalized Muscle Length to find the percentage 
		# above and below the optimal length.
		
		MaximumMuscleForce_FL = np.identity(4)*[force_length_curve(NormalizedMuscleLengths[0,i]) \
													+ passive_force_length_curve(NormalizedMuscleLengths[0,i],) \
														for i in range(4)]

		# Force-Velocity Considerations
		CurrentMuscleVelocity = -R.T*AngularVelocity
		NormalizedMuscleVelocity = [CurrentMuscleVelocity[i,0]/OptimalLengths.T[i,0] for i in range(4)]
		
		MaximumMuscleForce_FV = np.identity(4)* \
									[force_velocity_curve(NormalizedMuscleLengths[0,i],NormalizedMuscleVelocity[i]) \
										for i in range(4)]
		MaximumMuscleForce = (MaximumMuscleForce_FL+MaximumMuscleForce_FV)*[OptimalForces[0,i] for i in range(4)]
		return(MaximumMuscleForce)

def Hhfs(L,S,I):
    """Provides the I dot J matrix (magnetic dipole interaction)"""
    gS=int(2*S+1)
    Sx=jx(S)
    Sy=jy(S)
    Sz=jz(S)
    Si=identity(gS)

    gL=int(2*L+1)
    Lx=jx(L)
    Ly=jy(L)
    Lz=jz(L)
    Li=identity(gL)

    gJ=gL*gS
    Jx=kron(Lx,Si)+kron(Li,Sx)
    Jy=kron(Ly,Si)+kron(Li,Sy)
    Jz=kron(Lz,Si)+kron(Li,Sz)
    Ji=identity(gJ)
    J2=dot(Jx,Jx)+dot(Jy,Jy)+dot(Jz,Jz)

    gI=int(2*I+1)
    gF=gJ*gI
    Ix=jx(I)
    Iy=jy(I)
    Iz=jz(I)
    Ii=identity(gI)
    Fx=kron(Jx,Ii)+kron(Ji,Ix)
    Fy=kron(Jy,Ii)+kron(Ji,Iy)
    Fz=kron(Jz,Ii)+kron(Ji,Iz)
    Fi=identity(gF)
    F2=dot(Fx,Fx)+dot(Fy,Fy)+dot(Fz,Fz)
    Hhfs=0.5*(F2-I*(I+1)*Fi-kron(J2,Ii))
    return Hhfs

def rotation(ntheta):
    '''Find rotation matrix from axis-angle vector.'''
    theta = la.norm(ntheta)
    if theta == 0:
        return np.identity(3)
    Q = dual(ntheta / theta)
    return np.identity(3) - Q * math.sin(theta) + dot(Q, Q) * (1 - math.cos(theta))

 def on_key_down(self, keysym):
     # Control simulation
     if keysym.sym == sdl.SDLK_q:
         self.exit_()
     if keysym.sym == sdl.SDLK_r:
         self.reset()
     if keysym.sym == sdl.SDLK_p:
         self.runSimulation = not self.runSimulation
         print '-'*40
         print "Run Simulation: ", self.runSimulation
     if keysym.sym == sdl.SDLK_h:
         self.showHelp()
     if keysym.sym == sdl.SDLK_g:
         if self.worldView == 0:
             self.worldView = 1
             self.renderer.setProjection(self.sim.width, self.sim.height)
             self.renderer.matModelView = np.identity(4,'f')
             self.renderer.matView = np.identity(4,'f')
         elif self.worldView == 1:
             self.worldView = 0
             width, height = self.userController.agent.fov
             self.renderer.setProjection(width, height)
         self.renderer.worldView = self.worldView
     if keysym.sym == sdl.SDLK_f:
         self.frameStats.show()
     if keysym.sym == sdl.SDLK_i:
         # use for tests
         pass
     # Control agent
     self.userController.keyDown(keysym.sym)

def transform(args, workspace_lh, workspace_rh, nsubjs):

    transform_lh = np.zeros((args.nvoxel, args.nfeature, nsubjs))
    transform_rh = np.zeros((args.nvoxel, args.nfeature, nsubjs))

    if args.align_algo in ["ha", "srm_noneprob"]:
        transform_lh = workspace_lh["R"]
        transform_rh = workspace_rh["R"]
    elif args.align_algo in ["srm", "pca", "ica"]:
        bW_lh = workspace_lh["bW"]
        bW_rh = workspace_rh["bW"]
        for m in range(nsubjs):
            transform_lh[:, :, m] = bW_lh[m * args.nvoxel : (m + 1) * args.nvoxel, :]
            transform_rh[:, :, m] = bW_rh[m * args.nvoxel : (m + 1) * args.nvoxel, :]
    elif args.align_algo in ["ha_sm_retraction"]:
        bW_lh = workspace_lh["W"]
        bW_rh = workspace_rh["W"]
        for m in range(nsubjs):
            transform_lh[:, :, m] = bW_lh[:, :, m]
            transform_rh[:, :, m] = bW_rh[:, :, m]
    elif args.align_algo == "noalign":
        for m in range(nsubjs):
            transform_lh[:, :, m] = np.identity(args.nvoxel)
            transform_rh[:, :, m] = np.identity(args.nvoxel)
    else:
        exit("alignment algo not recognized")

    return (transform_lh, transform_rh)

def test_cell_triclinic():
    while True:
        rvecs = np.random.uniform(-1, 1, (3, 3))
        if abs(np.linalg.det(rvecs)) > 0.1:
            break
    cell = Cell(rvecs)

    # Test attributes
    assert cell.nvec == 3
    assert (cell.rvecs == rvecs).all()
    assert abs(cell.volume - abs(np.linalg.det(rvecs))) < 1e-10
    assert abs(np.dot(cell.gvecs, cell.rvecs.transpose()) - np.identity(3)).max() < 1e-5
    assert abs(np.dot(cell.gvecs.transpose(), cell.rvecs) - np.identity(3)).max() < 1e-5
    cell2 = Cell(-cell.rvecs)
    assert abs(cell2.volume - abs(np.linalg.det(rvecs))) < 1e-10
    for i in xrange(3):
        assert cell.get_rlength(i) == cell.rlengths[i]
        assert cell.get_glength(i) == cell.glengths[i]
        assert cell.get_rspacing(i) == cell.rspacings[i]
        assert cell.get_gspacing(i) == cell.gspacings[i]
        assert abs(cell.get_rlength(i) - 1.0/cell.get_gspacing(i)) < 1e-10
        assert abs(cell.get_glength(i) - 1.0/cell.get_rspacing(i)) < 1e-10

    # Test methods (1)
    vec1 = np.array([10.0, 0.0, 5.0])*angstrom
    cell.mic(vec1)
    cell.add_rvec(vec1, np.array([1,2,3]))

    # Test methods (2)
    check_frac_cart(cell)
    check_g_lincomb_dot_rvecs(cell)

def get_bands(h,output_file="BANDS2D_",nindex=[-1,1],
               nk=50,nsuper=1,reciprocal=False,
               operator=None,k0=[0.,0.]):
  """ Calculate band structure"""
  if h.dimensionality!=2: raise  # continue if two dimensional
  hk_gen = h.get_hk_gen() # gets the function to generate h(k)
  kxs = np.linspace(-nsuper,nsuper,nk)+k0[0]  # generate kx
  kys = np.linspace(-nsuper,nsuper,nk)+k0[1]  # generate ky
  kdos = [] # empty list
  kxout = []
  kyout = []
  if reciprocal: R = h.geometry.get_k2K() # get matrix
  else:  R = np.matrix(np.identity(3)) # get identity
  # setup a reasonable value for delta
  # setup the operator
  if operator is None:
    operator = np.matrix(np.identity(h.intra.shape[0]))
  os.system("rm -f "+output_file+"*") # delete previous files
  fo = [open(output_file+"_"+str(i)+".OUT","w") for i in nindex] # files        
  for x in kxs:
    for y in kxs:
      r = np.matrix([x,y,0.]).T # real space vectors
      k = np.array((R*r).T)[0] # change of basis
      hk = hk_gen(k) # get hamiltonian
      evals = lg.eigvalsh(hk) # eigenvalues
      epos = sorted(evals[evals>0]) # positive energies
      eneg = -np.array(sorted(np.abs(evals[evals<0]))) # negative energies
      for (i,j) in zip(nindex,range(len(nindex))): # loop over bands
        fo[j].write(str(x)+"     "+str(y)+"   ")
        if i>0: fo[j].write(str(epos[i-1])+"\n")
        if i<0: fo[j].write(str(eneg[abs(i)-1])+"\n")
  [f.close() for f in fo] # close file

def test_cell_cubic():
    rvecs = np.array([[9.865, 0.0, 0.0], [0.0, 9.865, 0.0], [0.0, 0.0, 9.865]])*angstrom
    cell = Cell(rvecs)

    # Test attributes
    assert cell.nvec == 3
    assert (cell.rvecs == rvecs).all()
    assert (cell.rspacings == 9.865*angstrom).all()
    assert (cell.gspacings == 1/(9.865*angstrom)).all()
    assert abs(cell.volume - (9.865*angstrom)**3) < 1e-10
    assert abs(np.dot(cell.gvecs, cell.rvecs.transpose()) - np.identity(3)).max() < 1e-5
    assert abs(np.dot(cell.gvecs.transpose(), cell.rvecs) - np.identity(3)).max() < 1e-5
    cell2 = Cell(-cell.rvecs)
    assert abs(cell2.volume - (9.865*angstrom)**3) < 1e-10
    for i in xrange(3):
        assert cell.get_rlength(i) == cell.rlengths[i]
        assert cell.get_glength(i) == cell.glengths[i]
        assert cell.get_rspacing(i) == cell.rspacings[i]
        assert cell.get_gspacing(i) == cell.gspacings[i]
        assert abs(cell.get_rlength(i) - 1.0/cell.get_gspacing(i)) < 1e-10
        assert abs(cell.get_glength(i) - 1.0/cell.get_rspacing(i)) < 1e-10

    # Test methods (1)
    vec1 = np.array([10.0, 0.0, 5.0])*angstrom
    cell.mic(vec1)
    assert abs(vec1 - np.array([0.135, 0.0, -4.865])*angstrom).max() < 1e-10
    cell.add_rvec(vec1, np.array([1,2,3]))
    assert abs(vec1 - np.array([10.0, 19.73, 24.73])*angstrom).max() < 1e-10

    # Test methods (2)
    check_frac_cart(cell)
    check_g_lincomb_dot_rvecs(cell)

def identity_like_generalized(a):
    a = asarray(a)
    if a.ndim == 3:
        return np.array([identity(a.shape[-2]) for ax in a])
    elif a.ndim > 3:
        raise ValueError("Not implemented...")
    return identity(a.shape[0])

def ale3dStrainOutToV(vecds):
    #takes 5 components of evecd and the 6th component is lndetv


    """convert from vecds representation to symmetry matrix"""
    eps = num.zeros([3,3],dtype='float64')
    #Akk_by_3 = sqr3i * vecds[5] # -p
    a = num.exp(vecds[5])**(1./3.)# -p
    t1 = sqr2i*vecds[0]
    t2 = sqr6i*vecds[1]

    eps[0, 0] =    t1 - t2
    eps[1, 1] =   -t1 - t2
    eps[2, 2] = sqr2b3*vecds[1]
    eps[1, 0] = vecds[2] * sqr2i
    eps[2, 0] = vecds[3] * sqr2i
    eps[2, 1] = vecds[4] * sqr2i

    eps[0, 1] = eps[1, 0]
    eps[0, 2] = eps[2, 0]
    eps[1, 2] = eps[2, 1]

    epstar=eps/a

    V=(num.identity(3)+epstar)*a
    Vinv=(num.identity(3)-epstar)/a

    return V,Vinv

    def _get_H(self, debug=False):
        """
        returns H_t as defined in algorithm 2
        
        Reference:
        https://en.wikipedia.org/wiki/Limited-memory_BFGS
        http://www.ccms.or.kr/data/pdfpaper/jcms21_1/21_1_117.pdf
        https://homes.cs.washington.edu/~galen/files/quasi-newton-notes.pdf
        """
        I = np.identity(len(self.w))
        
        if min(len(self.s), len(self.y)) == 0:
                print "Warning: No second order information used!"
                return I
            
        assert len(self.s) > 0, "s cannot be empty."
        assert len(self.s) == len(self.y), "s and y must have same length"
        assert self.s[0].shape == self.y[0].shape, \
            "s and y must have same shape"
        assert abs(self.y[-1]).sum() != 0, "latest y entry cannot be 0!"
        assert 1/np.inner(self.y[-1], self.s[-1]) != 0, "!"

        I = np.identity(len(self.s[0]))
        H = np.dot((np.inner(self.s[-1], self.y[-1]) / np.inner(self.y[-1],
                   self.y[-1])), I)

        for (s_j, y_j) in itertools.izip(self.s, self.y):
            rho = 1.0/np.inner(y_j, s_j)
            V = I - np.multiply(rho, np.outer(s_j, y_j))
            H = (V).dot(H).dot(V.T)
            H += np.multiply(rho, np.outer(s_j, s_j))

        return H

def test_symm_algorithm_equivalence():
    """Test different stabilization methods in the computation of modes,
    in the presence and/or absence of the discrete symmetries."""
    np.random.seed(400)
    for n in (12, 20, 40):
        for sym in kwant.rmt.sym_list:
            # Random onsite and hopping matrices in symmetry class
            h_cell = kwant.rmt.gaussian(n, sym)
            # Hopping is an offdiagonal block of a Hamiltonian. We rescale it
            # to ensure that there are modes at the Fermi level.
            h_hop = 10 * kwant.rmt.gaussian(2*n, sym)[:n, n:]

            if kwant.rmt.p(sym):
                p_mat = np.array(kwant.rmt.h_p_matrix[sym])
                p_mat = np.kron(np.identity(n // len(p_mat)), p_mat)
            else:
                p_mat = None

            if kwant.rmt.t(sym):
                t_mat = np.array(kwant.rmt.h_t_matrix[sym])
                t_mat = np.kron(np.identity(n // len(t_mat)), t_mat)
            else:
                t_mat = None

            if kwant.rmt.c(sym):
                c_mat = np.kron(np.identity(n // 2), np.diag([1, -1]))
            else:
                c_mat = None

            check_equivalence(h_cell, h_hop, n, sym=sym, particle_hole=p_mat,
                              chiral=c_mat, time_reversal=t_mat)

def Hfs(L,S,I):
    """Provides the L dot S matrix (fine structure)"""
    gS=int(2*S+1) #number of mS values
    Sx=jx(S)
    Sy=jy(S)
    Sz=jz(S)
    Si=identity(gS)

    gL=int(2*L+1)
    Lx=jx(L)
    Ly=jy(L)
    Lz=jz(L)
    Li=identity(gL)

    gJ=gL*gS
    Jx=kron(Lx,Si)+kron(Li,Sx)
    Jy=kron(Ly,Si)+kron(Li,Sy)
    Jz=kron(Lz,Si)+kron(Li,Sz)
    J2=dot(Jx,Jx)+dot(Jy,Jy)+dot(Jz,Jz)

    gI=int(2*I+1)
    Ii=identity(gI)
    gF=gJ*gI
    Fi=identity(gF)
    Hfs=0.5*(kron(J2,Ii)-L*(L+1)*Fi-S*(S+1)*Fi) # fine structure in m_L,m_S,m_I basis
    return Hfs

def get_H(ph_count1, at_count1, wc1, wa1, g1, ph_count2, at_count2, wc2, wa2, g2, m = 0.4, RWA=True):
	#------------------------------------------------------------------------------------------------------------------
	# get_H_err(ph_count, at_count, wc, wa, g, RWA)
	#------------------------------------------------------------------------------------------------------------------
	adiag1 = np.sqrt(np.arange(1, ph_count1+1))
	
	across1 = np.diagflat(adiag1, -1)
	a1 = np.diagflat(adiag1, 1)
	acrossa1 = np.dot(across1, a1)

	adiag2 = np.sqrt(np.arange(1, ph_count2+1))
	
	across2 = np.diagflat(adiag2, -1)
	a2 = np.diagflat(adiag2, 1)
	acrossa2 = np.dot(across2, a2)
	#------------------------------------------------------------------------------------------------------------------
	sigmadiag = [1]

	sigmacross = np.diagflat(sigmadiag, -1)
	sigma = np.diagflat(sigmadiag, 1)
	sigmacrosssigma = np.dot(sigmacross, sigma)
	#------------------------------------------------------------------------------------------------------------------
	ph1_dim = ph_count1+1
	I_ph1 = np.identity(ph1_dim)

	at1_dim = pow(2, at_count1)
	I_at1 = np.identity(at1_dim)
	
	ph2_dim = ph_count2+1
	I_ph2 = np.identity(ph2_dim)

	at2_dim = pow(2, at_count2)
	I_at2 = np.identity(at2_dim)
	#------------------------------------------------------------------------------------------------------------------
	if RWA:
		h1 = H1.get_H_RWA(ph_count1, at_count1, wc1, wa1, g1, RWA)
		h2 = H1.get_H_RWA(ph_count2, at_count2, wc2, wa2, g2, RWA)
	else:
		h1 = H1.get_H_EXACT(ph_count1, at_count1, wc1, wa1, g1)
		h2 = H1.get_H_EXACT(ph_count2, at_count2, wc2, wa2, g2)
		
	H = np.kron(h1, np.identity(h2.shape[0])) + np.kron(np.identity(h1.shape[0]), h2) 

	H1_m = np.kron(across1, np.identity(at1_dim))
	H1_m = np.kron(H1_m, a2)
	H1_m = np.kron(H1_m, np.identity(at2_dim))
	
	H2_m = np.kron(a1, np.identity(at1_dim))
	H2_m = np.kron(H2_m, across2)
	H2_m = np.kron(H2_m, np.identity(at2_dim))
	
	H_m = m * (H1_m + H2_m)
	#------------------------------------------------------------------------------------------------------------------
	H = np.matrix(H + H_m)
	H_size = np.shape(H)
	
	# print('H:\n', H, '\n')
	# # # print('H_size:', H_size, '\n')
	#------------------------------------------------------------------------------------------------------------------
	return H

    def __init__(self, population, sigma, **params):
        self.parents = population
        self.dim = len(self.parents[0])

        # Selection
        self.mu = params.get("mu", len(self.parents))
        self.lambda_ = params.get("lambda_", 1)

        # Step size control
        self.d = params.get("d", 1.0 + self.dim / 2.0)
        self.ptarg = params.get("ptarg", 1.0 / (5.0 + 0.5))
        self.cp = params.get("cp", self.ptarg / (2.0 + self.ptarg))

        # Covariance matrix adaptation
        self.cc = params.get("cc", 2.0 / (self.dim + 2.0))
        self.ccov = params.get("ccov", 2.0 / (self.dim ** 2 + 6.0))
        self.pthresh = params.get("pthresh", 0.44)

        # Internal parameters associated to the mu parent
        self.sigmas = [sigma] * len(population)
        # Lower Cholesky matrix (Sampling matrix)
        self.A = [numpy.identity(self.dim) for _ in range(len(population))]
        # Inverse Cholesky matrix (Used in the update of A)
        self.invCholesky = [numpy.identity(self.dim) for _ in range(len(population))]
        self.pc = [numpy.zeros(self.dim) for _ in range(len(population))]
        self.psucc = [self.ptarg] * len(population)

        self.indicator = params.get("indicator", tools.hypervolume)

def matrix_exp(M, ntaylor=20, nsquare=10):
    """Computes the exponential of a square matrix via a Taylor series.

    Calculates the matrix exponential by first calculating exp(M/(2**nsquare)),
    then squaring the result the appropriate number of times.

    Args:
       M: Matrix to be exponentiated.
       ntaylor: Optional integer giving the number of terms in the Taylor series.
          Defaults to 15.
       nsquare: Optional integer giving how many times the original matrix will
          be halved. Defaults to 15.

    Returns:
       The matrix exponential of M.
    """

    n = M.shape[1]
    tc = np.zeros(ntaylor + 1)
    tc[0] = 1.0
    for i in range(ntaylor):
        tc[i + 1] = tc[i] / (i + 1)

    SM = np.copy(M) / 2.0**nsquare

    EM = np.identity(n, float) * tc[ntaylor]
    for i in range(ntaylor - 1, -1, -1):
        EM = np.dot(SM, EM)
        EM += np.identity(n) * tc[i]

    for i in range(nsquare):
        EM = np.dot(EM, EM)
    return EM