Precision of Float64 in Array

question
embedding

#1

Always in the study of my wrapper for ACME.jl, I need to convert 1d and 2d array of float64 from julia (libjulia/c-api) to armadillo C++ (arma::vec / arma::mat).

In julia console, when I dump the array I need I get :

julia> model.c
1×2 Array{Float64,2}:
 -9.4e-8  0.0

This is my code to convert the array to armadillo matrix (assuming that type is always double precision aka Float64), the code is not optimized, more theorical :

static inline arma::mat make_arma_mat (jl_value_t* array_ptr)
{
    arma::mat container;

    if (jl_is_array(array_ptr))
    {
        jl_array_t* array_jl = (jl_array_t*)array_ptr;

        // Get array pointer
        double* raw = (double*)jl_array_ptr(array_jl);

        // Get number of dimensions
        int ndims = jl_array_ndims(array_jl);
        if (ndims != 2)
            throw std::exception("make_arma_mat: container is not a 2d-array");

        // Get the size of the i-th dim
        size_t nrows = jl_array_dim(array_jl, 0);
        size_t ncols = jl_array_dim(array_jl, 1);

        container = arma::mat(nrows, ncols);

        // Fill array with data
        for (size_t c = 0; c<ncols; c++)
        {
            for (size_t r = 0; r<nrows; r++)
            {
                double value = raw[r+nrows*c];
                container(r,c) = value;
            }
        }
    }

    return container;
}

The code work fine for me, I get a correct matrix except with the values precision.

This is the resulted armadillo matrix dump :

c
  0        0

I don’t understand why ? A simple test show that the libjulia float64 can be box and unbox with good resolution :

jl_init (NULL);

double cpp_64 = 1.8e-15;
jl_value_t* arg = jl_box_float64(cpp_64);
double julia_64 = jl_unbox_float64(arg);

std::cout << julia_64 << std::endl;

Output :

    1.8e-15

The array ‘eltype’ value is 8, I cast the raw array ptr to ‘double’, maybe this is not the right way to get the full precision?

Thank you very much for your answer.


#2

This looks like a problem with the precision of your computations, not the Float64 type?


#3

I just get a raw pointer on a double array (see my code on the first post) :

double* raw = (double*)jl_array_ptr(array_jl);

And the resulting double value (without any computation) is wrong :

double value = raw[r+nrows*c];

I don’t use the Float64 type directly.

Maybe it’s not the best way to get array values?

Note : the problem come with negative exponants values :

-9.4e-8
7.2e-14
1.0e-15

Become 0.0000000000 but others values seems to be correct.


#4

Are you sure this isn’t a printing precision issue? See here.


#5

This is not a armadillo precision issue because the problem come inside the MSVC debugger when I breakpoint to the double value.

double i_need = -9.4e-8;
double i_get = raw[r+nrows*c];
container(r,c) = i_get;

The debugger show me:

i_need -9.3999999999999995e-08
i_get   0.00000000000000000

The problem come not only when I get double by casting a Float64 array, but even when I dump some expr (expression) and need to unbox Float64. The following code just dump an expression recursively with Symbol, Int32, Float64 …

static inline void dump_expr (jl_expr_t* expr, std::string indent = "")
{
    if (jl_is_expr(expr))
    {
        std::cout << indent << jl_symbol_name(expr->head) << std::endl;
	size_t nargs = jl_expr_nargs(expr);
	for (int i = 0; i < nargs; ++i)
        {
            jl_value_t* arg = jl_exprarg(expr, i);
            if (jl_is_expr(arg))
            {
                jl_expr_t* subexpr = (jl_expr_t*)jl_exprarg(expr, i);
                if (subexpr != nullptr)
                    dump_expr(subexpr, indent + "    ");
            }
            else
            if (jl_is_symbol(arg))
            {
                const char* subexpr = jl_symbol_name((jl_sym_t*)arg);
                std::cout << indent << "Symbol: " << subexpr << std::endl;
            }
            else
            if (jl_is_int32(arg))
            {
                int subexpr = jl_unbox_int32(arg);
                std::cout << indent << "Int32: " << subexpr << std::endl;
            }
            else
            if (jl_is_float64(arg))
            {
                double subexpr = jl_unbox_float64(arg);
                std::cout << indent << "Float64: " << subexpr << std::endl;
            }
            else
            {
                std::cout << indent << std::string("NEW TYPE = ") << jl_typeof_str(arg) << std::endl;
            }
        }
    }
}

For example, I need to get 2 expression lines with Float64, this is a julia console dump :

res[1 + 1] = 1.8e-15 * (ex - 1) - i # C:\Users\maxprod\.julia\v0.5\ACME\src\elements.jl, line 239:
J[1 + 1,2 + 1] = 7.2e-14ex # C:\Users\maxprod\.julia\v0.5\ACME\src\elements.jl, line 240:

My libjulia/c-api dump is :

block
    line
    Int32: 238
    Symbol: C:\Users\maxprod\.julia\v0.5\ACME\src\elements.jl
    =
        ref
        Symbol: res
            call
            Symbol: +
            Int32: 1
            Int32: 1
        call
        Symbol: -
            call
            Symbol: *
            Float64: 0
                call
                Symbol: -
                Symbol: ex
                Int32: 1
        Symbol: i
    line
    Int32: 239
    Symbol: C:\Users\maxprod\.julia\v0.5\ACME\src\elements.jl
    =
        ref
        Symbol: J
            call
            Symbol: +
            Int32: 1
            Int32: 1
            call
            Symbol: +
            Int32: 2
            Int32: 1
        call
        Symbol: *
        Float64: 0
        Symbol: ex

The two Float64 are 0 value where the values must be, respectively:

Float64 1.8e-15
Float64 7.2e-14

So, this is not really a problem of array or cast of array, but (in my project ofcourse) a global precision problem when the Float64 use negative exponants because in this dump I use correctly the jl_unbox_float64 method.

And more, when the Float64 is a positive non decimal value, I get the correct value:

Symbol: ex
    call
    Symbol: exp
        call
        Symbol: *
        Symbol: v
        Float64: 40

Only the ex = exp(v * 40.0) part of the line dumped with the julia console :

let v = q[2 + 1], i = q[2 + 2], ex = exp(v * 40.0) # C:\Users\maxprod\.julia\v0.5\ACME\src\elements.jl, line 238:

Really, i’m dubious!


#6

It is difficult to figure out what is going on because your examples are not self-contained and don’t show how you are calling the code from Julia.

That said, here is a minimal example which works as expected. Hopefully it will help you figure out the problem.

test.c

#include "julia.h"

void test(jl_value_t* array_jl) {
    jl_array_t* ar = (jl_array_t*)array_jl;
    size_t nrows = jl_array_dim(array_jl, 0);
    size_t ncols = jl_array_dim(array_jl, 1);

    double* raw = (double*)jl_array_data(array_jl);
    for (size_t c = 0; c<ncols; c++)
    {
        for (size_t r = 0; r<nrows; r++)
        {
            double value = raw[r+nrows*c];
            printf("%.*f\n", 20, value);
        }
    }
}

compile with: cc -g -I$JL/src -I$JL/src/support -I$JL/usr/include test.c -otest -shared -ljulia -L$JL/usr/lib

test:

julia> a = rand(2,2); a[1] = 1.8e-15; a[2] = -9.4e-8; a[3] = 7.2e-14;

julia> ccall((:test, :test), Void, (Any,), a)
0.00000000000000180000
-0.00000009400000000000
0.00000000000007200000
0.08850464391479762050

#7

Hi ihnorton. Thank you very much for your answer. You’re right and you’re code too. But when I use extended precision printf like your example, I get :

i_need : -0.00000009400000000000
i_get  :  0.00000000000000000000

I’m sure that’s not a bug of julia ofcourse. I’ve tested some precision tips by the libjulia, no problem. But i’m sure at 99% that I use the same call in the libjulia (c-api) and the monitoring console (julia.exe). That’s why i’m dubious about having slight difference only in float64 with this kind of precision.

I will search in my code a little more, and will be back with some details of my julia call to the ACME.jl package if I don’t find solution. The only point that can generate problem is the sample period (1/sample rate) parameter. All other calls are copy/paste of the ACME example code.

One more time, I thank you very much for you help.