Class: Dobjects::Dvector

Inherits:
Object
  • Object
show all
Includes:
Enumerable
Defined in:
ext/Dobjects/Dvector/dvector.c,
lib/Dobjects/Dvector_extras.rb,
ext/Dobjects/Dvector/dvector.c,
ext/Tioga/FigureMaker/figures.c

Overview

Dvectors are a specialized implementation of one-dimensional arrays of double precision floating point numbers. They are intended for use in applications needing efficient processing of large vectors of numeric data. Essentially any of the operations you might do with a 1D Ruby Array of numbers can also be done with a Dvector. Just like Arrays, Dvector indexing starts at 0. A negative index is assumed to be relative to the end of the vector—that is, an index of -1 indicates the last element of the vector, -2 is the next to last element in the vector, and so on. Element reference and element assignment are the same as for Arrays, allowing for (start,length) or (range) as well as for simple indexing. All of the other Array operations that make sense for a 1D array of numbers are also provided for Dvectors. For example, you can “fetch” or “fill” with Dvectors, but there is no “assoc” or “compact” for them since one looks for arrays as elements and the other looks for nil elements, neither of which are found in Dvectors.

In addition to the usual Array methods, there are a variety of others that operate on the contents of the entire vector without the use of explicit iterators. These routines are crucial for efficient processing of large vectors. For example, “vec.sqrt” will create a vector of square roots of entries in vec more efficiently than the semantically equivalent form “vec.collect { |x| sqrt(x) }”.

All of the numeric methods also have ‘bang’ versions that modify the contents of the vector in place; for example, “vec.sqrt!” is more efficient than “vec = vec.sqrt”. By providing implicit iterators and in-place modification, Dvectors make it possible to operate on large vectors of doubles in Ruby at speeds closely approaching a C implementation.

As a final example, for diehards only, say we have large vectors holding values of abundances of total helium (xhe), singly ionized helium (xhe1), and doubly ionized helium (xhe2). We’re missing the values for neutral helium, but that is just what’s left of total helium after you subtract the two ionized forms, so it is easy to compute it. Finally, we need to calculate the log of the abundance of neutral helium and store it in another vector (log_xhe0). If we don’t care about creating work for the garbage collector, we can simply do this

log_xhe0 = (xhe - xhe1 - xhe2).log10

This works, but it creates multiple temporary vectors for intermediate results. If, like me, you’re compulsive about efficiency, you can do the whole thing with no garbage created at all:

log_xhe0 = Dvector.new
log_xhe0.replace(xhe).sub!(xhe1).sub!(xhe2).log10!

This copies xhe to log_xhe0, subtracts xhe1 and xhe2 from log_xhe0 in place, and then takes the log, also in place. It’s not pretty, but it is efficient – use if needed.

Please report problems with the Dvector extension to the tioga-users at rubyforge.org mailing list. [Note: for N-dimensional arrays or arrays of complex numbers or integers as well as doubles, along with a variety of matrix operations, check out the NArray extension.]

Dvector now also prides itselfs with a _dump and a _load function, which means you can Marshal them.

Constant Summary collapse

FANCY_READ_DEFAULTS =

Dvector.fancy_read’s defaults options. See that function for more details

{ 
  'sep' => /\s+/,
  'comments' => /^\s*\#/,
  'skip_first' => 0,
  'index_col' => false,
  'headers' => nil, # not used for now.
  'default' => 0.0/0.0, # defaults to NaN
  'initial_size' => 5001,
  'remove_space' => true ,# removes spaces at the beginning of the lines
  'last_col' => -1,       # Read all columns
  'text_columns' => [],   # Not a single column is text
}
WRITE_DEFAULTS = 
{
  'sep' => "\t",
  'write_mode' => "a",
}

Class Method Summary collapse

Instance Method Summary collapse

Constructor Details

#initialize ⇒ Object

Class Method Details

.[](*args) ⇒ Object

Returns a new Dvector populated with the given objects.

Dvector.[]( 1, 2, 3, 4 )      -> a_dvector
Dvector[ 1, 2, 3, 4 ]         -> a_dvector



146
147
148
149
150
151
152
153
154
 
# File 'ext/Dobjects/Dvector/dvector.c', line 146

static VALUE dvector_s_create(int argc, VALUE *argv, VALUE klass) {
   VALUE ary = make_new_dvector(klass, argc, argc);
   Dvector *d = Get_Dvector(ary);
   if (argc < 0) {
      rb_raise(rb_eArgError, "negative number of arguments");
   }
   ary_to_dvector(argv, argc, d->ptr);
   return ary;
}

._load ⇒ Object

.compute_formula(formula, a, modules = []) ⇒ Object

This function is a rudimentary formula computing stuff. Give it a text formula and an array of Dvectors (a), and it returns a Dvector with the result. The formula should contain the following;

column

represents the current element of the n th Dvector of the array

This is just a try, and should be implemented in C rather than in Ruby. But if you’re looking for simplicity, here you go ;-) !

modules are the modules you would wish the evaluator to include. This feature enables one to make sure custom functions are included



226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
 
# File 'lib/Dobjects/Dvector_extras.rb', line 226

def Dvector.compute_formula(formula, a, modules = []) # :doc:

  evaluator = MathEvaluator.new(formula, "column", modules)
  # if we reach this place, it means that there a no big syntax errors ;-)
  
  # we now need to inspect the array given, and make sure that there is
  # and transform it into a clean stuff (an array with only Dvectors and
  # nil elements).
  
  target = []
  last = nil
  a.each { |elem| 
    if elem.is_a? Dvector
      target << elem
      last = elem
    else
      target << nil
    end
  }
  
  raise "No Dvector found" unless last
  
  # we check all the vectors have the same length
  target.each {|x| 
    if x && x.length != last.length
      raise "Dvectors should have all the same length !" 
    end
  }
  
  res = Dvector.new
  
  last.each_index { |i|
    args = target.collect { |val| 
      if val
        val[i]
      else 
        nil
      end
    }
    # we add the index at the beginning:
    #         args.unshift(i) 
    # Commented out for simplicity
    
    # then we call the block:
    elem = evaluator.compute(args)
    res << elem
  }
  
  return res
end

.create_pm_cubic_interpolant ⇒ Object

.create_spline_interpolant ⇒ Object

.fancy_read(stream, cols = nil, opts = {}) ⇒ Object

This function is a wrapper for #fast_fancy_read that reflects the look-and-feel of #old_fancy_read

Raises:

  • (ArgumentError) 


98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
 
# File 'lib/Dobjects/Dvector_extras.rb', line 98

def Dvector.fancy_read(stream, cols = nil, opts = {}) # :doc:
  o = FANCY_READ_DEFAULTS.dup
  o.update(opts)

  if stream.is_a?(String)
    stream = File.open(stream)
  end
  raise ArgumentError.new("'stream' should have a gets method") unless 
    stream.respond_to? :gets
  
  o['sep'] = Regexp.new(o['sep']) unless o['sep'].is_a? Regexp
  
  res = Dvector.fast_fancy_read(stream, o)

  # Adding the index columns if necessary
  if o["index_col"] 
    res.unshift(Dvector.new(res[0].length) { |i| i})
  end

  if cols
    return cols.map {|i| res[i] }
  else
    return res
  end
end

.fast_fancy_read(stream, options) ⇒ Object

:call-seq: Dvector.fast_fancy_read(stream, options) => Array_of_Dvectors

Reads data from an IO stream (or anything that supports a gets method) and separate it into columns of data according to the options, a hash holding the following elements (compulsory, but you can use FANCY_READ_DEFAULTS): ‘sep’: a regular expression that separate the entries ‘comments’: any line matching this will be skipped ‘skip_first’: skips that many lines before reading anything ‘index_col’: if true, the first column returned contains the

number of the line read

‘remove_space’: whether to remove spaces at the beginning of a line. ‘comment_out’: this should be an array into which the comments

will be dumped one by one.

‘default’: what to put when nothing was found but a number must be used ‘last_col’: when this is specified, it represents the last column which

is read and parsed as numbers (0-based, so the 3rd column is 2).
The remaining is returned as text in one n+1 column

‘text_columns’: if provided, it is an array of integers containing the

0-based indices of columns to be parsed as text rather than numbers.

In addition to these options that control the output, here are a few others to tune memory allocation; these can strongly improve the performance (or make it worse if you wish):

‘initial_size’: the initial size of the memory buffers: if there

are not more lines than that, no additional memory allocation/copy
occurs.



5666
5667
5668
5669
5670
5671
5672
5673
5674
5675
5676
5677
5678
5679
5680
5681
5682
5683
5684
5685
5686
5687
5688
5689
5690
5691
5692
5693
5694
5695
5696
5697
5698
5699
5700
5701
5702
5703
5704
5705
5706
5707
5708
5709
5710
5711
5712
5713
5714
5715
5716
5717
5718
5719
5720
5721
5722
5723
5724
5725
5726
5727
5728
5729
5730
5731
5732
5733
5734
5735
5736
5737
5738
5739
5740
5741
5742
5743
5744
5745
5746
5747
5748
5749
5750
5751
5752
5753
5754
5755
5756
5757
5758
5759
5760
5761
5762
5763
5764
5765
5766
5767
5768
5769
5770
5771
5772
5773
5774
5775
5776
5777
5778
5779
5780
5781
5782
5783
5784
5785
5786
5787
5788
5789
5790
5791
5792
5793
5794
5795
5796
5797
5798
5799
5800
5801
5802
5803
5804
5805
5806
5807
5808
5809
5810
5811
5812
5813
5814
5815
5816
5817
5818
5819
5820
5821
5822
5823
5824
5825
5826
5827
5828
5829
5830
5831
5832
5833
5834
5835
5836
5837
5838
5839
5840
5841
5842
5843
5844
5845
5846
5847
5848
5849
5850
5851
5852
5853
5854
5855
5856
5857
5858
5859
5860
5861
5862
5863
5864
5865
5866
5867
5868
5869
5870
5871
5872
5873
5874
5875
5876
5877
5878
5879
5880
5881
5882
5883
5884
5885
5886
5887
5888
5889
 
# File 'ext/Dobjects/Dvector/dvector.c', line 5666

static VALUE dvector_fast_fancy_read(VALUE self, VALUE stream, VALUE options)
{
  /* First, we read up options: */
  double def = rb_num2dbl(rb_hash_aref(options, 
				       rb_str_new2("default")));
  int remove_space = RTEST(rb_hash_aref(options, 
					rb_str_new2("remove_space")));
//  int index_col = RTEST(rb_hash_aref(options, 
//				     rb_str_new2("index_col")));
  long skip_first = FIX2LONG(rb_hash_aref(options, 
					  rb_str_new2("skip_first")));
  VALUE sep = rb_hash_aref(options, rb_str_new2("sep"));
  VALUE comments = rb_hash_aref(options, rb_str_new2("comments"));
  VALUE comment_out = rb_hash_aref(options, rb_str_new2("comment_out"));

  /* Elements after that many columns  */
  VALUE lc = rb_hash_aref(options, rb_str_new2("last_col"));
  long last_col = RTEST(lc) ? FIX2LONG(lc) : -1;
  VALUE text_columns = rb_hash_aref(options, rb_str_new2("text_columns"));

  /* Then, some various variables: */
  VALUE line;

  ID chomp_id = rb_intern("chomp!");
  ID gets_id = rb_intern("gets");
  ID max_id = rb_intern("max");
  ID size_id = rb_intern("size");
  
  /* We compute the maximum number of text columns */
  long last_text_col = last_col+1;
  VALUE mx = RTEST(text_columns) ? rb_funcall(text_columns, max_id, 0) : Qnil;
  if(RTEST(mx) && last_text_col < 0) { /* Only taking the max into
                                          account if the last col
                                          stuff is not on */
    long d = FIX2LONG(mx);
    last_text_col = d;
  }


  /* array of Ruby arrays containing the text objects of interest */
  VALUE * text_cols = NULL;

  /*
    Handling of text columns.

    The number and position of text columns has to be known in
    advance. For each of those, the value of text_columns isn't Qnil,
    and the corresponding column is NULL.

   */
  if(last_text_col >= 0) {
    text_cols = ALLOC_N(VALUE, last_text_col + 1);
    int i;
    for(i = 0; i < last_text_col + 1; i++)
      text_cols[i] = Qnil;
    if(last_col >= 0) {
      text_cols[last_col+1] = marked_array();
    }
    if(RTEST(mx)) {
      /* Todo */
      int sz = 
#ifdef RARRAY_LENINT
      RARRAY_LENINT(text_columns);
#else
      RARRAY_LEN(text_columns);
#endif
      int i;
      for(i = 0; i <  sz; i++) {
        long idx = FIX2LONG(rb_ary_entry(text_columns, i));
        if(idx >= 0 && (last_col < 0 || idx < last_col)) {
          text_cols[idx] = marked_array();
        }
      }
    }
  }




  long line_number = 0;

  /* 
     Now come the fun part - rudimentary vectors management

     TODO: if the stream provides functionality to get its total size,
     it could be interesting to estimate the total number of lines
     based on some small heuristics
   */
  int nb_vectors = 0;		/* The number of vectors currently created */
  int current_size = 10;	/* The number of slots available */
  double ** vectors = ALLOC_N(double *, current_size);
  long index = 0;		/* The current index in the vectors */
  /* The size available in the vectors */
  int allocated_size = 
    FIX2LONG(rb_hash_aref(options, rb_str_new2("initial_size"))); 


  int i;
  for(i = 0; i < current_size; i++)
    vectors[i] = NULL;

  /* The return value */
  VALUE ary;

  /* We use a real gets so we can also use StringIO, for instance */
  while(RTEST(line = rb_funcall(stream, gets_id, 0))) {
    VALUE pre, post, match;
    const char * line_ptr;
    int col = 0;
    line_number++;
    /* Whether we should skip the line... */
    if(skip_first >= line_number)
      continue;

    /* We check for a blank line using isspace: */
    line_ptr = StringValueCStr(line);
    while(line_ptr && *line_ptr) {
      if(! isspace(*line_ptr))
	break;
      line_ptr++;
    }
    if(! *line_ptr)
      continue;			/* We found a blank line  */
    if(remove_space)		/* We replace the contents of the line  */
      line = rb_str_new2(line_ptr);

    /* ... or a comment line */
    if(RTEST(comments) && RTEST(rb_reg_match(comments, line))) {
      if(RTEST(comment_out))
	rb_ary_push(comment_out, line);
      continue;
    }

    /* Then, we remove the newline: */
    post = line;
    rb_funcall(post, chomp_id, 0);

    /* We iterate over the different portions between
       matches
    */
    while(RTEST(post)) {
      const char * a;
      char * b;
      if(RTEST(rb_reg_match(sep, post))) {
	match = rb_gv_get("$~");
	pre = rb_reg_match_pre(match);
	post = rb_reg_match_post(match);
      }
      else {
	pre = post;
	post = Qnil;
      }
      if(text_cols && col <= last_text_col && RTEST(text_cols[col])) {
        rb_ary_push(text_cols[col], pre);
        if(col >= nb_vectors) {
          nb_vectors ++;
          if(col < current_size)
            vectors[col] = NULL;
        }
      }
      else {
        a = StringValueCStr(pre);
        double c = strtod(a, &b);
        if(b == a) 
          c = def;
        if(col >= nb_vectors) {
          /* We need to create a new vector */
          if(col >= current_size) { /* Increase the available size */
            current_size = col + 5;
            REALLOC_N(vectors, double * , current_size);
          }
          for(; nb_vectors <= col; nb_vectors++)
            vectors[nb_vectors] = NULL; /* default to NULL */
          
          double * vals = vectors[col] = ALLOC_N(double, allocated_size);
          /* Filling it with the default value */
          for(i = 0; i < index; i++) {
            vals[i] = def;
          }
        }
        vectors[col][index] = c;
      }
      col++;
      if(last_col >= 0 && col > last_col) {
        rb_ary_push(text_cols[last_col + 1], post);
        nb_vectors = col + 1;
        col++;
        break;
      }
    }
    /* Now, we finish the line */
    for(; col < nb_vectors; col++) {
      if(text_cols && col <= last_text_col && RTEST(text_cols[col]))
        rb_ary_push(text_cols[col], Qnil);
      else
        vectors[col][index] = def;
    }
    index++;
    /* Now, we reallocate memory if necessary */
    if(index >= allocated_size) {
      allocated_size *= 2;	/* We double the size */
      for(col = 0; col < nb_vectors; col++) {
        if(col < current_size && vectors[col])
          REALLOC_N(vectors[col], double, allocated_size);
      }
    }
  }
  /* Now, we make up the array */
  ary = rb_ary_new();
  for(i = 0; i < nb_vectors; i++) {
    /* We create a vector */
    if(text_cols && i <= last_text_col && RTEST(text_cols[i]))
      rb_ary_store(ary, i, text_cols[i]);
    else {
      rb_ary_store(ary, i, make_dvector_from_data(cDvector, index, vectors[i]));
      /* And free the memory */
      free(vectors[i]);
    }
  }
  free(vectors);
  if(text_cols)
    free(text_cols);
  return ary;
}

.from_na(na) ⇒ Object

thanks to Dave MacMahon for from_na and to_na Create a Dvector from an NArray.



67
68
69
 
# File 'lib/Dobjects/Dvector_extras.rb', line 67

def Dvector.from_na(na)
  _load([1, na.length, na.to_s].pack('CIa*'))
end

.is_a_dvector ⇒ Object

.linear_interpolate ⇒ Object

.max_of_many ⇒ Object

.min_of_many ⇒ Object

.old_fancy_read(stream, cols = nil, opts = {}) ⇒ Object

This function reads in stream (can an IO object or a String, in which case it represents the name of a file to be opened) the columns specified by cols and returns them. column 0 is the first column. If cols is nil, then fancy_read attempts to find all the columns in the file, while filling absent data with NaN.

opts is a hash for tuning the behavior of the reading. It can hold the following keys:

‘sep’

the record separator

‘comments’

a regular expression matching comment lines

‘skip_first’

how many lines to skip at the beginning of the file,

‘default’

the value taken for missing elements

‘index_col’

if set to true, the first column contains the indices of the elements (starting from 0 for first and so on)

Raises:

  • (ArgumentError) 


140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
 
# File 'lib/Dobjects/Dvector_extras.rb', line 140

def Dvector.old_fancy_read(stream, cols = nil, opts = {}) # :doc:
  # first, we turn the stream into a real IO stream
  if stream.is_a?(String)
    stream = File.open(stream)
  end
  raise ArgumentError.new("'stream' should have a gets method") unless 
    stream.respond_to? :gets
  
  # we take default options and override them with opts
  o = FANCY_READ_DEFAULTS.merge(opts)
  
  # strip off the first lines.
  while o["skip_first"] > 0
    stream.gets
    o["skip_first"] -= 1
  end

  # then, parsing the lines. We store the results in an array. We read up
  # all columns, regardless of what is asked (it doesn't slow that much
  # the process -- or does it ?)
  
  columns = []
  line_number = 0 # the number of the significant lines read so far
  
  while line = stream.gets
    next if line =~ o["comments"]
    next if line =~ /^\s*$/ # skip empty lines
    if o["remove_space"]
      line.gsub!(/^\s+/,'')
    end
  
    elements = line.split(o["sep"])
    # now, the fun: the actual reading.
    # we first turn this elements into floats:
    numbers = elements.collect do |s| 
      begin 
        a = Float(s)
      rescue
        o["default"]
      end
    end

    if numbers.size < columns.size 
      # we pad it with default values
      while numbers.size < columns.size
        numbers << o["default"]
      end
    else
      # in that case, we need to create new Dvectors to match the
      # size of numbers
      while columns.size < numbers.size
        columns << Dvector.new(line_number, o["default"]) 
      end
    end
    # now, we should have the same number of elements both
    # in numbers and in columns
    columns.size.times do |i|
      columns[i] << numbers[i]
    end
    # and it's done ;-) !

    line_number += 1
  end
  # Adding the index columns if necessary
  if o["index_col"] 
    columns.unshift(Dvector.new(columns[0].length) { |i| i})
  end

  return columns unless cols
  return cols.collect { |i| 
    columns[i]
  }
end

.pm_cubic_interpolate ⇒ Object

.read ⇒ Object

.read_columns ⇒ Object

.read_row ⇒ Object

.read_rows ⇒ Object

.spline_interpolate ⇒ Object

.write(file, cols, options = {}) ⇒ Object

Writes an array of Dvectors into a text file



283
284
285
286
287
288
289
290
291
292
 
# File 'lib/Dobjects/Dvector_extras.rb', line 283

def Dvector.write(file, cols, options = {})
  ops = WRITE_DEFAULTS.update(options)
  if file.is_a?(String)
    file = File.open(file, ops["write_mode"])
  end 
  nb = cols.map {|d| d.size}.max # The number of lines
  nb.times do |i|
    file.puts(cols.map {|d| d[i].to_s }.join(ops["sep"]))
  end
end

Instance Method Details

#<< ⇒ Object

#<=> ⇒ Object

#[] ⇒ Object Also known as: slice

#[]= ⇒ Object

#_dump ⇒ Object

dvector marshalling

#abs ⇒ Object

numeric methods

#abs! ⇒ Object

#acos ⇒ Object

standard math functions

#acos! ⇒ Object

#acosh ⇒ Object

#acosh! ⇒ Object

#add ⇒ Object Also known as: +, plus

#add! ⇒ Object Also known as: plus!

#as_exponent_of ⇒ Object

#as_exponent_of! ⇒ Object

#asin ⇒ Object

#asin! ⇒ Object

#asinh ⇒ Object

#asinh! ⇒ Object

#at ⇒ Object

#atan ⇒ Object

#atan! ⇒ Object

#atan2 ⇒ Object

#atan2! ⇒ Object

#atanh ⇒ Object

#atanh! ⇒ Object

#bounds ⇒ Object

Returns the boundaries of a Dvector, that is [min, max]. It ignores

NaN and will complain if the Dvector contains only NaN.

 v = Dvector[0.0/0.0, 0.0/0.0, 1,2,4,5,9,0.0/0.0,0.1]
 v.bounds -> [0.1, 9]



5539
5540
5541
5542
5543
5544
5545
5546
5547
5548
5549
5550
5551
5552
5553
5554
5555
5556
5557
5558
5559
5560
5561
5562
5563
5564
5565
5566
5567
5568
5569
5570
5571
5572
 
# File 'ext/Dobjects/Dvector/dvector.c', line 5539

static VALUE dvector_bounds(VALUE self)
{
  double min, max;
  VALUE ret;
  long len;
  double * data = Dvector_Data_for_Read(self, &len);
  /* skip all NaNs at the beginning */
  while(len-- > 0)
    if(!isnan(*data++))
       break;
  if(len>=0)
    {
      min = max = *(data-1);
      while(len-- > 0)
	{
	  if(! isnan(*data))
	    {
	      if(*data < min)
		min = *data;
	      if(*data > max)
		max = *data;
	    }
	  data++;
	}
      ret = rb_ary_new2(2);
      rb_ary_store(ret, 0, rb_float_new(min));
      rb_ary_store(ret, 1, rb_float_new(max));
      return ret;
    }
  else
    rb_raise(rb_eRuntimeError, 
	     "bounds called on an array containing only NaN");
  return Qnil;
}

#ceil ⇒ Object

#ceil! ⇒ Object

#clean? ⇒ Boolean

Returns:

  • (Boolean)

#clear ⇒ Object

#collect ⇒ Object Also known as: map

#collect! ⇒ Object Also known as: map!

#collect2 ⇒ Object Also known as: map2

#collect2! ⇒ Object Also known as: map2!

#concat ⇒ Object

#convolve(kernel, middle) ⇒ Object

:call-seq:

vector.convolve(kernel, middle)

convolve applies a simple convolution to the vector using kernel centered at the point middle. (0 is the leftmost point of the kernel).



5582
5583
5584
5585
5586
5587
5588
5589
5590
5591
5592
5593
5594
5595
5596
5597
5598
5599
5600
5601
5602
5603
5604
5605
5606
5607
5608
5609
5610
5611
5612
5613
5614
5615
5616
5617
5618
5619
5620
5621
5622
5623
5624
 
# File 'ext/Dobjects/Dvector/dvector.c', line 5582

static VALUE dvector_convolve(VALUE self, VALUE kernel, VALUE middle)
{
  long len;
  const double * values = Dvector_Data_for_Read(self, &len);
  VALUE retval = dvector_new2(len,len);
  double * ret = Dvector_Data_for_Write(retval,NULL);
  long kernel_len;
  const double * ker = Dvector_Data_for_Read(kernel, &kernel_len);
  /* I guess */
  long mid = NUM2LONG(middle);
  if(mid > kernel_len)
    rb_raise(rb_eArgError, "middle should be within kernel's range");
  else
    {
      long i,j,k;
      for(i = 0; i < len; i++) 
	{
	  double sum = 0, k_sum = 0;
	  for(j = 0; j < kernel_len; j++) 
	    {
	      /* check that we are within the vector */
	      k = i - mid + j; 	/* The current index inside the vector */
	      /* This code is equivalent to saying that the vector is
		 prolongated until infinity with values at the boundaries
		 -> no, obnoxious, I think. Simply don't take care
		 of these points
		 -> yes, finally ?
	      */
	      if( k < 0)
/* 		continue; */
		k = 0;
	      if( k >= len)
/* 		continue; */
		k = len - 1;
	      sum += ker[j] * values[k];
	      k_sum += ker[j];
	    }
	  sum/= k_sum;
	  ret[i] = sum;
	}
    }
  return retval;
}

#cos ⇒ Object

#cos! ⇒ Object

#cosh ⇒ Object

#cosh! ⇒ Object

#delete ⇒ Object

#delete_at ⇒ Object

#delete_if ⇒ Object

#dirty= ⇒ Object

#dirty? ⇒ Boolean Also known as: dirty

dirty flag related methods

Returns:

  • (Boolean)

#div ⇒ Object Also known as: /

#div! ⇒ Object

#dot ⇒ Object

#dup ⇒ Object

#each ⇒ Object

#each2 ⇒ Object

#each2_with_index ⇒ Object

#each3 ⇒ Object

#each3_with_index ⇒ Object

#each_index ⇒ Object

#each_with_index ⇒ Object

#empty? ⇒ Boolean

Returns:

  • (Boolean)

#eql? ⇒ Boolean Also known as: ==

Returns:

  • (Boolean)

#exp ⇒ Object

#exp! ⇒ Object

#exp10 ⇒ Object

#exp10! ⇒ Object

#extrema(*args) ⇒ Object

Returns a list of local extrema of the vector, organized thus:

[ [:min, idmin1], [:max, idmax1], ...]

The values are pushed in the order in which they are found. It works thus: it scans the vector and looks around the current point in a given window. If the current point is the maximum or the minimum, it is considered as a local maximum/minimum. Control over which extrema are included is given to the used through threshold mechanisms.

The options hash controls how the peaks are detected: window: the number of elements on which we look on

both sides (default 5, ie the local maximum is over 11 points)

threshold: the minimum amplitude the extrema must have to

be considered (default 0)

dthreshold: how much over/under the average an extremum must be

(default 0)

or: whether the threshold and dthreshold tests are both

necessary or if only one is (default false: both tests are
necessary)

Note:* beware of NANs ! They will screw up peak detection, as they are neither bigger nor smaller than anything…



5916
5917
5918
5919
5920
5921
5922
5923
5924
5925
5926
5927
5928
5929
5930
5931
5932
5933
5934
5935
5936
5937
5938
5939
5940
5941
5942
5943
5944
5945
5946
5947
5948
5949
5950
5951
5952
5953
5954
5955
5956
5957
5958
5959
5960
5961
5962
5963
5964
5965
5966
5967
5968
5969
5970
5971
5972
5973
5974
5975
5976
5977
5978
5979
5980
5981
5982
5983
5984
5985
5986
5987
5988
5989
5990
5991
5992
5993
5994
5995
5996
5997
5998
5999
6000
6001
6002
6003
6004
6005
6006
6007
6008
6009
6010
 
# File 'ext/Dobjects/Dvector/dvector.c', line 5916

static VALUE dvector_extrema(int argc, VALUE *argv, VALUE self)
{
  long window = 5;
  double threshold = 0;
  double dthreshold = 0;
  int inclusive = 1;
  
  if(argc == 1) {
    VALUE t;
    t = rb_hash_aref(argv[0], rb_str_new2("window"));
    if(RTEST(t)) {
      window = FIX2LONG(t);
    }
    t = rb_hash_aref(argv[0], rb_str_new2("threshold"));
    if(RTEST(t)) {
      threshold = rb_num2dbl(t);
    }
    t = rb_hash_aref(argv[0], rb_str_new2("dthreshold"));
    if(RTEST(t)) {
      dthreshold = rb_num2dbl(t);
    }
    
    t = rb_hash_aref(argv[0], rb_str_new2("or"));
    inclusive = ! RTEST(t);
  } else if(argc > 1)
    rb_raise(rb_eArgError, "Dvector.extrema only takes 0 or 1 argument");

  /* Handling of the vector */
  long len, i,j;
  double * data = Dvector_Data_for_Read(self, &len);
  VALUE s_min = ID2SYM(rb_intern("min"));
  VALUE s_max = ID2SYM(rb_intern("max"));

  

  VALUE ret = rb_ary_new();
		       
  for(i = 0; i < len; i++) {

    /* This is stupid and will need decent optimization when I have
       time */
    long first = i > window ? i - window : 0;
    double cur_min = data[first];
    long cur_min_idx = first;
    double cur_max = data[first];
    long cur_max_idx = first;
    double average = 0;
    long nb = 0;
    
    for(j = first; (j < i+window) && (j < len); j++,nb++) {
      average += data[j];
      if(data[j] <= cur_min) {
	cur_min = data[j];
	cur_min_idx = j;
      }
      if(data[j] >= cur_max) {
	cur_max = data[j];
	cur_max_idx = j;
      }
    }
    average /= nb;

    if(cur_min_idx == i) {
      /* This is a potential minimum */
      if((inclusive && 
	  (fabs(cur_min) >= threshold) && 
	  (fabs(cur_min - average) >= dthreshold))
	 || (!inclusive && 
	     ((fabs(cur_min) >= threshold) ||
	      (fabs(cur_min - average) >= dthreshold))
	     )) {
	VALUE min = rb_ary_new();
	rb_ary_push(min, s_min);
	rb_ary_push(min, LONG2FIX(i));
	rb_ary_push(ret, min);
      }
    }
    else if(cur_max_idx == i) {
      /* A potential maximum */
      if((inclusive && 
	  (fabs(cur_max) >= threshold) && 
	  (fabs(cur_max - average) >= dthreshold))
	 || (!inclusive && 
	     ((fabs(cur_max) >= threshold) ||
	      (fabs(cur_max - average) >= dthreshold))
	     )) {
	VALUE max = rb_ary_new();
	rb_ary_push(max, s_max);
	rb_ary_push(max, LONG2FIX(i));
	rb_ary_push(ret, max);
      }
    }
  }
  return ret;
}

#fetch ⇒ Object

#fft! ⇒ Object

Performs an in-place Fourier transform of the vector. The results is stored in the so-called “half-complex” format (see www.fftw.org/fftw3_doc/The-Halfcomplex_002dformat-DFT.html for more information).



6023
6024
6025
6026
6027
6028
6029
6030
6031
6032
 
# File 'ext/Dobjects/Dvector/dvector.c', line 6023

static VALUE dvector_fft(VALUE self)
{
  long len;
  double * values = Dvector_Data_for_Write(self, &len);
  fftw_plan plan = fftw_plan_r2r_1d(len, values, values,
				    FFTW_R2HC, FFTW_ESTIMATE);
  fftw_execute(plan);
  fftw_destroy_plan(plan);
  return self;
}

#fft_conj! ⇒ Object

Converts the FFTed data in the complex conjugate



6096
6097
6098
6099
6100
6101
6102
6103
6104
6105
6106
 
# File 'ext/Dobjects/Dvector/dvector.c', line 6096

static VALUE dvector_fft_conj(VALUE self)
{
  long len;
  double * v1 = Dvector_Data_for_Write(self, &len);
  double * img;
  long i;
  for(i = 1, img = v1 + len-1; i < (len+1)/2;
      i++, img--)
    *img = -*img;
  return self;
}

#fft_mul!(m) ⇒ Object

Multiplies the FFTed data held in the vector by another vector. The behaviour depends on the size of the target vector:

if it is the same size, it is assumed to be FFTed data if it is the same size of a power spectrum, then it is assumed that it

is multiplication by real values

anything else won’t make this function happy.

As a side note, if you only want multiplication by a scalar, the standard #mul! should be what you look for.



6121
6122
6123
6124
6125
6126
6127
6128
6129
6130
6131
6132
6133
6134
6135
6136
6137
6138
6139
6140
6141
6142
6143
6144
6145
6146
6147
6148
6149
6150
6151
6152
6153
6154
6155
6156
6157
6158
6159
6160
6161
6162
6163
6164
6165
6166
6167
6168
 
# File 'ext/Dobjects/Dvector/dvector.c', line 6121

static VALUE dvector_fft_mul(VALUE self, VALUE m)
{
  long len;
  double * v1 = Dvector_Data_for_Write(self, &len);
  long len2;
  const double * v2 = Dvector_Data_for_Write(m, &len2);
  if(len2 == len) {		/* Full complex multiplication */
    const double * m_img;
    const double * m_real;
    double * v_img;
    double * v_real;
    long i;
    /* First, special cases */
    v1[0] *= v2[0];
    if(len % 2 == 0)
      v1[len/2] *= v2[len/2];
    
    for(i = 1, m_real = v2 + 1, m_img = v2 + len-1,
	  v_real = v1 + 1, v_img = v1 + len-1; i < (len+1)/2;
	i++, m_real++, v_real++, m_img--, v_img--) {
      double r = *m_real * *v_real - *m_img * *v_img;
      *v_img = *m_real * *v_img + *v_real * *m_img;
      *v_real = r;
    }
    return self;
  }
  else if(len2 == len/2+1) {		/* Complex * real*/
    const double * val;
    double * v_img;
    double * v_real;
    long i;
    /* First, special cases */
    v1[0] *= v2[0];
    if(len % 2 == 0)
      v1[len/2] *= v2[len/2];
    
    for(i = 1, val = v2 + 1,
	  v_real = v1 + 1, v_img = v1 + len-1; i < (len+1)/2;
	i++, val++, v_real++, v_img--) {
      *v_real *= *val;
      *v_img *= *val;
    }
    return self;
  }
  else {
    rb_raise(rb_eArgError, "incorrect Dvector size for fft_mul!");
  }
}

#fft_spectrum ⇒ Object

Returns the power spectra of the ffted data (ie the square of the norm of the complex fourier coefficients.

The returned value is a new Dvector of size about two times smaller than the original (precisely size/2 + 1)

For some reasons, convolutions don’t work for now.



6068
6069
6070
6071
6072
6073
6074
6075
6076
6077
6078
6079
6080
6081
6082
6083
6084
6085
6086
6087
6088
6089
6090
6091
 
# File 'ext/Dobjects/Dvector/dvector.c', line 6068

static VALUE dvector_fft_spectrum(VALUE self)
{
  long len;
  const double * values = Dvector_Data_for_Read(self, &len);
  /* First compute the size of the target: */
  long target_size = len/2+1;
  long i;
  VALUE retval = dvector_new2(target_size,target_size);
  double * ret = Dvector_Data_for_Write(retval,NULL);

  /* Pointer to real and imaginary parts */
  const double * real;
  const double * img;
  ret[0] = values[0] * values[0];


  /* The Nyquist frequency */
  if(len % 2 == 0)
    ret[target_size - 1] = values[target_size-1] * values[target_size-1];
  for(i = 1, real = values + 1, img = values + len-1; i < len/2;
      i++, real++, img--)
    ret[i] = *real * *real + *img * *img;
  return retval;
}

#fill ⇒ Object

#first ⇒ Object

#floor ⇒ Object

#floor! ⇒ Object

#freeze ⇒ Object

#frozen? ⇒ Boolean

Returns:

  • (Boolean)

#include? ⇒ Boolean

Returns:

  • (Boolean)

#index ⇒ Object

#initialize_copy ⇒ Object

#insert ⇒ Object

#inv ⇒ Object

#inv! ⇒ Object

#join ⇒ Object

#last ⇒ Object

#length ⇒ Object Also known as: size, nitems

#log ⇒ Object

#log! ⇒ Object

#log10 ⇒ Object

#log10! ⇒ Object

#make_bezier_control_points_for_cubic_in_x ⇒ Object

#max ⇒ Object

#max_lt ⇒ Object

#min ⇒ Object

#min_gt ⇒ Object

#modulo ⇒ Object Also known as: mod, %

#modulo! ⇒ Object Also known as: mod!

#mul ⇒ Object Also known as: *, times

#mul! ⇒ Object Also known as: times!

#neg ⇒ Object Also known as: -@

nonstandard math functions

#neg! ⇒ Object

#pop ⇒ Object

#pow ⇒ Object Also known as: raised_to, **

#pow! ⇒ Object Also known as: raised_to!

#prune ⇒ Object

#prune! ⇒ Object

#push ⇒ Object

#reject ⇒ Object

#reject! ⇒ Object

#remainder ⇒ Object

#remainder! ⇒ Object

#replace ⇒ Object

#resize ⇒ Object

#reverse ⇒ Object

#reverse! ⇒ Object

#reverse_each ⇒ Object

#reverse_each2 ⇒ Object

#reverse_each2_with_index ⇒ Object

#reverse_each3 ⇒ Object

#reverse_each3_with_index ⇒ Object

#reverse_each_index ⇒ Object

#reverse_each_with_index ⇒ Object

#rfft! ⇒ Object

Performs a reverse in-place Fourier transform of the vector. The original data must have been stored in the so called “half-complex” format (see #fft!).



6040
6041
6042
6043
6044
6045
6046
6047
6048
6049
 
# File 'ext/Dobjects/Dvector/dvector.c', line 6040

static VALUE dvector_rfft(VALUE self)
{
  long len;
  double * values = Dvector_Data_for_Write(self, &len);
  fftw_plan plan = fftw_plan_r2r_1d(len, values, values,
				    FFTW_HC2R, FFTW_ESTIMATE);
  fftw_execute(plan);
  fftw_destroy_plan(plan);
  return self;
}

#rindex ⇒ Object

#round ⇒ Object

#round! ⇒ Object

#safe_acos ⇒ Object

#safe_acos! ⇒ Object

#safe_asin ⇒ Object

#safe_asin! ⇒ Object

#safe_inv ⇒ Object

#safe_inv! ⇒ Object

#safe_log ⇒ Object

#safe_log! ⇒ Object

#safe_log10 ⇒ Object

#safe_log10! ⇒ Object

#safe_sqrt ⇒ Object

#safe_sqrt! ⇒ Object

#select ⇒ Object

#set ⇒ Object

#shift ⇒ Object

#sin ⇒ Object

#sin! ⇒ Object

#sinh ⇒ Object

#sinh! ⇒ Object

#slice! ⇒ Object

#sort ⇒ Object

#sort! ⇒ Object

#sqrt ⇒ Object

#sqrt! ⇒ Object

#sub ⇒ Object Also known as: -, minus

#sub! ⇒ Object Also known as: minus!

#sum ⇒ Object

#tan ⇒ Object

#tan! ⇒ Object

#tanh ⇒ Object

#tanh! ⇒ Object

#to_a ⇒ Object Also known as: to_ary

#to_dvector ⇒ Object 



77
78
79
 
# File 'lib/Dobjects/Dvector_extras.rb', line 77

def to_dvector
  self
end

#to_na ⇒ Object

Create an NArray with the same length and contents as self.



72
73
74
 
# File 'lib/Dobjects/Dvector_extras.rb', line 72

def to_na
  ::NArray.to_na(_dump(nil)[5..-1], ::NArray::DFLOAT)
end

#to_s ⇒ Object Also known as: inspect

#tridag ⇒ Object

#trim ⇒ Object

#trim! ⇒ Object

#uniq ⇒ Object

#uniq! ⇒ Object

#unshift ⇒ Object

#values_at ⇒ Object

#vector_length ⇒ Object

#where_closest ⇒ Object Also known as: where_first_closest

#where_first_eq ⇒ Object Also known as: where_eq

#where_first_ge ⇒ Object Also known as: where_ge

#where_first_gt ⇒ Object Also known as: where_gt

#where_first_le ⇒ Object Also known as: where_le

#where_first_lt ⇒ Object Also known as: where_lt

#where_first_ne ⇒ Object Also known as: where_ne

#where_last_closest ⇒ Object

#where_last_eq ⇒ Object

#where_last_ge ⇒ Object

#where_last_gt ⇒ Object

#where_last_le ⇒ Object

#where_last_lt ⇒ Object

#where_last_max ⇒ Object

#where_last_min ⇒ Object

#where_last_ne ⇒ Object

#where_max ⇒ Object Also known as: where_first_max

#where_min ⇒ Object Also known as: where_first_min