def spiral_sphere(N, Om=2*pi, b=array((0, 0, 1))):
    """
    Internal helper function for the raycasting that returns an array of
    unit vectors (N, 3) giving equally distributed directions on a part of
    sphere given by the center direction b and the solid angle Om
    """
    # first produce 'equally' distributed directions in spherical coords
    o = 4*pi/Om
    h = -1+ 2*arange(N)/(N*o-1.)
    theta = arccos(h)
    phi = zeros_like(theta)
    for i, hk in enumerate(h[1:]):
        phi[i+1] = phi[i]+3.6/sqrt(N*o*(1-hk*hk)) % (2*pi)
    # translate to cartesian coords
    xyz = vstack((sin(theta) * cos(phi), sin(theta) * sin(phi), cos(theta)))
    # mirror everything on a plane so that b points into the center
    a = xyz[:, 0]
    b = b/norm(b)
    ab = (a-b)[:, newaxis]
    if norm(ab)<1e-10:
        return xyz
    # this is the Householder matrix for mirroring
    H = identity(3)-dot(ab, ab.T)/dot(ab.T, a)
    # actual mirroring
    return dot(H, xyz)

 def j(self, x, u):
     j = super(CtmAdmm, self).j(x, u)
     x_strip = self.strip_x(x)
     u_strip = self.strip_u(u)
     ramp_strip = self.strip_ramp(u)
     j += tog.t * self.yl.dot(x_strip) + self.beta / 2.0 * norm(x_strip - self.xbarl)
     j += tog.t * self.yr.dot(u_strip) + self.beta / 2.0 * norm(u_strip - self.xbarr)
     j += tog.t * self.yramp.dot(ramp_strip) + self.beta / 2.0 * norm(ramp_strip - self.xbarramp)
     return j

 def test_3(self):
     kernel=GaussianKernel(sigma=10)
     X=randn(3000,10)
     K_chol, I, R, W=incomplete_cholesky(X, kernel, eta=0.001)
     K=kernel.kernel(X)
     
     self.assertEqual(shape(K_chol), (len(I), (len(I))))
     self.assertEqual(shape(R), (len(I), (len(X))))
     self.assertEqual(shape(W), (len(I), (len(X))))
     
     self.assertLessEqual(norm(K-R.T.dot(R)), 1)
     self.assertLessEqual(norm(K-W.T.dot(K_chol.dot(W))), 1)

def create_square_from_two_points(a_point,b_point):
    """   
    Let a and b two points, this function will return a square (polygon)
    in which the a_point is at one corner and in the other extreme corner a point in the direction of b.
    
    Parameters
    ----------
        
        a_point : geometry (shapely Point)
        b_point : geometry (shapely Point) the direction is what matters
    
    """

        
    d = create_rectangle_from_two_points(a_point, b_point)
    
    apt = d['a']
    bpt = d['b']
    appt = d['a_p']
    bppt = d['b_p']

    a = asarray(apt)
    b = asarray(bpt)
    bp = asarray(bppt)
    ap = asarray(appt)
    
    a_m_bp = a - bp
    a_m_b = a - b
    a_m_ap = a - ap
    
    n_a_m_b = norm(a_m_b)
    n_a_m_ap = norm(a_m_ap)
    n_a_m_bp = norm(a_m_bp)
    
    if n_a_m_b < 0 :
        sig = -1.0
    elif n_a_m_b > 0 :
        sig = 1.0
    else:
        logger.error("The points are equal. It's not possible to generate an area with only one point.")
    
    
    p = (sig * ( n_a_m_ap / n_a_m_bp ) ) * a_m_bp

    pp = p + (a_m_ap) + a
    
    #n_chiqui = min((n_a_m_bp,n_a_m_ap))
    
    pp_pt = Point(pp)
    a_point = Point(a_point)
    new_dic = create_rectangle_from_two_points(a_point, pp_pt)
    return new_dic

    def _compute_colors(self, array_x, array_y):

        # on calcule le maximum absolu de toutes valeurs pour former un carré
        abs_maximum = max([max(map(abs,array_x)), max(map(abs,array_y))])
        diagonal_length = norm(array([abs_maximum, abs_maximum])) # longueur de la projection
        diag = array([diagonal_length, diagonal_length])
        anti_diag = array([-diagonal_length, diagonal_length])

        # on instancie le gradient de couleur sur le modèle de couleur du centre
        linear_normalizer = mpl.colors.Normalize(vmin=-abs_maximum, vmax=abs_maximum)
        log_normalizer = mpl.colors.SymLogNorm(abs_maximum/5, vmin=-abs_maximum, vmax=abs_maximum)
        r_to_b_gradient = cm.ScalarMappable(norm=linear_normalizer, cmap=redtoblue)

        # on calcule le produit scalaire de chaque valeur avec la diagonale
        # ensuite, on calcule la couleur à partir de la valeur de la projection sur la diagonale
        hex_color_values = []
        for i, x in enumerate(array_x):
            # on calcule les produits scalaire du point avec la diagonale et l'antidiagonale
            scal_p_diag = dot(array([array_x[i], array_y[i]]), diag) / diagonal_length
            scal_p_antidiag = dot(array([array_x[i], array_y[i]]), anti_diag) / diagonal_length

            #on calcule le gradient de couleur sur la diagonale
            on_diag_color = colorConverter.to_rgb(r_to_b_gradient.to_rgba(scal_p_diag))
            # puis on utilise cette couleur (en rgb) pour définir un gradient, dont la valeur sera estimée
            # sur l'antidiagonale
            on_diag_gradient = make_white_gradient(on_diag_color, log_normalizer)
            final_color = on_diag_gradient.to_rgba(scal_p_antidiag)

            #on traduit en HEX
            hex_color_values.append(rgb2hex(colorConverter.to_rgb(final_color)))

        return hex_color_values, abs_maximum

    def test_mode_newton_2d(self):
        X = asarray([-1, 1])
        X = reshape(X, (len(X), 1))
        y = asarray([+1 if x >= 0 else -1 for x in X])
        covariance = SquaredExponentialCovariance(sigma=1, scale=1)
        likelihood = LogitLikelihood()
        gp = GaussianProcess(y, X, covariance, likelihood)
        laplace = LaplaceApproximation(gp, newton_start=asarray([3, 3]))
        
        f_mode, _, steps = laplace.find_mode_newton(return_full=True)
        F = linspace(-10, 10, 20)
        D = zeros((len(F), len(F)))
        for i in range(len(F)):
            for j in range(len(F)):
                f = asarray([F[i], F[j]])
                D[i, j] = gp.log_posterior_unnormalised(f)
           
        idx = unravel_index(D.argmax(), D.shape)
        empirical_max = asarray([F[idx[0]], F[idx[1]]])
        
        pcolor(F, F, D)
        hold(True)
        plot(steps[:, 0], steps[:, 1])
        plot(f_mode[1], f_mode[0], 'mo', markersize=10)
        hold(False)
        colorbar()
        clf()
#        show()
           
        self.assertLessEqual(norm(empirical_max - f_mode), 1)

def on_mouse_move(event):
	if event.press_event is None:
		return

	modifiers = event.modifiers
	pos = event.press_event.pos
	if is_in_view(pos, view1.camera):
		if modifiers is not ():
			if 1 in event.buttons and modifiers[0].name=='Control':
				# Translate: camera._scene_transform.imap(event.pos)
				p1 = np.array(pos)[:2]
				p2 = np.array(event.last_event.pos)[:2]
				p1 = p1 - view1.pos
				p2 = p2 - view1.pos
				# print p1,p2
				p1s = view1.camera._scene_transform.imap(p1)[:2]
				p2s = view1.camera._scene_transform.imap(p2)[:2]
				print p1s, p2s
				pos_ = np.vstack((p2s,p1s))
				# print pos_
				measure_line.set_data(pos=pos_)
				measure_line.visible = True
				d_pixel = norm(pos_[1,:]-pos_[0,:])
				d_um = d_pixel*get_mpp(_id)
				print 'distance =',d_um
				measure_text.visible = True
				measure_text.text = '%.2f um' % d_um
				measure_text.pos = pos_[1,:]
				measure_text.pos[0] -= 10
				event.handled = True

    def integrate_gyro( self ):
        w_gyr = array( self.gyr.get_gyr(  ) ) \
                * self.deg_to_rad;
        g_acc = array( self.acc.get_acc(  ) );
        dg    = cross( self.g, w_gyr );

        is_to_suppress_acc                 \
            = norm( g_acc ) - self.g0_norm \
            > self.d_gravity_threashold;
        #
        
        for i in range( 3 ):
            if is_to_suppress_acc:
                self.kalman[ i ].y_innov_modulate = 0.;
                print 'Suppressed acc';
            else:
                self.kalman[ i ].y_innov_modulate = 1.;
            #
            res = self.kalman[ i ] \
                  ( dg[ i ], g_acc[ i ], self.dt );
            if not isnan( res ):
                self.g[ i ] = res;
            #
        #

        self.g = self.normalize( self.g );
        return;

    def test_adapt_does_nothing(self):
        dimension = 3
        ps = rand(dimension)
        distribution = Bernoulli(ps)
        kernel = HypercubeKernel(1.)
        Z = zeros((2, distribution.dimension))
        threshold = 0.5
        spread = .5
        sampler = DiscreteKameleon(distribution, kernel, Z, threshold, spread)
        
        # serialise, call adapt, load, compare
        f = NamedTemporaryFile()
        dump(sampler, f)
        f.seek(0)
        sampler_copy = load(f)
        f.close()
        
        sampler.adapt(None, None)
        
        # rough check for equality, dont do a proper one here
        self.assertEqual(type(sampler_copy.kernel), type(sampler.kernel))
        self.assertEqual(sampler_copy.kernel.gamma, sampler.kernel.gamma)
        
        self.assertEqual(type(sampler_copy.distribution), type(sampler.distribution))
        self.assertEqual(sampler_copy.distribution.dimension, sampler.distribution.dimension)
        
        self.assertEqual(type(sampler_copy.Z), type(sampler.Z))
        self.assertEqual(sampler_copy.Z.shape, sampler.Z.shape)
        self.assertAlmostEqual(norm(sampler_copy.Z - sampler.Z), 0)
        
        self.assertEqual(sampler_copy.spread, sampler.spread)

        # this is none, so just compare
        self.assertEqual(sampler.Q, sampler_copy.Q)

    def v( self, xx):
        """
        Provides the flow field as a function of the location. This is
        implemented here only for the component in the direction of :attr:`flow`;
        entrainment components are set to zero.

        Parameters
        ----------
        xx : array of floats of shape (3, )
            Location in the fluid for which to provide the data.

        Returns
        -------
        tuple with two elements
            The first element in the tuple is the velocity vector and the
            second is the Jacobian of the velocity vector field, both at the
            given location.
        """
        # TODO: better to make sure that self.flow and self.plane are indeed unit vectors before
        # normalize
        flow = self.flow/norm(self.flow)
        plane = self.plane/norm(self.plane)
        # additional axes of global co-ordinate system
        yy = -cross(flow,plane)
        zz = cross(flow,yy)
        # distance from slot exit plane
        xx1 = xx-self.origin
        # local co-ordinate system 
        x = dot(flow,xx1)
        y = dot(yy,xx1)
        x1 = 0.109*x
        h1 = abs(y)+sqrt(pi)*0.5*x1-0.5*self.B
        if h1 < 0.0:
            # core jet
            Ux = self.v0
            Udx = 0
            Udy = 0
        else:
            # shear layer
            Ux = self.v0*exp(-h1*h1/(2*x1*x1))
            Udx = (h1*h1/(x*x1*x1)-sqrt(pi)*0.5*h1/(x*x1))*Ux
            Udy = -sign(y)*h1*Ux/(x1*x1)
        # Jacobi matrix
        dU = array(((Udx,0,0),(Udy,0,0),(0,0,0))).T
        # rotation matrix
        R = array((flow,yy,zz)).T
        return dot(R,array((Ux,0,0))), dot(dot(R,dU),R.T)

	def move_axis( self , dist , axis ) :
		''' translate cursor in self.node space '''

		axis = np.array( axis )
		axis = axis / la.norm( axis )
		axis*= dist

		return self.move_vec( axis )

def dibujarpf(bdcha,cadena,pasos=1,ppo=0, fin = np.inf,multi = False):
    '''dibuja a partir de datos recogidos en arrays de datos tipo barco y cadena
    ver el sistema de preparar matrices para guardar datos dibuja también
    un buque fondeado buque = [eslora, manga, posición, orientación,de     sup'''
    
    
    
    if fin > bdcha.shape[0]-1:
        fin = bdcha.shape[0]-1
        
    for i in range(ppo,fin,pasos):
        if multi:
            figure()
        
        cms = np.array([cadena[i,0,:],cadena[i,1,:]])
        para = np.array([cadena[i,-2,:],cadena[i,-1,:]])
        pl.plot(cms[0,:],cms[1,:],'bo')
        pl.hold(True)
        barrasi = cms + cadena[-1,1,0] * para
        barrasd = cms - cadena[-1,1,0] * para
        pl.plot([barrasi[0,:],barrasd[0,:]],[barrasi[1,:],barrasd[1,:]],'k')
                        

        
        

        vertices = np.array([[-bdcha[-1,6]/2.,-0.25*bdcha[-1,6]/2],\
        [-bdcha[-1,6]/2.,0.25*bdcha[-1,6]/2],\
        [-0.25*bdcha[-1,6]/2,0.35*bdcha[-1,6]/2],[bdcha[-1,6]/2.,0],\
        [-0.25*bdcha[-1,6]/2,-0.35*bdcha[-1,6]/2],[-bdcha[-1,6]/2.,\
        -0.25*bdcha[-1,6]/2]])        
        rot = np.array([[np.cos(bdcha[i,6]),- np.sin(bdcha[i,6])],[np.sin(bdcha[i,6]),\
        np.cos(bdcha[i,6])]])       
        vertrot = np.array([np.dot(rot,j) for j in vertices]) + [bdcha[i,0],bdcha[i,1]]
        codes = [Path.MOVETO,Path.LINETO,Path.CURVE3,Path.CURVE3,Path.CURVE3,\
        Path.CURVE3]
     
        pathd = Path(vertrot,codes)
        patchd = patches.PathPatch(pathd,facecolor = 'green') #'green'
        pl.gca().add_patch(patchd)
         #######################dibujar cable de arrastre derecha#######################
        rot = np.array([[np.cos(bdcha[i,6]),- np.sin(bdcha[i,6])],[np.sin(bdcha[i,6]),\
        np.cos(bdcha[i,6])]])  
        popad =  np.dot(rot, np.array([-bdcha[-1,6]/2.,0])) + [bdcha[i,0],bdcha[i,1]]  
        tipd = - para[:,-1] * cadena[-1,1,0] + cms[:,-1]
        
        
        distd = norm(popad - tipd)
        dd = distd/cadena[-1,3,1]
        #print dd
        if dd > 1: dd = 1
        r = bezier_cvr.bezier4p([[tipd[0]],[tipd[1]]],[[popad[0]],[popad[1]]],1,1,1.5,\
        (1-dd) * bdcha[i,6]\
         +dd * np.arctan2(popad[1]-tipd[1],popad[0] - tipd[0]),\
        (1-dd) * np.arctan2(-para[0,0],-para[0,1])\
         +dd * np.arctan2(popad[1]-tipd[1],popad[0] - tipd[0]),\
        100)
        bezier_cvr.pintar_bezier(r[0],color = 'b')

 def test_log_pdf_1d_2n(self):
     mu = asarray([0], dtype=numpy.bool8)
     spread = rand()
     dist = DiscreteRandomWalkProposal(mu, spread)
     X = asarray([[1], [0]], dtype=numpy.bool8)
     log_liks = dist.log_pdf(X)
     expected = asarray([log(1.), -inf])
     self.assertAlmostEqual(norm(log_liks[0] - expected[0]), 0)
     self.assertEqual(log_liks[1], expected[1])

 def test_kernel_X_two_points_fixed(self):
     gamma = .2
     k = HypercubeKernel(gamma)
     X = asarray([[1, 0], [1, 1]], dtype=numpy.bool8)
     K = zeros((2, 2))
     for i in range(2):
         for j in range(2):
             dist = sum(X[i] != X[j])
             K[i, j] = tanh(gamma) ** dist
     self.assertAlmostEqual(norm(K - k.kernel(X)), 0)

def testFitNorm():
    X = coo_matrix((ones(4),([0, 1, 2, 2], [1, 1, 0, 1])), shape=(3, 3), dtype=np.uint8).tolil()
    A = np.array([[0.9, 0.1],
         [0.8, 0.2],
         [0.1, 0.9]])
    R = np.array([[0.9, 0.1],
         [0.1, 0.9]])
    expectedNorm = norm(X - dot(A,dot(R, A.T)))**2
    assert_almost_equal(fitNorm(X, A, R), expectedNorm)  
    assert_almost_equal(fitNormWithoutNormX(X, A, R) + squareFrobeniusNormOfSparse(X), expectedNorm)

def varify_gradient_decent(lambda_num):
    input_layer_size = 3
    hidden_layer_size = 5
    num_labels = 3
    m = 5
    input = generate_debug_input(m, input_layer_size - 1)
    print(input)
    output = 1 + np.mod(np.arange(m).reshape(1, m), num_labels).T
    print(output)
    theta1 = rand_initialize_weights(input_layer_size, hidden_layer_size)
    theta2 = rand_initialize_weights(hidden_layer_size, num_labels)
    nn_param = unroll_params(theta1, theta2)

    cost_func = lambda param: cost_function(param, input_layer_size, hidden_layer_size, input, output, lambda_num,
                                            num_labels)
    (cost, grad) = cost_func(nn_param)
    num_grad = compute_numerical_gradient(cost_func, nn_param)
    diff = linalg.norm(num_grad - grad, 2) / linalg.norm(num_grad + grad,2)
    print("The relative difference will be small (less than 1e-9)", diff)

 def test_log_lik_vector_multiple2(self):
     n=100
     y=randint(0,2,n)*2-1
     F=randn(10,n)
     
     lik=LogitLikelihood()
     multiples=lik.log_lik_vector_multiple(y, F)
     singles=asarray([lik.log_lik_vector(y, f) for f in F])
     
     self.assertLessEqual(norm(singles-multiples), 1e-10)

 def test_log_lik_vector_multiple1(self):
     n=100
     y=randint(0,2,n)*2-1
     f=randn(n)
     
     lik=LogitLikelihood()
     multiple=lik.log_lik_vector_multiple(y, f.reshape(1,n))
     single=lik.log_lik_vector(y, f)
     
     self.assertLessEqual(norm(single-multiple), 1e-10)

 def rknewton(self,f,t,x0,tol=1e-5, nmax=50, nmax_gss=100):
     F = lambda x,x0:  Matrix(lsolve(-self.J(t,x),f(x0)))
     #F = lambda x,x0:  Matrix(lsolve(-J(x),f(x0)))
     #F = lambda x,x0: -inv(J(x))*f(x0)
     if norm(f(x0),1) == 0:
         return [x0,f(x0), 0]
     else:
         for n in range(1,nmax+1):
             x = self.RK.RKX(lambda t,Y:F(Y,x0), 0, x0,1,1)[:,1]
             if norm(f(x),2) > norm(f(x0),2):
                 s = x - x0
                 f2 = lambda Y: (f(Y).T*f(Y))[0]
                 f2_= lambda alpha:f2(x0 + alpha*s)
                 alpha = gnewton(f2_,0,min(1.0*tol/norm(s,2),1e-2),nmax_gss)[0]
                 x = x0 + alpha*s
             if norm(x-x0, 1)<tol:
                 break
             else:
                 x0 = x
     return [x,f(x), n]

def create_matrix(m, k, p, N, iteration_number, starting_conditions):
    x = [random.choice([-1, 1]) for _ in range(0, N)]
    A = zeros(shape=(N, N))
    for i in range(0, N):
        for j in range(0, N):
            if i == j:
                A[i][j] = k
            elif j > i:
                A[i][j] = (-1) ** (j + 1) * (m / (j + 1))
            elif j == (i - 1):
                A[i][j] = m / (i + 1)
    x_copy = x
    b = dot(A, x)
    D = diag(A)
    R = A - diagflat(D)
    x = starting_conditions
    x_norm = p + 1
    i = 0
    B = R/D
    e_vals, e_vect = linalg.eig(B)
    print(";".join((str(N), str(max(abs(e_vals))))))
    # print "results for ||x(i+1) - x(i): "
    while (x_norm >= p) or (i > iteration_number):
        prev_x = x
        x = (b - dot(R, x)) / D
        x_norm = linalg.norm(x - prev_x, inf)  # norma typu max po kolumnach
        i += 1

    # print ";".join((str(N), str("%.8f" % p), str(i), str("%.15f" % linalg.norm(x_copy - x)), str("%.15f" % linalg.norm(x_copy - x, inf)))) + ";",
    x = x_copy
    b = dot(A, x)
    D = diag(A)
    R = A - diagflat(D)
    x = starting_conditions
    b_norm = p + 1
    i = 0
    # print "results for ||Ax(i) -b ||"
    while (b_norm >= p) or (i > iteration_number):
        x = (b - dot(R, x)) / D
        b_norm = linalg.norm(dot(A, x) - b, inf)
        i += 1