ext/numru/gphys/dim_op.c in gphys-1.5.2 vs ext/numru/gphys/dim_op.c in gphys-1.5.3
- old
+ new
@@ -29,10 +29,13 @@
#define ID3T(i,j,k) ((i) + (j)*n0 + (k)*n0*(n1-wlen+1))
#define ID3o(i,j,k) ((i) + (j)*n0 + (k)*n0*n1o)
#define ID3z(i,j,k) ((i) + (j)*n0 + (k)*n0*nz)
#define ID3c(i,j,k) ((i) + (j)*n0 + (k)*n0*nc)
#define ID2(i,j) ((i) + (j)*n0)
+#define ID5(i,j,k,l,m) ((i) + (j)*ni + (k)*ni*nd1 + (l)*ni*nd1*no1 + (m)*ni*nd1*no1*nd2)
+#define NMax(i,j) ( (i) >= (j) ? (i) : (j) )
+#define NMin(i,j) ( (i) <= (j) ? (i) : (j) )
static int
convol_result_length(na_shape_t n1, na_shape_t wlen, int ibc)
{
int n1o;
@@ -408,11 +411,220 @@
}
return vo;
}
+/* Running tile (2D) mean for data with missing.
+ (if no missing, you can just run the 1D running mean twice.)
+ */
static void
+running_mean_2D_miss(zi, ni,nd1,no1,nd2,no2,
+ with_w1, w1, wlen1, ibc1, with_w2, w2, wlen2, ibc2,
+ zmiss, nminvalid, zo)
+ // IN
+ double *zi;
+ na_shape_t ni,nd1,no1,nd2,no2; /* treated as 5D array (inner => outer) */
+ int with_w1, with_w2;
+ double *w1, *w2;
+ na_shape_t wlen1, wlen2;
+ int ibc1, ibc2;
+ double zmiss;
+ int nminvalid;
+ // OUT
+ double *zo; // must have been set to 0 everywhere
+{
+ na_shape_t i,j1,k1,j2,k2,m1,m2, jj1, jj2;
+ na_shape_t mstr1, mend1, mstr2, mend2;
+ na_shape_t wlhs1, wlhs2;
+ double wsum, wv;
+ int nm;
+
+ wlhs1 = (wlen1 - 1)/2; // e.g. 7->3, 6->2
+ wlhs2 = (wlen2 - 1)/2; // e.g. 7->3, 6->2
+
+ // < weighted running mean >
+
+ if ( !( ibc1 == BC_SIMPLE || ibc1 == BC_CYCLIC ) ) {
+ rb_raise(rb_eArgError,"Undefined boundary condision(1st D) (%d)", ibc1);
+ }
+ if ( !( ibc2 == BC_SIMPLE || ibc2 == BC_CYCLIC ) ) {
+ rb_raise(rb_eArgError,"Undefined boundary condision(2nd D) (%d)", ibc2);
+ }
+
+ for (k2=0; k2<no2; k2++) {
+ for (j2=0; j2<nd2; j2++) {
+ for (k1=0; k1<no1; k1++) {
+ for (j1=0; j1<nd1; j1++) {
+ for (i=0; i<ni; i++) {
+ wsum = 0.0;
+ nm = 0;
+ if (ibc1==BC_SIMPLE) {
+ mstr1 = NMax(0,wlhs1-j1);
+ mend1 = NMin(wlen1,nd1+wlhs1-j1);
+ } else { /* BC_CYCLIC */
+ mstr1 = 0;
+ mend1 = wlen1;
+ }
+ if (ibc2==BC_SIMPLE) {
+ mstr2 = NMax(0,wlhs2-j2);
+ mend2 = NMin(wlen2,nd2+wlhs2-j2);
+ } else { /* BC_CYCLIC */
+ mstr2 = 0;
+ mend2 = wlen2;
+ }
+ for (m2=mstr2; m2<mend2; m2++) {
+ for (m1=mstr1; m1<mend1; m1++) {
+ jj2 = j2-wlhs2+m2;
+ jj1 = j1-wlhs1+m1;
+ if (ibc1==BC_CYCLIC) {
+ if ( jj1 < 0 ) {
+ jj1 += nd1;
+ } else {
+ jj1 %= nd1;
+ }
+ }
+ if (ibc2==BC_CYCLIC) {
+ if ( jj2 < 0 ) {
+ jj2 += nd2;
+ } else {
+ jj2 %= nd2;
+ }
+ }
+ if ( zi[ID5(i,jj1,k1,jj2,k2)] != zmiss ) {
+ wv = 1.0;
+ if(with_w1){ wv *= w1[m1]; }
+ if(with_w2){ wv *= w2[m2]; }
+ zo[ID5(i,j1,k1,j2,k2)] +=
+ wv*zi[ID5(i,jj1,k1,jj2,k2)];
+ wsum += wv;
+ nm++;
+ }
+ }
+ }
+ if (nm >= nminvalid) {
+ zo[ID5(i,j1,k1,j2,k2)] /= wsum;
+ } else {
+ zo[ID5(i,j1,k1,j2,k2)] = zmiss;
+ }
+ }
+ }
+ }
+ }
+ }
+}
+
+// driver of running_mean_2D_miss
+// (no-missing version is not needed)
+//
+static VALUE
+running_mean_2D(obj, vi, dim1, len_or_wgt1, bc1, dim2, len_or_wgt2, bc2,
+ missv, nminvalid)
+ VALUE obj;
+ VALUE vi; // input NArray
+ VALUE dim1; // Integer
+ VALUE len_or_wgt1; // weight 1D NArray or Integer (length of uniform wgt)
+ VALUE bc1 ; // boundary condition (Integer)
+ VALUE dim2; // Integer
+ VALUE len_or_wgt2; // weight 1D NArray or Integer (length of uniform wgt)
+ VALUE bc2 ; // boundary condition (Integer)
+ VALUE missv; // missing value in vi
+ VALUE nminvalid; // miminum count of non-missing vales
+{
+ struct NARRAY *na;
+ VALUE wgt;
+ na_shape_t *shi, wlen1, wlen2;
+ int with_w1, with_w2;
+ na_shape_t ni,nd1,no1,nd2,no2;
+ double *zi, *w1, *w2, zmiss;
+ int d1, d2, i, rank, ibc1, ibc2;
+ VALUE vo;
+ na_shape_t *sho, n1o;
+ int nminvalid_; // miminum count of non-missing vales
+ double *zo;
+
+ // < process arguments >
+
+ if (!IsNArray(vi)) { rb_raise(rb_eArgError,"1st arg must be a NArray"); }
+ vi = na_cast_object(vi, NA_DFLOAT);
+ rank = NA_RANK(vi);
+ if (rank < 2) { rb_raise(rb_eArgError,"rank 1st arg must be >= 2"); }
+ zi = NA_PTR_TYPE(vi, double *);
+ GetNArray(vi, na);
+ shi = na->shape;
+
+ d1 = NUM2INT( dim1 );
+ d2 = NUM2INT( dim2 );
+ nminvalid_ = NUM2INT( nminvalid );
+
+ if (d1 >= d2) { rb_raise(rb_eArgError,"d1 < d2 is required"); }
+ if (d2 >= rank) { rb_raise(rb_eArgError,"d2 must be < rank"); }
+
+ with_w1 = IsNArray(len_or_wgt1);
+ if (with_w1) {
+ wgt = len_or_wgt1;
+ if (NA_RANK(wgt)!=1){rb_raise(rb_eArgError,"wgt1 must be a 1D NArray");}
+ wgt = na_cast_object(wgt, NA_DFLOAT);
+ wlen1 = NA_TOTAL(wgt);
+ w1 = NA_PTR_TYPE(wgt, double *);
+ } else {
+ wlen1 = NUM2INT( len_or_wgt1 );
+ w1 == NULL; /* will not be used */
+ }
+
+ with_w2 = IsNArray(len_or_wgt2);
+ if (with_w2) {
+ wgt = len_or_wgt2;
+ if (NA_RANK(wgt)!=1){rb_raise(rb_eArgError,"wgt2 must be a 1D NArray");}
+ wgt = na_cast_object(wgt, NA_DFLOAT);
+ wlen2 = NA_TOTAL(wgt);
+ w2 = NA_PTR_TYPE(wgt, double *);
+ } else {
+ wlen2 = NUM2INT( len_or_wgt2 );
+ w2 == NULL; /* will not be used */
+ }
+
+ ibc1 = NUM2INT( bc1 );
+ ibc2 = NUM2INT( bc2 );
+
+ zmiss = NUM2DBL(missv);
+
+
+ // < shape as 5D >
+
+ nd1 = shi[d1]; // length of the dim to filter
+ if (wlen1 >= nd1) {rb_raise(rb_eArgError,"filter len >= len of the dim1");}
+ nd2 = shi[d2]; // length of the dim to filter
+ if (wlen2 >= nd1) {rb_raise(rb_eArgError,"filter len >= len of the dim2");}
+
+ ni = no1 = no2 = 1;
+ for (i=0; i<d1; i++) {
+ ni *= shi[i]; // total lengths of dims before d
+ }
+ for (i=d1+1; i<d2; i++) {
+ no1 *= shi[i]; // total lengths of dims after d
+ }
+ for (i=d2+1; i<rank; i++) {
+ no2 *= shi[i]; // total lengths of dims after d
+ }
+
+ // < initialize the NArray to ruturn >
+
+ sho = shi; // Limited to shape conserved cases
+ vo = na_make_object(NA_DFLOAT, rank, sho, cNArray);
+ GetNArray(vo, na);
+ na_clear_data(na);
+ zo = NA_PTR_TYPE(vo, double *);
+
+ // < do the job >
+ running_mean_2D_miss(zi, ni,nd1,no1,nd2,no2,
+ with_w1, w1, wlen1, ibc1, with_w2, w2, wlen2, ibc2,
+ zmiss, nminvalid_, zo);
+
+ return vo;
+}
+
+static void
bin_mean_nomiss(zi,n0,n1,n2, len, zo)
// IN
double *zi;
na_shape_t n0;
na_shape_t n1;
@@ -1207,13 +1419,18 @@
// rb_require("narray"); // it does not work ??
mNumRu = rb_define_module("NumRu");
cGPhys = rb_define_class_under(mNumRu, "GPhys", rb_cObject);
rb_define_private_method(cGPhys, "c_running_mean", running_mean, -1);
+ rb_define_private_method(cGPhys, "c_running_mean_2D", running_mean_2D, 9);
+ rb_define_singleton_method(cGPhys, "c_running_mean", running_mean, -1);
+ rb_define_singleton_method(cGPhys, "c_running_mean_2D", running_mean_2D, 9);
cVArray = rb_define_class_under(mNumRu, "VArray", rb_cObject);
rb_define_private_method(cVArray, "c_bin_mean", bin_mean, -1);
rb_define_private_method(cVArray, "c_bin_sum", bin_sum, -1);
+ rb_define_singleton_method(cGPhys, "c_bin_mean", bin_mean, -1);
+ rb_define_singleton_method(cGPhys, "c_bin_sum", bin_sum, -1);
rb_define_private_method(cVArray, "c_cum_sum", cum_sum, 2);
rb_define_singleton_method(cGPhys, "c_cell_integ_irreg",
cell_integ_irreg, 6);
rb_define_singleton_method(cGPhys, "c_cum_integ_irreg", cum_integ_irreg, 6);