30
|
1 |
#======================================================================
|
|
2 |
# . N M I N T E R P . L I N E A R
|
|
3 |
# doc: Wed Nov 22 21:01:09 2000
|
39
|
4 |
# dlm: Fri Aug 30 13:26:58 2019
|
30
|
5 |
# (c) 2000 A.M. Thurnherr
|
39
|
6 |
# uE-Info: 101 9 NIL 0 0 72 0 2 4 NIL ofnI
|
30
|
7 |
#======================================================================
|
|
8 |
|
|
9 |
# linear interpolation for [fillgaps]
|
|
10 |
|
|
11 |
# HISTORY:
|
|
12 |
# Apr 6, 2006: - created from [.interp.linear]
|
|
13 |
# Jul 28, 2006: - added xf to ISInit() args
|
|
14 |
# Jul 30, 2006: - added max-gap support
|
|
15 |
# Aug 8, 2008: - documentation
|
|
16 |
# Sep 23, 2011: - added support for %RECNO
|
|
17 |
# - removed xv from interp() args
|
|
18 |
# Oct 19, 2011: - code did not work for %RECNO and trailing missing vals
|
|
19 |
# Jan 23, 2014: - implemented huge speedup for sparse files
|
|
20 |
|
|
21 |
# NOTES:
|
|
22 |
# - in contrast to the [.interp.*] routines, the [.nminterp.*] routines:
|
|
23 |
# 1) do not assume that x values increase monotonically
|
|
24 |
# 2) can handle nans
|
|
25 |
# 3) cannot be used with multiple buffers
|
|
26 |
# 4) use %RECNO when xfnr==-1
|
|
27 |
|
|
28 |
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
|
|
29 |
#
|
|
30 |
# THE FOLLOWING VARIABLES MUST BE SET GLOBALLY (i.e. during loading)
|
|
31 |
#
|
|
32 |
# $ISOpts string of allowed options
|
|
33 |
# $ISOptsUsage usage information string for options
|
|
34 |
#
|
|
35 |
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
|
|
36 |
|
|
37 |
$ISOpts = "";
|
|
38 |
$ISOptsUsage = "";
|
|
39 |
|
|
40 |
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
|
|
41 |
#
|
|
42 |
# &ISUsage() mangle parameters (options, really)
|
|
43 |
#
|
|
44 |
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
|
|
45 |
|
|
46 |
sub ISUsage() {}
|
|
47 |
|
|
48 |
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
|
|
49 |
#
|
|
50 |
# &ISInit(f,xf) init interpolation of field f
|
|
51 |
# f field number
|
|
52 |
# xf x field number or -1 for %RECNO
|
|
53 |
# <ret val> none
|
|
54 |
#
|
|
55 |
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
|
|
56 |
|
|
57 |
sub ISInit($$) {}
|
|
58 |
|
|
59 |
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
|
|
60 |
#
|
|
61 |
# &interpolate(xf,mg,f,lvr,cr) interpolate field f
|
|
62 |
# xf x field or -1 for %RECNO
|
|
63 |
# mg max allowed x gap
|
|
64 |
# f field number to interpolate
|
|
65 |
# lvr last record with valid y val
|
|
66 |
# cr current record
|
|
67 |
# <ret val> interpolated value
|
|
68 |
#
|
|
69 |
# NB:
|
|
70 |
# - handle f == xf
|
|
71 |
# - return undef if interpolation cannot be carried out
|
|
72 |
# - x VALUES MAY NOT BE MONOTONIC
|
|
73 |
#
|
|
74 |
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
|
|
75 |
|
|
76 |
{ my(@nvr); # static: next valid record (per field)
|
|
77 |
|
|
78 |
sub interpolate($$$$$$)
|
|
79 |
{
|
|
80 |
my($xf,$mg,$f,$lvr,$cr) = @_;
|
|
81 |
|
|
82 |
return $ants_[$cr][$xf] if ($xf == $f);
|
|
83 |
return $ants_[$lvr][$f] # edge value is ok
|
|
84 |
if ($xf>=0 && equal($ants_[$cr][$xf],$ants_[$lvr][$xf])) ||
|
|
85 |
($xf<0 && $cr==$lvr);
|
|
86 |
|
|
87 |
unless ($nvr[$f] > $cr) {
|
|
88 |
$nvr[$f] = $cr + 1;
|
|
89 |
while ($nvr[$f] <= $#ants_ && !numberp($ants_[$nvr[$f]][$f])) { $nvr[$f]++; }
|
|
90 |
}
|
|
91 |
|
|
92 |
return undef if ($nvr[$f] > $#ants_);
|
|
93 |
return undef
|
|
94 |
if ($xf>=0 && abs($ants_[$nvr[$f]][$xf]-$ants_[$lvr][$xf])>$mg);
|
|
95 |
return $ants_[$nvr[$f]][$f] # edge value is ok
|
|
96 |
if ($xf>=0 && equal($ants_[$cr][$xf],$ants_[$nvr[$f]][$xf])) ||
|
|
97 |
($xf<0 && $cr==$nvr[$f]);
|
|
98 |
|
39
|
99 |
if ($xf < 0) {
|
|
100 |
die("$nvr[$f] - $lvr") unless ($nvr[$f] - $lvr > 0);
|
|
101 |
}
|
30
|
102 |
my($sc) = ($xf >= 0)
|
|
103 |
? ($ants_[$cr][$xf] - $ants_[$lvr][$xf]) / ($ants_[$nvr[$f]][$xf] - $ants_[$lvr][$xf])
|
|
104 |
: ($cr - $lvr) / ($nvr[$f] - $lvr);
|
|
105 |
return $ants_[$lvr][$f] + $sc * ($ants_[$nvr[$f]][$f] - $ants_[$lvr][$f]);
|
|
106 |
}
|
|
107 |
} # static scope
|
|
108 |
|
|
109 |
#----------------------------------------------------------------------
|
|
110 |
|
|
111 |
1;
|