def reflect(self, ray, intercept):
     """
     Takes a Ray object and its intercept with the optical element
     as arguments and returns the new normalised direction vector of 
     the reflected ray.
     
     If the ray does not intersect the optical element, it returns
     the original direction of the vector.
     """
     if intercept is None:
         k_refracted = ray.k()
         return k_refracted
     
     if self._C != 0:
         n_unnorm = intercept - self._O
         n = n_unnorm/(np.linalg.norm(n_unnorm))
         cos_dot_product = - np.inner(n, ray.k())
         if cos_dot_product < 0:
             n = -n
     else:
         n = np.array([0,0,1.])
         cos_dot_product = - np.inner(n, ray.k())
         if cos_dot_product < 0:
             n = -n
         
     #Vector form of reflection.
     n_dot_k = np.inner(n, ray.k())
     k_reflected_unnorm = ray.k() - 2*n_dot_k*n
     k_reflected = k_reflected_unnorm/(np.linalg.norm(k_reflected_unnorm))
     return k_reflected

    def __init__(self,reciprocals,path):
        '''
        Constructor.

        Parameters
        ----------
        reciprocals : iterable of 1d ndarray
            The translation vectors of the reciprocal lattice.
        path : str
            The str-formed path.
        '''
        path=path.replace(' ', '')
        assert path[0] in KMap.database and path[1]==':'
        space,path,database,reciprocals=path[0],path[2:].split(','),KMap.database[path[0]],np.asarray(reciprocals)
        if space=='L':
            assert len(reciprocals)==1
        elif space=='S':
            assert len(reciprocals)==2
            inner=np.inner(reciprocals[0],reciprocals[1])/nl.norm(reciprocals[0])/nl.norm(reciprocals[1])
            assert np.abs(inner)<RZERO
        elif space=='H':
            assert len(reciprocals)==2
            inner=np.inner(reciprocals[0],reciprocals[1])/nl.norm(reciprocals[0])/nl.norm(reciprocals[1])
            assert np.abs(np.abs(inner)-0.5)<RZERO
            if np.abs(inner+0.5)<RZERO: reciprocals[1]=-reciprocals[1]
        for segment in path:
            segment=segment.split('-')
            assert len(segment)==2
            self.append([reciprocals.T.dot(database[segment[0]]),reciprocals.T.dot(database[segment[1]])])

    def _parallax_correction(self, image_num, s1):
        """ 
        Compute the parallax correction, which accounts for the shift in the position 
        between where photons are detected and recorded due to non-negligible detector
        thickness. Below code follows the framework outlined by cxtbx/dxtbx:
        https://github.com/cctbx/cctbx_project/blob/\
        b0460954eac07a3a3639dbe7efd942c21e12ebc1/dxtbx/model/parallax_correction.h
        
        """

        u_fast = self.system['f']/np.linalg.norm(self.system['f'])
        u_slow = self.system['s']/np.linalg.norm(self.system['s'])

        normal = np.cross(u_fast, u_slow)
        dist = np.dot(self.system['p'][self.system['img2batch'][image_num]], normal)
        if dist < 0:
            normal = -1*normal

        cos_t = np.inner(s1.copy(), normal)
        attenuation = 1.0/self.system['mu'] - (self.system['t0']/cos_t + 1.0/self.system['mu'])\
            *np.exp(-self.system['mu']*self.system['t0']/cos_t)

        xcorr = attenuation*np.inner(s1.copy(), u_fast)
        ycorr = attenuation*np.inner(s1.copy(), u_slow)

        corrections = np.zeros(s1.shape)
        corrections[:,0], corrections[:,1] = xcorr, ycorr

        return corrections

def _batch_omp_step(G, alpha_0, m, eps_0=None, eps=None):
    idx = []
    L = np.ones((1,1))
    alpha = alpha_0
    eps_curr = eps_0
    delta = 0
    it = 0
    if eps == None:
        stopping_condition = lambda: it == m
    else:
        stopping_condition = lambda: eps_curr <=eps

    while not stopping_condition():
        lam = np.abs(alpha).argmax()
        if len(idx) > 0:
            w = linalg.solve_triangular(L, G[idx, lam],
                                        lower = True, unit_diagonal=True)
            L = np.r_[np.c_[L, np.zeros(len(L))],
                      np.atleast_2d(np.append(w, np.sqrt(1-np.inner(w,w))))]
        idx.append(lam)
        it +=1
        Ltc = linalg.solve_triangular(L, alpha_0[idx], lower=True)
        gamma = linalg.solve_triangular(L, Ltc, trans=1, lower=True)
        beta = np.dot(G[:, idx], gamma)
        alpha = alpha_0 - beta
        if eps != None:
            eps_curr += delta
            delta = np.inner(gamma, beta[idx])
            eps_curr -= delta
    return gamma, idx

def calculatePhi(X, B, Pi, n, epsilon=1e-4, C=None):
    """
    Calculate the matrix for multiplicative update

    Parameters
    ----------
    X       : the observed tensor
    B       : the factor matrix associated with mode n
    Pi      : the product of all matrices but the n-th from above
    n       : the mode that we are trying to solve the subproblem for
    epsilon : the 
    C       : the augmented / non-augmented tensor (\alpha u \Psi or B \Phi) in sparse form
    """
    Phi = None
    if X.__class__ == sptensor.sptensor:
        Phi = -np.ones((X.shape[n], B.shape[1]))
        xsubs = X.subs[:,n]
        if C !=  None:
            v = np.sum(np.multiply(B[xsubs,:], Pi) + C, axis=1)
        else:
            v = np.sum(np.multiply(B[xsubs,:], Pi), axis=1)
        wvals = X.vals.flatten() / v
        for r in range(B.shape[1]):
            Phi[:,r] = accumarray.accum_np(xsubs, np.multiply(wvals, Pi[:,r]), size=X.shape[n])
    else:
        Xn = tenmat.tenmat(X,[n])
        V = np.inner(B,Pi)
        W = Xn.data / np.maximum(V, epsilon)
        Phi = np.inner(W, Pi.transpose())  
    return Phi

 def calculate_alpha(self, error, jacobi):
     """
     :param error: vector which contains the error
     :param jacobi: jacobian matrix
     """
     JJTe = np.dot(np.dot(jacobi, np.transpose(jacobi)), error)
     return float(np.inner(error, JJTe)) / float(np.inner(JJTe, JJTe))

    def intersects(self, seg):
        '''Calculate intersection point bewteen self and another
        segment seg.
        Return True is two segments overlap, False if two segments do not
        intersect, intersection point if two segments have intersection.
        Read documents for mechanism.'''
        p = np.array(self.left.coor)
        r = np.array(self.right.coor) - p
        q = np.array(seg.left.coor)
        s = np.array(seg.right.coor) - q

        if np.cross(r, s) == 0:
            if np.cross(q - p, r) != 0:
                return False
            elif 0 <= np.inner(q - p, r) <= np.inner(r, r) or \
                0 <= np.inner(p - q, s) <= np.inner(s, s):
                return True
            else:
                return False
        else:
            # np.array() and np.cross() does not guarantee float number.
            btm = float(np.cross(r, s))
            t = np.cross(q - p, s) / btm
            u = np.cross(q - p, r) / btm
            if 0 < t < 1 and 0 < u < 1:
                coor = tuple(p + t * r)
                intersection = Point(coor)
                return intersection
            else:
                return False

def cholesky_inplace(A):
    for i in xrange(A.shape[0]):
        A[i,i+1:] = 0.
        for j in range(i):
            A[i,j] = (A[i,j] - np.inner(A[i,:j],A[j,:j])) / A[j,j]
        sl = A[i,:i]
        A[i,i] = sqrt(A[i,i] - np.inner(sl, sl))

    def ttv(self, v, dims):
        """ 
        Computes the product of the Kruskal tensor with the column vector along
        specified dimensions.
        
        Parameters
        ----------
        v - column vector 
        dims - dimensions to multiply the product

        Returns
        -------
        out : 
        """
        (dims,vidx) = tools.tt_dimscheck(dims, self.ndims(), len(v));
        remdims = np.setdiff1d(range(self.ndims()), dims);
        
        ## Collapse dimensions that are being multiplied out
        newlmbda = self.lmbda;
        for i in range(self.ndims()):
            newlmbda = np.inner(np.inner(self.U[dims[i]], v[vidx[i]]));
        
        if len(remdims) == 0:
            return np.sum(newlmbda);
        
        return ktensor(newlmbda, self.u[remdims]);

def lanczos(psi0, H, N=200, stabilize=False):
    """Perform a Lanczos iteration building the tridiagonal matrix T and ONB of the Krylov space."""
    if psi0.ndim != 1:
        raise ValueError("psi0 should be a vector, "
                         "i.e., a numpy array with a single dimension of len 2**L")
    if H.shape[1] != psi0.shape[0]:
        raise ValueError("Shape of H doesn't match len of psi0.")
    psi0 = psi0/np.linalg.norm(psi0)
    vecs = [psi0]
    T = np.zeros((N, N))
    psi = H @ psi0  # @ means matrix multiplication
    # and works both for numpy arrays and scipy.sparse.csr_matrix
    alpha = T[0, 0] = np.inner(psi0.conj(), psi).real
    psi = psi - alpha* vecs[-1]
    for i in range(1, N):
        beta = np.linalg.norm(psi)
        if beta  < 1.e-13:
            print("Lanczos terminated early after i={i:d} steps:"
                  "full Krylov space built".format(i=i))
            T = T[:i, :i]
            break
        psi /= beta
        # note: mathematically, psi should be orthogonal to all other states in `vecs`
        if stabilize:
            for vec in vecs:
                psi -= vec * np.inner(vec.conj(), psi)
            psi /= np.linalg.norm(psi)
        vecs.append(psi)
        psi = H @ psi - beta * vecs[-2]
        alpha = np.inner(vecs[-1].conj(), psi).real
        psi = psi - alpha * vecs[-1]
        T[i, i] = alpha
        T[i-1, i] = T[i, i-1] = beta
    return T, vecs

 def filterbydistancefromplane(self, terncoordlist, compvert0, compvert1, compvert2, critdist, withintriangle=True,  affine=False, invlogic=False, returnall=False): #not sure fi affine transformation makes sense here but haven't through through it
     xyzarr=numpy.array(self.toCart(terncoordlist, affine=affine)).T
     xyz0=numpy.array(self.toCart([compvert0], affine=affine)).T[0]
     xyz1=numpy.array(self.toCart([compvert1], affine=affine)).T[0]
     xyz2=numpy.array(self.toCart([compvert2], affine=affine)).T[0]
     nhat=numpy.cross(xyz1-xyz0, xyz2-xyz1)
     nhat/=numpy.linalg.norm(nhat)
     
     distfromplane=numpy.array([numpy.abs(numpy.linalg.norm(numpy.inner(nhat, xyz1-xyz))) for xyz in xyzarr])
     intriangle=None
     if returnall or withintriangle:
         xphat=xyz1-xyz0
         xphat/=numpy.linalg.norm(xphat)
         yphat=numpy.cross(nhat, xphat)
         xyparr=numpy.array([[numpy.inner(xyz-xyz0, xphat), numpy.inner(xyz-xyz0, yphat)] for xyz in xyzarr]) #this makes xyz0 the origin, x axis points to xyz1
         xyp_verts=numpy.array([[numpy.inner(xyz-xyz0, xphat), numpy.inner(xyz-xyz0, yphat)] for xyz in [xyz0, xyz1, xyz2]])
     if withintriangle:
         intriangle=numpy.array([self.point_wrt_polygon(xyp, xyp_verts) for xyp in xyparr])
         if invlogic:
             inds=numpy.where(numpy.logical_not((distfromplane<=critdist) & (intriangle==1)))[0]
         else:
             inds=numpy.where((distfromplane<=critdist) & (intriangle==1))[0]
     else:
         if invlogic:
             inds=numpy.where(numpy.logical_not(distfromplane<=critdist))[0]
         else:
             inds=numpy.where(distfromplane<=critdist)[0]
     if returnall:
         if intriangle is None:
             intriangle=numpy.array([self.point_wrt_polygon(xyp, xyp_verts) for xyp in xyparr])
         return inds, distfromplane, xyparr, xyp_verts,intriangle#xyparr is array of x,y projections into the selected plan with xyz0 as the origin
     else:
         return inds

def dft(t,y,f):
	"""
	It's written to be as fast as possible.
	I've used numpy ufuncs. If you don't use them, this function becomes
	EXTREMELY slow!
	Maybe I'll have to write it in C later but I'm satisfied now.
	"""
	
	#  t  = numpy.delete( gv.t, gv.delete ) # gv.t without deleted points
	nt = float(len(t))
	nf=len(f)

	dft = np.zeros(nf, dtype=float) # Allocate memory
	w = 2*np.pi*np.array(f)                  # Angular frequencies
  
	# Transform y values subtracted from mean
	#dy = []
	#for i in range(0,ny):
	#y = np.delete( Cy[i], delete ) # Cy without the deleted points
	mean = np.mean( y )          # Scalar
	dy = y-mean             # list of lists dim: [ny,nt]

	for k in range(0,nf):
		C2 = np.cos( w[k]*t )           # Array with nt positions
		S2 = np.sin( w[k]*t )           # Array with nt positions
		#for i in range(0,ny):
		sum1  = np.inner(dy,C2)    # Inner vector product between dy and C2
		sum2  = np.inner(dy,S2)    # Inner vector product between dy and C2
		dft[k] = 2.0*np.hypot(sum1,sum2)/nt # Scalar
	return dft

def get_uv_coords(array, station1, station2, ra, dec, ha):

    # Convert hour angle, ra and dec to radians
    ha %= 24.0
    ha *= pi / 12.0
    center = array.arrxyz
    ra *= pi / 180.0
    dec *= pi / 180.0
    # Same for lat/long of array center
    latitude = vlti.latitude * pi / 180.0
    longitude = vlti.longitude * pi / 180.0

    # Calculate an "east" unit vector
    east = np.array([center[1]*center[2],-center[0]*center[2],0])
    if center[2] > 0: east *= -1
    east /= np.sqrt(np.sum(east**2))
    # Calculate a "north" unit vector
    north = np.cross(center, east)
    north /= np.sqrt(np.sum(north**2))

    up = center / np.sqrt(np.sum(center**2))

    # Eq. 3 of Segransan 2007
    B = np.array([np.inner(s2.staxyz - s1.staxyz, north),
                  np.inner(s2.staxyz - s1.staxyz, east),
                  np.inner(s2.staxyz - s1.staxyz, up)])

    # From Lawson's book
    u = B[1] * cos(ha) - B[0] * sin(latitude) * sin(ha)
    v = B[1] * sin(dec) * sin(ha) + B[0] * (sin(latitude) * sin(dec) * cos(ha) + cos(latitude) * cos(dec))

    return u, v

def compute_voroni_area_of_triangle(w1,w2,l1,l2):
    """
    computes the part of the triangle for the Voroni area
    descriped in [1]

            x_i
          +
          |\
          | \
          |  \
        l1|   \l2
          |    \
          |     \
          |w1  w2\
          +-------+
        x_j  l3  x_{j+1}

    """
    # Check if triangle is obtuse in x_i
    if w1+w2 < pi/2.0:
        # Then get Area(T)/2
        return norm(cross(l1,l2))/4.0
    
    if w1 > pi/2.0 or w2 > pi/2:
        # Then get Area(T)/4
        return norm(cross(l1,l2))/8.0
            

    #Else use formula on page 9 in [1]
    return ((1/tan(w1))*inner(l2,l2) + (1/tan(w2))*inner(l1,l1))/8.0

    def _simplex(self, x, c, B, Binv):
        iteration = 1
        z = np.inner(c,x) # first computation of z
        self.log("starting with z = {}".format(z),2)
        while iteration <= self.max_it : # avoid infinite loop
            q , rq = self._blandRule(c, B, Binv) # get q, entring basic variable
            if q == -1 : # => optimal
                self.log("simplex it {:2d}: Optimal solution found, z* = {:.2f}".format(iteration, z),2)
                return True, x, B, np.inner(c,x) # x is optimal with basis B, recompute z for precision
            # compute directions u = -dq = Binv*A[q]
            u = np.dot(Binv, self.A[:,q])
            # select B(p) leaving variable
            theta = INF # infinite at first
            p = -1
            for i in range(len(u)):
                if u[i] > 0 : # for positive components (negative for dq)
                    theta_i = x[B[i]]/u[i]
                    if theta_i < theta :
                        theta = theta_i
                        p = i
            if theta == INF: # => all directions are non-positive => -d >= 0
                self.log("simplex it {:2d}: Infinite ray found, unbounded problem".format(iteration),2)
                return False, x, B, -INF # infinite direction, z = -infinite
            # compute new z
            z += theta*rq
            self.log("simplex it {:2d}: B({:d}) = {:2d} <-> {:2d}, theta = {:.3f}, \
z = {:.2f}".format(iteration, p, B[p], q, theta,z),2)
            # compute new feasible solution
            x[B] -= theta*u # move along direction
            x[q] = theta
            x[B[p]] = 0 # put a real 0 to fix epsilons
            B[p] = q  # replace basic variable in basis
            self._recomputeBinv(p, u, Binv) # compute new Binv after changing column p of Binv by u
            iteration += 1

def test_inner_scalar_and_matrix_of_objects():
    # Ticket #4482
    # 2018-04-29: moved here from core.tests.test_multiarray
    arr = np.matrix([1, 2], dtype=object)
    desired = np.matrix([[3, 6]], dtype=object)
    assert_equal(np.inner(arr, 3), desired)
    assert_equal(np.inner(3, arr), desired)

def algo_backpropagation(X,eta,seuil,sigma,nu,W,Y,T):
	res=[]
	iter=0
	for i in xrange(0,len(X)):
		sortie=calcul_couche(X[i],W)
		res.append(np.inner(nu,sortie))
	error_list=[0]*len(Y)
	for i in xrange(0,len(Y)):
		error_list[i]=res[i]-Y[i]
	error_list=np.array(error_list)
	error=(1+0.0)/len(X)*linalg.norm(error_list)**2
	while seuil<=error and iter<=T:
		for i in xrange(0,len(X)): 	
			coeff=back_prop(W,nu,X[i],Y[i],X,0.01)
			for k in xrange(0,W.shape[0]):
				ref=coeff[k+1]
				ref.shape=(1,len(ref))
				W[k,:]=ref
			ref_0=coeff[0]
			ref_0=np.array(ref_0)
			nu_res=[0]*3
			for i in xrange(0,3):
				nu_res[i]=ref_0[i][0]
			nu=nu_res
		res=[]
		for i in xrange(0,len(X)):
			sortie=calcul_couche(X[i],W)
			res.append(np.inner(nu,sortie))
		error_list=[0]*len(Y)
		for i in xrange(0,len(Y)):
			error_list[i]=res[i]-Y[i]
		error_list=np.array(error_list)
		error=(1+0.0)/len(X)*linalg.norm(error_list)**2
		iter=iter+1
		print error

    def inverse_kinematics(self, effector_name, transform):
        '''solve the inverse kinematics

        :param str effector_name: name of end effector, e.g. LLeg, RLeg
        :param transform: 4x4 transform matrix
        :return: list of joint angles
        '''
        
        fitting_joints = {key : self.perception.joint[key] for key in CHAINS[effector_name]} #get (joint->angle ) for all effector joints
        
        effector_chain = CHAINS[effector_name] #end effector for the given effectr chain
        end_effector = effector_chain[-1] #end_effector
        
        Te = self.extract_values(transform)
        theta = np.zeros(len(fitting_joints)) #initial return value
        
                    
        while True: #break loop if angles are satisfiable
        
            self.forward_kinematics(fitting_joints) #execute forward kinematics
            T = self.transforms #result of forward kinematics for effector joints
            error = Te - self.extract_values(T[end_effector]) #calculate distance between target and actual result
            
            jacobi = self.calculate_jacobian_matrix(Te, T, effector_chain) #get jacobian matrix
            
            alpha = self.calculate_alpha(error, jacobi) #get scalar for jacobian tranpose method
            
            d_theta = alpha * np.dot(jacobi.T, error)#jacobian transpose method
            theta += d_theta #update theta value
            self.update_joint_angles(fitting_joints, theta) #update joint angles for next iteration
            print np.inner(error, error)
            if np.inner(error, error) < EPSILON: #break if scalar product of error is lower than tolerated error
                break
            
        return theta

def directionalDisplacement(_firstPoints, _secondPoints, _direction):
    #Given a direction to displace in, determine the displacement needed to get it out
    first_dots = numpy.inner(_firstPoints, _direction)
    second_dots = numpy.inner(_secondPoints, _direction)
    projected_displacement = max(second_dots)-min(first_dots)
    norm_sqr = 1.0 if _direction == [0, 0] else _direction[0]*_direction[0]+_direction[1]*_direction[1]
    return [projected_displacement/norm_sqr*_direction[0], projected_displacement/norm_sqr*_direction[1]]

 def intersect(self, p,d,v0,v1,v2):
     e1 = self.vectorSubtract(v1, v0)
     e2 = self.vectorSubtract(v2, v0)
 
     h = numpy.cross(d, e2)
     a = numpy.inner(e1, h)
 
     if a > -0.00001 and a < 0.00001:
         return False
 
     f = 1.0 / a
     s = self.vectorSubtract(p,v0)
     u = f * numpy.inner(s,h)
 
     if u < 0.0 or u > 1.0:
         return False
 
     q = numpy.cross(s, e1)
     v = f * numpy.inner(d, q)
 
     if v < 0.0 or u + v > 1.0:
         return False
 
     t = f * numpy.inner(e2, q)
 
     if t > 0.00001:
         return True
     else:
         return False