! ************************************************************************************************** ! ____________________ ___ ________ __ ! / ____/ _/_ __/ __ \/ | / ____/ //_/ ! / /_ / / / / / /_/ / /| |/ / / ,< ! / __/ _/ / / / / ____/ ___ / /___/ /| | ! /_/ /___/ /_/ /_/ /_/ |_\____/_/ |_| ! ! A Curve Fitting Package ! ! Refactored by Federico Perini, 10/6/2022 ! Based on the netlib library by Paul Dierckx ! ! References : ! - C. De Boor, "On calculating with b-splines", J Approx Theory 6 (1972) 50-62 ! - M. G. Cox, "The numerical evaluation of b-splines", J Inst Maths Applics 10 (1972) 134-149 ! - P. Dierckx, "Curve and surface fitting with splines", Monographs on numerical analysis, ! Oxford university press, 1993. ! ! ************************************************************************************************** module fitpack use fitpack_core implicit none private ! Public interface public :: RKIND public :: fitpack_curve public :: fitpack_surface ! Max fitting degree integer, parameter :: MAX_K = 5 !> A public type describing a surface fitter z = s(x,y) type :: fitpack_surface !> The data points integer :: m = 0 real(RKIND), allocatable :: x(:),y(:),z(:) !> Spline degree integer :: order(2) = 3 !> Interval boundaries real(RKIND) :: left(2),right(2) ! Node weights real(RKIND), allocatable :: w(:) ! Estimated and actual number of knots and their allocations integer :: nest(2) = 0 integer :: nmax = 0 integer :: lwrk1 = 0, lwrk2 = 0, liwrk = 0 integer, allocatable :: iwrk(:) real(RKIND), allocatable :: wrk1(:),wrk2(:) ! Curve fit smoothing parameter (fit vs. points MSE) real(RKIND) :: smoothing = 1000.d0 ! Actual curve MSE real(RKIND) :: fp = zero ! Curve extrapolation behavior integer :: bc = OUTSIDE_NEAREST_BND ! Knots integer :: knots(2) = 0 real(RKIND), allocatable :: t(:,:) ! Knot locations (:,1)=x; (:,2)=y ! Spline coefficients [knots-order-1] real(RKIND), allocatable :: c(:) ! Runtime flag integer :: iopt = 0 contains !> Clean memory procedure :: destroy => surf_destroy !> Set new points procedure :: new_points => surf_new_points !> Generate new fit procedure :: new_fit => surf_new_fit !> Generate/update fitting curve, with optional smoothing procedure :: fit => surface_fit_automatic_knots end type fitpack_surface !> A public type describing a curve fitter y = c(x) type :: fitpack_curve !> The data points integer :: m = 0 real(RKIND), allocatable :: x(:),y(:) !> Spline degree integer :: order = 3 !> Interval boundaries real(RKIND) :: xleft,xright ! Node weights real(RKIND), allocatable :: sp(:),w(:) ! Estimated and actual number of knots and their allocations integer :: nest = 0 integer :: lwrk = 0 integer, allocatable :: iwrk(:) real(RKIND), allocatable :: wrk(:) ! Curve fit smoothing parameter (fit vs. points MSE) real(RKIND) :: smoothing = 1000.0_RKIND ! Actual curve MSE real(RKIND) :: fp = zero ! Curve extrapolation behavior integer :: bc = OUTSIDE_NEAREST_BND ! Knots integer :: knots = 0 real(RKIND), allocatable :: t(:) ! Knot location ! Spline coefficients [knots-order-1] real(RKIND), allocatable :: c(:) ! Runtime flag integer :: iopt = 0 contains !> Clean memory procedure :: destroy !> Set new points procedure :: new_points !> Generate new fit procedure :: new_fit !> Generate/update fitting curve, with optional smoothing procedure :: fit => curve_fit_automatic_knots procedure :: interpolate => interpolating_curve !> Evaluate curve at given coordinates procedure, private :: curve_eval_one procedure, private :: curve_eval_many generic :: eval => curve_eval_one,curve_eval_many !> Evaluate derivative at given coordinates procedure, private :: curve_derivative procedure, private :: curve_derivatives generic :: dfdx => curve_derivative,curve_derivatives !> Properties: MSE procedure, non_overridable :: mse => curve_error end type fitpack_curve ! Default constructor interface fitpack_curve module procedure new_from_points end interface fitpack_curve interface fitpack_surface module procedure surf_new_from_points end interface fitpack_surface contains ! A default constructor type(fitpack_curve) function new_from_points(x,y,w,ierr) result(this) real(RKIND), intent(in) :: x(:),y(size(x)) real(RKIND), optional, intent(in) :: w(size(x)) ! node weights integer, optional, intent(out) :: ierr integer :: ierr0 ierr0 = this%new_fit(x,y,w) ! Error handling call fitpack_error_handling(ierr0,ierr,'new curve fit') end function new_from_points ! Fit a new curve integer function new_fit(this,x,y,w,smoothing) class(fitpack_curve), intent(inout) :: this real(RKIND), intent(in) :: x(:),y(size(x)) real(RKIND), optional, intent(in) :: w(size(x)) ! node weights real(RKIND), optional, intent(in) :: smoothing call this%new_points(x,y,w) new_fit = this%fit(smoothing) end function new_fit elemental subroutine destroy(this) class(fitpack_curve), intent(inout) :: this integer :: ierr this%m = 0 deallocate(this%x,stat=ierr) deallocate(this%y,stat=ierr) deallocate(this%w,stat=ierr) deallocate(this%sp,stat=ierr) deallocate(this%iwrk,stat=ierr) deallocate(this% wrk,stat=ierr) deallocate(this%t,stat=ierr) deallocate(this%c,stat=ierr) this%xleft = zero this%xright = zero this%smoothing = 1000.0_RKIND this%order = 3 this%iopt = 0 this%nest = 0 this%lwrk = 0 this%knots = 0 this%fp = 0.0_RKIND this%bc = OUTSIDE_NEAREST_BND end subroutine destroy subroutine new_points(this,x,y,w) class(fitpack_curve), intent(inout) :: this real(RKIND), intent(in) :: x(:),y(size(x)) real(RKIND), optional, intent(in) :: w(size(x)) ! node weights integer, allocatable :: isort(:) integer, parameter :: SAFE = 10 associate(m=>this%m,nest=>this%nest,lwrk=>this%lwrk) call this%destroy() m = size(x) ! Ensure x are sorted isort = RKIND_argsort(x) allocate(this%x,source=x(isort)) allocate(this%y,source=y(isort)) ! set up uniform weights if (present(w)) then allocate(this%w(m),source=w(isort)) else allocate(this%w(m),source=1.0_RKIND) endif allocate(this%sp(m),source=0.0_RKIND) ! Setup boundaries this%xleft = minval(x) this%xright = maxval(x) ! Reset run flag this%iopt = 0 ! Reset estimated knots nest = max(SAFE*2*MAX_K+2,ceiling(1.4*this%m)) allocate(this%iwrk(nest),this%t(nest),this%c(nest)) ! Setup working space. lwrk = (m*(MAX_K+1)+nest*(7+3*MAX_K)) allocate(this%wrk(lwrk),source=0.0_RKIND) endassociate end subroutine new_points real(RKIND) function curve_eval_one(this,x,ierr) result(y) class(fitpack_curve), intent(inout) :: this real(RKIND), intent(in) :: x ! Evaluation point integer, optional, intent(out) :: ierr ! Optional error flag real(RKIND) :: y1(1) y1 = curve_eval_many(this,[x],ierr) y = y1(1) end function curve_eval_one ! Curve evaluation driver function curve_eval_many(this,x,ierr) result(y) class(fitpack_curve), intent(inout) :: this real(RKIND), intent(in) :: x(:) ! Evaluation points integer, optional, intent(out) :: ierr ! Optional error flag real(RKIND) :: y(size(x)) integer :: npts,ier npts = size(x) ! subroutine splev evaluates in a number of points x(i),i=1,2,...,m ! a spline s(x) of degree k, given in its b-spline representation. ! ! calling sequence: ! call splev(t,n,c,k,x,y,m,e,ier) ! ! input parameters: ! t : array,length n, which contains the position of the knots. ! n : integer, giving the total number of knots of s(x). ! c : array,length n, which contains the b-spline coefficients. ! k : integer, giving the degree of s(x). ! x : array,length m, which contains the points where s(x) must ! be evaluated. ! m : integer, giving the number of points where s(x) must be ! evaluated. ! e : integer, if 0 the spline is extrapolated from the end ! spans for points not in the support, if 1 the spline ! evaluates to zero for those points, if 2 ier is set to ! 1 and the subroutine returns, and if 3 the spline evaluates ! to the value of the nearest boundary point. call splev(t=this%t,& ! the position of the knots n=this%knots,& ! total number of knots of s(x) c=this%c,& ! the b-spline coefficients k=this%order,& ! the degree of s(x) x=x,m=npts,& ! the points where s(x) must be evaluated. y=y, & ! the predictions e=this%bc, & ! What to do outside mapped knot range ier=ier) call fitpack_error_handling(ier,ierr,'evaluate 1d spline') end function curve_eval_many ! Interpolating curve integer function interpolating_curve(this) result(ierr) class(fitpack_curve), intent(inout) :: this ! Set zero smoothing ierr = curve_fit_automatic_knots(this,zero) end function interpolating_curve ! Curve fitting driver: automatic number of knots integer function curve_fit_automatic_knots(this,smoothing) result(ierr) class(fitpack_curve), intent(inout) :: this real(RKIND), optional, intent(in) :: smoothing integer :: loop,nit real(RKIND), parameter :: smoothing_trajectory(*) = [1000.d0,60.d0,30.d0] real(RKIND), dimension(size(smoothing_trajectory)) :: smooth_now if (present(smoothing)) then smooth_now = smoothing nit = 1 else smooth_now = smoothing_trajectory nit = size(smoothing_trajectory) end if ! First iteration lets solver decide knots this%iopt = 0 do loop=1,nit ! Set current smoothing this%smoothing = smooth_now(loop) ! Call curvfit call curfit(this%iopt, & ! option this%m,this%x,this%y,this%w, & ! points this%xleft,this%xright, & ! x range this%order,this%smoothing, & ! spline accuracy this%nest,this%knots,this%t, & ! spline output this%c,this%fp, & ! spline output this%wrk,this%lwrk,this%iwrk,& ! memory ierr) ! Error flag end do end function curve_fit_automatic_knots ! Return fitting MSE elemental real(RKIND) function curve_error(this) class(fitpack_curve), intent(in) :: this curve_error = this%fp end function curve_error ! Utilities: argsort ! Return indices of sorted array pure function RKIND_argsort(list) result(ilist) real(RKIND), dimension(:), intent(in) :: list integer(RSIZE), dimension(size(list)) :: ilist real(RKIND), allocatable :: copy(:) integer(RSIZE) :: i ! Prepare data allocate(copy(size(list)),source=list) forall(i=1:size(list,kind=RSIZE)) ilist(i) = i ! Perform sort call RKIND_quicksort_andlist(copy,ilist) deallocate(copy) end function RKIND_argsort ! Quicksort pure recursive subroutine RKIND_quicksort_andlist(list,ilist,down) real(RKIND), dimension(:), intent(inout) :: list integer(RSIZE), dimension(size(list)), intent(inout) :: ilist logical, optional, intent(in) :: down integer(RSIZE) :: i, j, n real(RKIND) :: chosen integer(RSIZE), parameter :: SMALL_SIZE = 8 logical :: descending descending = .false.; if (present(down)) descending = down n = size(list,kind=RSIZE) choose_sorting_algorithm: if (n <= SMALL_SIZE) then ! Use interchange sort for small lists do i = 1, n - 1 do j = i + 1, n if (toBeSwapped(list(i),list(j),.false.)) then call swap_data(list(i),list(j)) call swap_size(ilist(i),ilist(j)) end if end do end do else ! Use partition (quick) sort if the list is big chosen = list(int(n/2)) i = 0 j = n + 1 scan_lists: do ! Scan list from left end ! until element >= chosen is found scan_from_left: do i = i + 1 if (toBeSwapped(list(i),chosen,.true.) .or. i==n) exit scan_from_left end do scan_from_left ! Scan list from right end ! until element <= chosen is found scan_from_right: do j = j - 1 if (toBeSwapped(chosen,list(j),.true.) .or. j==1) exit scan_from_right end do scan_from_right swap_element: if (i < j) then ! Swap two out of place elements call swap_data(list(i),list(j)) call swap_size(ilist(i),ilist(j)) else if (i == j) then i = i + 1 exit else exit endif swap_element end do scan_lists if (1 < j) call RKIND_quicksort_andlist(list(:j),ilist(:j),down) if (i < n) call RKIND_quicksort_andlist(list(i:),ilist(i:),down) end if choose_sorting_algorithm ! test for small array contains elemental logical function toBeSwapped(a,b,orEqual) real(RKIND), intent(in) :: a,b logical, intent(in) :: orEqual toBeSwapped = merge(is_before(a,b),is_after(a,b),descending) if (orEqual .and. a==b) toBeSwapped = .true. end function toBeSwapped end subroutine RKIND_quicksort_andlist elemental subroutine swap_data(a,b) real(RKIND), intent(inout) :: a, b real(RKIND) :: tmp tmp = a a = b b = tmp return end subroutine swap_data elemental subroutine swap_size(a,b) integer(RSIZE), intent(inout) :: a, b integer(RSIZE) :: tmp tmp = a a = b b = tmp return end subroutine swap_size elemental logical function is_before(a,b) real(RKIND), intent(in) :: a,b is_before = a<b end function is_before elemental logical function is_after(a,b) real(RKIND), intent(in) :: a,b is_after = a>b end function is_after elemental logical function is_ge(a,b) real(RKIND), intent(in) :: a,b is_ge = a>=b end function is_ge elemental logical function is_le(a,b) real(RKIND), intent(in) :: a,b is_le = a<=b end function is_le ! Flow control: on output flag present, return it; ! otherwise, halt on error subroutine fitpack_error_handling(ierr,ierr_out,whereAt) integer, intent(in) :: ierr integer, optional, intent(out) :: ierr_out character(*), intent(in) :: whereAt if (present(ierr_out)) then ierr_out = ierr elseif (.not.FITPACK_SUCCESS(ierr)) then print *, '[fitpack] at '//trim(whereAt)//' failed with error '//FITPACK_MESSAGE(ierr) stop ierr end if end subroutine fitpack_error_handling !> Evaluate k-th derivative of the curve at points x !> Use 1st derivative if order not present function curve_derivatives(this, x, order, ierr) result(ddx) class(fitpack_curve), intent(inout) :: this real(RKIND), intent(in) :: x(:) ! Evaluation point (scalar) integer, optional, intent(in) :: order ! Derivative order. Default 1 integer, optional, intent(out) :: ierr ! Optional error flag real(RKIND), dimension(size(x)) :: ddx integer :: ddx_order,m,ierr0 if (present(order)) then ddx_order = max(0,order) else ddx_order = 1 end if ierr0 = FITPACK_OK m = size(x); if (m<=0) goto 1 ! subroutine splder evaluates in a number of points x(i),i=1,2,...,m the derivative of ! order nu of a spline s(x) of degree k, given in its b-spline representation. call splder(this%t, & ! Position of the knots this%knots, & ! Number of knots this%c, & ! spline coefficients this%order, & ! spline degree ddx_order, & ! derivative order (0<=order<=this%order) x,y,m,e,wrk,ier) x, & ! Array of points where this should be evaluated ddx, & ! Evaluated derivatives m, & ! Number of input points this%bc, & ! Extrapolation behavior this%wrk, & ! Temporary working space ierr0) ! Output flag 1 call fitpack_error_handling(ierr0,ierr,'evaluate derivative') end function curve_derivatives !> Evaluate k-th derivative of the curve at points x !> Use 1st derivative if order not present real(RKIND) function curve_derivative(this, x, order, ierr) result(ddx) class(fitpack_curve), intent(inout) :: this real(RKIND), intent(in) :: x ! Evaluation point (scalar) integer, optional, intent(in) :: order ! Derivative order. Default 1 integer, optional, intent(out) :: ierr ! Optional error flag real(RKIND) :: ddxa(1) ! Use the array-based wrapper ddxa = curve_derivatives(this,[x],order,ierr) ddx = ddxa(1) end function curve_derivative ! Fit a surface z = s(x,y) defined on a meshgrid: x[1:n], y[1:m] integer function surface_fit_automatic_knots(this,smoothing,order) result(ierr) class(fitpack_surface), intent(inout) :: this real(RKIND), optional, intent(in) :: smoothing integer, optional, intent(in) :: order integer :: loop,nit real(RKIND), parameter :: smoothing_trajectory(*) = [1000.d0,60.d0,30.d0] real(RKIND), dimension(size(smoothing_trajectory)) :: smooth_now if (present(smoothing)) then smooth_now = smoothing nit = 1 else smooth_now = smoothing_trajectory nit = size(smoothing_trajectory) end if ! User may want to change the order for both x and y if (present(order)) this%order = order ! First iteration lets solver decide knots this%iopt = 0 do loop=1,nit ! Set current smoothing this%smoothing = smooth_now(loop) ! Call curvfit call surfit(this%iopt, & ! [-1]=lsq on given knots; [0,1]=smoothing spline this%m,this%x,this%y,this%z, & ! points and their coordinates this%w, & ! weights this%left(1),this%right(1), & ! x range this%left(2),this%right(2), & ! y range this%order(1),this%order(2), & ! [1:5] x,y spline order. Recommended: bicubic (x=y=3) this%smoothing, & ! spline accuracy (iopt>=0) this%nest(1),this%nest(2), & ! estimated number of knots and storage nxest >= 2*(kx+1), nyest >= 2*(ky+1) this%nmax, & ! nmax>=max(nxest,nyest) size of tx,ty epsilon(zero), & ! numeric limits of the current kind this%knots(1),this%t(:,1), & ! x knots (out) this%knots(2),this%t(:,2), & ! y knots (out) this%c,this%fp, & ! spline output. size(c)>=(nxest-kx-1)*(nyest-ky-1) this%wrk1,this%lwrk1, & ! memory this%wrk2,this%lwrk2, & ! memory this%iwrk,this%liwrk, & ! memory ierr) ! Error flag end do end function surface_fit_automatic_knots elemental subroutine surf_destroy(this) class(fitpack_surface), intent(inout) :: this integer :: ierr this%m = 0 deallocate(this%x,stat=ierr) deallocate(this%y,stat=ierr) deallocate(this%z,stat=ierr) deallocate(this%w,stat=ierr) deallocate(this%iwrk,stat=ierr) deallocate(this%wrk1,stat=ierr) deallocate(this%wrk2,stat=ierr) deallocate(this%t,stat=ierr) deallocate(this%c,stat=ierr) this%left = zero this%right = zero this%smoothing = 1000.0_RKIND this%order = 3 this%iopt = 0 this%nest = 0 this%nmax = 0 this%lwrk1 = 0 this%lwrk2 = 0 this%liwrk = 0 this%knots = 0 this%fp = zero this%bc = OUTSIDE_NEAREST_BND end subroutine surf_destroy subroutine surf_new_points(this,x,y,z,w) class(fitpack_surface), intent(inout) :: this real(RKIND), intent(in) :: x(:),y(size(x)),z(size(x)) real(RKIND), optional, intent(in) :: w(size(x)) ! node weights integer :: clen,uv(2),km,bxy(2),b1,b2 integer, parameter :: SAFE = 2 associate(m=>this%m,nest=>this%nest,nmax=>this%nmax,order=>this%order) call this%destroy() m = size(x) ! Ensure x are sorted allocate(this%x,source=x) allocate(this%y,source=y) allocate(this%z,source=z) ! set up uniform weights if (present(w)) then allocate(this%w(m),source=w) else allocate(this%w(m),source=one) endif ! Setup boundaries this%left(1) = minval(x,1) this%left(2) = minval(y,1) this%right(1) = maxval(x,1) this%right(2) = maxval(y,1) ! Reset run flag this%iopt = 0 ! Knot space: overestimate (2*order+1 => order+m+1) nest = SAFE*(order + m + 1) nmax = maxval(nest) nest = nmax allocate(this%t(nmax,2),source=zero) ! Spline coefficients clen = product(nest-order-1) allocate(this%c(clen),source=zero) this%fp = zero ! Working space this%liwrk = m+product(nest-2*order-1) allocate(this%iwrk(this%liwrk),source=0) ! wrk1 uv = nest-order-1 km = maxval(order,1)+1 bxy(1) = order(1)*uv(2)+order(2)+1 bxy(2) = order(2)*uv(1)+order(1)+1 if (bxy(1)<=bxy(2)) then b1 = bxy(1); b2 = b1 + uv(2) - order(2) else b1 = bxy(2); b2 = b1 + uv(1) - order(1) end if this%lwrk1 = product(uv)*(2+b1+b2)+2*(sum(uv)+km*(m+nmax)+nmax-sum(order))+b2+1 allocate(this%wrk1(this%lwrk1),source=zero) ! wrk2 this%lwrk2 = product(uv)*(b2+1)+b2 allocate(this%wrk2(this%lwrk2),source=zero) endassociate end subroutine surf_new_points ! A default constructor type(fitpack_surface) function surf_new_from_points(x,y,z,w,ierr) result(this) real(RKIND), intent(in) :: x(:),y(size(x)),z(size(x)) real(RKIND), optional, intent(in) :: w(size(x)) ! node weights integer, optional, intent(out) :: ierr integer :: ierr0 ierr0 = this%new_fit(x,y,z,w) ! Error handling call fitpack_error_handling(ierr0,ierr,'new surface fit') end function surf_new_from_points ! Fit a new curve integer function surf_new_fit(this,x,y,z,w,smoothing,order) class(fitpack_surface), intent(inout) :: this real(RKIND), intent(in) :: x(:),y(size(x)),z(size(x)) real(RKIND), optional, intent(in) :: w(size(x)) ! node weights real(RKIND), optional, intent(in) :: smoothing integer , optional, intent(in) :: order call this%new_points(x,y,z,w) surf_new_fit = this%fit(smoothing,order) end function surf_new_fit end module fitpack