ext/gsl_narray.c in gsl-1.14.5 vs ext/gsl_narray.c in gsl-1.14.6

- old
+ new

@@ -16,10 +16,11 @@ static VALUE rb_gsl_na_to_gsl_vector_int_method(VALUE na); static VALUE rb_gsl_vector_int_to_na(VALUE obj); static VALUE rb_gsl_na_to_gsl_vector_int(VALUE obj, VALUE na); /* GSL::Vector -> NArray */ + static VALUE rb_gsl_vector_to_narray(VALUE obj, VALUE klass) { gsl_vector *v = NULL; VALUE nary; int shape[1]; @@ -37,27 +38,72 @@ } } return nary; } +static VALUE rb_gsl_vector_complex_to_narray(VALUE obj, VALUE klass) +{ + gsl_vector_complex *v = NULL; + VALUE nary; + int shape[1]; + Data_Get_Struct(obj, gsl_vector_complex, v); + shape[0] = v->size; + nary = na_make_object(NA_DCOMPLEX, 1, shape, klass); + if (v->stride == 1) { + memcpy(NA_PTR_TYPE(nary,double*), v->data, shape[0]*2*sizeof(double)); + } else { + int i; + struct NARRAY *na; + GetNArray(nary, na); + for(i=0; i < 2*v->size; i++) { + (NA_PTR_TYPE(nary,gsl_complex*))[i] = gsl_vector_complex_get(v, i); + } + } + return nary; +} + static VALUE rb_gsl_vector_to_na(VALUE obj) { - return rb_gsl_vector_to_narray(obj, cNArray); + VALUE na = Qnil; + + if(VECTOR_P(obj)) + na = rb_gsl_vector_to_narray(obj, cNArray); + else if(VECTOR_COMPLEX_P(obj)) + na = rb_gsl_vector_complex_to_narray(obj, cNArray); + else + rb_raise(rb_eRuntimeError, "unexpected type '%s'", + rb_obj_classname(obj)); + + return na; } +static VALUE rb_gsl_vector_complex_to_na(VALUE obj) +{ + return rb_gsl_vector_complex_to_narray(obj, cNArray); +} + static VALUE rb_gsl_vector_to_nvector(VALUE obj) { return rb_gsl_vector_to_narray(obj, cNVector); } +static VALUE rb_gsl_vector_complex_to_nvector(VALUE obj) +{ + return rb_gsl_vector_complex_to_narray(obj, cNVector); +} + +/* GSL::Vector -> NArray view */ + static struct NARRAY* rb_gsl_na_view_alloc(int rank, int total, int type) { struct NARRAY *na; na = (struct NARRAY*) malloc(sizeof(struct NARRAY)); na->rank = rank; na->total = total; na->type = type; + // TODO Set na->ref to a one element NArray of type NAObject that contains + // the GSL::Vector being referenced. na->ref = Qtrue; /* to initialize */ na->shape = (int *) malloc(sizeof(int)*rank); return na; } @@ -68,10 +114,11 @@ } static VALUE rb_gsl_vector_to_narray_ref(VALUE obj, VALUE klass) { gsl_vector *v = NULL; + gsl_vector_complex *vc = NULL; gsl_vector_int *vi = NULL; VALUE nary; struct NARRAY *na; if (VECTOR_P(obj)) { Data_Get_Struct(obj, gsl_vector, v); @@ -87,12 +134,21 @@ rb_raise(rb_eRuntimeError, "Cannot make a reference obj: stride!=1"); } na = rb_gsl_na_view_alloc(1, vi->size, NA_LINT); na->shape[0] = vi->size; na->ptr = (char *) vi->data; + } else if (VECTOR_COMPLEX_P(obj)) { + Data_Get_Struct(obj, gsl_vector_complex, vc); + if (vc->stride != 1) { + rb_raise(rb_eRuntimeError, "Cannot make a reference obj: stride!=1"); + } + na = rb_gsl_na_view_alloc(1, vc->size, NA_DCOMPLEX); + na->shape[0] = vc->size; + na->ptr = (char *) vc->data; } else { - rb_raise(rb_eRuntimeError, "???"); + rb_raise(rb_eRuntimeError, "cannot convert %s to NArray reference object", + rb_obj_classname(obj)); } nary = Data_Wrap_Struct(klass, 0, rb_gsl_na_view_free, na); return nary; } @@ -139,95 +195,160 @@ /* singleton method */ static VALUE rb_gsl_na_to_gsl_vector(VALUE obj, VALUE na) { return Data_Wrap_Struct(cgsl_vector, 0, gsl_vector_free, - na_to_gv(na)); + na_to_gv(na)); } static VALUE rb_gsl_na_to_gsl_vector_view(VALUE obj, VALUE na) { return Data_Wrap_Struct(cgsl_vector_view, 0, gsl_vector_view_free, - na_to_gv_view(na)); + na_to_gv_view(na)); } +static VALUE rb_gsl_na_to_gsl_vector_complex(VALUE obj, VALUE na) +{ + return Data_Wrap_Struct(cgsl_vector_complex, 0, gsl_vector_complex_free, + na_to_gv_complex(na)); +} + +static VALUE rb_gsl_na_to_gsl_vector_complex_view(VALUE obj, VALUE na) +{ + return Data_Wrap_Struct(cgsl_vector_complex_view, 0, gsl_vector_complex_view_free, + na_to_gv_complex_view(na)); +} + static VALUE rb_gsl_na_to_gsl_vector_int(VALUE obj, VALUE na) { return Data_Wrap_Struct(cgsl_vector_int, 0, gsl_vector_int_free, - na_to_gv_int(na)); + na_to_gv_int(na)); } static VALUE rb_gsl_na_to_gsl_vector_int_view(VALUE obj, VALUE na) { return Data_Wrap_Struct(cgsl_vector_int_view, 0, rb_gsl_vector_int_view_free, - na_to_gv_int_view(na)); + na_to_gv_int_view(na)); } static VALUE rb_gsl_na_to_gsl_vector_method(VALUE na) { - return Data_Wrap_Struct(cgsl_vector, 0, gsl_vector_free, - na_to_gv(na)); + VALUE v; + + if(NA_TYPE(na) == NA_SCOMPLEX || NA_TYPE(na) == NA_DCOMPLEX) + v = Data_Wrap_Struct(cgsl_vector_complex, 0, gsl_vector_complex_free, + na_to_gv_complex(na)); + else + v = Data_Wrap_Struct(cgsl_vector, 0, gsl_vector_free, + na_to_gv(na)); + + return v; } static VALUE rb_gsl_na_to_gsl_vector_view_method(VALUE na) { - return Data_Wrap_Struct(cgsl_vector_view, 0, gsl_vector_view_free, - na_to_gv_view(na)); + VALUE v; + + if(NA_TYPE(na) == NA_SCOMPLEX || NA_TYPE(na) == NA_DCOMPLEX) + v = Data_Wrap_Struct(cgsl_vector_complex_view, 0, gsl_vector_complex_view_free, + na_to_gv_complex_view(na)); + else + v = Data_Wrap_Struct(cgsl_vector_view, 0, gsl_vector_view_free, + na_to_gv_view(na)); + + return v; } static VALUE rb_gsl_na_to_gsl_vector_int_method(VALUE na) { return Data_Wrap_Struct(cgsl_vector_int, 0, gsl_vector_int_free, - na_to_gv_int(na)); + na_to_gv_int(na)); } static VALUE rb_gsl_na_to_gsl_vector_int_view_method(VALUE na) { return Data_Wrap_Struct(cgsl_vector_int_view, 0, rb_gsl_vector_int_view_free, - na_to_gv_int_view(na)); + na_to_gv_int_view(na)); } gsl_vector* na_to_gv(VALUE na) { gsl_vector *v = NULL; - VALUE nary; + VALUE nary = na; v = gsl_vector_alloc(NA_TOTAL(na)); - nary = na_change_type(na, NA_DFLOAT); - memcpy(v->data, NA_PTR_TYPE(nary,double*), v->size*sizeof(double)); + if(NA_TYPE(na) != NA_DFLOAT) { + nary = na_change_type(na, NA_DFLOAT); + } + memcpy(v->data, NA_PTR_TYPE(na,double*), v->size*sizeof(double)); return v; } gsl_vector_view* na_to_gv_view(VALUE na) { gsl_vector_view *v = NULL; - VALUE ary2; + + // Raise exception if na's type is not NA_DFLOAT. + if(NA_TYPE(na) != NA_DFLOAT) + rb_raise(rb_eTypeError, "GSL::Vector::View requires NArray be DFLOAT"); + v = gsl_vector_view_alloc(); - ary2 = na_change_type(na, NA_DFLOAT); - v->vector.data = NA_PTR_TYPE(ary2,double*); + v->vector.data = NA_PTR_TYPE(na,double*); v->vector.size = NA_TOTAL(na); v->vector.stride = 1; v->vector.owner = 0; return v; } +gsl_vector_complex* na_to_gv_complex(VALUE na) +{ + gsl_vector_complex *v = NULL; + VALUE nary = na; + v = gsl_vector_complex_alloc(NA_TOTAL(na)); + if(NA_TYPE(na) != NA_DCOMPLEX) { + nary = na_change_type(na, NA_DCOMPLEX); + } + memcpy(v->data, NA_PTR_TYPE(na,gsl_complex*), v->size*sizeof(gsl_complex)); + return v; +} + +gsl_vector_complex_view* na_to_gv_complex_view(VALUE na) +{ + gsl_vector_complex_view *v = NULL; + + // Raise exception if na's type is not NA_DCOMPLEX + if(NA_TYPE(na) != NA_DCOMPLEX) + rb_raise(rb_eTypeError, "GSL::Vector::Complex::View requires NArray be DCOMPLEX"); + + v = gsl_vector_complex_view_alloc(); + v->vector.data = NA_PTR_TYPE(na,double*); + v->vector.size = NA_TOTAL(na); + v->vector.stride = 1; + v->vector.owner = 0; + return v; +} + gsl_vector_int* na_to_gv_int(VALUE na) { gsl_vector_int *v = NULL; - VALUE nary; + VALUE nary = na; v = gsl_vector_int_alloc(NA_TOTAL(na)); - nary = na_change_type(na, NA_LINT); + if(NA_TYPE(na) != NA_LINT) { + nary = na_change_type(na, NA_LINT); + } memcpy(v->data, NA_PTR_TYPE(nary,int*), v->size*sizeof(int)); return v; } gsl_vector_int_view* na_to_gv_int_view(VALUE na) { gsl_vector_int_view *v = NULL; - VALUE ary2; + + // Raise exception if na's type is not NA_LINT + if(NA_TYPE(na) != NA_LINT) + rb_raise(rb_eTypeError, "GSL::Vector::Int::View requires NArray be LINT"); v = rb_gsl_vector_int_view_alloc(NA_TOTAL(na)); - ary2 = na_change_type(na, NA_LINT); - v->vector.data = NA_PTR_TYPE(ary2,int*); + v->vector.data = NA_PTR_TYPE(na,int*); v->vector.size = NA_TOTAL(na); v->vector.stride = 1; v->vector.owner = 0; return v; } @@ -242,11 +363,11 @@ shape[0] = m->size2; shape[1] = m->size1; nary = na_make_object(NA_DFLOAT, 2, shape, klass); for (i = 0; i < shape[1]; i++) { memcpy(NA_PTR_TYPE(nary,double*)+(i*shape[0]), m->data+(i*m->tda), - shape[0]*sizeof(double)); + shape[0]*sizeof(double)); } return nary; } static VALUE rb_gsl_matrix_to_na(VALUE obj, VALUE klass) @@ -269,11 +390,11 @@ shape[0] = m->size2; shape[1] = m->size1; nary = na_make_object(NA_LINT, 2, shape, klass); for (i = 0; i < shape[1]; i++) { memcpy(NA_PTR_TYPE(nary,int*)+(i*shape[0]), m->data+(i*m->tda), - shape[0]*sizeof(int)); + shape[0]*sizeof(int)); } return nary; } static VALUE rb_gsl_matrix_int_to_na(VALUE obj) @@ -412,10 +533,14 @@ gsl_matrix_view* na_to_gm_view(VALUE nna) { gsl_matrix_view *m = NULL; VALUE ary2; struct NARRAY *na = NULL; + + // Raise exception if nna's type is not NA_DFLOAT + if(NA_TYPE(nna) != NA_DFLOAT) + rb_raise(rb_eTypeError, "GSL::Matrix::View requires NArray be DFLOAT"); GetNArray(nna, na); m = gsl_matrix_view_alloc(); ary2 = na_change_type(nna, NA_DFLOAT); m->matrix.data = NA_PTR_TYPE(ary2,double*); m->matrix.size1 = na->shape[1]; @@ -440,10 +565,14 @@ gsl_matrix_int_view* na_to_gm_int_view(VALUE nna) { gsl_matrix_int_view *m = NULL; VALUE ary2; struct NARRAY *na = NULL; + + // Raise exception if nna's type is not NA_LINT + if(NA_TYPE(nna) != NA_LINT) + rb_raise(rb_eTypeError, "GSL::Matrix::Int::View requires NArray be LINT"); GetNArray(nna, na); m = rb_gsl_matrix_int_view_alloc(na->shape[1], na->shape[0]); ary2 = na_change_type(nna, NA_LINT); m->matrix.data = NA_PTR_TYPE(ary2,int*); m->matrix.size1 = na->shape[1]; @@ -487,18 +616,18 @@ h = gsl_histogram_alloc(n); for (i = 0; i <= n; i++) h->range[i] = NUM2DBL(rb_ary_entry(argv[0], i)); break; default: if (VECTOR_P(argv[0])) { - Data_Get_Struct(argv[0], gsl_vector, ranges); - n = ranges->size - 1; - h = gsl_histogram_alloc(n); - gsl_histogram_set_ranges(h, ranges->data, ranges->size); + Data_Get_Struct(argv[0], gsl_vector, ranges); + n = ranges->size - 1; + h = gsl_histogram_alloc(n); + gsl_histogram_set_ranges(h, ranges->data, ranges->size); } else if (NA_IsNArray(argv[0])) { - ptr_range = get_vector_ptr(argv[0], &stride, &n); - h = gsl_histogram_alloc(n); - gsl_histogram_set_ranges(h, ptr_range, n); + ptr_range = get_vector_ptr(argv[0], &stride, &n); + h = gsl_histogram_alloc(n); + gsl_histogram_set_ranges(h, ptr_range, n); } break; } break; case 2: @@ -508,11 +637,11 @@ min = NUM2DBL(rb_ary_entry(argv[1], 0)); max = NUM2DBL(rb_ary_entry(argv[1], 1)); break; default: rb_raise(rb_eTypeError, "wrong argument type %s (Array expected)", - rb_class2name(CLASS_OF(argv[1]))); + rb_class2name(CLASS_OF(argv[1]))); break; } h = gsl_histogram_alloc(n); gsl_histogram_set_ranges_uniform(h, min, max); break; @@ -554,19 +683,44 @@ rb_define_singleton_method(cgsl_vector, "to_gslv_view", rb_gsl_na_to_gsl_vector_view, 1); rb_define_singleton_method(cgsl_vector, "to_gv_view", rb_gsl_na_to_gsl_vector_view, 1); rb_define_singleton_method(cgsl_vector, "na_to_gslv_view", rb_gsl_na_to_gsl_vector_view, 1); - rb_define_singleton_method(cgsl_vector, "na_to_gv_view", rb_gsl_na_to_gsl_vector, 1); + rb_define_singleton_method(cgsl_vector, "na_to_gv_view", rb_gsl_na_to_gsl_vector_view, 1); rb_define_method(cNArray, "to_gv", rb_gsl_na_to_gsl_vector_method, 0); rb_define_alias(cNArray, "to_gslv", "to_gv"); rb_define_method(cNArray, "to_gv_view", rb_gsl_na_to_gsl_vector_view_method, 0); rb_define_alias(cNArray, "to_gslv_view", "to_gv_view"); rb_define_alias(cNArray, "to_gv2", "to_gv_view"); rb_define_alias(cNArray, "to_gv_ref", "to_gv_view"); /*****/ + rb_define_method(cgsl_vector_complex, "to_na", rb_gsl_vector_complex_to_na, 0); + rb_define_alias(cgsl_vector_complex, "to_narray", "to_na"); + + rb_define_method(cgsl_vector_complex, "to_na2", rb_gsl_vector_to_na_ref, 0); + rb_define_alias(cgsl_vector_complex, "to_narray_ref", "to_na2"); + rb_define_alias(cgsl_vector_complex, "to_na_ref", "to_na2"); + + rb_define_method(cgsl_vector_complex, "to_nv", rb_gsl_vector_complex_to_nvector, 0); + rb_define_alias(cgsl_vector_complex, "to_nvector", "to_nv"); + + rb_define_method(cgsl_vector_complex, "to_nv2", rb_gsl_vector_to_nvector_ref, 0); + rb_define_alias(cgsl_vector_complex, "to_nv_ref", "to_nv2"); + rb_define_alias(cgsl_vector_complex, "to_nvector_ref", "to_nv2"); + + rb_define_singleton_method(cgsl_vector_complex, "to_gslv", rb_gsl_na_to_gsl_vector_complex, 1); + rb_define_singleton_method(cgsl_vector_complex, "to_gv", rb_gsl_na_to_gsl_vector_complex, 1); + rb_define_singleton_method(cgsl_vector_complex, "na_to_gslv", rb_gsl_na_to_gsl_vector_complex, 1); + rb_define_singleton_method(cgsl_vector_complex, "na_to_gv", rb_gsl_na_to_gsl_vector_complex, 1); + + rb_define_singleton_method(cgsl_vector_complex, "to_gslv_view", rb_gsl_na_to_gsl_vector_complex_view, 1); + rb_define_singleton_method(cgsl_vector_complex, "to_gv_view", rb_gsl_na_to_gsl_vector_complex_view, 1); + rb_define_singleton_method(cgsl_vector_complex, "na_to_gslv_view", rb_gsl_na_to_gsl_vector_complex_view, 1); + rb_define_singleton_method(cgsl_vector_complex, "na_to_gv_view", rb_gsl_na_to_gsl_vector_complex_view, 1); + + /*****/ rb_define_method(cgsl_vector_int, "to_na", rb_gsl_vector_int_to_na, 0); rb_define_alias(cgsl_vector_int, "to_narray", "to_na"); rb_define_method(cgsl_vector_int, "to_nv", rb_gsl_vector_int_to_nvector, 0); rb_define_alias(cgsl_vector_int, "to_nvector", "to_nv"); rb_define_method(cgsl_vector_int, "to_na2", rb_gsl_vector_to_na_ref, 0); @@ -582,11 +736,11 @@ rb_define_singleton_method(cgsl_vector_int, "na_to_gv", rb_gsl_na_to_gsl_vector_int, 1); rb_define_singleton_method(cgsl_vector_int, "to_gslv_view", rb_gsl_na_to_gsl_vector_int_view, 1); rb_define_singleton_method(cgsl_vector_int, "to_gv_view", rb_gsl_na_to_gsl_vector_int_view, 1); - rb_define_singleton_method(cgsl_vector_int, "na_to_gslv_view", rb_gsl_na_to_gsl_vector_int, 1); + rb_define_singleton_method(cgsl_vector_int, "na_to_gslv_view", rb_gsl_na_to_gsl_vector_int_view, 1); rb_define_singleton_method(cgsl_vector_int, "na_to_gv_view", rb_gsl_na_to_gsl_vector_int_view, 1); rb_define_method(cNArray, "to_gv_int", rb_gsl_na_to_gsl_vector_int_method, 0); rb_define_alias(cNArray, "to_gslv_int", "to_gv_int"); rb_define_method(cNArray, "to_gv_int_view", rb_gsl_na_to_gsl_vector_int_view_method, 0); @@ -617,9 +771,12 @@ rb_define_method(cNArray, "to_gslm", rb_gsl_na_to_gsl_matrix_method, 0); rb_define_alias(cNArray, "to_gm", "to_gslm"); rb_define_method(cNArray, "to_gslm_view", rb_gsl_na_to_gsl_matrix_view_method, 0); rb_define_alias(cNArray, "to_gm_view", "to_gslm_view"); + + /*****/ + // TODO Complex matrix /*****/ rb_define_method(cgsl_matrix_int, "to_na", rb_gsl_matrix_int_to_na, 0); rb_define_alias(cgsl_matrix_int, "to_narray", "to_na"); rb_define_method(cgsl_matrix_int, "to_nm", rb_gsl_matrix_int_to_nmatrix, 0);