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   601 use Text::NumericData::App;
  1         3  
  1         41  
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   9 use strict;
  1         3  
  1         1869  
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 108 my $class = shift;
29 1         10 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         20 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 23 my $self = shift;
56 7         23 my $p = $self->{param};
57              
58             return $self->error("Append mode and dualgrid don't mix.")
59 7 50 66     58 if($p->{append} and $p->{dualgrid});
60              
61 7 50       23 if(@{$self->{argv}})
  7         48  
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       54 unless(valid_column([$p->{xcol}]));
76             return $self->error("Bad y column index/indices!")
77 7 50       32 unless(valid_column($p->{ycol}));
78              
79             # Translate to plain 0-based indices.
80 7         35 $self->{xi} = $p->{xcol}-1;
81 7         23 $self->{yi} = [@{$p->{ycol}}];
  7         38  
82 7         19 for(@{$self->{yi}}){ --$_; }
  7         31  
  8         27  
83              
84 7         30 return 0;
85             }
86              
87              
88             sub process_file
89             {
90 7     7 0 19 my $self = shift;
91 7         30 my $p = $self->{param};
92 7         18 my $txd = $self->{txd};
93              
94             # sort out titles
95 7         35 my @cidx = ($p->{xcol}-1);
96 7         34 my $cols = $txd->columns();
97 7 50       30 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         20 for my $y (@{$p->{ycol}})
  7         30  
106             {
107 8         52 push(@derivtitles, derivtitle($txd->{titles}[$y-1], $y));
108 8         43 push(@cidx, $y-1);
109             }
110              
111 7 100       22 if(@{$txd->{titles}})
  7         33  
112             {
113 4 100       21 if($p->{append})
114             {
115 2 50       16 $txd->{titles}[$cols-1] = '' unless defined $txd->{titles}[$cols-1];
116 2         6 push(@{$txd->{titles}}, @derivtitles);
  2         13  
117             }
118             else
119             {
120 2         14 my $xt = $txd->{titles}[$p->{xcol}-1];
121 2         8 @{$txd->{titles}} = ($xt, @derivtitles);
  2         13  
122             }
123             }
124              
125 7 100       18 if(@{$txd->{raw_header}}){ $txd->write_header($self->{out}); }
  7         29  
  4         61  
126 7 100       18 if(@{$txd->{titles}}){ print {$self->{out}} ${$txd->title_line()}; }
  7         31  
  4         12  
  4         15  
  4         27  
127              
128             # compute derivatives, after sorting data
129 7 50       85 $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       67 : central_diff_samegrid($txd->{data}, $self->{xi}, $self->{yi});
133 7 100       35 if($p->{append})
134             {
135 3         12 for(my $i=0; $i<=$#{$diffdata}; ++$i)
  33         102  
136             {
137 30         61 shift(@{$diffdata->[$i]});
  30         68  
138 30         77 push(@{$txd->{data}[$i]}, @{$diffdata->[$i]});
  30         70  
  30         96  
139             }
140             }
141             else
142             {
143 4         33 $txd->{data} = $diffdata;
144             }
145 7         48 $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 62 my ($cols, $max) = @_;
152 14 50       38 return 0 if (grep {$_ != int($_) or $_ < 1} @{$cols});
  15 50       152  
  14         43  
153 14 50 33     60 return 0 if (defined $max and grep {$_ > $max} @{$cols});
  0         0  
  0         0  
154 14         68 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 12 my ($data, $xi, $yi) = @_;
163 2         6 my @diff;
164 2         6 for(my $i=1; $i<=$#{$data}; ++$i)
  20         70  
165             {
166 18         41 my $prev = $data->[$i-1];
167 18         38 my $next = $data->[$i];
168             # treat empty lines as separations
169 18 50 33     34 if(not @{$prev} or not @{$next})
  18         52  
  18         67  
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         53 my $dx = $next->[$xi] - $prev->[$xi];
177 18 50       62 next if not $dx > 1e-12; # some safety precaution, better tunable?
178 18         52 my $invdx = 1./$dx;
179 18         76 my @point = ( 0.5*($prev->[$xi]+$next->[$xi]) );
180 18         40 for(@{$yi})
  18         50  
181             {
182 27         129 push(@point, $invdx*($next->[$_]-$prev->[$_]));
183             }
184 18         71 push(@diff, \@point);
185             }
186 2         9 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 28 my ($data, $xi, $yi) = @_;
195 5         14 my @diff;
196             # Cannot do anything with a single value ... well, could call it a zero.
197 5 50       11 return \@diff if(@{$data} < 2);
  5         23  
198              
199 5         19 for(my $i=0; $i<=$#{$data}; ++$i)
  55         198  
200             {
201 50         118 my $this = $data->[$i];
202 50 100       140 my $prev = $i > 0 ? $data->[$i-1] : undef;
203 50 100       95 my $next = $i < $#{$data} ? $data->[$i+1] : undef;
  50         149  
204 50         184 my @point = ( $this->[$xi] );
205             # empty line, treat separate pieces
206 50         132 for ($prev, $this, $next)
207             {
208 150 50 66     389 $_ = undef if(defined $_ and not @{$_});
  140         530  
209             }
210 50 50       137 unless(defined $this)
211             {
212 0         0 push(@diff, []);
213 0         0 next;
214             }
215 50 100       179 my $v = defined $prev ? $this->[$xi] - $prev->[$xi] : undef;
216 50 100       180 my $w = defined $next ? $next->[$xi] - $this->[$xi] : undef;
217              
218 50 100       162 if(not defined $prev)
    100          
219             {
220             # f'(x) = (f(x+w)-f(x)) / w
221 5         20 my $invdx = 1./$w;
222 5         14 for(@{$yi})
  5         18  
223             {
224 5         37 push(@point, $invdx*($next->[$_]-$this->[$_]));
225             }
226             }
227             elsif(not defined $next)
228             {
229             # f'(x) = (f(x+w)-f(x)) / w
230 5         20 my $invdx = 1./$v;
231 5         12 for(@{$yi})
  5         18  
232             {
233 5         25 push(@point, $invdx*($this->[$_]-$prev->[$_]));
234             }
235             }
236             else
237             {
238 40         110 my $wbyv = $w/$v;
239 40         144 my $invdx = 1./($w*(1+$wbyv));
240 40         169 my $thisweight = 1-$wbyv**2;
241 40         123 my $prevweight = $wbyv**2;
242 40         79 for(@{$yi})
  40         106  
243             {
244 40         293 push(@point, $invdx*($next->[$_] - $thisweight*$this->[$_] - $prevweight*$prev->[$_]));
245             }
246             }
247 50         216 push(@diff, \@point);
248             }
249 5         22 return \@diff;
250             }
251              
252             sub derivtitle
253             {
254 8     8 0 29 my $ot = shift;
255 8         26 my $col = shift;
256 8 100       63 return 'd('.(defined $ot ? $ot : "col$col").')';
257             }