Function compareJacobianValue

test/testingutilities.cpp:317–439

Compares the jacobian of the approximant to the central difference of the approximant function value */

/* Compares the jacobian of the approximant to the central difference of the approximant function value */
void compareJacobianValue(TestFunction *exact,
                          std::function<Function *(const DataTable &table)> approx_gen_func,
                          size_t numSamplePoints, size_t numEvalPoints,
                          double one_eps, double two_eps, double inf_eps)
{
    auto dim = exact->getNumVariables();

    auto samplePoints = linspace(dim, -5, 5, std::pow(numSamplePoints, 1.0/dim));
    auto evalPoints = linspace(dim, -4.95, 4.95, std::pow(numEvalPoints, 1.0/dim));

    DataTable table = sample(exact, samplePoints);

    Function *approx = approx_gen_func(table);

    INFO("Approximant: " << approx->getDescription());
    INFO("Function: " << exact->getFunctionStr());

    DenseVector oneNormVec(evalPoints.size());
    DenseVector twoNormVec(evalPoints.size());
    DenseVector infNormVec(evalPoints.size());

    double maxOneNormError = 0.0;
    double maxTwoNormError = 0.0;
    double maxInfNormError = 0.0;

    DenseVector maxOneNormErrorPoint(dim);
    maxOneNormErrorPoint.fill(0.0);
    DenseVector maxTwoNormErrorPoint(dim);
    maxTwoNormErrorPoint.fill(0.0);
    DenseVector maxInfNormErrorPoint(dim);
    maxInfNormErrorPoint.fill(0.0);

    int i = 0;
    for (auto &point : evalPoints)
    {
        DenseVector x = vecToDense(point);

        // Compare the central difference to the approximated jacobian
        DenseMatrix exactValue = approx->centralDifference(x);
        DenseMatrix approxValue = approx->evalJacobian(x);

        DenseVector error = DenseVector::Zero(exactValue.cols());
        for (size_t j = 0; j < (size_t) error.size(); ++j)
        {
            error(j) = getError(exactValue(j), approxValue(j));
        }

        oneNormVec(i) = getOneNorm(error) / error.size(); // "Average"
        twoNormVec(i) = getTwoNorm(error);
        infNormVec(i) = getInfNorm(error);

        if (oneNormVec(i) > maxOneNormError)
        {
            maxOneNormError = oneNormVec(i);
            maxOneNormErrorPoint = x;
        }
        if (twoNormVec(i) > maxTwoNormError)
        {