AtomicFunctions

This section describes how to use an atomic function (an inner reusable model) in another model (the outer model).

Both of these models are will be added to the same dynamic library.

The Models

First the model inner and outer models need to be defined:

std::vector<ADCGD> innerModel(const std::vector<ADCGD>& x) {
    std::vector<ADCGD> y(MInner);
    for (size_t i = 0; i < x.size(); i++) {
        y[i] = x[i] * (i + 2);
    }
    return y;
}

std::vector<ADCGD> outerModel(const std::vector<ADCGD>& x, atomic_base<CGD>& atomicfun) {
    std::vector<ADCGD> z(MOuter);

    for (size_t k = 0; k < 2; k++) {
        std::vector<ADCGD> y(MInner);
        std::vector<ADCGD> xInner(x.begin() + k, x.begin() + k + NInner);
        atomicfun(xInner, y);

        for (size_t i = 0; i < y.size(); i++) {
            z[i + k * NInner] = y[i];
        }
    }
    z[4] = x[1] * x[1];
    z[5] = x[0] * x[2] + 1;

    return z;
}

Tape the Models Using CppAD

As part of the process to generate the source code, an ADFun needs to be generated for each of these models.

This involves:

defining the independent variables
calling the model
creating an ADFun and for the dependent variables provided by the models

std::unique_ptr<ADFun<CGD>> tapeInnerModel() {
    std::vector<ADCGD> ax(NInner);
    CppAD::Independent(ax);

    std::vector<ADCGD> ay = innerModel(ax);

    std::unique_ptr<ADFun<CGD>> funInner(new ADFun<CGD>());
    funInner->Dependent(ay);

    return funInner;
}

std::unique_ptr<ADFun<CGD>> tapeOuterModel(CGAtomicFunBridge<double>& atomicFun) {
    std::vector<ADCGD> ax(NOuter);
    CppAD::Independent(ax);

    std::vector<ADCGD> az = outerModel(ax, atomicFun);

    std::unique_ptr<ADFun<CGD>> funOuter(new ADFun<CGD>());
    funOuter->Dependent(az);

    return funOuter;
}

Create a Dynamic Library

Then, these ADFuns are used to create a dynamic library.

It is essential to configure the generation of source code for the Forward/Reverse modes for the inner model.

std::unique_ptr<ModelCSourceGen<double>> prepareInnerModelCompilation(ADFun<CGD>& funInner) {
    auto cSourceInner = std::make_unique<ModelCSourceGen<double>> (funInner, "my_inner_model");

    cSourceInner->setCreateForwardZero(true);
    cSourceInner->setCreateForwardOne(true);
    cSourceInner->setCreateReverseOne(true);
    cSourceInner->setCreateReverseTwo(true);

    return cSourceInner;
}

std::unique_ptr<ModelCSourceGen<double>> prepareOuterModelCompilation(ADFun<CGD>& funOuter) {
    auto cSourceOuter = std::make_unique<ModelCSourceGen<double>>(funOuter, "my_outer_model");

    cSourceOuter->setCreateForwardZero(true);
    cSourceOuter->setCreateSparseJacobian(true);
    cSourceOuter->setCreateSparseHessian(true);

    return cSourceOuter;
}

void compileLibrary() {
    auto funInner = tapeInnerModel();
    CGAtomicFunBridge<double> atomicFun("my_inner_model", *funInner, true);
    auto funOuter = tapeOuterModel(atomicFun);

    auto cSourceInner = prepareInnerModelCompilation(*funInner);
    auto cSourceOuter = prepareOuterModelCompilation(*funOuter);

    ModelLibraryCSourceGen<double> libSourceGen(*cSourceInner, *cSourceOuter);

    GccCompiler<double> compiler;
    // the following options are used so that it is easier to debug generated sources
    compiler.setCompileFlags({"-O0", "-g", "-ggdb", "-D_FORTIFY_SOURCE=2"});
    compiler.setSaveToDiskFirst(true);

    DynamicModelLibraryProcessor<double> p(libSourceGen, LIBRARY_NAME);
    bool loadLib = false;
    p.createDynamicLibrary(compiler, loadLib);
}

Evaluate Derivatives

The dynamic library can now be used to compute derivatives.

Before the outer model can be called, the inner model must be configured as an atomic function by calling outerModel->addExternalModel(). This is required since the atomic function could have been created in a different dynamic library.

void useLibrary() {
    LinuxDynamicLib<double> dynamicLib(LIBRARY_NAME_EXT);
    std::unique_ptr<GenericModel<double>> outerModel = dynamicLib.model("my_outer_model");
    std::unique_ptr<GenericModel<double>> innerModel = dynamicLib.model("my_inner_model");

    outerModel->addExternalModel(*innerModel);

    std::vector<double> xv{1.1, 2.5, 3.5};
    std::vector<double> jac;
    std::vector<size_t> row, col;
    outerModel->SparseJacobian(xv, jac, row, col);

    // print out the result
    for (size_t i = 0; i < jac.size(); ++i)
        std::cout << "dz[" << row[i] << "]/dx[" << col[i] << "] = " << jac[i] << ";" << std::endl;
}

Altogether

The complete example is shown below:

#include <vector>
#include <cppad/cg.hpp>

using namespace CppAD;
using namespace CppAD::cg;

using CGD = CppAD::cg::CG<double>;
using ADCGD = CppAD::AD<CGD>;

const int MOuter = 6;
const int NOuter = 3;
const int MInner = 2;
const int NInner = 2;
const std::string LIBRARY_NAME = "./two_model_library";
const std::string LIBRARY_NAME_EXT = LIBRARY_NAME + system::SystemInfo<>::DYNAMIC_LIB_EXTENSION;


std::vector<ADCGD> innerModel(const std::vector<ADCGD>& x) {
    std::vector<ADCGD> y(MInner);
    for (size_t i = 0; i < x.size(); i++) {
        y[i] = x[i] * (i + 2);
    }
    return y;
}

std::vector<ADCGD> outerModel(const std::vector<ADCGD>& x, atomic_base<CGD>& atomicfun) {
    std::vector<ADCGD> z(MOuter);

    for (size_t k = 0; k < 2; k++) {
        std::vector<ADCGD> y(MInner);
        std::vector<ADCGD> xInner(x.begin() + k, x.begin() + k + NInner);
        atomicfun(xInner, y);

        for (size_t i = 0; i < y.size(); i++) {
            z[i + k * NInner] = y[i];
        }
    }
    z[4] = x[1] * x[1];
    z[5] = x[0] * x[2] + 1;

    return z;
}

std::unique_ptr<ADFun<CGD>> tapeInnerModel() {
    std::vector<ADCGD> ax(NInner);
    CppAD::Independent(ax);

    std::vector<ADCGD> ay = innerModel(ax);

    std::unique_ptr<ADFun<CGD>> funInner(new ADFun<CGD>());
    funInner->Dependent(ay);

    return funInner;
}

std::unique_ptr<ADFun<CGD>> tapeOuterModel(CGAtomicFunBridge<double>& atomicFun) {
    std::vector<ADCGD> ax(NOuter);
    CppAD::Independent(ax);

    std::vector<ADCGD> az = outerModel(ax, atomicFun);

    std::unique_ptr<ADFun<CGD>> funOuter(new ADFun<CGD>());
    funOuter->Dependent(az);

    return funOuter;
}

std::unique_ptr<ModelCSourceGen<double>> prepareInnerModelCompilation(ADFun<CGD>& funInner) {
    auto cSourceInner = std::make_unique<ModelCSourceGen<double>> (funInner, "my_inner_model");

    cSourceInner->setCreateForwardZero(true);
    cSourceInner->setCreateForwardOne(true);
    cSourceInner->setCreateReverseOne(true);
    cSourceInner->setCreateReverseTwo(true);

    return cSourceInner;
}

std::unique_ptr<ModelCSourceGen<double>> prepareOuterModelCompilation(ADFun<CGD>& funOuter) {
    auto cSourceOuter = std::make_unique<ModelCSourceGen<double>>(funOuter, "my_outer_model");

    cSourceOuter->setCreateForwardZero(true);
    cSourceOuter->setCreateSparseJacobian(true);
    cSourceOuter->setCreateSparseHessian(true);

    return cSourceOuter;
}

void compileLibrary() {
    auto funInner = tapeInnerModel();
    CGAtomicFunBridge<double> atomicFun("my_inner_model", *funInner, true);
    auto funOuter = tapeOuterModel(atomicFun);

    auto cSourceInner = prepareInnerModelCompilation(*funInner);
    auto cSourceOuter = prepareOuterModelCompilation(*funOuter);

    ModelLibraryCSourceGen<double> libSourceGen(*cSourceInner, *cSourceOuter);

    GccCompiler<double> compiler;
    // the following options are used so that it is easier to debug generated sources
    compiler.setCompileFlags({"-O0", "-g", "-ggdb", "-D_FORTIFY_SOURCE=2"});
    compiler.setSaveToDiskFirst(true);

    DynamicModelLibraryProcessor<double> p(libSourceGen, LIBRARY_NAME);
    bool loadLib = false;
    p.createDynamicLibrary(compiler, loadLib);
}

void useLibrary() {
    LinuxDynamicLib<double> dynamicLib(LIBRARY_NAME_EXT);
    std::unique_ptr<GenericModel<double>> outerModel = dynamicLib.model("my_outer_model");
    std::unique_ptr<GenericModel<double>> innerModel = dynamicLib.model("my_inner_model");

    outerModel->addExternalModel(*innerModel);

    std::vector<double> xv{1.1, 2.5, 3.5};
    std::vector<double> jac;
    std::vector<size_t> row, col;
    outerModel->SparseJacobian(xv, jac, row, col);

    // print out the result
    for (size_t i = 0; i < jac.size(); ++i)
        std::cout << "dz[" << row[i] << "]/dx[" << col[i] << "] = " << jac[i] << ";" << std::endl;
}

int main() {
    if (!system::isFile(LIBRARY_NAME_EXT)) {
        std::cout << "Creating a new library" << std::endl;
        compileLibrary();
    } else {
        std::cout << "Reusing existing library" << std::endl;
    }

    useLibrary();
}

Provide feedback

Saved searches