File Coverage

blib/lib/Text/NumericData/App/txdderiv.pm
Criterion Covered Total %
statement 125 145 86.2
branch 41 60 68.3
condition 6 12 50.0
subroutine 9 9 100.0
pod 0 7 0.0
total 181 233 77.6


line stmt bran cond sub pod time code
1             package Text::NumericData::App::txdderiv;
2              
3 1     1   506 use Text::NumericData::App;
  1         3  
  1         30  
4             # Not using Math::Derivative, central difference mid-interval suits me better.
5             # Actually changing my mind there ... computing on the points again.
6             # You got the dualgrid option now.
7              
8 1     1   5 use strict;
  1         2  
  1         1583  
9              
10             # This is just a placeholder because of a past build system bug.
11             # The one and only version for Text::NumericData is kept in
12             # the Text::NumericData module itself.
13             our $VERSION = '1';
14             $VERSION = eval $VERSION;
15              
16             #the infostring says it all
17             my $infostring = "compute derivative of data sets using central differences
18              
19             Usage:
20             $0 [parameters] [ycol|xcol ycol [ycol2 ...]] < data.dat
21              
22             You can provide the x and y column(s) to work on on the command line as plain numbers of as values for the --xcol and --ycol parameters. The latter are overridden by the former.";
23              
24             our @ISA = ('Text::NumericData::App');
25              
26             sub new
27             {
28 1     1 0 98 my $class = shift;
29 1         8 my @pars =
30             (
31             'xcol',1,'x','abscissa column'
32             ,'ycol',[2],'y','ordinate column, or columns as array'
33             ,'dualgrid',0,'','define the derivative on points between the initial sampling points, thusly creating a new grid; otherwise, the central differences are computed on each input point with its neighbours'
34             ,'append',0,'','do not replace data with the derivatives, append them as additional columns instead (only works if dual-grid mode is off)'
35             ,'sort',1,'','sort the data according to x column (usually desirable unless you got a specific order wich should be monotonic)'
36             );
37 1         14 return $class->SUPER::new
38             ({
39             parconf=>
40             {
41             info=>$infostring # default version
42             # default author
43             # default copyright
44             }
45             ,pardef=>\@pars
46             ,filemode=>1
47             ,pipemode=>1
48             ,pipe_init=>\&prepare
49             ,pipe_file=>\&process_file
50             });
51             }
52              
53             sub prepare
54             {
55 7     7 0 14 my $self = shift;
56 7         11 my $p = $self->{param};
57              
58             return $self->error("Append mode and dualgrid don't mix.")
59 7 50 66     23 if($p->{append} and $p->{dualgrid});
60              
61 7 50       14 if(@{$self->{argv}})
  7         15  
62             {
63 0 0       0 if(@{$self->{argv}} == 1)
  0         0  
64             {
65 0         0 $p->{ycol} = [shift @{$self->{argv}}];
  0         0  
66             }
67             else
68             {
69 0         0 $p->{xcol} = shift @{$self->{argv}};
  0         0  
70 0         0 @{$p->{ycol}} = @{$self->{argv}};
  0         0  
  0         0  
71             }
72             }
73              
74             return $self->error("Bad x column index!")
75 7 50       20 unless(valid_column([$p->{xcol}]));
76             return $self->error("Bad y column index/indices!")
77 7 50       18 unless(valid_column($p->{ycol}));
78              
79             # Translate to plain 0-based indices.
80 7         15 $self->{xi} = $p->{xcol}-1;
81 7         9 $self->{yi} = [@{$p->{ycol}}];
  7         21  
82 7         11 for(@{$self->{yi}}){ --$_; }
  7         19  
  8         29  
83              
84 7         17 return 0;
85             }
86              
87              
88             sub process_file
89             {
90 7     7 0 13 my $self = shift;
91 7         83 my $p = $self->{param};
92 7         15 my $txd = $self->{txd};
93              
94             # sort out titles
95 7         23 my @cidx = ($p->{xcol}-1);
96 7         24 my $cols = $txd->columns();
97 7 50       19 unless($cols > 0)
98             {
99 0         0 print STDERR "No data columns?\n";
100 0         0 $txd->write_all($self->{out});
101 0         0 return;
102             }
103              
104 7         14 my @derivtitles;
105 7         9 for my $y (@{$p->{ycol}})
  7         23  
106             {
107 8         26 push(@derivtitles, derivtitle($txd->{titles}[$y-1], $y));
108 8         25 push(@cidx, $y-1);
109             }
110              
111 7 100       9 if(@{$txd->{titles}})
  7         18  
112             {
113 4 100       12 if($p->{append})
114             {
115 2 50       9 $txd->{titles}[$cols-1] = '' unless defined $txd->{titles}[$cols-1];
116 2         4 push(@{$txd->{titles}}, @derivtitles);
  2         7  
117             }
118             else
119             {
120 2         5 my $xt = $txd->{titles}[$p->{xcol}-1];
121 2         4 @{$txd->{titles}} = ($xt, @derivtitles);
  2         6  
122             }
123             }
124              
125 7 100       13 if(@{$txd->{raw_header}}){ $txd->write_header($self->{out}); }
  7         18  
  4         13  
126 7 100       11 if(@{$txd->{titles}}){ print {$self->{out}} ${$txd->title_line()}; }
  7         27  
  4         6  
  4         8  
  4         11  
127              
128             # compute derivatives, after sorting data
129 7 50       47 $txd->sort_data([$p->{xcol}-1], [0]) if $p->{sort};
130             my $diffdata = $p->{dualgrid}
131             ? central_diff_dualgrid($txd->{data}, $self->{xi}, $self->{yi})
132 7 100       36 : central_diff_samegrid($txd->{data}, $self->{xi}, $self->{yi});
133 7 100       19 if($p->{append})
134             {
135 3         6 for(my $i=0; $i<=$#{$diffdata}; ++$i)
  33         62  
136             {
137 30         38 shift(@{$diffdata->[$i]});
  30         46  
138 30         40 push(@{$txd->{data}[$i]}, @{$diffdata->[$i]});
  30         47  
  30         48  
139             }
140             }
141             else
142             {
143 4         18 $txd->{data} = $diffdata;
144             }
145 7         26 $txd->write_data($self->{out});
146             }
147              
148             # Perhaps a general helper in TextData for validation of input.
149             sub valid_column
150             {
151 14     14 0 31 my ($cols, $max) = @_;
152 14 50       20 return 0 if (grep {$_ != int($_) or $_ < 1} @{$cols});
  15 50       70  
  14         23  
153 14 50 33     35 return 0 if (defined $max and grep {$_ > $max} @{$cols});
  0         0  
  0         0  
154 14         37 return 1;
155             }
156              
157             # local functions, not methods
158              
159             # f'(x) = f(x+h)-f(x-h)/2h + O(h^2)
160             sub central_diff_dualgrid
161             {
162 2     2 0 8 my ($data, $xi, $yi) = @_;
163 2         4 my @diff;
164 2         5 for(my $i=1; $i<=$#{$data}; ++$i)
  20         44  
165             {
166 18         31 my $prev = $data->[$i-1];
167 18         27 my $next = $data->[$i];
168             # treat empty lines as separations
169 18 50 33     19 if(not @{$prev} or not @{$next})
  18         40  
  18         44  
170             {
171 0         0 push(@diff, []);
172             # Really skip it if nothing follows.
173 0 0       0 ++$i unless not @{$next};
  0         0  
174 0         0 next;
175             }
176 18         32 my $dx = $next->[$xi] - $prev->[$xi];
177 18 50       30 next if not $dx > 1e-12; # some safety precaution, better tunable?
178 18         31 my $invdx = 1./$dx;
179 18         38 my @point = ( 0.5*($prev->[$xi]+$next->[$xi]) );
180 18         22 for(@{$yi})
  18         34  
181             {
182 27         59 push(@point, $invdx*($next->[$_]-$prev->[$_]));
183             }
184 18         38 push(@diff, \@point);
185             }
186 2         12 return \@diff;
187             }
188              
189             # f'(x) = (v^2 f(x+w) - (v^2-w^2) f(x) - w^2 f(x-v)) / (wv^2+vw^2) + O(v^3w^2/(wv^2+vw^2)) + O(w^3v^2/(wv^2+vw^2))
190             # = (f(x+w) - (1-(w/v)^2) f(x) - (w/v)^2 f(x-v)) / w(1+w/v) + O(v^3w^2/(wv^2+vw^2)) + O(w^3v^2/(wv^2+vw^2))
191             # you better have some distance between the points, no precaution here
192             sub central_diff_samegrid
193             {
194 5     5 0 60 my ($data, $xi, $yi) = @_;
195 5         9 my @diff;
196             # Cannot do anything with a single value ... well, could call it a zero.
197 5 50       8 return \@diff if(@{$data} < 2);
  5         15  
198              
199 5         12 for(my $i=0; $i<=$#{$data}; ++$i)
  55         107  
200             {
201 50         78 my $this = $data->[$i];
202 50 100       90 my $prev = $i > 0 ? $data->[$i-1] : undef;
203 50 100       68 my $next = $i < $#{$data} ? $data->[$i+1] : undef;
  50         95  
204 50         114 my @point = ( $this->[$xi] );
205             # empty line, treat separate pieces
206 50         81 for ($prev, $this, $next)
207             {
208 150 50 66     274 $_ = undef if(defined $_ and not @{$_});
  140         337  
209             }
210 50 50       98 unless(defined $this)
211             {
212 0         0 push(@diff, []);
213 0         0 next;
214             }
215 50 100       99 my $v = defined $prev ? $this->[$xi] - $prev->[$xi] : undef;
216 50 100       91 my $w = defined $next ? $next->[$xi] - $this->[$xi] : undef;
217              
218 50 100       93 if(not defined $prev)
    100          
219             {
220             # f'(x) = (f(x+w)-f(x)) / w
221 5         11 my $invdx = 1./$w;
222 5         8 for(@{$yi})
  5         12  
223             {
224 5         20 push(@point, $invdx*($next->[$_]-$this->[$_]));
225             }
226             }
227             elsif(not defined $next)
228             {
229             # f'(x) = (f(x+w)-f(x)) / w
230 5         8 my $invdx = 1./$v;
231 5         12 for(@{$yi})
  5         7  
232             {
233 5         12 push(@point, $invdx*($this->[$_]-$prev->[$_]));
234             }
235             }
236             else
237             {
238 40         63 my $wbyv = $w/$v;
239 40         66 my $invdx = 1./($w*(1+$wbyv));
240 40         71 my $thisweight = 1-$wbyv**2;
241 40         52 my $prevweight = $wbyv**2;
242 40         51 for(@{$yi})
  40         66  
243             {
244 40         102 push(@point, $invdx*($next->[$_] - $thisweight*$this->[$_] - $prevweight*$prev->[$_]));
245             }
246             }
247 50         115 push(@diff, \@point);
248             }
249 5         11 return \@diff;
250             }
251              
252             sub derivtitle
253             {
254 8     8 0 19 my $ot = shift;
255 8         11 my $col = shift;
256 8 100       30 return 'd('.(defined $ot ? $ot : "col$col").')';
257             }