alglib_impl::spline1dbuildcubic fails

jambulat · **Joined:** Mon Nov 22, 2010 11:19 am **Posts:** 10

Hello!
First I would like to say a great thank to the author of this nice library. I used Maple for calculations and because of its low speed decided to translate my code in C++ using the alglib library. It seems the speed-up will be about 100 times!

I didn't know what namespace should I use and set my choice on alglib_impl as it's done in test_c.cpp example file. Later I was sorry about it because of facing anywhere unnecessary, as it seems, parameter "_state". But now I see that it was right choice because I need to calculate complex integrals depending on parameters and it's more convenient to do it directly in autogkiteration cycle.

And now about my current problem. It appeared when I tried to use alglib_impl::spline1dbuildcubic function. And I mention again that before it I successfully used alglib_impl:: function for integration without any fail. Here a code example demonstrating crash on the alglib_impl::spline1dbuildcubic function:

Code:

#include <stdafx.h>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <interpolation.h>

using namespace alglib_impl;


int main(int argc, char **argv)
{
    ae_state state;
   ae_frame _frame_block;
   alglib::real_1d_array x = "[-1.0,-0.5,0.0,+0.5,+1.0]";
   alglib::real_1d_array y = "[+1.0,0.25,0.0,0.25,+1.0]";
    double t = 0.25;
    double v;
    spline1dinterpolant s;

   ae_frame_make(&state, &_frame_block);
   //ae_vector_init(&x, 0, DT_REAL, &state, ae_true);
    //ae_vector_init(&y, 0, DT_REAL, &state, ae_true);
   _spline1dinterpolant_init(&s, &state, ae_true);

   spline1dbuildcubic(x.c_ptr(), y.c_ptr(), x.length(), 2, 0.0, 2, 0.0, &s, &state);
    
   v = spline1dcalc(&s, t, &state);
    printf("%.3f\n", double(v)); // EXPECTED: 0.0580
   getchar();
    return 0;
}

this is the place rising system exception:

Code:

void ae_db_free(ae_dyn_block *block)
{
    if( block->ptr!=NULL )
        ((ae_deallocator)block->deallocator)(block->ptr);
    block->ptr = NULL;
    block->deallocator = ae_free;
}

Could you please anybody say me how to make this code to work properly? Maybe "state" parameter first must be initialized somehow?

PS: OS Windows7, compiler&linker by MSWS9, CPU: 64 bit, multicore.

jambulat · **Joined:** Mon Nov 22, 2010 11:19 am **Posts:** 10

Sorry, I really forgot some initialization procedures:

This works OK:

Code:

#include <stdafx.h>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <interpolation.h>

using namespace alglib_impl;


int main(int argc, char **argv)
{
    ae_state state;
   ae_frame _frame_block;
   alglib::real_1d_array x = "[-1.0,-0.5,0.0,+0.5,+1.0]";
   alglib::real_1d_array y = "[+1.0,0.25,0.0,0.25,+1.0]";
    double t = 0.25;
    double v;
    spline1dinterpolant s;

   ae_frame_make(&state, &_frame_block);
   //ae_vector_init(&x, 0, DT_REAL, &state, ae_true);
    //ae_vector_init(&y, 0, DT_REAL, &state, ae_true);
   _spline1dinterpolant_init(&s, &state, ae_true);

   spline1dbuildcubic(x.c_ptr(), y.c_ptr(), x.length(), 2, 0.0, 2, 0.0, &s, &state);
    
   v = spline1dcalc(&s, t, &state);
    printf("%.3f\n", double(v)); // EXPECTED: 0.0580
   getchar();
    return 0;
}

Sergey.Bochkanov · **Joined:** Fri May 07, 2010 7:06 am **Posts:** 906

Yes, ae_state needs initialization. This structure is used to automatically allocate/deallocate arrays, exception handling, to store environment state and properties, etc.

alglib_impl namespace is not intended to be used directly, it is better to use alglib namespace. Performance penalty for transition between C++ interface and C cores is not very high (especially when you calculate numerical integrals), and alglib namespace is more user-friendly.

However, if you want, you can use internal core directly. But I recommend you to study implementation of corresponding C++ functions before starting to work with core. For example, you've forgot to clear state after you used it (in the current version it is OK, but in future versions of ALGLIB it can lead to memory leaks).

Sergey.Bochkanov · **Joined:** Fri May 07, 2010 7:06 am **Posts:** 906

And, forgot to tell, functions from alglib namespace are more user-friendly and support more features - like automatic determination of input sizes.

jambulat · **Joined:** Mon Nov 22, 2010 11:19 am **Posts:** 10

Sergey.Bochkanov wrote:

For example, you've forgot to clear state after you used it (in the current version it is OK, but in future versions of ALGLIB it can lead to memory leaks).

Sergey, thank you for your answer! And could you please tell me what function clears a state and when it should be called?

Today I faced another problem, now with the spline1dbuildlinear function on such x and y vectors:

Code:

i=0, x=-1.000000e+001, y=4.506061e-009
i=1, x=-9.300000e+000, y=1.820550e-008
i=2, x=-8.600000e+000, y=7.352574e-008
i=3, x=-7.900000e+000, y=2.968071e-007
i=4, x=-7.200000e+000, y=1.197468e-006
i=5, x=-6.500000e+000, y=4.827799e-006
i=6, x=-5.800000e+000, y=1.944683e-005
i=7, x=-5.100000e+000, y=7.824295e-005
i=8, x=-4.400000e+000, y=3.143099e-004
i=9, x=-3.700000e+000, y=1.259774e-003
i=10, x=-3.000000e+000, y=5.031704e-003
i=11, x=-2.880000e+000, y=6.376730e-003
i=12, x=-2.760000e+000, y=8.079761e-003
i=13, x=-2.640000e+000, y=1.023544e-002
i=14, x=-2.520000e+000, y=1.296314e-002
i=15, x=-2.400000e+000, y=1.641328e-002
i=16, x=-2.280000e+000, y=2.077518e-002
i=17, x=-2.160000e+000, y=2.628679e-002
i=18, x=-2.040000e+000, y=3.324660e-002
i=19, x=-1.920000e+000, y=4.202824e-002
i=20, x=-1.800000e+000, y=5.309789e-002
i=21, x=-1.680000e+000, y=6.703490e-002
i=22, x=-1.560000e+000, y=8.455517e-002
i=23, x=-1.440000e+000, y=1.065362e-001
i=24, x=-1.320000e+000, y=1.340411e-001
i=25, x=-1.200000e+000, y=1.683336e-001
i=26, x=-1.080000e+000, y=2.108710e-001
i=27, x=-9.600000e-001, y=2.632419e-001
i=28, x=-8.400000e-001, y=3.269837e-001
i=29, x=-7.200000e-001, y=4.031410e-001
i=30, x=-6.000000e-001, y=4.912729e-001
i=31, x=-4.800000e-001, y=5.873412e-001
i=32, x=-3.600000e-001, y=6.795966e-001
i=33, x=-2.400000e-001, y=7.423721e-001
i=34, x=-1.200000e-001, y=7.341857e-001
i=35, x=0.000000e+000, y=6.210876e-001
i=36, x=1.200000e-001, y=4.258490e-001
i=37, x=2.400000e-001, y=2.228184e-001
i=38, x=3.600000e-001, y=6.719326e-002
i=39, x=4.800000e-001, y=-3.072114e-002
i=40, x=6.000000e-001, y=-8.318193e-002
i=41, x=7.200000e-001, y=-1.056627e-001
i=42, x=8.400000e-001, y=-1.102537e-001
i=43, x=9.600000e-001, y=-1.050420e-001
i=44, x=1.080000e+000, y=-9.505269e-002
i=45, x=1.200000e+000, y=-8.325916e-002
i=46, x=1.320000e+000, y=-7.133737e-002
i=47, x=1.440000e+000, y=-6.016893e-002
i=48, x=1.560000e+000, y=-5.016170e-002
i=49, x=1.680000e+000, y=-4.144936e-002
i=50, x=1.800000e+000, y=-3.401360e-002
i=51, x=1.920000e+000, y=-2.775799e-002
i=52, x=2.040000e+000, y=-2.255166e-002
i=53, x=2.160000e+000, y=-1.825448e-002
i=54, x=2.280000e+000, y=-1.473085e-002
i=55, x=2.400000e+000, y=-1.185671e-002
i=56, x=2.520000e+000, y=-9.522383e-003
i=57, x=2.640000e+000, y=-7.633225e-003
i=58, x=2.760000e+000, y=-6.108885e-003
i=59, x=2.880000e+000, y=-4.882010e-003
i=60, x=3.000000e+000, y=-3.896677e-003
i=61, x=3.700000e+000, y=-1.027049e-003
i=62, x=4.400000e+000, y=-2.652145e-004
i=63, x=5.100000e+000, y=-6.766669e-005
i=64, x=5.800000e+000, y=-1.713168e-005
i=65, x=6.500000e+000, y=-4.314561e-006
i=66, x=7.200000e+000, y=-1.082521e-006
i=67, x=7.900000e+000, y=-2.708450e-007
i=68, x=8.600000e+000, y=-6.762023e-008
i=69, x=9.300000e+000, y=-1.685407e-008
i=70, x=1.000000e+001, y=-4.195189e-009

spline1dbuildlinear and spline1dbuildakima work with them OK, but spline1dbuildcubic rises exception:

ae_assert(isfinitevector(d, n, _state), "Spline1DBuildHermite: D contains infinite or NAN values!", _state);

through this code:

Code:

 

void spline1dbuildhermite(/* Real    */ ae_vector* x,.........................................   
     /*
     * check and sort points
     */
    ae_assert(isfinitevector(x, n, _state), "Spline1DBuildHermite: X contains infinite or NAN values!", _state);
    ae_assert(isfinitevector(y, n, _state), "Spline1DBuildHermite: Y contains infinite or NAN values!", _state);
    ae_assert(isfinitevector(d, n, _state), "Spline1DBuildHermite: D contains infinite or NAN values!", _state);
    heapsortdpoints(x, y, d, n, _state);
    ae_assert(aredistinct(x, n, _state), "Spline1DBuildHermite: at least two consequent points are too close!", _state);
    
    

My code generating this is following:

Code:

   for(i=0;i<y.length();i++)
   {
      y(i) = _ZXcorr(x(i),k*d,_state);
      printf("i=%d, x=%e, y=%e\n",i,x(i),y(i));
   }

   spline1dbuildcubic(x.c_ptr(), y.c_ptr(), n, 2, 0.0, 2, 0.0, &ZXcorr_ipt, _state);
   //spline1dbuildlinear(x.c_ptr(), y.c_ptr(), n, &ZXcorr_ipt, _state);

Sergey.Bochkanov · **Joined:** Fri May 07, 2010 7:06 am **Posts:** 906

ALGLIB state is cleared with alglib_impl::ae_state_clear(&_alglib_env_state) call.

I've checked your data (not code because you've omitted a lot of initialization code) and everything works correctly. Maybe you've forgot to initialize something or somehow accidentally corrupted memory. Again - why don't you try to work with C++ interface instead of C layer? It is more error-prone, more convenient, and only marginally slower. And if you need very efficient calculation of spline values in many points, why don't you use gridconv() family of functions?

jambulat · **Joined:** Mon Nov 22, 2010 11:19 am **Posts:** 10

I use C-interface because I find it more convenient for computing multiple integrals with parameters, it's no need to define new function with one parameter, as it is realized in alglib:: namespace functions, instead of it I can do all I need directly in iteration cycle.

Now spline1dbuildakima also falls under some iteration (I call in a cycle a function that builds splines many times). Data as follows:

Code:

i=0, x=-1.000000e+001, y=5.118948e-009
i=1, x=-9.300000e+000, y=2.078599e-008
i=2, x=-8.600000e+000, y=8.442384e-008
i=3, x=-7.900000e+000, y=3.429986e-007
i=4, x=-7.200000e+000, y=1.394093e-006
i=5, x=-6.500000e+000, y=5.669162e-006
i=6, x=-5.800000e+000, y=2.307045e-005
i=7, x=-5.100000e+000, y=9.397858e-005
i=8, x=-4.400000e+000, y=3.833889e-004
i=9, x=-3.700000e+000, y=1.567595e-003
i=10, x=-3.000000e+000, y=6.433627e-003
i=11, x=-2.880000e+000, y=8.200172e-003
i=12, x=-2.760000e+000, y=1.045403e-002
i=13, x=-2.640000e+000, y=1.333059e-002
i=14, x=-2.520000e+000, y=1.700332e-002
i=15, x=-2.400000e+000, y=2.169466e-002
i=16, x=-2.280000e+000, y=2.769020e-002
i=17, x=-2.160000e+000, y=3.535704e-002
i=18, x=-2.040000e+000, y=4.516793e-002
i=19, x=-1.920000e+000, y=5.773275e-002
i=20, x=-1.800000e+000, y=7.384033e-002
i=21, x=-1.680000e+000, y=9.451391e-002
i=22, x=-1.560000e+000, y=1.210859e-001
i=23, x=-1.440000e+000, y=1.552995e-001
i=24, x=-1.320000e+000, y=1.994497e-001
i=25, x=-1.200000e+000, y=2.565825e-001
i=26, x=-1.080000e+000, y=3.307846e-001
i=27, x=-9.600000e-001, y=4.276208e-001
i=28, x=-8.400000e-001, y=5.548246e-001
i=29, x=-7.200000e-001, y=7.234596e-001
i=30, x=-6.000000e-001, y=9.500280e-001
i=31, x=-4.800000e-001, y=1.260695e+000
i=32, x=-3.600000e-001, y=1.700832e+000
i=33, x=-2.400000e-001, y=2.359456e+000
i=34, x=-1.200000e-001, y=3.420713e+000
i=35, x=0.000000e+000, y=3.964653e+000
i=36, x=1.200000e-001, y=1.045948e+000
i=37, x=2.400000e-001, y=2.475926e-002
i=38, x=3.600000e-001, y=-2.795873e-001
i=39, x=4.800000e-001, y=-3.577587e-001
i=40, x=6.000000e-001, y=-3.530987e-001
i=41, x=7.200000e-001, y=-3.179900e-001
i=42, x=8.400000e-001, y=-2.739382e-001
i=43, x=9.600000e-001, y=-2.301128e-001
i=44, x=1.080000e+000, y=-1.902663e-001
i=45, x=1.200000e+000, y=-1.556587e-001
i=46, x=1.320000e+000, y=-1.263953e-001
i=47, x=1.440000e+000, y=-1.020708e-001
i=48, x=1.560000e+000, y=-8.208536e-002
i=49, x=1.680000e+000, y=-6.580045e-002
i=50, x=1.800000e+000, y=-5.261168e-002
i=51, x=1.920000e+000, y=-4.197982e-002
i=52, x=2.040000e+000, y=-3.344002e-002
i=53, x=2.160000e+000, y=-2.660021e-002
i=54, x=2.280000e+000, y=-2.113460e-002
i=55, x=2.400000e+000, y=-1.677533e-002
i=56, x=2.520000e+000, y=-1.330390e-002
i=57, x=2.640000e+000, y=-1.054310e-002
i=58, x=2.760000e+000, y=-8.349900e-003
i=59, x=2.880000e+000, y=-6.609255e-003
i=60, x=3.000000e+000, y=-5.228911e-003
i=61, x=3.700000e+000, y=-1.323616e-003
i=62, x=4.400000e+000, y=-3.323377e-004
i=63, x=5.100000e+000, y=-8.304436e-005
i=64, x=5.800000e+000, y=-2.068720e-005
i=65, x=6.500000e+000, y=-5.142581e-006
i=66, x=7.200000e+000, y=-1.276469e-006
i=67, x=7.900000e+000, y=-3.164890e-007
i=68, x=8.600000e+000, y=-7.840453e-008
i=69, x=9.300000e+000, y=-1.941059e-008
i=70, x=1.000000e+001, y=-4.802969e-009

exception:
ae_assert(isfinitevector(y, n, _state), "Spline1DBuildAkima: Y contains infinite or NAN values!", _state);

About optimization with gridconv() functions I'll think when the all remaining code works correctly))

I can send you my whole project in Visual Studio... If you have some time to look through it, maybe we could find where corruption happening...

OK, first I try to translate whole project in C++ interface (Is this the same as using the functions of the alglib:: namespace?)...

Sergey.Bochkanov · **Joined:** Fri May 07, 2010 7:06 am **Posts:** 906

Quote:

I use C-interface because I find it more convenient for computing multiple integrals with parameters, it's no need to define new function with one parameter, as it is realized in alglib:: namespace functions, instead of it I can do all I need directly in iteration cycle.

You can do exactly same thing by calling alglib::autogkiteration and working with alglib::autogkstate (if you use autogk subpackage). It is better than working with C interface, and more reliable.

You can post your code here (as zip archive), I'll look at it

jambulat · **Joined:** Mon Nov 22, 2010 11:19 am **Posts:** 10

Have removed all unnecessary code as:

Code:

    
_autogkstate_init(&state, _state, ae_true);
_autogkreport_init(&rep, _state, ae_true);
autogksmooth(0, r_x_lim, &state, _state);

... And now all works perfectly, without any crash! I treated the file test_c.cpp as an example that I should imitate... That's why I've wrote several hundreds lines of code using alglib_impl:: namespace functions.

Once again, Sergey, thank you for your support and for your excellent library!

forum.alglib.net

Forum rules

alglib_impl::spline1dbuildcubic fails

Who is online