def test_pseudo_inverse():
    "Tests of the pseudo-inverse"

    # Numerical version of the pseudo-inverse:
    pinv_num = core.wrap_array_func(numpy.linalg.pinv)

    ##########
    # Full rank rectangular matrix:
    m = unumpy.matrix([[ufloat((10, 1)), -3.1],
                       [0, ufloat((3, 0))],
                       [1, -3.1]])

    # Numerical and package (analytical) pseudo-inverses: they must be
    # the same:
    rcond = 1e-8  # Test of the second argument to pinv()
    m_pinv_num = pinv_num(m, rcond)
    m_pinv_package = core._pinv(m, rcond)
    assert matrices_close(m_pinv_num, m_pinv_package)

    ##########
    # Example with a non-full rank rectangular matrix:
    vector = [ufloat((10, 1)), -3.1, 11]
    m = unumpy.matrix([vector, vector])
    m_pinv_num = pinv_num(m, rcond)
    m_pinv_package = core._pinv(m, rcond)
    assert matrices_close(m_pinv_num, m_pinv_package)

    ##########
    # Example with a non-full-rank square matrix:
    m = unumpy.matrix([[ufloat((10, 1)), 0], [3, 0]])
    m_pinv_num = pinv_num(m, rcond)
    m_pinv_package = core._pinv(m, rcond)
    assert matrices_close(m_pinv_num, m_pinv_package)

def test_list_pseudo_inverse():
    "Test of the pseudo-inverse"

    x = ufloat((1, 0.1))
    y = ufloat((2, 0.1))
    mat = unumpy.matrix([[x, x], [y, 0]])

    # Internal consistency: the inverse and the pseudo-inverse yield
    # the same result on square matrices:
    assert matrices_close(mat.I, unumpy.ulinalg.pinv(mat), 1e-4)
    assert matrices_close(
        unumpy.ulinalg.inv(mat),
        # Support for the optional pinv argument is
        # tested:
        unumpy.ulinalg.pinv(mat, 1e-15),
        1e-4,
    )

    # Non-square matrices:
    x = ufloat((1, 0.1))
    y = ufloat((2, 0.1))
    mat1 = unumpy.matrix([[x, y]])  # "Long" matrix
    mat2 = unumpy.matrix([[x, y], [1, 3 + x], [y, 2 * x]])  # "Tall" matrix

    # Internal consistency:
    assert matrices_close(mat1.I, unumpy.ulinalg.pinv(mat1, 1e-10))
    assert matrices_close(mat2.I, unumpy.ulinalg.pinv(mat2, 1e-8))

def inertia_components(jay, beta):
    '''Returns the 2D orthogonal inertia tensor.

    When at least three moments of inertia and their axes orientations are
    known relative to a common inertial frame of a planar object, the orthoganl
    moments of inertia relative the frame are computed.

    Parameters
    ----------
    jay : ndarray, shape(n,)
        An array of at least three moments of inertia. (n >= 3)
    beta : ndarray, shape(n,)
        An array of orientation angles corresponding to the moments of inertia
        in jay.

    Returns
    -------
    eye : ndarray, shape(3,)
        Ixx, Ixz, Izz

    '''
    sb = unumpy.sin(beta)
    cb = unumpy.cos(beta)
    betaMat = unumpy.matrix(np.vstack((cb**2, -2 * sb * cb, sb**2)).T)
    eye = np.squeeze(np.asarray(np.dot(betaMat.I, jay)))
    return eye

def test_wrap_array_func():
    '''
    Test of numpy.wrap_array_func(), with optional arguments and
    keyword arguments.
    '''

    # Function that works with numbers with uncertainties in mat (if
    # mat is an uncertainties.unumpy.matrix):
    def f_unc(mat, *args, **kwargs):
        return mat.I + args[0]*kwargs['factor']

    # Test with optional arguments and keyword arguments:
    def f(mat, *args, **kwargs):
        # This function is wrapped: it should only be called with pure
        # numbers:
        assert not any(isinstance(v, uncert_core.UFloat) for v in mat.flat)
        return f_unc(mat, *args, **kwargs)


    # Wrapped function:
    f_wrapped = core.wrap_array_func(f)

    ##########
    # Full rank rectangular matrix:
    m = unumpy.matrix([[ufloat(10, 1), -3.1],
                       [0, ufloat(3, 0)],
                       [1, -3.1]])

    # Numerical and package (analytical) pseudo-inverses: they must be
    # the same:
    m_f_wrapped = f_wrapped(m, 2, factor=10)
    m_f_unc = f_unc(m, 2, factor=10)

    assert arrays_close(m_f_wrapped, m_f_unc)

def test_component_extraction():
    "Extracting the nominal values and standard deviations from an array"

    arr = unumpy.uarray(([1, 2], [0.1, 0.2]))

    assert numpy.all(unumpy.nominal_values(arr) == [1, 2])
    assert numpy.all(unumpy.std_devs(arr) == [0.1, 0.2])

    # unumpy matrices, in addition, should have nominal_values that
    # are simply numpy matrices (not unumpy ones, because they have no
    # uncertainties):
    mat = unumpy.matrix(arr)
    assert numpy.all(unumpy.nominal_values(mat) == [1, 2])
    assert numpy.all(unumpy.std_devs(mat) == [0.1, 0.2])
    assert type(unumpy.nominal_values(mat)) == numpy.matrix

def test_matrix():
    "Matrices of numbers with uncertainties"
    # Matrix inversion:

    # Matrix with a mix of Variable objects and regular
    # Python numbers:

    m = unumpy.matrix([[ufloat((10, 1)), -3.1],
                       [0, ufloat((3, 0))]])
    m_nominal_values = unumpy.nominal_values(m)

    # Test of the nominal_value attribute:
    assert numpy.all(m_nominal_values == m.nominal_values)

    assert type(m[0, 0]) == uncertainties.Variable

def test_matrix():
    "Matrices of numbers with uncertainties"
    # Matrix inversion:

    # Matrix with a mix of Variable objects and regular
    # Python numbers:

    m = unumpy.matrix([[ufloat(10, 1), -3.1],
                       [0, ufloat(3, 0)]])
    m_nominal_values = unumpy.nominal_values(m)

    # Test of the nominal_value attribute:
    assert numpy.all(m_nominal_values == m.nominal_values)

    assert type(m[0, 0]) == uncert_core.Variable

    # Test of scalar multiplication, both sides:
    3*m
    m*3

def center_of_mass(slopes, intercepts):
    '''Returns the center of mass relative to the slopes and intercepts
    coordinate system.

    Parameters
    ----------
    slopes : ndarray, shape(n,)
        The slope of every line used to calculate the center of mass.
    intercepts : ndarray, shape(n,)
        The intercept of every line used to calculate the center of mass.

    Returns
    -------
    x : float
        The abscissa of the center of mass.
    y : float
        The ordinate of the center of mass.

    '''
    num = range(len(slopes))
    allComb = cartesian((num, num))
    comb = []
    # remove doubles
    for row in allComb:
        if row[0] != row[1]:
            comb.append(row)
    comb = np.array(comb)

    # initialize the matrix to store the line intersections
    lineX = np.zeros((len(comb), 2), dtype='object')
    # for each line intersection...
    for j, row in enumerate(comb):
        sl = np.array([slopes[row[0]], slopes[row[1]]])
        a = unumpy.matrix(np.vstack((-sl, np.ones((2)))).T)
        b = np.array([intercepts[row[0]], intercepts[row[1]]])
        lineX[j] = np.dot(a.I, b)
    com = np.mean(lineX, axis=0)

    return com[0], com[1]

def test_list_inverse():
    "Test of the inversion of a square matrix"

    mat_list = [[1, 1], [1, 0]]

    # numpy.linalg.inv(mat_list) does calculate the inverse even
    # though mat_list is a list of lists (and not a matrix).  Can
    # ulinalg do the same?  Here is a test:
    mat_list_inv = unumpy.ulinalg.inv(mat_list)

    # More type testing:
    mat_matrix = numpy.asmatrix(mat_list)
    assert isinstance(unumpy.ulinalg.inv(mat_matrix),
                      type(numpy.linalg.inv(mat_matrix)))

    # unumpy.ulinalg should behave in the same way as numpy.linalg,
    # with respect to types:
    mat_list_inv_numpy = numpy.linalg.inv(mat_list)
    assert type(mat_list_inv) == type(mat_list_inv_numpy)

    # The resulting matrix does not have to be a matrix that can
    # handle uncertainties, because the input matrix does not have
    # uncertainties:
    assert not isinstance(mat_list_inv, unumpy.matrix)

    # Individual element check:
    assert isinstance(mat_list_inv[1,1], float)    
    assert mat_list_inv[1,1] == -1

    x = ufloat((1, 0.1))
    y = ufloat((2, 0.1))
    mat = unumpy.matrix([[x, x], [y, 0]])
    
    # Internal consistency: ulinalg.inv() must coincide with the
    # unumpy.matrix inverse, for square matrices (.I is the
    # pseudo-inverse, for non-square matrices, but inv() is not).
    assert matrices_close(unumpy.ulinalg.inv(mat), mat.I)

def test_inverse():
    "Tests of the matrix inverse"

    m = unumpy.matrix([[ufloat((10, 1)), -3.1],
                       [0, ufloat((3, 0))]])
    m_nominal_values = unumpy.nominal_values(m)

    # "Regular" inverse matrix, when uncertainties are not taken
    # into account:
    m_no_uncert_inv = m_nominal_values.I

    # The matrix inversion should not yield numbers with uncertainties:
    assert m_no_uncert_inv.dtype == numpy.dtype(float)

    # Inverse with uncertainties:
    m_inv_uncert = m.I  # AffineScalarFunc elements
    # The inverse contains uncertainties: it must support custom
    # operations on matrices with uncertainties:
    assert isinstance(m_inv_uncert, unumpy.matrix)
    assert type(m_inv_uncert[0, 0]) == uncertainties.AffineScalarFunc

    # Checks of the numerical values: the diagonal elements of the
    # inverse should be the inverses of the diagonal elements of
    # m (because we started with a triangular matrix):
    assert _numbers_close(1/m_nominal_values[0, 0],
                          m_inv_uncert[0, 0].nominal_value), "Wrong value"

    assert _numbers_close(1/m_nominal_values[1, 1],
                          m_inv_uncert[1, 1].nominal_value), "Wrong value"


    ####################

    # Checks of the covariances between elements:
    x = ufloat((10, 1))
    m = unumpy.matrix([[x, x],
                       [0, 3+2*x]])

    m_inverse = m.I

    # Check of the properties of the inverse:
    m_double_inverse = m_inverse.I
    # The initial matrix should be recovered, including its
    # derivatives, which define covariances:
    assert _numbers_close(m_double_inverse[0, 0].nominal_value,
                          m[0, 0].nominal_value)
    assert _numbers_close(m_double_inverse[0, 0].std_dev(),
                          m[0, 0].std_dev())

    assert matrices_close(m_double_inverse, m)

    # Partial test:
    assert _derivatives_close(m_double_inverse[0, 0], m[0, 0])
    assert _derivatives_close(m_double_inverse[1, 1], m[1, 1])

    ####################

    # Tests of covariances during the inversion:

    # There are correlations if both the next two derivatives are
    # not zero:
    assert m_inverse[0, 0].derivatives[x]
    assert m_inverse[0, 1].derivatives[x]

    # Correlations between m and m_inverse should create a perfect
    # inversion:
    assert matrices_close(m * m_inverse,  numpy.eye(m.shape[0]))

U_CU_C, U_CU_H, U_CU_M, U_AL_C, U_AL_H, U_AL_M = np.loadtxt("Messdaten/Messung_Material.txt",
                                                            unpack=True)

# Erstellen der Fehlerbehaftete Spannungen
uU_CU_C = unp.uarray(U_CU_C, len(U_CU_C)*[U_ERR])
uU_CU_H = unp.uarray(U_CU_H, len(U_CU_H)*[U_ERR])
uU_CU_M = unp.uarray(U_CU_M, len(U_CU_M)*[U_ERR])

uU_AL_C = unp.uarray(U_AL_C, len(U_AL_C)*[U_ERR])
uU_AL_H = unp.uarray(U_AL_H, len(U_AL_H)*[U_ERR])
uU_AL_M = unp.uarray(U_AL_M, len(U_AL_M)*[U_ERR])

# Erstellen von 3x3 Matrizen für die Werte von Cu und Al
# Spalte entspricht einer Versuchsreihe, Zeile: Entspricht einer Größe

uU_CU = unp.matrix([uU_CU_C, uU_CU_H, uU_CU_M])
uU_AL = unp.matrix([uU_AL_C, uU_AL_H, uU_AL_M])


# Umrechnung der Spannungen in Temperaturen

uT_CU = TensToTemp(uU_CU)
uT_AL = TensToTemp(uU_AL)

TEST = True
if TEST:
    for i in range(3):
        for j in range(3):
            uT_CU[i,j] = ufloat(noms(uT_CU[i,j]), stds(uT_CU[i,j]))

if TEST:

d6 = genfromtxt(f6, delimiter=',')

ar8 = unumpy.uarray(d6[:,0],d6[:,1])
ar08 = unumpy.uarray(d1[:,0],d1[:,1])
ar008 = unumpy.uarray(d2[:,0],d2[:,1])
arabs = unumpy.uarray(d3[:,0],d3[:,1])
aragua = unumpy.uarray(d4[:,0],d4[:,1])
arVBC = unumpy.uarray(d5[:,0],d4[:,1])


vol = (ar8 + ar08 + ar008 + aragua + arVBC)
mol = (ar8 * 0.8 + ar08 *0.08 + ar008 *0.008)

conc = mol/vol

mT = unumpy.matrix(conc)
filename = "resconc.csv"
numpy.savetxt(filename, mT, fmt='%r', delimiter='\n')

XB = arabs/(0.757)

mT = unumpy.matrix(XB)
filename = "resXB.csv"
numpy.savetxt(filename, mT, fmt='%r', delimiter='\n')

XHB = 1 - XB

mT = unumpy.matrix(XHB)
filename = "resXHB.csv"
numpy.savetxt(filename, mT, fmt='%r', delimiter='\n')