Navigation Menu

Skip to content

Commit

Permalink
added repeated sparse jacobian driver
Browse files Browse the repository at this point in the history 
fixed some input checks, they now convert to arrays of dtype=float
  • Loading branch information
@b45ch1
b45ch1 committed May 8, 2009
1 parent 00a9712 commit d8fd666
Show file tree
Hide file tree
Showing 6 changed files with 307 additions and 35 deletions.
46 changes: 45 additions & 1 deletion SConstruct.EXAMPLE
View file
Expand Up @@ -2,11 +2,13 @@ import distutils.sysconfig
import os
import numpy

# 1: BUILD ADOL-C FUNCTIONALITY WITHOUT SPARSE DRIVERS

adolc_include_path = os.getcwd() + '/adolc-2.0.0'
adolc_library_path = os.getcwd() + '/adolc-2.0.0/adolc/.libs'

LIBS = ['adolc',
'boost_python'
'boost_python',
]
LIBPATH = [
adolc_library_path,
Expand All @@ -28,3 +30,45 @@ Default('.')
adolc = env.SharedLibrary(target='_adolc', source=['py_adolc.cpp', 'num_util.cpp'])
#env.Install("./release/adolc", adolc)


# 2: BUILD SPARSE SUPPORT
LIBS = ['adolc',
'boost_python',
'colpack',
]
LIBPATH = [
adolc_library_path,
'/data/walter/opt_software/ColPack/',
]
INCLUDEPATH = [
adolc_include_path,
adolc_include_path + '/adolc/sparse',
adolc_include_path + '/colpack',
'/usr/include/python2.5'
]

env = Environment(
CPPPATH=[distutils.sysconfig.get_python_inc(),numpy.get_include()] + INCLUDEPATH,
CXXFLAGS="-ftemplate-depth-100 -DBOOST_PYTHON_DYNAMIC_LIB -O2",
LIBPATH=["/usr/lib/python2.5/config"] + LIBPATH,
LIBS= LIBS,
RPATH = LIBPATH, #include information where shared libraries can be found to avoid errors like: "ImportError: libboost_python-gcc42-mt-1_34_1.so.1.34.1: cannot open shared object file: No such file or directory"
SHLIBPREFIX="", #gets rid of lib prefix
)
Default('.')


colpack = env.SharedLibrary(target='_colpack',
source=['py_colpack.cpp',
'num_util.cpp',
adolc_include_path +'/adolc/sparse/sparsedrivers.cpp',
adolc_include_path +'/adolc/sparse/sparse_fo_rev.cpp',
adolc_include_path +'/adolc/int_forward_s.c',
adolc_include_path +'/adolc/int_forward_t.c',
adolc_include_path +'/adolc/int_reverse_s.c',
adolc_include_path +'/adolc/int_reverse_t.c',
adolc_include_path +'/adolc/nonl_ind_forward_s.c',
adolc_include_path +'/adolc/nonl_ind_forward_t.c',
adolc_include_path +'/adolc/indopro_forward_s.c',
adolc_include_path +'/adolc/indopro_forward_t.c',
])
70 changes: 36 additions & 34 deletions adolc.py
View file
@@ -1,6 +1,8 @@
# -*- coding: utf-8 -*-
import numpy
import _adolc
import sparse

from _adolc import *

__doc__ = """
Expand Down Expand Up @@ -42,7 +44,7 @@ def adouble(x):
elif isinstance(x,_adolc.adouble):
return _adolc.adouble(x)
else:
x = numpy.asarray(x,dtype=float)
x = numpy.asarray(x, dtype=float)
shp = numpy.shape(x)
xr = numpy.ravel(x)
axr = numpy.array([_adolc.adouble(xr[n]) for n in range(len(xr))])
Expand Down Expand Up @@ -130,31 +132,31 @@ def function(tape_tag,x):
evaluate the function f(x) recorded on tape with index tape_tag
"""
assert type(tape_tag) == int
x = numpy.asarray(x)
x = numpy.asarray(x, dtype=float)
return _adolc.function(tape_tag,x)

def gradient(tape_tag,x):
"""
evaluate the gradient g = f'(x), f:R^N -> R
"""
assert type(tape_tag) == int
x = numpy.asarray(x)
x = numpy.asarray(x, dtype=float)
return _adolc.gradient(tape_tag,x)

def hessian(tape_tag,x):
"""
evaluate the hessian H = f\"(x), f:R^N -> R"
"""
assert type(tape_tag) == int
x = numpy.asarray(x)
x = numpy.asarray(x, dtype=float)
return _adolc.hessian(tape_tag,x)

def jacobian(tape_tag,x):
"""
evaluate the jacobian J = F'(x), F:R^N -> R^M
"""
assert type(tape_tag) == int
x = numpy.asarray(x)
x = numpy.asarray(x, dtype=float)
return _adolc.jacobian(tape_tag,x)

def vec_jac(tape_tag,x,u, repeat):
Expand All @@ -163,26 +165,26 @@ def vec_jac(tape_tag,x,u, repeat):
"""
assert type(repeat) == bool or type(repeat) == int
assert type(tape_tag) == int
x = numpy.asarray(x)
u = numpy.asarray(u)
x = numpy.asarray(x, dtype=float)
u = numpy.asarray(u,dtype=float)
return _adolc.vec_jac(tape_tag,x,u, repeat)

def jac_vec(tape_tag,x,v):
"""
evaluate F'(x)v, F:R^N -> R^M
"""
assert type(tape_tag) == int
x = numpy.asarray(x)
v = numpy.asarray(v)
x = numpy.asarray(x, dtype=float)
v = numpy.asarray(v, dtype=float)
return _adolc.jac_vec(tape_tag,x,v)

def hess_vec(tape_tag,x,v):
"""
evaluate f''(x)v, f:R^N -> R
"""
assert type(tape_tag) == int
x = numpy.asarray(x)
v = numpy.asarray(v)
x = numpy.asarray(x, dtype=float)
v = numpy.asarray(v, dtype=float)
return _adolc.hess_vec(tape_tag,x,v)


Expand All @@ -191,9 +193,9 @@ def lagra_hess_vec(tape_tag,x,u,v):
evaluate u^T F''(x)v, F:R^N -> R^M
"""
assert type(tape_tag) == int
y = numpy.asarray(x)
u = numpy.asarray(u)
v = numpy.asarray(v)
x = numpy.asarray(x, dtype=float)
u = numpy.asarray(u, dtype=float)
v = numpy.asarray(v, dtype=float)
return _adolc.lagra_hess_vec(tape_tag,x,u,v)

def zos_forward(tape_tag,x,keep):
Expand All @@ -206,10 +208,10 @@ def zos_forward(tape_tag,x,keep):
"""
assert type(tape_tag) == int
assert type(keep) == int
x = numpy.asarray(x)
return _adolc.zos_forward(tape_tag,x,keep)
x = numpy.asarray(x, dtype=float)
return _adolc.zos_forward(tape_tag, x, keep)

def fos_forward(tape_tag,x,v,keep):
def fos_forward(tape_tag, x, v, keep):
"""
first order scalar forward:
(y,w) = fos_forward(tape_tag, x, v, keep)
Expand All @@ -222,8 +224,8 @@ def fos_forward(tape_tag,x,v,keep):
"""
assert type(tape_tag) == int
assert type(keep) == int
x = numpy.asarray(x)
v = numpy.asarray(v)
x = numpy.asarray(x, dtype=float)
v = numpy.asarray(v, dtype=float)
return _adolc.fos_forward(tape_tag,x,v,keep)

def fov_forward(tape_tag,x,V):
Expand All @@ -237,8 +239,8 @@ def fov_forward(tape_tag,x,V):
"""
assert type(tape_tag) == int
assert numpy.ndim(V) == 2
x = numpy.asarray(x)
V = numpy.asarray(V)
x = numpy.asarray(x, dtype=float)
V = numpy.asarray(V, dtype=float)
return _adolc.fov_forward(tape_tag,x,V)

def hos_forward(tape_tag, x, V, keep):
Expand All @@ -256,8 +258,8 @@ def hos_forward(tape_tag, x, V, keep):
assert type(tape_tag) == int
assert numpy.ndim(V) == 2
assert type(keep) == int
x = numpy.asarray(x)
V = numpy.asarray(V)
x = numpy.asarray(x, dtype=float)
V = numpy.asarray(V, dtype=float)
return _adolc.hos_forward(tape_tag, x, V, keep)

def hov_forward(tape_tag, x, V):
Expand All @@ -274,8 +276,8 @@ def hov_forward(tape_tag, x, V):
"""
assert type(tape_tag) == int
assert numpy.ndim(V) == 3
x = numpy.asarray(x)
V = numpy.asarray(V)
x = numpy.asarray(x, dtype=float)
V = numpy.asarray(V, dtype=float)
return _adolc.hov_forward(tape_tag, x, V)


Expand All @@ -292,8 +294,8 @@ def hov_wk_forward(tape_tag, x, V, keep):
assert type(tape_tag) == int
assert type(keep) == int
assert numpy.ndim(V) == 3
x = numpy.asarray(x)
V = numpy.asarray(V)
x = numpy.asarray(x, dtype=float)
V = numpy.asarray(V, dtype=float)
return _adolc.hov_wk_forward(tape_tag, x, V, keep)


Expand All @@ -310,7 +312,7 @@ def fos_reverse(tape_tag, u):
"""
assert type(tape_tag) == int
assert numpy.ndim(u) == 1
u = numpy.asarray(u)
u = numpy.asarray(u,dtype=float)
return _adolc.fos_reverse(tape_tag, u)


Expand All @@ -326,7 +328,7 @@ def fov_reverse(tape_tag, U):
"""
assert type(tape_tag) == int
assert numpy.ndim(U) == 2
U = numpy.asarray(U)
U = numpy.asarray(U, dtype=float)
return _adolc.fov_reverse(tape_tag, U)


Expand All @@ -344,7 +346,7 @@ def hos_reverse(tape_tag, D, u):
assert type(tape_tag) == int
assert type(D) == int
assert numpy.ndim(u) == 1
u = numpy.asarray(u)
u = numpy.asarray(u, dtype=float)
return _adolc.hos_reverse(tape_tag, D, u)


Expand All @@ -363,7 +365,7 @@ def hov_reverse(tape_tag, D, U):
assert type(tape_tag) == int
assert type(D) == int
assert numpy.ndim(U) == 2
U = numpy.asarray(U)
U = numpy.asarray(U,dtype=float)
return _adolc.hov_reverse(tape_tag, D, U)


Expand All @@ -381,7 +383,7 @@ def hov_ti_reverse(tape_tag, U):
"""
assert type(tape_tag) == int
assert numpy.ndim(U) == 3
U = numpy.asarray(U)
U = numpy.asarray(U,dtype=float)
return _adolc.hov_ti_reverse(tape_tag, U)

def tape_to_latex(tape_tag,x,y):
Expand All @@ -394,8 +396,8 @@ def tape_to_latex(tape_tag,x,y):
"""
assert type(tape_tag) == int
x = numpy.asarray(x)
y = numpy.asarray(y)
x = numpy.asarray(x, dtype=float)
y = numpy.asarray(y, dtype=float)

return _adolc.tape_to_latex(tape_tag, x, y)

Expand Down
93 changes: 93 additions & 0 deletions py_colpack.cpp
View file
@@ -0,0 +1,93 @@
#include "py_colpack.hpp"

bp::list wrapped_jac_pat(short tape_tag, bpn::array &bpn_x,bpn::array &bpn_options){
int tape_stats[STAT_SIZE];
tapestats(tape_tag, tape_stats);
npy_intp N = tape_stats[NUM_INDEPENDENTS];
npy_intp M = tape_stats[NUM_DEPENDENTS];

double* x = (double*) nu::data(bpn_x);
npy_intp* options = (npy_intp*) nu::data(bpn_options);
unsigned int* JP[M];

jac_pat(tape_tag, M, N, x, JP, options);

bp::list ret_JP(M);

for(int m = 0; m != M; ++m){
ret_JP.append(bp::list(JP[m][0]));
for(int c = 1; c <= JP[m][0]; ++c){
ret_JP[m][c-1] = JP[m][c];
}
}

return ret_JP;
}


bp::list wrapped_sparse_jac_no_repeat(short tape_tag, bpn::array &bpn_x, bpn::array &bpn_options){
int tape_stats[STAT_SIZE];
tapestats(tape_tag, tape_stats);
npy_intp N = tape_stats[NUM_INDEPENDENTS];
npy_intp M = tape_stats[NUM_DEPENDENTS];

double* x = (double*) nu::data(bpn_x);
npy_intp* options = (npy_intp*) nu::data(bpn_options);

npy_intp nnz=-1;
size_t *rind;
size_t *cind;
double *values;
sparse_jac(tape_tag, M, N, 0, x, &nnz, &rind, &cind, &values, options);

bp::object bp_rind ( bp::handle<>(PyArray_SimpleNewFromData(1, &nnz, PyArray_INT, (char*) rind )));
bp::object bp_cind ( bp::handle<>(PyArray_SimpleNewFromData(1, &nnz, PyArray_INT, (char*) cind )));
bp::object bp_values ( bp::handle<>(PyArray_SimpleNewFromData(1, &nnz, PyArray_DOUBLE, (char*) values )));

bpn::array ret_rind = boost::python::extract<boost::python::numeric::array>(bp_rind);
bpn::array ret_cind = boost::python::extract<boost::python::numeric::array>(bp_cind);
bpn::array ret_values = boost::python::extract<boost::python::numeric::array>(bp_values);

bp::list retvals;
retvals.append(nnz);
retvals.append(ret_rind);
retvals.append(ret_cind);
retvals.append(ret_values);

return retvals;

}


bp::list wrapped_sparse_jac_repeat(short tape_tag, bpn::array &bpn_x, npy_intp nnz, bpn::array &bpn_rind, bpn::array &bpn_cind, bpn::array &bpn_values){
int tape_stats[STAT_SIZE];
tapestats(tape_tag, tape_stats);
npy_intp N = tape_stats[NUM_INDEPENDENTS];
npy_intp M = tape_stats[NUM_DEPENDENTS];

double* x = (double*) nu::data(bpn_x);
size_t* rind = (size_t*) nu::data(bpn_rind);
size_t* cind = (size_t*) nu::data(bpn_cind);
double *values = (double*) nu::data(bpn_values);
npy_intp options[4]={0,0,0,0};

sparse_jac(tape_tag, M, N, 1, x, &nnz, &rind, &cind, &values, options);

bp::object bp_rind ( bp::handle<>(PyArray_SimpleNewFromData(1, &nnz, PyArray_INT, (char*) rind )));
bp::object bp_cind ( bp::handle<>(PyArray_SimpleNewFromData(1, &nnz, PyArray_INT, (char*) cind )));
bp::object bp_values ( bp::handle<>(PyArray_SimpleNewFromData(1, &nnz, PyArray_DOUBLE, (char*) values )));

bpn::array ret_rind = boost::python::extract<boost::python::numeric::array>(bp_rind);
bpn::array ret_cind = boost::python::extract<boost::python::numeric::array>(bp_cind);
bpn::array ret_values = boost::python::extract<boost::python::numeric::array>(bp_values);

bp::list retvals;
retvals.append(nnz);
retvals.append(ret_rind);
retvals.append(ret_cind);
retvals.append(ret_values);

return retvals;

}

0 comments on commit d8fd666

Please sign in to comment.