Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with HTTPS or Subversion.

Download ZIP
Browse files

condassign should now work also with adubs

  • Loading branch information...
commit 7a12389d9318fce678ea2b78007ba40ed6541101 1 parent 6e7e1dd
@b45ch1 authored
View
32 adolc/src/py_adolc.cpp
@@ -195,7 +195,7 @@ void c_wrapped_hov_forward (short tape_tag, int M, int N, int D, int P, bpn::ar
for(int n = 0; n != N; ++n){
V[n] = &V1[ n * P];
}
-
+
double* W_data = (double*) nu::data(bpn_W);
double** W[M];
double* W1[M*P];
@@ -222,7 +222,7 @@ void c_wrapped_hov_wk_forward (short tape_tag, int M, int N, int D, int keep, in
for(int n = 0; n != N; ++n){
V[n] = &V1[ n * P];
}
-
+
double* W_data = (double*) nu::data(bpn_W);
double** W[M];
double* W1[M*P];
@@ -249,7 +249,7 @@ void c_wrapped_fov_reverse (short tape_tag, int M, int N, int Q, bpn::array &bp
for(int q = 0; q != Q; ++q){
U[q] = &U_data[M * q];
}
-
+
double* Z_data = (double*) nu::data(bpn_Z);
double* Z[Q];
for(int q = 0; q != Q; ++q){
@@ -280,7 +280,7 @@ void c_wrapped_hos_ti_reverse (short tape_tag, int M, int N, int D, bpn::array
}
hos_ti_reverse(tape_tag, M, N, D, U, Z);
}
-
+
void c_wrapped_hov_reverse (short tape_tag, int M, int N, int D, int Q, bpn::array &bpn_U, bpn::array &bpn_Z, bpn::array &bpn_nz){
double* U_data = (double*) nu::data(bpn_U);
@@ -314,11 +314,11 @@ void c_wrapped_hov_ti_reverse (short tape_tag, int M, int N, int D, int Q, bpn::
double* U_data = (double*) nu::data(bpn_U);
double** U[Q];
double* U1[Q*M];
-
+
for( int qn = 0; qn != Q*M; ++qn){
U1[qn] = &U_data[qn * (D+1)];
- }
-
+ }
+
for(int q = 0; q != Q; ++q){
U[q] = &U1[ q * M];
}
@@ -327,15 +327,15 @@ void c_wrapped_hov_ti_reverse (short tape_tag, int M, int N, int D, int Q, bpn::
double* Z_data = (double*) nu::data(bpn_Z);
double** Z[Q];
double* Z1[Q*N];
-
+
for( int qn = 0; qn != Q*N; ++qn){
Z1[qn] = &Z_data[qn * (D+1)];
- }
-
+ }
+
for(int q = 0; q != Q; ++q){
Z[q] = &Z1[ q * N];
}
-
+
/* nz is (Q,N) matrix */
short* nz_data = (short*) nu::data(bpn_nz);
short* nz[Q];
@@ -350,27 +350,27 @@ void c_wrapped_hov_ti_reverse (short tape_tag, int M, int N, int D, int Q, bpn::
void c_wrapped_hos_ov_reverse (short tape_tag, int M, int N, int D, int P, bpn::array &bpn_U, bpn::array &bpn_Z){
/* this function is experimental and likely not to work ... */
-
+
/* U is (M,D+1) array */
double* U_data = (double*) nu::data(bpn_U);
double* U[M];
for(int m = 0; m != M; ++m){
U[m] = &U_data[ m * (D+1)];
}
-
+
/* Z is (N, P, D+1) array???? */
double* Z_data = (double*) nu::data(bpn_Z);
double** Z[N];
double* Z1[N*P];
-
+
for( int np = 0; np != N*P; ++np){
Z1[np] = &Z_data[np * (D+1)];
}
-
+
for(int n = 0; n != N; ++n){
Z[n] = &Z1[ n * P];
}
-
+
hos_ov_reverse(tape_tag, M, N, D, P, U, Z);
}
View
32 adolc/src/py_adolc.hpp
@@ -207,13 +207,13 @@ BOOST_PYTHON_MODULE(_adolc)
using namespace boost::python;
import_array(); /* some kind of hack to get numpy working */
bpn::array::set_module_and_type("numpy", "ndarray"); /* some kind of hack to get numpy working */
-
+
scope().attr("__doc__") ="unused: moved docstring to adolc.py";
-
+
def("get_size_of_short", get_size_of_short);
def("get_size_of_int", get_size_of_int);
def("get_size_of_long", get_size_of_long);
-
+
def("trace_on",trace_on_default_argument);
def("trace_off",trace_off_default_argument);
@@ -245,7 +245,7 @@ BOOST_PYTHON_MODULE(_adolc)
def("hov_reverse", &c_wrapped_hov_reverse);
def("hov_ti_reverse", &c_wrapped_hov_ti_reverse);
def("hos_ov_reverse", &c_wrapped_hos_ov_reverse);
-
+
def("depends_on", &depends_on);
def("tape_to_latex", py_tape_doc);
@@ -254,7 +254,7 @@ BOOST_PYTHON_MODULE(_adolc)
def("exp", adub_exp_badouble, return_value_policy<manage_new_object>() );
def("log", adub_log_badouble, return_value_policy<manage_new_object>() );
def("sin", adub_sin_badouble, return_value_policy<manage_new_object>() );
- def("cos", adub_cos_badouble, return_value_policy<manage_new_object>() );
+ def("cos", adub_cos_badouble, return_value_policy<manage_new_object>() );
def("tan", adub_tan_badouble, return_value_policy<manage_new_object>() );
def("asin", adub_asin_badouble, return_value_policy<manage_new_object>() );
def("acos", adub_acos_badouble, return_value_policy<manage_new_object>() );
@@ -266,7 +266,7 @@ BOOST_PYTHON_MODULE(_adolc)
def("ceil", adub_ceil_badouble, return_value_policy<manage_new_object>() );
def("floor", adub_floor_badouble, return_value_policy<manage_new_object>() );
def("log10", adub_log10_badouble, return_value_policy<manage_new_object>() );
-
+
def("condassign", &wrapped_condassign_double_if);
def("condassign", &wrapped_condassign_double_if_else);
def("condassign", &wrapped_condassign_adouble_if);
@@ -277,7 +277,7 @@ BOOST_PYTHON_MODULE(_adolc)
.add_property("val", &badouble::value)
.add_property("loc", &badouble::loc)
-
+
.def("is_independent", &badouble::operator<<=, return_internal_reference<>())
.def("__ilshift__", operator_eq_double, return_internal_reference<>())
.def("__ilshift__", operator_eq_badouble, return_internal_reference<>())
@@ -290,34 +290,34 @@ BOOST_PYTHON_MODULE(_adolc)
.def(self -= double() )
.def(self *= double() )
.def(self /= double() )
-
+
.def(self += int() )
.def(self -= int() )
.def(self *= int() )
- .def(self /= int() )
+ .def(self /= int() )
.def(self += self )
.def(self -= self )
.def(self *= self )
.def(self /= self )
-
+
.def(self < double() )
.def(self <= double() )
.def(self > double() )
.def(self >= double() )
-
+
.def(self < int() )
.def(self <= int() )
.def(self > int() )
- .def(self >= int() )
-
-
+ .def(self >= int() )
+
+
.def(self < self )
.def(self <= self )
.def(self > self )
.def(self >= self )
-
+
// .def(-self) using this unary operator somehow screws up LATER computations, i.e. the operator works correctly, but subsequent calculations screw up!!
// .def(+self)
@@ -342,7 +342,7 @@ BOOST_PYTHON_MODULE(_adolc)
.def("__mul__", adub_mul_badouble_double, return_value_policy<manage_new_object>())
.def("__div__", adub_div_badouble_double, return_value_policy<manage_new_object>())
.def("__truediv__", adub_div_badouble_double, return_value_policy<manage_new_object>())
-
+
.def("__pow__", adub_pow_badouble_double, return_value_policy<manage_new_object>())
.def("__pow__", adouble_pow_badouble_badouble, return_value_policy<manage_new_object>())
.def("__rpow__",adouble_pow_double_badouble, return_value_policy<manage_new_object>())
View
312 adolc/tests/test_wrapped_functions.py
@@ -12,211 +12,246 @@ def test_constructors(self):
a = adouble(13.);
b = adouble(5)
c = adouble(a)
-
+
assert a.val == 13.
assert b.val == 5
assert c.val == 13.
-
+
def test_unary_operators(self):
a = adouble(1.)
b = -a
assert b.val == -1.
assert a.val == 1.
-
+
print type(b)
print type(a)
-
+
def test_conditional_operators(self):
ax = adouble(2.)
ay = adouble(1.)
-
+
assert ax <= 2
assert ax <= 2.
assert not ax < 2
assert not ax < 2.
-
-
+
+
assert ax >= 2
assert ax >= 2.
assert not ax > 2
assert not ax > 2.
-
+
assert ax > ay
assert ax >= ay
assert not ax < ay
assert not ax <= ay
-
-
+
+
def test_radd(self):
a = adouble(1.)
b = a + 2.
c = a + 2.
d = 2.+ a
-
+
assert a.val == 1.
-
+
def test_add(self):
a = adouble(1.)
b = a + 2.
c = a + 2
d = 2.+ a
e = 2 + a
-
+
assert b.val == 3.
assert c.val == 3.
assert d.val == 3.
assert e.val == 3.
-
+
def test_sub(self):
a = adouble(1.)
b = a - 2.
c = 2.- a
-
+
assert b.val == -1.
assert c.val == 1.
-
+
def test_mul(self):
a = adouble(1.5)
b = a * 2.
c = 2.* a
-
+
assert b.val == 3.
assert c.val == 3.
-
+
def test_div(self):
a = adouble(3.)
b = a/2.
c = 2./a
-
+
assert b.val == 3./2.
assert c.val == 2./3.
-
+
def test_truediv(self):
x=1
y=2
ax=adouble(x)
ay=adouble(y)
-
+
z= x.__truediv__(y)
az1=ax.__truediv__(y)
az2=ay.__rtruediv__(x)
az3=ax.__truediv__(ay)
-
+
assert_almost_equal(az1.val, z)
assert_almost_equal(az2.val, z)
assert_almost_equal(az3.val, z)
-
+
def test_pow(self):
r = 5
x = 3.
y = 2.
ax = adouble(x)
ay = adouble(y)
-
+
az1 = ax**ay
az2 = ax**r
az3 = r**ax
-
+
assert_almost_equal(az1.val, x**y)
assert_almost_equal(az2.val, x**r)
assert_almost_equal(az3.val, r**x)
-
+
def test_hyperbolic_functions(self):
x = 3.
ax = adouble(x)
-
+
ash = numpy.sinh(ax)
ach = numpy.cosh(ax)
ath = numpy.tanh(ax)
-
+
assert_almost_equal(ash.val, numpy.sinh(x))
assert_almost_equal(ach.val, numpy.cosh(x))
assert_almost_equal(ath.val, numpy.tanh(x))
-
-
-
+
+
+
def test_fabs(self):
x = 3.
xs = numpy.array([1.,2.,3.])
ax = adouble(x)
axs = adouble(xs)
-
+
aabs = numpy.fabs(ax)
aabss = numpy.fabs(axs)
-
+
assert_almost_equal(aabs.val, numpy.fabs(x))
-
+
def test_abs(self):
x = 3.
xs = numpy.array([1.,2.,3.])
ax = adouble(x)
axs = adouble(xs)
-
+
afabs = abs(ax)
afabss = abs(axs)
-
+
assert_almost_equal(afabs.val, abs(x))
-
+
def test_numpyabs(self):
x = 3.
xs = numpy.array([1.,2.,3.])
ax = adouble(x)
axs = adouble(xs)
-
+
afabs = numpy.abs(ax)
afabss = numpy.abs(axs)
-
+
assert_almost_equal(afabs.val, numpy.abs(x))
-
+
#test_expression('fabs (a) : ', lambda x: numpy.fabs (x), a, a.val)
-
-
+
+
def test_double_condassign_if(self):
x = 3.
y = 4.
cond = 1.
-
+
x = condassign(x,cond,y)
- assert x == 4.
-
+ print x
+ assert_almost_equal(x,4.)
+
x = 3.
y = 4.
cond = -1.
x = condassign(x,cond,y)
print x
- assert x == 3.
-
+ assert_almost_equal(x,3.)
+
def test_double_condassign_if_else(self):
x = 3.
y = 4.
z = 5.
cond = 1.
-
+
x = condassign(x,cond,y,z)
assert x == 4.
-
+
x = 3.
y = 4.
z = 5.
cond = -1.
-
+
x = condassign(x,cond,y,z)
assert x == 5
-
-
+
+
def test_adouble_condassign_if(self):
x = adouble(3.)
y = adouble(4.)
cond = adouble(1.)
-
+
x = condassign(x,cond,y)
- assert x.val == 4.
-
+ print x
+ assert_almost_equal(x.val, 4.)
+
x = adouble(3.)
y = adouble(4.)
cond = adouble(-3.)
x = condassign(x,cond,y)
- assert x.val == 3.
+ print x
+ assert_almost_equal(x.val, 3.)
+
+
+ def test_xuchen_condassign(self):
+ """
+ see https://github.com/b45ch1/pyadolc/issues/12
+ """
+
+ def f(x):
+ a = x + 1
+ b = 2*x
+ x = condassign(x, b-a, b)
+ return x
+
+
+ def g(x):
+ a = x + 1
+ b = 2*x
+ if b-a > 0:
+ x = b
+ return x
+
+
+ trace_on(0)
+ ax = adouble(2.0)
+ independent(ax)
+ ay = f(ax)
+ dependent(ay)
+ trace_off()
+
+ assert_array_almost_equal(g(-1.), zos_forward(0, -1., keep=0)[0])
+ assert_array_almost_equal(g(1.), zos_forward(0, 1., keep=0)[0])
+ assert_array_almost_equal(g(2.), zos_forward(0, 2., keep=0)[0])
def test_adouble_condassign_if_else(self):
@@ -224,62 +259,65 @@ def test_adouble_condassign_if_else(self):
y = adouble(4.)
z = adouble(5.)
cond = adouble(1.)
-
+
x = condassign(x,cond,y,z)
- assert x.val == 4.
-
+ print x
+ assert_almost_equal(x.val, 4.)
+
x = adouble(3.)
y = adouble(4.)
z = adouble(5.)
cond = adouble(-3.)
-
+
x = condassign(x,cond,y,z)
- assert x.val == 5
-
+ print x
+ assert_almost_equal(x.val, 5.)
+
+
class LowLevelFunctionsTests ( TestCase ):
-
+
def test_independent(self):
# 0D
ax = adouble(1)
bx = independent(ax)
assert ax == bx
-
+
# 1D
N = 10
ax = numpy.array([adouble(n) for n in range(N)])
bx = independent(ax)
assert numpy.prod( ax == bx )
-
+
# 2D
N = 2; M=3
ax = numpy.array([[adouble(n+m) for n in range(N)] for m in range(M)])
bx = independent(ax)
assert numpy.prod( ax == bx )
-
+
def test_dependent(self):
# 0D
ax = adouble(1)
bx = dependent(ax)
assert ax == bx
-
+
# 1D
N = 10
ax = numpy.array([adouble(n) for n in range(N)])
bx = dependent(ax)
assert numpy.prod( ax == bx )
-
+
# 2D
N = 2; M=3
ax = numpy.array([[adouble(n+m) for n in range(N)] for m in range(M)])
bx = dependent(ax)
assert numpy.prod( ax == bx )
-
-
+
+
def test_hos_forward_with_keep_then_hos_ti_reverse(self):
"""compute the first columnt of the hessian of f = x_1 x_2 x_3"""
def f(x):
return x[0]*x[1]*x[2]
-
+
#tape f
ax = numpy.array([adouble(0.) for i in range(3)])
trace_on(0)
@@ -288,28 +326,28 @@ def f(x):
ay = f(ax)
dependent(ay)
trace_off()
-
+
x = numpy.array([3.,5.,7.])
V = numpy.zeros((3,1))
V[0,0]=1
-
+
(y,W) = hos_forward(0,x,V,2)
assert y[0] == 105.
assert W[0] == 35.
-
+
U = numpy.zeros((1,2), dtype=float)
U[0,0] = 1.
-
+
Z = hos_ti_reverse(0,U)
assert numpy.prod( Z[:,0] == numpy.array([35., 21., 15.]))
assert numpy.prod( Z[:,1] == numpy.array([0., 7., 5.]))
-
-
+
+
def test_hov_ti_reverse(self):
"""compute the first columnt of the hessian of f = x_1 x_2 x_3"""
def f(x):
return x[0]*x[1]*x[2]
-
+
#tape f
ax = numpy.array([adouble(0.) for i in range(3)])
trace_on(0)
@@ -318,29 +356,29 @@ def f(x):
ay = f(ax)
dependent(ay)
trace_off()
-
+
x = numpy.array([3.,5.,7.])
V = numpy.zeros((3,1))
V[0,0]=1
-
+
(y,W) = hos_forward(0,x,V,2)
assert y[0] == 105.
assert W[0] == 35.
-
+
U = numpy.zeros((1,1,2), dtype=float)
U[0,0,0] = 1.
-
+
Z = hov_ti_reverse(0,U)[0]
print Z[0,:,0]
assert numpy.prod( Z[0,:,0] == numpy.array([35., 21., 15.]))
assert numpy.prod( Z[0,:,1] == numpy.array([0., 7., 5.]))
-
-
+
+
def test_hov_wk_forward_with_keep_then_hos_ov_reverse(self):
"""compute the full hessian of f = x_1 x_2 x_3"""
def f(x):
return x[0]*x[1]*x[2]
-
+
#tape f
ax = numpy.array([adouble(0.) for i in range(3)])
trace_on(0)
@@ -349,27 +387,27 @@ def f(x):
ay = f(ax)
dependent(ay)
trace_off()
-
+
x = numpy.array([3.,5.,7.])
P = 3
V = numpy.zeros((3,P,1))
V[0,0,0] = 1.
V[1,1,0] = 1.
V[2,2,0] = 1.
-
+
(y,W) = hov_wk_forward(0,x,V,2)
assert_almost_equal(y[0],105.)
-
+
U = numpy.zeros((1,2), dtype=float)
U[0,0] = 1.
-
+
Z = hos_ov_reverse(0, P ,U)
-
+
H = numpy.array([[0, x[2], x[1]],[x[2], 0, x[0]], [x[1], x[0],0]],dtype=float)
assert_array_almost_equal(Z[:,:,1],H)
-
-
-
+
+
+
def test_simple_function(self):
def f(x):
y1 = 1./(numpy.fabs(x))
@@ -378,7 +416,7 @@ def f(x):
return y3
def g(x):
return -1./numpy.fabs(x)**2 + 5.
-
+
#tape f
trace_on(0)
x = 2.
@@ -388,15 +426,15 @@ def g(x):
depends_on(ay)
trace_off()
assert_array_almost_equal(g(x), gradient(0,numpy.array([x])))
-
+
def test_tape_to_latex(self):
N = 40
def scalar_f(x):
return 0.5*numpy.dot(x,x)
-
+
x = numpy.array([1.*n for n in range(N)])
ax = adouble(x)
-
+
trace_on(123)
independent(ax)
ay = scalar_f(ax)
@@ -410,8 +448,8 @@ def scalar_f(x):
os.chdir("/tmp")
os.system("pdflatex tape_123.tex ")
os.chdir(cwd)
-
-
+
+
class HighLevelFunctionsTests ( TestCase ):
"""
TESTING HIGH LEVEL CONVENICENCE FUNCTIONS (GRADIENT,HESSIAN, ETC..)
@@ -421,40 +459,40 @@ def test_function(self):
N = 10
def scalar_f(x):
return numpy.dot(x,x)
-
+
x = numpy.ones(N)
ax = adouble(x)
-
+
trace_on(0)
independent(ax)
ay = scalar_f(ax)
dependent(ay)
trace_off()
assert_almost_equal(scalar_f(x),function(0,x))
-
+
def test_gradient(self):
N = 10
def scalar_f(x):
return 0.5*numpy.dot(x,x)
-
+
x = numpy.array([1.*n for n in range(N)])
ax = adouble(x)
-
+
trace_on(0)
independent(ax)
ay = scalar_f(ax)
dependent(ay)
trace_off()
assert_array_almost_equal(x,gradient(0,x))
-
+
def test_hessian(self):
N = 10
def scalar_f(x):
return 0.5*numpy.dot(x,x)
-
+
x = numpy.array([1.*n for n in range(N)])
ax = adouble(x)
-
+
trace_on(0)
independent(ax)
ay = scalar_f(ax)
@@ -462,53 +500,53 @@ def scalar_f(x):
trace_off()
true_H = numpy.eye(N)
assert_array_almost_equal(true_H, hessian(0,x))
-
+
def test_jacobian(self):
N = 31 # dimension
M = 29 # codimension
A = numpy.array([[ 1./N +(n==m) for n in range(N)] for m in range(M)])
def vector_f(x):
return numpy.dot(A,x)
-
+
x = numpy.array([1.*n for n in range(N)])
ax = adouble(x)
-
+
trace_on(123)
independent(ax)
ay = vector_f(ax)
dependent(ay)
trace_off()
assert_array_almost_equal(A, jacobian(123,x))
-
+
def test_hess_vec(self):
N = 1132
def scalar_f(x):
return 0.5*numpy.dot(x,x)
-
+
x = numpy.array([1.*n for n in range(N)])
ax = adouble(x)
-
+
trace_on(0)
independent(ax)
ay = scalar_f(ax)
dependent(ay)
trace_off()
-
+
v = numpy.random.rand(N)
H = numpy.eye(N)
Hv = numpy.dot(H,v)
assert_array_almost_equal( Hv, hess_vec(0,x,v))
-
+
def test_vec_jac(self):
N = 3 # dimension
M = 2 # codimension
A = numpy.array([[ 1./N +(n==m) for n in range(N)] for m in range(M)])
def vector_f(x):
return numpy.dot(A,x)
-
+
x = numpy.array([1.*n for n in range(N)])
ax = adouble(x)
-
+
trace_on(1)
independent(ax)
ay = vector_f(ax)
@@ -517,18 +555,18 @@ def vector_f(x):
u = numpy.random.rand(M)
uJ = numpy.dot(u,A)
assert_array_almost_equal( uJ, vec_jac(1,x,u, 0))
-
-
+
+
def test_jac_vec(self):
N = 3 # dimension
M = 2 # codimension
A = numpy.array([[ 1./N +(n==m) for n in range(N)] for m in range(M)])
def vector_f(x):
return numpy.dot(A,x)
-
+
x = numpy.array([1.*n for n in range(N)])
ax = adouble(x)
-
+
trace_on(1)
independent(ax)
ay = vector_f(ax)
@@ -537,7 +575,7 @@ def vector_f(x):
v = numpy.random.rand(N)
Jv = numpy.dot(A,v)
assert_array_almost_equal( Jv, jac_vec(1,x,v) )
-
+
def test_lagra_hess_vec(self):
""" This test needs improvement: the result is always 0!!"""
N = 3 # dimension
@@ -545,10 +583,10 @@ def test_lagra_hess_vec(self):
A = numpy.array([[ 1./N +(n==m) for n in range(N)] for m in range(M)])
def vector_f(x):
return numpy.dot(A,x)
-
+
x = numpy.array([1.*n for n in range(N)])
ax = adouble(x)
-
+
trace_on(1)
independent(ax)
ay = vector_f(ax)
@@ -557,19 +595,19 @@ def vector_f(x):
u = numpy.random.rand(M)
v = numpy.random.rand(N)
assert_array_almost_equal(numpy.zeros(N,dtype=float), lagra_hess_vec(1,x,u,v) )
-
+
def test_repeated_taping(self):
R = 20 # number of repetitions of the taping
-
+
N = 3 # dimension
M = 2 # codimension
A = numpy.array([[ 1./N +(n==m) for n in range(N)] for m in range(M)])
def vector_f(x):
return numpy.dot(A,x)
-
+
x = numpy.array([1.*n for n in range(N)])
ax = adouble(x)
-
+
for r in range(R):
trace_on(1)
independent(ax)
@@ -579,7 +617,7 @@ def vector_f(x):
u = numpy.random.rand(M)
uJ = numpy.dot(u,A)
assert_array_almost_equal( uJ, vec_jac(1,x,u, 0))
-
+
for r in range(R):
trace_on(r)
independent(ax)
@@ -589,18 +627,18 @@ def vector_f(x):
u = numpy.random.rand(M)
uJ = numpy.dot(u,A)
assert_array_almost_equal( uJ, vec_jac(r,x,u, 0))
-
-
+
+
def test_hov_forward(self):
""" checks only first order"""
N = 3
P = 1
D = 1
epsilon1 = numpy.sqrt(10**-16)
-
+
def f(x):
return numpy.array([x[0]*x[1] + x[0]*x[2], x[1]*x[2]])
-
+
x = numpy.array([1.,2.,3.])
ax = adouble(x)
trace_on(1)
@@ -610,14 +648,14 @@ def f(x):
trace_off()
x = numpy.random.rand(N)
V = numpy.random.rand(N,P,D)
-
+
(y,W) = hov_forward(1, x, V)
-
+
W2 = (f(x+epsilon1*V[:,0,0]) - f(x))/epsilon1
W2 = W2.reshape((2,P,D))
-
+
assert_array_almost_equal(y, f(x))
- assert_array_almost_equal(W, W2)
+ assert_array_almost_equal(W, W2)
if __name__ == '__main__':
try:
View
207 adolc/wrapped_functions.py
@@ -12,7 +12,7 @@ def adouble(x):
return _adolc.adouble(float(x))
elif isinstance(x,_adolc.adouble) or isinstance(x,_adolc.adub):
return _adolc.adouble(x)
-
+
else:
x = numpy.ascontiguousarray(x, dtype=float)
shp = numpy.shape(x)
@@ -36,7 +36,7 @@ def iadouble(x):
is equivalent to:
ax = adouble(0.)
ax.is_independent(x)
-
+
"""
ax = adouble(0.)
@@ -62,14 +62,14 @@ def independent(ax):
xr = numpy.array([axr[n].val for n in range(N)])
map(_adolc.adouble.is_independent,axr,xr)
return ax
-
+
def dependent(ax):
"""
dependent(ax)
INPUT : numpy.array of type adouble
OUTPUT: ax
-
+
Mark ax as dependent.
"""
if isinstance(ax, _adolc.adouble) or isinstance(ax, _adolc.adub):
@@ -78,7 +78,7 @@ def dependent(ax):
else:
axr = numpy.ravel(ax)
N = numpy.size(axr)
-
+
for n in range(N):
if numpy.isscalar(axr[n]):
depends_on(adouble(axr[n]))
@@ -97,7 +97,7 @@ def trace_on(tape_tag):
def trace_off():
"""turn off tracing"""
return _adolc.trace_off()
-
+
def function(tape_tag,x):
@@ -141,12 +141,12 @@ def hessian(tape_tag, x, format='full' ):
x = numpy.ascontiguousarray(x, dtype=float)
assert numpy.size(x) == N
assert numpy.ndim(x) == 1
-
+
if format == 'full':
H = numpy.zeros((N,N), dtype=float)
ints = [i for i in range(N)]
_adolc.hessian(tape_tag, N, x, H)
- H[:] = H[:] + H.T
+ H[:] = H[:] + H.T
H[ints, ints] /= 2.
return H
else:
@@ -166,7 +166,7 @@ def jacobian(tape_tag,x):
J = numpy.zeros((M,N), dtype=float)
_adolc.jacobian(tape_tag, M, N, x, J)
return J
-
+
def vec_jac(tape_tag, x, u, repeat = False):
"""
@@ -225,7 +225,7 @@ def hess_vec(tape_tag, x, v):
def lagra_hess_vec(tape_tag, x, u, v):
"""
evaluate z = u^T F''(x)v, F:R^N -> R^M
-
+
v N-array
u M-array
z N-array
@@ -233,21 +233,21 @@ def lagra_hess_vec(tape_tag, x, u, v):
ts = tapestats(tape_tag)
N = ts['NUM_INDEPENDENTS']
M = ts['NUM_DEPENDENTS']
-
+
x = numpy.ascontiguousarray(x, dtype=float)
assert numpy.size(x) == N
assert numpy.ndim(x) == 1
-
+
v = numpy.ascontiguousarray(v, dtype=float)
assert numpy.size(v) == N
assert numpy.ndim(v) == 1
-
+
u = numpy.ascontiguousarray(u, dtype=float)
assert numpy.size(u) == M
assert numpy.ndim(u) == 1
-
+
z = numpy.zeros(N, dtype=float)
-
+
_adolc.lagra_hess_vec(tape_tag, M, N, x, v, u, z)
return z
@@ -263,14 +263,14 @@ def zos_forward(tape_tag,x,keep):
ts = tapestats(tape_tag)
N = ts['NUM_INDEPENDENTS']
M = ts['NUM_DEPENDENTS']
-
+
x = numpy.ascontiguousarray(x, dtype=float)
assert numpy.size(x) == N
assert numpy.ndim(x) == 1
-
+
y = numpy.zeros(M, dtype=float)
_adolc.zos_forward(tape_tag, M, N, keep, x, y)
-
+
return y
def fos_forward(tape_tag, x, v, keep):
@@ -287,14 +287,14 @@ def fos_forward(tape_tag, x, v, keep):
assert type(tape_tag) == int
assert type(keep) == int
ts = tapestats(tape_tag)
-
+
N = ts['NUM_INDEPENDENTS']
M = ts['NUM_DEPENDENTS']
-
+
x = numpy.ascontiguousarray(x, dtype=float)
assert numpy.size(x) == N
assert numpy.ndim(x) == 1
-
+
v = numpy.ascontiguousarray(v, dtype=float)
assert numpy.size(v) == N
assert numpy.ndim(v) == 1
@@ -302,7 +302,7 @@ def fos_forward(tape_tag, x, v, keep):
y = numpy.zeros(M, dtype=float)
w = numpy.zeros(M, dtype=float)
_adolc.fos_forward(tape_tag, M, N, keep, x, v, y, w)
-
+
return (y,w)
def fov_forward(tape_tag, x, V):
@@ -314,28 +314,28 @@ def fov_forward(tape_tag, x, V):
V is (N x P)-matrix. P directions
W is (M x P)-matrix. P directiona derivatives
"""
-
+
assert type(tape_tag) == int
assert type(keep) == int
ts = tapestats(tape_tag)
-
+
N = ts['NUM_INDEPENDENTS']
M = ts['NUM_DEPENDENTS']
-
+
x = numpy.ascontiguousarray(x, dtype=float)
assert numpy.size(x) == N
assert numpy.ndim(x) == 1
-
+
V = numpy.ascontiguousarray(V, dtype=float)
assert numpy.shape(V)[0] == N
assert numpy.ndim(V) == 2
P = numpy.shape(V)[1]
-
+
y = numpy.zeros(M, dtype=float)
W = numpy.zeros((M,P), dtype=float)
-
+
_adolc.fov_forward(tape_tag, M, N, P, x, V, y, W)
-
+
return (y,W)
def hos_forward(tape_tag, x, V, keep):
@@ -354,24 +354,24 @@ def hos_forward(tape_tag, x, V, keep):
assert type(tape_tag) == int
assert type(keep) == int
ts = tapestats(tape_tag)
-
+
N = ts['NUM_INDEPENDENTS']
M = ts['NUM_DEPENDENTS']
-
+
x = numpy.ascontiguousarray(x, dtype=float)
assert numpy.size(x) == N
assert numpy.ndim(x) == 1
-
+
V = numpy.ascontiguousarray(V, dtype=float)
assert numpy.shape(V)[0] == N
assert numpy.ndim(V) == 2
D = numpy.shape(V)[1]
-
+
y = numpy.zeros(M, dtype=float)
W = numpy.zeros((M,D), dtype=float)
_adolc.hos_forward(tape_tag, M, N, D, keep, x, V, y, W)
-
+
return (y,W)
def hov_forward(tape_tag, x, V):
@@ -384,27 +384,27 @@ def hov_forward(tape_tag, x, V):
V is (N x P x D)-matrix, P directions
W is (M x P x D)-matrix, P directional derivatives
"""
-
+
assert type(tape_tag) == int
ts = tapestats(tape_tag)
-
+
N = ts['NUM_INDEPENDENTS']
M = ts['NUM_DEPENDENTS']
-
+
x = numpy.ascontiguousarray(x, dtype=float)
assert numpy.size(x) == N
assert numpy.ndim(x) == 1
-
+
V = numpy.ascontiguousarray(V, dtype=float)
assert numpy.shape(V)[0] == N
assert numpy.ndim(V) == 3
P = numpy.shape(V)[1]
D = numpy.shape(V)[2]
-
-
+
+
y = numpy.zeros(M, dtype=float)
W = numpy.zeros((M,P,D), dtype=float)
-
+
_adolc.hov_forward(tape_tag, M, N, D, P, x, V, y, W)
return (y,W)
@@ -421,24 +421,24 @@ def hov_wk_forward(tape_tag, x, V, keep):
assert type(keep) == int
assert type(tape_tag) == int
ts = tapestats(tape_tag)
-
+
N = ts['NUM_INDEPENDENTS']
M = ts['NUM_DEPENDENTS']
-
+
x = numpy.ascontiguousarray(x, dtype=float)
assert numpy.size(x) == N
assert numpy.ndim(x) == 1
-
+
V = numpy.ascontiguousarray(V, dtype=float)
assert numpy.shape(V)[0] == N
assert numpy.ndim(V) == 3
P = numpy.shape(V)[1]
D = numpy.shape(V)[2]
-
-
+
+
y = numpy.zeros(M, dtype=float)
W = numpy.zeros((M,P,D), dtype=float)
-
+
_adolc.hov_wk_forward(tape_tag, M, N, D, keep, P, x, V, y, W)
return (y,W)
# raise NotImplementedError
@@ -455,14 +455,14 @@ def fos_reverse(tape_tag, u):
"""
assert type(tape_tag) == int
ts = tapestats(tape_tag)
-
+
N = ts['NUM_INDEPENDENTS']
M = ts['NUM_DEPENDENTS']
-
+
u = numpy.ascontiguousarray(u, dtype=float)
assert numpy.size(u) == M
assert numpy.ndim(u) == 1
-
+
z = numpy.zeros(N, dtype=float)
_adolc.fos_reverse(tape_tag, M, N, u, z)
return z
@@ -476,20 +476,20 @@ def fov_reverse(tape_tag, U):
U is (QxM)-matrix, Q adjoint directions
Z is (QxN)-matrix, adjoint directional derivative Z = U F'(x)
after calling zos_forward, fos_forward or hos_forward with keep = 1
-
+
"""
-
+
assert type(tape_tag) == int
ts = tapestats(tape_tag)
-
+
N = ts['NUM_INDEPENDENTS']
M = ts['NUM_DEPENDENTS']
-
+
U = numpy.ascontiguousarray(U, dtype=float)
assert numpy.ndim(U) == 2
assert numpy.shape(U)[1] == M
Q = numpy.shape(U)[0]
-
+
Z = numpy.zeros((Q,N), dtype=float)
_adolc.fov_reverse(tape_tag, M, N, Q, U, Z)
return Z
@@ -503,16 +503,16 @@ def hos_reverse(tape_tag, D, u):
u is M-vector, adjoint vector
Z is (N x D+1)-matrix, adjoint directional derivative Z = [u^T F'(x), u^T F\" v[:,0],U F\" v[:,1] + 0.5 u^T F^(3) v[:,0],...]
after calling fos_forward or hos_forward with keep = D+1
-
+
"""
assert type(tape_tag) == int
assert type(D) == int
-
+
ts = tapestats(tape_tag)
-
+
N = ts['NUM_INDEPENDENTS']
M = ts['NUM_DEPENDENTS']
-
+
u = numpy.ascontiguousarray(u, dtype=float)
assert numpy.size(u) == M
assert numpy.ndim(u) == 1
@@ -531,12 +531,12 @@ def hos_ti_reverse(tape_tag, U):
after calling fos_forward or hos_forward with keep = D+1
"""
assert type(tape_tag) == int
-
+
ts = tapestats(tape_tag)
-
+
N = ts['NUM_INDEPENDENTS']
M = ts['NUM_DEPENDENTS']
-
+
U = numpy.ascontiguousarray(U, dtype=float)
assert numpy.ndim(U) == 2
assert numpy.shape(U)[0] == M
@@ -550,7 +550,7 @@ def hos_ti_reverse(tape_tag, U):
def hov_reverse(tape_tag, D, U):
"""
this function is deprecated, use hov_ti_reverse instead!
-
+
higher order vector reverse:
(Z,nz) = hov_reverse(tape_tag, D, U)
F:R^N -> R^M
@@ -563,21 +563,21 @@ def hov_reverse(tape_tag, D, U):
"""
assert type(tape_tag) == int
assert type(D) == int
-
+
ts = tapestats(tape_tag)
N = ts['NUM_INDEPENDENTS']
M = ts['NUM_DEPENDENTS']
-
+
U = numpy.ascontiguousarray(U, dtype=float)
assert numpy.ndim(U) == 2
assert numpy.shape(U)[1] == M
Q = numpy.shape(U)[0]
-
+
Z = numpy.zeros((Q, N, D+1), dtype=float)
nz = numpy.zeros((Q,N), dtype=numpy.int16)
_adolc.hov_reverse(tape_tag, M, N, D, Q, U, Z, nz)
-
+
return (Z, nz)
@@ -594,18 +594,18 @@ def hov_ti_reverse(tape_tag, U):
after calling fos_forward or hos_forward with keep = D+1
"""
assert type(tape_tag) == int
-
+
ts = tapestats(tape_tag)
N = ts['NUM_INDEPENDENTS']
M = ts['NUM_DEPENDENTS']
-
+
U = numpy.ascontiguousarray(U, dtype=float)
assert numpy.ndim(U) == 3
assert numpy.shape(U)[1] == M
Q = numpy.shape(U)[0]
Dp1 = numpy.shape(U)[2]
D = Dp1 - 1
-
+
Z = numpy.zeros((Q, N, Dp1), dtype=float)
nz = numpy.zeros((Q,N), dtype=numpy.int16)
_adolc.hov_ti_reverse(tape_tag, M, N, D, Q, U, Z, nz)
@@ -624,22 +624,53 @@ def hos_ov_reverse(tape_tag, P, U):
"""
assert type(tape_tag) == int
assert type(P) == int
-
+
ts = tapestats(tape_tag)
-
+
N = ts['NUM_INDEPENDENTS']
M = ts['NUM_DEPENDENTS']
-
+
U = numpy.ascontiguousarray(U, dtype=float)
assert numpy.ndim(U) == 2
assert numpy.shape(U)[0] == M
Dp1 = numpy.shape(U)[1]
D = Dp1 - 1
-
+
Z = numpy.zeros((N, P, D+1), dtype=float)
_adolc.hos_ov_reverse(tape_tag, M, N, D, P, U, Z)
return Z
+def condassign(x, cond, y, z = None):
+ """
+
+ x = condassign(cond, y, z = None)
+
+ equivalent to:
+ if cond:
+ x = y
+
+ else:
+ if z != None:
+ x = z
+
+ """
+
+ def c(v):
+ if isinstance(v, _adolc.adub):
+ return _adolc.adouble(v)
+ return v
+
+ x = c(x)
+ cond = c(cond)
+ y = c(y)
+
+ if z == None:
+ return _adolc.condassign(x, cond, y)
+ else:
+ z = c(z)
+ return _adolc.condassign(x, cond, y, z)
+
+
def tape_to_latex(tape_tag,x,y):
"""
tape_to_latex(tape_tag,x,y)
@@ -656,20 +687,20 @@ def tape_to_latex(tape_tag,x,y):
def tapestats(tape_tag):
"""
returns a dictionary with information on the tape:
- NUM_INDEPENDENTS, /* # of independent variables */
- NUM_DEPENDENTS, /* # of dependent variables */
- NUM_MAX_LIVES, /* max # of live variables */
- TAY_STACK_SIZE, /* # of values in the taylor (value) stack */
- OP_BUFFER_SIZE, /* # of operations per buffer == OBUFSIZE (usrparms.h) */
- NUM_OPERATIONS, /* overall # of operations */
- OP_FILE_ACCESS, /* operations file written or not */
- NUM_LOCATIONS, /* overall # of locations */
- LOC_FILE_ACCESS, /* locations file written or not */
- NUM_VALUES, /* overall # of values */
- VAL_FILE_ACCESS, /* values file written or not */
- LOC_BUFFER_SIZE, /* # of locations per buffer == LBUFSIZE (usrparms.h) */
- VAL_BUFFER_SIZE, /* # of values per buffer == CBUFSIZE (usrparms.h) */
- TAY_BUFFER_SIZE, /* # of taylors per buffer <= TBUFSIZE (usrparms.h) */
+ NUM_INDEPENDENTS, /* # of independent variables */
+ NUM_DEPENDENTS, /* # of dependent variables */
+ NUM_MAX_LIVES, /* max # of live variables */
+ TAY_STACK_SIZE, /* # of values in the taylor (value) stack */
+ OP_BUFFER_SIZE, /* # of operations per buffer == OBUFSIZE (usrparms.h) */
+ NUM_OPERATIONS, /* overall # of operations */
+ OP_FILE_ACCESS, /* operations file written or not */
+ NUM_LOCATIONS, /* overall # of locations */
+ LOC_FILE_ACCESS, /* locations file written or not */
+ NUM_VALUES, /* overall # of values */
+ VAL_FILE_ACCESS, /* values file written or not */
+ LOC_BUFFER_SIZE, /* # of locations per buffer == LBUFSIZE (usrparms.h) */
+ VAL_BUFFER_SIZE, /* # of values per buffer == CBUFSIZE (usrparms.h) */
+ TAY_BUFFER_SIZE, /* # of taylors per buffer <= TBUFSIZE (usrparms.h) */
"""
assert type(tape_tag) == int
return _adolc.tapestats(tape_tag)
Please sign in to comment.
Something went wrong with that request. Please try again.