Python odl.FunctionSpace類代碼示例

def test_addition():
    fspace = FunctionSpace(odl.Interval(0, pi))
    sine = fspace.element(np.sin)
    cosine = fspace.element(np.cos)

    sum_func = sine + cosine

    for x in [0.0, 0.2, 1.0]:
        assert almost_equal(sum_func(x), np.sin(x) + np.cos(x))

def test_fspace_vector_eval_complex():
    rect, points, mg = _standard_setup_2d()

    fspace = FunctionSpace(rect, field=odl.ComplexNumbers())
    f_novec = fspace.element(cfunc_2d_novec, vectorized=False)
    f_vec_oop = fspace.element(cfunc_2d_vec_oop, vectorized=True)
    f_vec_ip = fspace.element(cfunc_2d_vec_ip, vectorized=True)
    f_vec_dual = fspace.element(cfunc_2d_vec_dual, vectorized=True)

    true_arr = cfunc_2d_vec_oop(points)
    true_mg = cfunc_2d_vec_oop(mg)

    # Out-of-place
    assert f_novec([0.5, 1.5]) == cfunc_2d_novec([0.5, 1.5])
    assert f_vec_oop([0.5, 1.5]) == cfunc_2d_novec([0.5, 1.5])
    assert all_equal(f_vec_oop(points), true_arr)
    assert all_equal(f_vec_oop(mg), true_mg)

    # In-place standard implementation
    out_arr = np.empty((5,), dtype='complex128')
    out_mg = np.empty((2, 3), dtype='complex128')

    f_vec_oop(points, out=out_arr)
    f_vec_oop(mg, out=out_mg)
    assert all_equal(out_arr, true_arr)
    assert all_equal(out_mg, true_mg)

    with pytest.raises(TypeError):  # ValueError: invalid vectorized input
        f_vec_oop(points[0])

    # In-place only
    out_arr = np.empty((5,), dtype='complex128')
    out_mg = np.empty((2, 3), dtype='complex128')

    f_vec_ip(points, out=out_arr)
    f_vec_ip(mg, out=out_mg)
    assert all_equal(out_arr, true_arr)
    assert all_equal(out_mg, true_mg)

    # Standard out-of-place evaluation
    assert f_vec_ip([0.5, 1.5]) == cfunc_2d_novec([0.5, 1.5])
    assert all_equal(f_vec_ip(points), true_arr)
    assert all_equal(f_vec_ip(mg), true_mg)

    # Dual use
    assert f_vec_dual([0.5, 1.5]) == cfunc_2d_novec([0.5, 1.5])
    assert all_equal(f_vec_dual(points), true_arr)
    assert all_equal(f_vec_dual(mg), true_mg)

    out_arr = np.empty((5,), dtype='complex128')
    out_mg = np.empty((2, 3), dtype='complex128')

    f_vec_dual(points, out=out_arr)
    f_vec_dual(mg, out=out_mg)
    assert all_equal(out_arr, true_arr)
    assert all_equal(out_mg, true_mg)

def test_fspace_vector_ufunc():
    intv = odl.Interval(0, 1)
    points = _points(intv, num=5)
    mg = _meshgrid(intv, shape=(5,))

    fspace = FunctionSpace(intv)
    f_vec = fspace.element(np.sin)

    assert f_vec(0.5) == np.sin(0.5)
    assert all_equal(f_vec(points), np.sin(points.squeeze()))
    assert all_equal(f_vec(mg), np.sin(mg[0]))

def test_interval(exponent):
    fspace = FunctionSpace(odl.Interval(0, pi))
    lpdiscr = odl.uniform_discr(fspace, 10, exponent=exponent)

    if exponent == float("inf"):
        sine = fspace.element(np.sin)
        discr_sine = lpdiscr.element(sine)
        assert discr_sine.norm() <= 1
    else:
        sine_p = fspace.element(lambda x: np.sin(x) ** (1 / exponent))
        discr_sine_p = lpdiscr.element(sine_p)
        assert almost_equal(discr_sine_p.norm(), 2 ** (1 / exponent), places=2)

def test_rectangle(exponent):
    fspace = FunctionSpace(odl.Rectangle((0, 0), (pi, 2 * pi)))
    n, m = 10, 10
    lpdiscr = odl.uniform_discr(fspace, (n, m), exponent=exponent)

    if exponent == float("inf"):
        sine2 = fspace.element(lambda x, y: np.sin(x) * np.sin(y))
        discr_sine = lpdiscr.element(sine2)
        assert discr_sine.norm() <= 1
    else:
        sine_p = fspace.element(lambda x, y: (np.sin(x) * np.sin(y)) ** (1 / exponent))
        discr_sine_p = lpdiscr.element(sine_p)
        assert almost_equal(discr_sine_p.norm(), 4 ** (1 / exponent), places=2)

def test_fspace_vector_equality():
    rect = odl.Rectangle([0, 0], [1, 2])
    fspace = FunctionSpace(rect)

    f_novec = fspace.element(func_2d_novec, vectorized=False)

    f_vec_oop = fspace.element(func_2d_vec_oop, vectorized=True)
    f_vec_oop_2 = fspace.element(func_2d_vec_oop, vectorized=True)

    f_vec_ip = fspace.element(func_2d_vec_ip, vectorized=True)
    f_vec_ip_2 = fspace.element(func_2d_vec_ip, vectorized=True)

    f_vec_dual = fspace.element(func_2d_vec_dual, vectorized=True)
    f_vec_dual_2 = fspace.element(func_2d_vec_dual, vectorized=True)

    assert f_novec == f_novec
    assert f_novec != f_vec_oop
    assert f_novec != f_vec_ip
    assert f_novec != f_vec_dual

    assert f_vec_oop == f_vec_oop
    assert f_vec_oop == f_vec_oop_2
    assert f_vec_oop != f_vec_ip
    assert f_vec_oop != f_vec_dual

    assert f_vec_ip == f_vec_ip
    assert f_vec_ip == f_vec_ip_2
    assert f_vec_ip != f_vec_dual

    assert f_vec_dual == f_vec_dual
    assert f_vec_dual == f_vec_dual_2

def test_fspace_vector_copy():
    fspace = FunctionSpace(odl.Interval(0, 1))

    f_novec = fspace.element(func_1d_oop, vectorized=False)
    f_vec_ip = fspace.element(func_1d_ip, vectorized=True)
    f_vec_dual = fspace.element(func_1d_dual, vectorized=True)

    f_out = f_novec.copy()
    assert f_out == f_novec

    f_out = f_vec_ip.copy()
    assert f_out == f_vec_ip

    f_out = f_vec_dual.copy()
    assert f_out == f_vec_dual

def test_fspace_out_dtype():
    rect = odl.Rectangle([0, 0], [3, 5])
    points = np.array([[0, 1], [0, 3], [3, 4], [2, 5]], dtype='int').T
    vec1 = np.array([0, 1, 3])[:, None]
    vec2 = np.array([1, 2, 4, 5])[None, :]
    mg = (vec1, vec2)

    true_arr = func_2d_vec_oop(points)
    true_mg = func_2d_vec_oop(mg)

    fspace = FunctionSpace(rect, out_dtype='int')
    f_vec = fspace.element(func_2d_vec_oop)

    assert all_equal(f_vec(points), true_arr)
    assert all_equal(f_vec(mg), true_mg)
    assert f_vec(points).dtype == np.dtype('int')
    assert f_vec(mg).dtype == np.dtype('int')

def test_fspace_one():
    rect, points, mg = _standard_setup_2d()

    # real
    fspace = FunctionSpace(rect)
    one_vec = fspace.one()

    assert one_vec([0.5, 1.5]) == 1.0
    assert all_equal(one_vec(points), np.ones(5, dtype=float))
    assert all_equal(one_vec(mg), np.ones((2, 3), dtype=float))

    # complex
    fspace = FunctionSpace(rect, field=odl.ComplexNumbers())
    one_vec = fspace.one()

    assert one_vec([0.5, 1.5]) == 1.0 + 1j * 0.0
    assert all_equal(one_vec(points), np.ones(5, dtype=complex))
    assert all_equal(one_vec(mg), np.ones((2, 3), dtype=complex))

def test_fspace_vector_arithmetic(variant, op):
    if variant == 'sv' and '=' in op:  # makes no sense, quit
        return

    # Setup
    rect, points, mg = _standard_setup_2d()
    point = points.T[0]

    fspace = FunctionSpace(rect)
    a = -1.5
    b = 2.0
    array_out = np.empty((5,), dtype=float)
    mg_out = np.empty((2, 3), dtype=float)

    # Initialize a bunch of elements
    f_novec = fspace.element(func_2d_novec, vectorized=False)
    f_vec = fspace.element(func_2d_vec_dual, vectorized=True)
    g_novec = fspace.element(other_func_2d_novec, vectorized=False)
    g_vec = fspace.element(other_func_2d_vec_dual, vectorized=True)

    out_novec = fspace.element(vectorized=False)
    out_vec = fspace.element(vectorized=True)

    if variant[0] == 'v':
        true_l_novec = func_2d_novec(point)
        true_l_arr = func_2d_vec_oop(points)
        true_l_mg = func_2d_vec_oop(mg)

        test_l_novec = f_novec
        test_l_vec = f_vec
    else:  # 's'
        true_l_novec = true_l_arr = true_l_mg = a
        test_l_novec = test_l_vec = a

    if variant[1] == 'v':
        true_r_novec = other_func_2d_novec(point)
        true_r_arr = other_func_2d_vec_oop(points)
        true_r_mg = other_func_2d_vec_oop(mg)

        test_r_novec = g_novec
        test_r_vec = g_vec
    else:  # 's'
        true_r_novec = true_r_arr = true_r_mg = b
        test_r_novec = test_r_vec = b

    true_novec = _op(true_l_novec, op, true_r_novec)
    true_arr = _op(true_l_arr, op, true_r_arr)
    true_mg = _op(true_l_mg, op, true_r_mg)

    out_novec = _op(test_l_novec, op, test_r_novec)
    out_vec = _op(test_l_vec, op, test_r_vec)

    assert almost_equal(out_novec(point), true_novec)
    assert all_equal(out_vec(points), true_arr)
    out_vec(points, out=array_out)
    assert all_equal(array_out, true_arr)
    assert all_equal(out_vec(mg), true_mg)
    out_vec(mg, out=mg_out)
    assert all_equal(mg_out, true_mg)

def test_fspace_power(power):
    rect, points, mg = _standard_setup_2d()
    point = points.T[0]
    out_arr = np.empty(5)
    out_mg = np.empty((2, 3))

    fspace = FunctionSpace(rect)

    # Not vectorized
    true_novec = func_2d_novec(point) ** power

    f_novec = fspace.element(func_2d_novec, vectorized=False)
    pow_novec = f_novec ** power
    assert almost_equal(pow_novec(point), true_novec)

    pow_novec = f_novec.copy()
    pow_novec **= power

    assert almost_equal(pow_novec(point), true_novec)

    # Vectorized
    true_arr = func_2d_vec_oop(points) ** power
    true_mg = func_2d_vec_oop(mg) ** power

    f_vec = fspace.element(func_2d_vec_dual, vectorized=True)
    pow_vec = f_vec ** power

    assert all_almost_equal(pow_vec(points), true_arr)
    assert all_almost_equal(pow_vec(mg), true_mg)

    pow_vec = f_vec.copy()
    pow_vec **= power

    assert all_almost_equal(pow_vec(points), true_arr)
    assert all_almost_equal(pow_vec(mg), true_mg)

    pow_vec(points, out=out_arr)
    pow_vec(mg, out=out_mg)

    assert all_almost_equal(out_arr, true_arr)
    assert all_almost_equal(out_mg, true_mg)

def test_fspace_vector_real_imag():
    rect, _, mg = _standard_setup_2d()
    cspace = FunctionSpace(rect, field=odl.ComplexNumbers())
    f = cspace.element(cfunc_2d_vec_oop)

    # real / imag on complex functions
    assert all_equal(f.real(mg), cfunc_2d_vec_oop(mg).real)
    assert all_equal(f.imag(mg), cfunc_2d_vec_oop(mg).imag)
    out_mg = np.empty((2, 3))
    f.real(mg, out=out_mg)
    assert all_equal(out_mg, cfunc_2d_vec_oop(mg).real)
    f.imag(mg, out=out_mg)
    assert all_equal(out_mg, cfunc_2d_vec_oop(mg).imag)

    # real / imag on real functions, should be the function itself / zero
    rspace = FunctionSpace(rect)
    f = rspace.element(func_2d_vec_oop)
    assert all_equal(f.real(mg), f(mg))
    assert all_equal(f.imag(mg), rspace.zero()(mg))

    # Complex conjugate
    f = cspace.element(cfunc_2d_vec_oop)
    fbar = f.conj()
    assert all_equal(fbar(mg), cfunc_2d_vec_oop(mg).conj())
    out_mg = np.empty((2, 3), dtype='complex128')
    fbar(mg, out=out_mg)
    assert all_equal(out_mg, cfunc_2d_vec_oop(mg).conj())

def test_fspace_astype():

    rspace = FunctionSpace(odl.Interval(0, 1))
    cspace = FunctionSpace(odl.Interval(0, 1), field=odl.ComplexNumbers())
    rspace_s = FunctionSpace(odl.Interval(0, 1), out_dtype='float32')
    cspace_s = FunctionSpace(odl.Interval(0, 1), out_dtype='complex64')

    assert rspace.astype('complex64') == cspace_s
    assert rspace.astype('complex128') == cspace
    assert rspace.astype('complex128') is rspace._complex_space
    assert rspace.astype('float32') == rspace_s
    assert rspace.astype('float64') is rspace._real_space

    assert cspace.astype('float32') == rspace_s
    assert cspace.astype('float64') == rspace
    assert cspace.astype('float64') is cspace._real_space
    assert cspace.astype('complex64') == cspace_s
    assert cspace.astype('complex128') is cspace._complex_space

def test_fspace_vector_assign():
    fspace = FunctionSpace(odl.IntervalProd(0, 1))

    f_novec = fspace.element(func_1d_oop, vectorized=False)
    f_vec_ip = fspace.element(func_1d_ip, vectorized=True)
    f_vec_dual = fspace.element(func_1d_dual, vectorized=True)

    f_out = fspace.element()
    f_out.assign(f_novec)
    assert f_out == f_novec

    f_out = fspace.element()
    f_out.assign(f_vec_ip)
    assert f_out == f_vec_ip

    f_out = fspace.element()
    f_out.assign(f_vec_dual)
    assert f_out == f_vec_dual

def test_fspace_lincomb(a, b):
    rect, points, mg = _standard_setup_2d()
    point = points.T[0]

    fspace = FunctionSpace(rect)

    # Note: Special cases and alignment are tested later in the magic methods

    # Not vectorized
    true_novec = a * func_2d_novec(point) + b * other_func_2d_novec(point)
    f_novec = fspace.element(func_2d_novec, vectorized=False)
    g_novec = fspace.element(other_func_2d_novec, vectorized=False)
    out_novec = fspace.element(vectorized=False)
    fspace.lincomb(a, f_novec, b, g_novec, out_novec)
    assert almost_equal(out_novec(point), true_novec)

    # Vectorized
    true_arr = (a * func_2d_vec_oop(points) +
                b * other_func_2d_vec_oop(points))
    true_mg = (a * func_2d_vec_oop(mg) + b * other_func_2d_vec_oop(mg))

    # Out-of-place
    f_vec_oop = fspace.element(func_2d_vec_oop, vectorized=True)
    g_vec_oop = fspace.element(other_func_2d_vec_oop, vectorized=True)
    out_vec = fspace.element()
    fspace.lincomb(a, f_vec_oop, b, g_vec_oop, out_vec)

    assert all_equal(out_vec(points), true_arr)
    assert all_equal(out_vec(mg), true_mg)
    assert almost_equal(out_vec(point), true_novec)
    out_arr = np.empty((5,), dtype=float)
    out_mg = np.empty((2, 3), dtype=float)
    out_vec(points, out=out_arr)
    out_vec(mg, out=out_mg)
    assert all_equal(out_arr, true_arr)
    assert all_equal(out_mg, true_mg)

    # In-place
    f_vec_ip = fspace.element(func_2d_vec_ip, vectorized=True)
    g_vec_ip = fspace.element(other_func_2d_vec_ip, vectorized=True)
    out_vec = fspace.element()
    fspace.lincomb(a, f_vec_ip, b, g_vec_ip, out_vec)

    assert all_equal(out_vec(points), true_arr)
    assert all_equal(out_vec(mg), true_mg)
    assert almost_equal(out_vec(point), true_novec)
    out_arr = np.empty((5,), dtype=float)
    out_mg = np.empty((2, 3), dtype=float)
    out_vec(points, out=out_arr)
    out_vec(mg, out=out_mg)
    assert all_equal(out_arr, true_arr)
    assert all_equal(out_mg, true_mg)

    # Dual use
    f_vec_dual = fspace.element(func_2d_vec_dual, vectorized=True)
    g_vec_dual = fspace.element(other_func_2d_vec_dual, vectorized=True)
    out_vec = fspace.element()
    fspace.lincomb(a, f_vec_dual, b, g_vec_dual, out_vec)

    assert all_equal(out_vec(points), true_arr)
    assert all_equal(out_vec(mg), true_mg)
    assert almost_equal(out_vec(point), true_novec)
    out_arr = np.empty((5,), dtype=float)
    out_mg = np.empty((2, 3), dtype=float)
    out_vec(points, out=out_arr)
    out_vec(mg, out=out_mg)
    assert all_equal(out_arr, true_arr)
    assert all_equal(out_mg, true_mg)

    # Mix of vectorized and non-vectorized -> manual vectorization
    fspace.lincomb(a, f_vec_dual, b, g_novec, out_vec)
    assert all_equal(out_vec(points), true_arr)
    assert all_equal(out_vec(mg), true_mg)