1# interpolate.tcl -- 2# 3# Package for interpolation methods (one- and two-dimensional) 4# 5# Remarks: 6# None of the methods deal gracefully with missing values 7# 8# To do: 9# Add B-splines as methods 10# For spatial interpolation in two dimensions also quadrant method? 11# Method for destroying a table 12# Proper documentation 13# Proper test cases 14# 15# version 0.1: initial implementation, january 2003 16# version 0.2: added linear and Lagrange interpolation, straightforward 17# spatial interpolation, april 2004 18# version 0.3: added Neville algorithm. 19# version 1.0: added cubic splines, september 2004 20# 21# Copyright (c) 2004 by Arjen Markus. All rights reserved. 22# Copyright (c) 2004 by Kevin B. Kenny. All rights reserved. 23# 24# See the file "license.terms" for information on usage and redistribution 25# of this file, and for a DISCLAIMER OF ALL WARRANTIES. 26# 27# RCS: @(#) $Id: interpolate.tcl,v 1.10 2009/10/22 18:19:52 arjenmarkus Exp $ 28# 29#---------------------------------------------------------------------- 30 31package require Tcl 8.4 32package require struct::matrix 33 34# ::math::interpolate -- 35# Namespace holding the procedures and variables 36# 37 38namespace eval ::math::interpolate { 39 variable search_radius {} 40 variable inv_dist_pow 2 41 42 namespace export interp-1d-table interp-table interp-linear \ 43 interp-lagrange 44 namespace export neville 45} 46 47# defineTable -- 48# Define a two-dimensional table of data 49# 50# Arguments: 51# name Name of the table to be created 52# cols Names of the columns (for convenience and for counting) 53# values List of values to fill the table with (must be sorted 54# w.r.t. first column or first column and first row) 55# 56# Results: 57# Name of the new command 58# 59# Side effects: 60# Creates a new command, which is used in subsequent calls 61# 62proc ::math::interpolate::defineTable { name cols values } { 63 64 set table ::math::interpolate::__$name 65 ::struct::matrix $table 66 67 $table add columns [llength $cols] 68 $table add row 69 $table set row 0 $cols 70 71 set row 1 72 set first 0 73 set nocols [llength $cols] 74 set novals [llength $values] 75 while { $first < $novals } { 76 set last [expr {$first+$nocols-1}] 77 $table add row 78 $table set row $row [lrange $values $first $last] 79 80 incr first $nocols 81 incr row 82 } 83 84 return $table 85} 86 87# inter-1d-table -- 88# Interpolate in a one-dimensional table 89# (first column is independent variable, all others dependent) 90# 91# Arguments: 92# table Name of the table 93# xval Value of the independent variable 94# 95# Results: 96# List of interpolated values, including the x-variable 97# 98proc ::math::interpolate::interp-1d-table { table xval } { 99 100 # 101 # Search for the records that enclose the x-value 102 # 103 set xvalues [lrange [$table get column 0] 1 end] 104 105 foreach {row row2} [FindEnclosingEntries $xval $xvalues] break 106 107 set prev_values [$table get row $row] 108 set next_values [$table get row $row2] 109 110 set xprev [lindex $prev_values 0] 111 set xnext [lindex $next_values 0] 112 113 if { $row == $row2 } { 114 return [concat $xval [lrange $prev_values 1 end]] 115 } else { 116 set wprev [expr {($xnext-$xval)/($xnext-$xprev)}] 117 set wnext [expr {1.0-$wprev}] 118 set results {} 119 foreach vprev $prev_values vnext $next_values { 120 set vint [expr {$vprev*$wprev+$vnext*$wnext}] 121 lappend results $vint 122 } 123 return $results 124 } 125} 126 127# interp-table -- 128# Interpolate in a two-dimensional table 129# (first column and first row are independent variables) 130# 131# Arguments: 132# table Name of the table 133# xval Value of the independent row-variable 134# yval Value of the independent column-variable 135# 136# Results: 137# Interpolated value 138# 139# Note: 140# Use bilinear interpolation 141# 142proc ::math::interpolate::interp-table { table xval yval } { 143 144 # 145 # Search for the records that enclose the x-value 146 # 147 set xvalues [lrange [$table get column 0] 2 end] 148 149 foreach {row row2} [FindEnclosingEntries $xval $xvalues] break 150 incr row 151 incr row2 152 153 # 154 # Search for the columns that enclose the y-value 155 # 156 set yvalues [lrange [$table get row 1] 1 end] 157 158 foreach {col col2} [FindEnclosingEntries $yval $yvalues] break 159 160 set yvalues [concat "." $yvalues] ;# Prepend a dummy column! 161 162 set prev_values [$table get row $row] 163 set next_values [$table get row $row2] 164 165 set x1 [lindex $prev_values 0] 166 set x2 [lindex $next_values 0] 167 set y1 [lindex $yvalues $col] 168 set y2 [lindex $yvalues $col2] 169 170 set v11 [lindex $prev_values $col] 171 set v12 [lindex $prev_values $col2] 172 set v21 [lindex $next_values $col] 173 set v22 [lindex $next_values $col2] 174 175 # 176 # value = v0 + a*(x-x1) + b*(y-y1) + c*(x-x1)*(y-y1) 177 # if x == x1 and y == y1: value = v11 178 # if x == x1 and y == y2: value = v12 179 # if x == x2 and y == y1: value = v21 180 # if x == x2 and y == y2: value = v22 181 # 182 set a 0.0 183 if { $x1 != $x2 } { 184 set a [expr {($v21-$v11)/($x2-$x1)}] 185 } 186 set b 0.0 187 if { $y1 != $y2 } { 188 set b [expr {($v12-$v11)/($y2-$y1)}] 189 } 190 set c 0.0 191 if { $x1 != $x2 && $y1 != $y2 } { 192 set c [expr {($v11+$v22-$v12-$v21)/($x2-$x1)/($y2-$y1)}] 193 } 194 195 set result \ 196 [expr {$v11+$a*($xval-$x1)+$b*($yval-$y1)+$c*($xval-$x1)*($yval-$y1)}] 197 198 return $result 199} 200 201# FindEnclosingEntries -- 202# Search within a sorted list 203# 204# Arguments: 205# val Value to be searched 206# values List of values to be examined 207# 208# Results: 209# Returns a list of the previous and next indices 210# 211proc FindEnclosingEntries { val values } { 212 set found 0 213 set row2 1 214 foreach v $values { 215 if { $val <= $v } { 216 set row [expr {$row2-1}] 217 set found 1 218 break 219 } 220 incr row2 221 } 222 223 # 224 # Border cases: extrapolation needed 225 # 226 if { ! $found } { 227 incr row2 -1 228 set row $row2 229 } 230 if { $row == 0 } { 231 set row $row2 232 } 233 234 return [list $row $row2] 235} 236 237# interp-linear -- 238# Use linear interpolation 239# 240# Arguments: 241# xyvalues List of x/y values to be interpolated 242# xval x-value for which a value is sought 243# 244# Results: 245# Estimated value at $xval 246# 247# Note: 248# The list xyvalues must be sorted w.r.t. the x-value 249# 250proc ::math::interpolate::interp-linear { xyvalues xval } { 251 # 252 # Border cases first 253 # 254 if { [lindex $xyvalues 0] > $xval } { 255 return [lindex $xyvalues 1] 256 } 257 if { [lindex $xyvalues end-1] < $xval } { 258 return [lindex $xyvalues end] 259 } 260 261 # 262 # The ordinary case 263 # 264 set idxx -2 265 set idxy -1 266 foreach { x y } $xyvalues { 267 if { $xval < $x } { 268 break 269 } 270 incr idxx 2 271 incr idxy 2 272 } 273 274 set x2 [lindex $xyvalues $idxx] 275 set y2 [lindex $xyvalues $idxy] 276 277 if { $x2 != $x } { 278 set yval [expr {$y+($y2-$y)*($xval-$x)/($x2-$x)}] 279 } else { 280 set yval $y 281 } 282 return $yval 283} 284 285# interp-lagrange -- 286# Use the Lagrange interpolation method 287# 288# Arguments: 289# xyvalues List of x/y values to be interpolated 290# xval x-value for which a value is sought 291# 292# Results: 293# Estimated value at $xval 294# 295# Note: 296# The list xyvalues must be sorted w.r.t. the x-value 297# Furthermore the Lagrange method is not a very practical 298# method, as potentially the errors are unbounded 299# 300proc ::math::interpolate::interp-lagrange { xyvalues xval } { 301 # 302 # Border case: xval equals one of the "nodes" 303 # 304 foreach { x y } $xyvalues { 305 if { $x == $xval } { 306 return $y 307 } 308 } 309 310 # 311 # Ordinary case 312 # 313 set nonodes2 [llength $xyvalues] 314 315 set yval 0.0 316 317 for { set i 0 } { $i < $nonodes2 } { incr i 2 } { 318 set idxn 0 319 set xn [lindex $xyvalues $i] 320 set yn [lindex $xyvalues [expr {$i+1}]] 321 322 foreach { x y } $xyvalues { 323 if { $idxn != $i } { 324 set yn [expr {$yn*($x-$xval)/($x-$xn)}] 325 } 326 incr idxn 2 327 } 328 329 set yval [expr {$yval+$yn}] 330 } 331 332 return $yval 333} 334 335# interp-spatial -- 336# Use a straightforward interpolation method with weights as 337# function of the inverse distance to interpolate in 2D and N-D 338# space 339# 340# Arguments: 341# xyvalues List of coordinates and values at these coordinates 342# coord List of coordinates for which a value is sought 343# 344# Results: 345# Estimated value(s) at $coord 346# 347# Note: 348# The list xyvalues is a list of lists: 349# { {x1 y1 z1 {v11 v12 v13 v14} 350# {x2 y2 z2 {v21 v22 v23 v24} 351# ... 352# } 353# The last element of each inner list is either a single number 354# or a list in itself. In the latter case the return value is 355# a list with the same number of elements. 356# 357# The method is influenced by the search radius and the 358# power of the inverse distance 359# 360proc ::math::interpolate::interp-spatial { xyvalues coord } { 361 variable search_radius 362 variable inv_dist_pow 363 364 set result {} 365 foreach v [lindex [lindex $xyvalues 0] end] { 366 lappend result 0.0 367 } 368 369 set total_weight 0.0 370 371 if { $search_radius != {} } { 372 set max_radius2 [expr {$search_radius*$search_radius}] 373 } else { 374 set max_radius2 {} 375 } 376 377 foreach point $xyvalues { 378 set dist 0.0 379 foreach c [lrange $point 0 end-1] cc $coord { 380 set dist [expr {$dist+($c-$cc)*($c-$cc)}] 381 } 382 if { $max_radius2 == {} || $dist <= $max_radius2 } { 383 if { $inv_dist_pow == 1 } { 384 set dist [expr {sqrt($dist)}] 385 } 386 set total_weight [expr {$total_weight+1.0/$dist}] 387 388 set idx 0 389 foreach v [lindex $point end] r $result { 390 lset result $idx [expr {$r+$v/$dist}] 391 incr idx 392 } 393 } 394 } 395 396 if { $total_weight == 0.0 } { 397 set idx 0 398 foreach r $result { 399 lset result $idx {} 400 incr idx 401 } 402 } else { 403 set idx 0 404 foreach r $result { 405 lset result $idx [expr {$r/$total_weight}] 406 incr idx 407 } 408 } 409 410 return $result 411} 412 413# interp-spatial-params -- 414# Set the parameters for spatial interpolation 415# 416# Arguments: 417# max_search Search radius (if none: use {} or "") 418# power Power for the inverse distance (1 or 2, defaults to 2) 419# 420# Results: 421# None 422# 423proc ::math::interpolate::interp-spatial-params { max_search {power 2} } { 424 variable search_radius 425 variable inv_dist_pow 426 427 set search_radius $max_search 428 if { $power == 1 } { 429 set inv_dist_pow 1 430 } else { 431 set inv_dist_pow 2 432 } 433} 434 435#---------------------------------------------------------------------- 436# 437# neville -- 438# 439# Interpolate a function between tabulated points using Neville's 440# algorithm. 441# 442# Parameters: 443# xtable - Table of abscissae. 444# ytable - Table of ordinates. Must be a list of the same 445# length as 'xtable.' 446# x - Abscissa for which the function value is desired. 447# 448# Results: 449# Returns a two-element list. The first element is the 450# requested ordinate. The second element is a rough estimate 451# of the absolute error, that is, the magnitude of the first 452# neglected term of a power series. 453# 454# Side effects: 455# None. 456# 457#---------------------------------------------------------------------- 458 459proc ::math::interpolate::neville { xtable ytable x } { 460 461 set n [llength $xtable] 462 463 # Initialization. Set c and d to the ordinates, and set ns to the 464 # index of the nearest abscissa. Set y to the zero-order approximation 465 # of the nearest ordinate, and dif to the difference between x 466 # and the nearest tabulated abscissa. 467 468 set c [list] 469 set d [list] 470 set i 0 471 set ns 0 472 set dif [expr { abs( $x - [lindex $xtable 0] ) }] 473 set y [lindex $ytable 0] 474 foreach xi $xtable yi $ytable { 475 set dift [expr { abs ( $x - $xi ) }] 476 if { $dift < $dif } { 477 set ns $i 478 set y $yi 479 set dif $dift 480 } 481 lappend c $yi 482 lappend d $yi 483 incr i 484 } 485 486 # Compute successively higher-degree approximations to the fit 487 # function by using the recurrence: 488 # d_m[i] = ( c_{m-1}[i+1] - d{m-1}[i] ) * (x[i+m]-x) / 489 # (x[i] - x[i+m]) 490 # c_m[i] = ( c_{m-1}[i+1] - d{m-1}[i] ) * (x[i]-x) / 491 # (x[i] - x[i+m]) 492 493 for { set m 1 } { $m < $n } { incr m } { 494 for { set i 0 } { $i < $n - $m } { set i $ip1 } { 495 set ip1 [expr { $i + 1 }] 496 set ipm [expr { $i + $m }] 497 set ho [expr { [lindex $xtable $i] - $x }] 498 set hp [expr { [lindex $xtable $ipm] - $x }] 499 set w [expr { [lindex $c $ip1] - [lindex $d $i] }] 500 set q [expr { $w / ( $ho - $hp ) }] 501 lset d $i [expr { $hp * $q }] 502 lset c $i [expr { $ho * $q }] 503 } 504 505 # Take the straighest path possible through the tableau of c 506 # and d approximations back to the tabulated value 507 if { 2 * $ns < $n - $m } { 508 set dy [lindex $c $ns] 509 } else { 510 incr ns -1 511 set dy [lindex $d $ns] 512 } 513 set y [expr { $y + $dy }] 514 } 515 516 # Return the approximation and the highest-order correction term. 517 518 return [list $y [expr { abs($dy) }]] 519} 520 521# prepare-cubic-splines -- 522# Prepare interpolation based on cubic splines 523# 524# Arguments: 525# xcoord The x-coordinates 526# ycoord Y-values for these x-coordinates 527# Result: 528# Intermediate parameters describing the spline function, 529# to be used in the second step, interp-cubic-splines. 530# Note: 531# Implicitly it is assumed that the function decribed by xcoord 532# and ycoord has a second derivative 0 at the end points. 533# To minimise the work if more than one value is needed, the 534# algorithm is divided in two steps 535# (Derived from the routine SPLINT in Davis and Rabinowitz: 536# Methods for Numerical Integration, AP, 1984) 537# 538proc ::math::interpolate::prepare-cubic-splines {xcoord ycoord} { 539 540 if { [llength $xcoord] < 3 } { 541 return -code error "At least three points are required" 542 } 543 if { [llength $xcoord] != [llength $ycoord] } { 544 return -code error "Equal number of x and y values required" 545 } 546 547 set m2 [expr {[llength $xcoord]-1}] 548 549 set s 0.0 550 set h {} 551 set c {} 552 for { set i 0 } { $i < $m2 } { incr i } { 553 set ip1 [expr {$i+1}] 554 set h1 [expr {[lindex $xcoord $ip1]-[lindex $xcoord $i]}] 555 lappend h $h1 556 if { $h1 <= 0.0 } { 557 return -code error "X values must be strictly ascending" 558 } 559 set r [expr {([lindex $ycoord $ip1]-[lindex $ycoord $i])/$h1}] 560 lappend c [expr {$r-$s}] 561 set s $r 562 } 563 set s 0.0 564 set r 0.0 565 set t {--} 566 lset c 0 0.0 567 568 for { set i 1 } { $i < $m2 } { incr i } { 569 set ip1 [expr {$i+1}] 570 set im1 [expr {$i-1}] 571 set y2 [expr {[lindex $c $i]+$r*[lindex $c $im1]}] 572 set t1 [expr {2.0*([lindex $xcoord $im1]-[lindex $xcoord $ip1])-$r*$s}] 573 set s [lindex $h $i] 574 set r [expr {$s/$t1}] 575 lset c $i $y2 576 lappend t $t1 577 } 578 lappend c 0.0 579 580 for { set j 1 } { $j < $m2 } { incr j } { 581 set i [expr {$m2-$j}] 582 set ip1 [expr {$i+1}] 583 set h1 [lindex $h $i] 584 set yp1 [lindex $c $ip1] 585 set y1 [lindex $c $i] 586 set t1 [lindex $t $i] 587 lset c $i [expr {($h1*$yp1-$y1)/$t1}] 588 } 589 590 set b {} 591 set d {} 592 for { set i 0 } { $i < $m2 } { incr i } { 593 set ip1 [expr {$i+1}] 594 set s [lindex $h $i] 595 set yp1 [lindex $c $ip1] 596 set y1 [lindex $c $i] 597 set r [expr {$yp1-$y1}] 598 lappend d [expr {$r/$s}] 599 set y1 [expr {3.0*$y1}] 600 lset c $i $y1 601 lappend b [expr {([lindex $ycoord $ip1]-[lindex $ycoord $i])/$s 602 -($y1+$r)*$s}] 603 } 604 605 lappend d 0.0 606 lappend b 0.0 607 608 return [list $d $c $b $ycoord $xcoord] 609} 610 611# interp-cubic-splines -- 612# Interpolate based on cubic splines 613# 614# Arguments: 615# coeffs Coefficients resulting from the preparation step 616# x The x-coordinate for which to estimate the value 617# Result: 618# Interpolated value at x 619# 620proc ::math::interpolate::interp-cubic-splines {coeffs x} { 621 foreach {dcoef ccoef bcoef acoef xcoord} $coeffs {break} 622 623 # 624 # Check the bounds - no extrapolation 625 # 626 if { $x < [lindex $xcoord 0] } {error "X value too small"} 627 if { $x > [lindex $xcoord end] } {error "X value too large"} 628 629 # 630 # Which interval? 631 # 632 set idx -1 633 foreach xv $xcoord { 634 if { $xv > $x } { 635 break 636 } 637 incr idx 638 } 639 640 set a [lindex $acoef $idx] 641 set b [lindex $bcoef $idx] 642 set c [lindex $ccoef $idx] 643 set d [lindex $dcoef $idx] 644 set dx [expr {$x-[lindex $xcoord $idx]}] 645 646 return [expr {(($d*$dx+$c)*$dx+$b)*$dx+$a}] 647} 648 649 650 651# 652# Announce our presence 653# 654package provide math::interpolate 1.0.3 655