Line data Source code
1 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3 : !!!! !!!!
4 : !!!! ParaMonte: Parallel Monte Carlo and Machine Learning Library. !!!!
5 : !!!! !!!!
6 : !!!! Copyright (C) 2012-present, The Computational Data Science Lab !!!!
7 : !!!! !!!!
8 : !!!! This file is part of the ParaMonte library. !!!!
9 : !!!! !!!!
10 : !!!! LICENSE !!!!
11 : !!!! !!!!
12 : !!!! https://github.com/cdslaborg/paramonte/blob/main/LICENSE.md !!!!
13 : !!!! !!!!
14 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
15 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
16 :
17 : !> \brief
18 : !> This include file contains procedure implementation of the generic interfaces of [pm_polation](@ref pm_polation).
19 : !>
20 : !> \finmain
21 : !>
22 : !> \author
23 : !> \AmirShahmoradi, Apr 21, 2017, 1:54 AM, Institute for Computational Engineering and Sciences (ICES), The University of Texas at Austin
24 :
25 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
26 :
27 : #define CHECK_CRDX_LEN(PROC) \
28 : CHECK_ASSERTION(__LINE__, size(crdx, 1, IK) == size(func, 1, IK), PROC//SK_": The condition `size(crdx) == size(func)` must hold. size(crdx), size(func) = "//getStr([size(crdx, 1, IK), size(func, 1, IK)]))
29 : #define CHECK_QBOUNDED(PROC) \
30 : CHECK_ASSERTION(__LINE__, minval(crdx, 1) <= queryx .and. queryx <= maxval(crdx, 1), PROC//SK_": The condition `minval(crdx) <= queryx .and. queryx <= maxval(crdx)` must hold. minval(crdx), queryx, maxval(crdx) = "//getStr([minval(crdx), queryx, maxval(crdx)]))
31 : ! Define the procedure name.
32 : #if getExtrap_ENABLED || setExtrap_ENABLED
33 : #define SET_POLATION setExtrap
34 : #define POLATION extrap
35 : #define POLSTR SK_"extrap"
36 : character(*, SK), parameter :: PROC_NAME = "@setExtrap()"
37 : #elif getInterp_ENABLED || setInterp_ENABLED
38 : #define SET_POLATION setInterp
39 : #define POLATION interp
40 : #define POLSTR SK_"interp"
41 : character(*, SK), parameter :: PROC_NAME = "@setInterp()"
42 : #else
43 : #error "Unrecognized interface."
44 : #endif
45 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
46 : #if (getExtrap_ENABLED || getInterp_ENABLED) && ND1_ENABLED
47 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
48 :
49 14 : call SET_POLATION(method, crdx, func, queryx, POLATION)
50 :
51 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
52 : #elif (setExtrap_ENABLED || setInterp_ENABLED) && ND1_ENABLED && QD1_ENABLED
53 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
54 :
55 : integer(IK) :: ique
56 760 : CHECK_ASSERTION(__LINE__, size(queryx, 1, IK) == size(POLATION, 1, IK), PROC_NAME//SK_": The condition `size(queryx) == size("//POLSTR//SK_")` must hold. shape(queryx), shape("//POLSTR//SK_") = "//getStr([shape(queryx, IK), shape(POLATION, IK)]))
57 : #if MNPLE_ENABLED
58 20 : CHECK_ASSERTION(__LINE__, size(queryx, 1, IK) == size(relerr, 1, IK), PROC_NAME//SK_": The condition `size(queryx) == size(relerr)` must hold. shape(queryx), shape(relerr) = "//getStr([shape(queryx, IK), shape(relerr, IK)]))
59 : #endif
60 3670 : do ique = 1, size(queryx, 1, IK)
61 : #if MNPLD_ENABLED || PWLN_ENABLED || MEAN_ENABLED || NEAR_ENABLED || NEXT_ENABLED || PREV_ENABLED
62 2700 : call SET_POLATION(method, crdx, func, queryx(ique), POLATION(ique))
63 : #elif MNPLE_ENABLED
64 970 : call SET_POLATION(method, crdx, func, queryx(ique), POLATION(ique), relerr(ique))
65 : #else
66 : #error "Unrecognized interface."
67 : #endif
68 : end do
69 :
70 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
71 : #elif (setExtrap_ENABLED || setInterp_ENABLED) && ND1_ENABLED && (MNPLD_ENABLED || MNPLE_ENABLED) && QD0_ENABLED
72 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
73 :
74 : #if MNPLD_ENABLED
75 : real(RKC) :: relerr
76 : #elif !MNPLE_ENABLED
77 : #error "Unrecognized interface."
78 : #endif
79 3952 : real(RKC) :: diff, difX, correctionC(size(crdx, kind = IK)), correctionD(size(crdx, kind = IK))
80 : integer(IK) :: i, level, ns
81 : integer(IK) :: lenx
82 : lenx = size(crdx, kind = IK)
83 5928 : CHECK_CRDX_LEN(PROC_NAME)
84 : #if setInterp_ENABLED
85 43470 : CHECK_QBOUNDED(PROC_NAME)
86 : #endif
87 : ns = 1_IK
88 : ! Find the nearest neighbor.
89 1976 : diff = abs(queryx - crdx(1))
90 19641 : do i = 1_IK, lenx
91 17665 : difX = abs(queryx - crdx(i))
92 17665 : if (difX < diff) then
93 : ns = i
94 69 : diff = difX
95 : endif
96 17665 : correctionC(i) = func(i)
97 19641 : correctionD(i) = func(i)
98 : end do
99 1976 : POLATION = func(ns)
100 1976 : ns = ns - 1_IK
101 17665 : do level = 1_IK, lenx - 1_IK
102 86098 : do i = 1_IK, lenx - level
103 70409 : diff = (correctionC(i + 1) - correctionD(i)) / (crdx(i) - crdx(i + level))
104 70409 : correctionD(i) = (crdx(i + level) - queryx) * diff
105 86098 : correctionC(i) = (crdx(i) - queryx) * diff
106 : end do
107 15689 : if (2 * ns < lenx - level)then
108 7732 : relerr = correctionC(ns + 1)
109 : else
110 7957 : relerr = correctionD(ns)
111 7957 : ns = ns - 1_IK
112 : endif
113 17665 : POLATION = POLATION + relerr
114 : end do
115 :
116 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
117 : #elif (setExtrap_ENABLED || setInterp_ENABLED) && ND1_ENABLED && PWLN_ENABLED && QD0_ENABLED
118 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
119 :
120 : integer(IK) :: ibin
121 : real(RKC) :: const, slope
122 747 : CHECK_CRDX_LEN(PROC_NAME)
123 249 : CHECK_ASSERTION(__LINE__, isAscendingAll(crdx), PROC_NAME//SK_": The condition `isAscendingAll(crdx)` must hold. crdx = "//getStr(crdx))
124 : #if setInterp_ENABLED
125 2970 : CHECK_QBOUNDED(PROC_NAME)
126 66 : ibin = getBin(crdx, queryx)
127 : #elif setExtrap_ENABLED
128 : ibin = size(crdx, 1, IK)
129 183 : if (crdx(ibin) < queryx) then
130 12 : ibin = ibin - 1_IK
131 171 : elseif (queryx < crdx(1)) then
132 : ibin = 1_IK
133 : else
134 133 : ibin = getBin(crdx, queryx)
135 : end if
136 : #else
137 : #error "Unrecognized interface."
138 : #endif
139 249 : if (crdx(ibin) /= queryx) then
140 210 : const = 1._RKC / (crdx(ibin) - crdx(ibin + 1))
141 210 : slope = -const * func(ibin + 1)
142 210 : const = +const * func(ibin)
143 210 : POLATION = const * (queryx - crdx(ibin + 1)) + slope * (queryx - crdx(ibin))
144 : else
145 39 : POLATION = func(ibin)
146 : end if
147 :
148 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
149 : #elif (setExtrap_ENABLED || setInterp_ENABLED) && ND1_ENABLED && (MEAN_ENABLED || NEAR_ENABLED || NEXT_ENABLED || PREV_ENABLED) && QD0_ENABLED
150 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
151 :
152 : integer(IK) :: ibin
153 4230 : CHECK_CRDX_LEN(PROC_NAME)
154 1410 : CHECK_ASSERTION(__LINE__, isAscending(crdx), PROC_NAME//SK_": The condition `isAscending(crdx)` must hold. crdx = "//getStr(crdx))
155 : #if setInterp_ENABLED
156 11880 : CHECK_QBOUNDED(PROC_NAME)
157 : #elif setExtrap_ENABLED
158 1146 : if (crdx(size(crdx, 1, IK)) < queryx) then
159 48 : POLATION = func(size(crdx, 1, IK))
160 48 : return
161 1098 : elseif (queryx < crdx(1)) then
162 202 : POLATION = func(1)
163 202 : return
164 : end if
165 : #else
166 : #error "Unrecognized interface."
167 : #endif
168 1160 : ibin = getBin(crdx, queryx)
169 1160 : if (crdx(ibin) /= queryx) then
170 : ! ibin cannot be the last element of `crdx`, unless `queryx` is out of the range covered by `crdx`.
171 : #if MEAN_ENABLED
172 516 : POLATION = .5_RKC * (func(ibin) + func(ibin + 1))
173 : #elif NEAR_ENABLED
174 150 : if (queryx - crdx(ibin) < crdx(ibin + 1) - queryx) then
175 72 : POLATION = func(ibin)
176 : else
177 78 : POLATION = func(ibin + 1)
178 : end if
179 : #elif NEXT_ENABLED
180 150 : POLATION = func(ibin + 1)
181 : #elif PREV_ENABLED
182 150 : POLATION = func(ibin)
183 : #else
184 : #error "Unrecognized interface."
185 : #endif
186 : else
187 194 : POLATION = func(ibin)
188 : end if
189 : #else
190 : !%%%%%%%%%%%%%%%%%%%%%%%%
191 : #error "Unrecognized interface."
192 : !%%%%%%%%%%%%%%%%%%%%%%%%
193 : #endif
194 : #undef SET_POLATION
195 : #undef POLATION
196 : #undef POLSTR
|
