gnuplot分段拟合颜色映射表达式

类别:    标签: gnuplot 颜色   阅读次数:   版权: (CC) BY-NC-SA

以前的时候讲过如何使用gnuplot拟合颜色映射解析表达式, 在拟合时使用了分段三次函数. 如果要处理的拟合数据不多, 分段的范围可以手动指定, 虽然显得很笨拙. 如果要拟合大量的数据, 那就需要自动处理, 手动指定拟合范围就不实际了. 最近需要拟合很多的数据, 就重新思考了一下这个问题, 完成了拟合需要的gnuplot脚本.

问题

给出一系列n个数据点, (x_i, y_i), i=1..n, 求y在整个x范围内的分段解析表达式.

解析表达式可以自定义, 一般可以使用多项式.

要求自动确定分段节点. 节点的选择标准是所得解析表达式在两个相邻节点间的平均误差不超过指定值.

得到的解析表达式可用于不方便使用查表或内插的情况. 如果给定的解析表达式是三次函数, 得到的节点可作为三次样条的节点.

算法

读入数据(x_i, y_i), i=1..n, 将x_1作为第一个节点, x_n作为第二个节点, 拟合, 计算平均拟合误差.

如果平均拟合误差大于指定值, 说明选区的范围过大, 需要减小范围. 因此将x_n的前一个点x_n-1作为新的节点, 重新拟合, 计算平均拟合误差.

重复上面的步骤, 直至所得平均拟合误差满足要求, 得到两个节点x_1, x_k.

将x_k作为新的第一个节点, 重复以上步骤, 直至所有x值都处理完毕.

之所以选择从使用x值反向递减的拟合策略, 而不使用x值正向递增的拟合粗略, 是为了得到的节点范围尽量大, 这样所得的分段函数使用更方便. 但这种策略的一个缺点在于无法准确地找到那些明显的不连续点.

代码

使用gnuplot自带的非线性拟合功能即可完成. 但鉴于这个过程中需要多次拟合循环, 而且要记录每次循环得到的节点, 所以需要使用数组功能. 只有gnuplot 5.2以上的版本才支持数组, 所以下面的代码不适用于低版本的gnuplot. 当然, 如果你非要将它用于低版本的gnuplot, 那你肯定能想到办法.

cmfit.gpl
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
# 载入数据集
load 'C:\Users\Jicun\Dropbox\!Script\gpl_cm.dat'

# 定义用到的拟合函数, 绘图函数
f(x)=a0+a1*x+a2*x**2+a3*x**3
F(x, beg, end, a0, a1, a2, a3)= beg<=x ? (x<end?a0+a1*x+a2*x**2+a3*x**3:1/0) : 1/0

set fit quiet # 屏幕不输出拟合过程

Nord=3        # 拟合阶数, 取决于所用拟合函数
maxSTD=1E-3   # 最大平均拟合偏差

col=2         # 要拟合的数据列
#cm='rainbow' # 使用的数据集
#cm='darkrainbow'
#cm='cetr2'
cm='viridis';
#cm='bgy'
#cm='parula'
#cm='avocado'
#cm='moreland';
#cm='thermometer'
#cm='redbluetones'
#cm='tmap'
#cm='ltmap'
#cm='cividis';

data = '$' . cm

# 不同数据集使用的数据点数, 间距可能不同
if(cm eq 'parula'       \
|| cm eq 'viridis'      \
|| cm eq 'moreland')    { N=1024; dt=1./(N-1); dx=1.; nev=30 }
if(cm eq 'cividis'      \
|| cm eq 'cetr2'        \
|| cm eq 'rainbow'      \
|| cm eq 'darkrainbow'  \
|| cm eq 'avocado'      \
|| cm eq 'bgy'          \
|| cm eq 'thermometer'  \
|| cm eq 'redbluetones' \
|| cm eq 'tmap'         \
|| cm eq 'ltmap')       { N=256;  dt=1./(N-1); dx=dt; nev=10 }

# 以下为所得拟合结果

set tit cm
set xl"x"
set yl "RGB"

#plot [] [0:1] $rainbow \
   u ($1*dx):2 ev nev w p pt 7 ps 2 lc 'red'  , \
'' u ($1*dx):3 ev nev w p pt 7 ps 2 lc 'green', \
'' u ($1*dx):4 ev nev w p pt 7 ps 2 lc 'blue' , \
F(x, 0.000000,     0.062745,     0.471375,    -2.554017,    -0.603914,     7.812982), \
F(x, 0.062745,     0.152941,     0.292460,     1.809704,   -31.728520,   117.994977), \
F(x, 0.152941,     0.615686,     0.321657,    -1.002195,     3.507752,    -1.465506), \
F(x, 0.615686,     0.807843,     1.060710,    -5.522637,    12.504767,    -7.328945), \
F(x, 0.807843,     0.933333,     0.431446,    -0.873251,     4.037986,    -2.774247), \
F(x, 0.933333,     1.000000,     5.221261,   -12.772499,    12.792531,    -4.383959), \
F(x, 0.000000,     0.062745,     0.108729,     0.152568,    -0.600651,     7.770915), \
F(x, 0.062745,     0.176471,     0.019507,     1.603318,    -1.113349,    12.925571), \
F(x, 0.176471,     0.364706,    -0.087439,     2.477820,     0.408927,    -4.700034), \
F(x, 0.364706,     0.705882,     0.144965,     1.791961,    -0.862001,    -0.770124), \
F(x, 0.705882,     0.784314,    16.562714,   -66.930100,    95.145955,   -45.549403), \
F(x, 0.784314,     0.878431,   -44.677970,   163.297275,  -194.114854,    75.927463), \
F(x, 0.878431,     1.000000,     2.436863,    -1.845566,    -0.624881,     0.164628), \
F(x, 0.000000,     0.105882,     0.526474,     2.243463,     0.792247,   -31.571500), \
F(x, 0.105882,     0.184314,     0.026914,    13.337161,   -80.796176,   171.577097), \
F(x, 0.184314,     0.254902,     0.913699,    -1.400599,     7.239935,   -13.759153), \
F(x, 0.254902,     0.556863,     0.442544,     4.247120,   -14.052824,    11.110080), \
F(x, 0.556863,     0.952941,     2.960527,    -9.677238,    11.755653,    -4.927506), \
F(x, 0.952941,     1.000000,     0.585773,    -0.453645,     0.000000,    -0.000000); exit

#plot [] [0:1] $darkrainbow \
   u ($1*dx):2 ev nev w p pt 7 ps 2 lc 'red'  , \
'' u ($1*dx):3 ev nev w p pt 7 ps 2 lc 'green', \
'' u ($1*dx):4 ev nev w p pt 7 ps 2 lc 'blue' , \
F(x, 0.000000,     0.309804,     0.236537,     0.246053,    -1.136582,     3.202751), \
F(x, 0.309804,     0.407843,    -1.510327,    13.561592,   -35.298055,    33.624234), \
F(x, 0.407843,     0.603922,    -0.599178,     2.859926,    -0.741459,    -0.178725), \
F(x, 0.603922,     0.701961,     3.197059,   -12.278383,    20.066775,   -10.373348), \
F(x, 0.701961,     0.901961,     0.468624,     2.127260,    -2.783342,     0.824166), \
F(x, 0.901961,     1.000000,     0.729870,     0.000000,    -0.000000,     0.000000), \
F(x, 0.000000,     0.109804,     0.339471,     0.168823,    -3.598735,    26.873972), \
F(x, 0.109804,     0.309804,     0.218689,     1.512175,    -3.197075,     3.440555), \
F(x, 0.309804,     0.458824,    -0.261191,     5.401505,   -14.030694,    14.113778), \
F(x, 0.458824,     0.603922,    -1.534715,    10.091878,   -15.985773,     9.257180), \
F(x, 0.603922,     0.701961,     4.759391,   -17.980549,    27.364646,   -14.145868), \
F(x, 0.701961,     0.815686,    18.275342,   -66.730000,    86.570221,   -38.633344), \
F(x, 0.815686,     0.901961,   -19.096521,    74.250435,   -90.514567,    35.431661), \
F(x, 0.901961,     1.000000,     0.239399,     0.000000,    -0.000000,     0.000000), \
F(x, 0.000000,     0.101961,     0.575370,    -0.214218,     1.381410,   -10.814899), \
F(x, 0.101961,     0.239216,     0.620844,     0.637717,   -16.547215,    37.635934), \
F(x, 0.239216,     0.329412,    -1.404892,    21.875973,   -87.291389,   110.089969), \
F(x, 0.329412,     0.443137,    -1.127449,    12.499131,   -36.589253,    34.785129), \
F(x, 0.443137,     0.717647,     0.971421,    -4.438096,     8.640928,    -5.173692), \
F(x, 0.717647,     0.847059,    12.089098,   -46.895840,    62.535294,   -27.917903), \
F(x, 0.847059,     0.929412,   -64.437519,   225.686918,  -261.540385,   100.686232), \
F(x, 0.929412,     1.000000,     0.230961,     0.000000,    -0.000000,     0.000000); exit

#plot [] [0:1] $cetr2 \
   u ($1*dx):2 ev nev w p pt 7 ps 2 lc 'red'  , \
'' u ($1*dx):3 ev nev w p pt 7 ps 2 lc 'green', \
'' u ($1*dx):4 ev nev w p pt 7 ps 2 lc 'blue' , \
F(x, 0.000000,     0.137255,    -0.000316,     0.038887,    -0.883031,     5.101004), \
F(x, 0.137255,     0.203922,    -4.023305,    61.355168,  -303.724536,   513.861584), \
F(x, 0.203922,     0.305882,    -1.622921,    20.396258,   -74.241922,    90.162344), \
F(x, 0.305882,     0.396078,     8.471779,   -73.154898,   211.089180,  -195.491355), \
F(x, 0.396078,     0.623529,    -0.543610,     2.677936,    -0.137169,    -0.445552), \
F(x, 0.623529,     0.749020,   -15.225713,    68.550324,   -96.645216,    45.458704), \
F(x, 0.749020,     1.000000,     2.274626,    -4.643168,     5.626198,    -2.267654), \
F(x, 0.000000,     0.207843,     0.204965,     2.709276,    -7.313826,     9.194078), \
F(x, 0.207843,     0.396078,     0.376606,     0.109729,     4.680670,    -7.355355), \
F(x, 0.396078,     0.619608,     0.865401,    -1.956191,     5.447395,    -3.937938), \
F(x, 0.619608,     0.760784,    -8.802577,    42.865907,   -62.175709,    29.088263), \
F(x, 0.760784,     0.992157,     7.266933,   -22.020637,    25.525012,   -10.575381), \
F(x, 0.992157,     1.000000,   -17.632920,    38.955075,   -21.133125,     0.000000), \
F(x, 0.000000,     0.203922,     0.963162,    -2.861736,     0.251083,     4.557713), \
F(x, 0.203922,     0.313725,     0.656233,     0.043612,    -5.790983,     0.106573), \
F(x, 0.313725,     0.400000,    17.430800,  -140.294221,   376.158694,  -334.690862), \
F(x, 0.400000,     0.631373,     0.328764,    -2.247241,     5.614561,    -3.859553), \
F(x, 0.631373,     0.968627,    -0.646066,     4.007488,    -5.858817,     2.485090), \
F(x, 0.968627,     1.000000,    -0.000000,     0.000000,    -0.000000,     0.000000); exit

plot [] [0:1] $viridis \
   u ($1*dx):2 ev nev w p pt 7 ps 2 lc 'red'  , \
'' u ($1*dx):3 ev nev w p pt 7 ps 2 lc 'green', \
'' u ($1*dx):4 ev nev w p pt 7 ps 2 lc 'blue' , \
F(x, 0.000000,     0.455523,     0.282113,     0.585552,    -3.199668,     3.273305), \
F(x, 0.455523,     0.700880,    -2.797136,    18.565436,   -37.531934,    24.521874), \
F(x, 0.700880,     1.000000,     4.563142,   -18.410284,    23.559766,    -8.759172), \
F(x, 0.000000,     0.074291,    -0.000438,    -0.351490,    48.354922,  -360.945985), \
F(x, 0.074291,     0.702835,    -0.014997,     1.564580,    -1.154700,     0.700928), \
F(x, 0.702835,     1.000000,    -0.960756,     4.427797,    -3.482215,     0.926468), \
F(x, 0.000000,     0.476051,     0.327782,     1.734561,    -4.300750,     3.500006), \
F(x, 0.476051,     0.932551,     0.218208,     1.243982,    -0.736753,    -0.849211), \
F(x, 0.932551,     1.000000,   348.076269, -1047.336499,  1048.282355,  -348.910600); exit

#plot [] [0:1] $parula \
   u ($1*dx):2 ev nev w p pt 7 ps 2 lc 'red'  , \
'' u ($1*dx):3 ev nev w p pt 7 ps 2 lc 'green', \
'' u ($1*dx):4 ev nev w p pt 7 ps 2 lc 'blue' , \
F(x, 0.000000,     0.115347,     0.209319,     0.021835,     8.049552,  -184.003036), \
F(x, 0.115347,     0.184751,     2.374031,   -44.761896,   276.704512,  -556.160051), \
F(x, 0.184751,     0.317693,    -0.399346,     3.459321,    -4.534020,    -6.611694), \
F(x, 0.317693,     0.447703,    -2.173981,    20.360889,   -61.942289,    62.029278), \
F(x, 0.447703,     0.675464,     4.243855,   -25.971352,    50.451225,   -29.381695), \
F(x, 0.675464,     0.840665,     1.758620,    -8.829890,    15.820527,    -7.604765), \
F(x, 0.840665,     1.000000,   -47.324025,   160.987718,  -178.222332,    65.537864), \
F(x, 0.000000,     0.122190,     0.166515,     1.467473,    -0.294203,    23.661735), \
F(x, 0.122190,     0.328446,     0.142739,     2.829204,    -8.770193,    13.750274), \
F(x, 0.328446,     0.695015,     0.308801,     0.881104,     0.574870,    -1.368778), \
F(x, 0.695015,     0.842620,   -14.989531,    63.658741,   -85.592272,    38.217482), \
F(x, 0.842620,     1.000000,   -14.773917,    50.485538,   -55.874747,    21.147001), \
F(x, 0.000000,     0.109482,     0.530296,     2.857383,     6.806557,   -41.144233), \
F(x, 0.109482,     0.232649,     0.465835,     7.313556,   -40.797955,    70.197547), \
F(x, 0.232649,     0.384164,     1.953887,   -11.423272,    39.003550,   -44.999622), \
F(x, 0.384164,     0.675464,    -0.247051,     7.832808,   -17.887612,    11.458347), \
F(x, 0.675464,     0.862170,     9.308603,   -35.135626,    47.027607,   -21.482656), \
F(x, 0.862170,     1.000000,    14.096229,   -43.850498,    47.004022,   -17.195887); exit

#plot [] [0:1] $bgy \
   u ($1*dx):2 ev nev w p pt 7 ps 2 lc 'red'  , \
'' u ($1*dx):3 ev nev w p pt 7 ps 2 lc 'green', \
'' u ($1*dx):4 ev nev w p pt 7 ps 2 lc 'blue' , \
F(x, 0.000000,     0.250980,     0.121646,    -0.172462,    -1.127437,     5.431302), \
F(x, 0.250980,     0.392157,    -0.131357,     2.516024,    -9.597266,    12.395703), \
F(x, 0.392157,     0.556863,    -1.492039,    10.422082,   -22.967599,    17.609067), \
F(x, 0.556863,     0.713725,    -5.895145,    28.778858,   -46.038523,    25.322913), \
F(x, 0.713725,     0.847059,    -8.586891,    32.582707,   -40.596624,    17.627591), \
F(x, 0.847059,     1.000000,    -1.142666,     2.057658,    -0.000297,     0.000113), \
F(x, 0.000000,     0.250980,     0.009478,     1.567212,     1.134523,    -5.465413), \
F(x, 0.250980,     0.392157,     0.263542,    -1.131873,     9.632994,   -12.442064), \
F(x, 0.392157,     0.592157,     0.200958,     0.462406,     1.933537,    -2.125990), \
F(x, 0.592157,     1.000000,    -0.091472,     2.264930,    -1.894472,     0.618395), \
F(x, 0.000000,     0.188235,     0.399324,     0.574092,     1.663684,    -7.302784), \
F(x, 0.188235,     0.407843,     0.584562,    -1.180634,     5.786781,    -7.789910), \
F(x, 0.407843,     1.000000,     0.380303,     1.004679,    -1.825543,     0.793271); exit

#plot [] [0:1] $avocado \
   u ($1*dx):2 ev nev w p pt 7 ps 2 lc 'red'  , \
'' u ($1*dx):3 ev nev w p pt 7 ps 2 lc 'green', \
'' u ($1*dx):4 ev nev w p pt 7 ps 2 lc 'blue' , \
F(x, 0.000000,     0.254902,    -0.000315,     0.019821,    -0.238792,     0.738267), \
F(x, 0.254902,     0.521569,    -0.384667,     1.975454,    -2.290351,     2.092677), \
F(x, 0.521569,     0.784314,     0.251672,    -2.065961,     5.852542,    -3.112441), \
F(x, 0.784314,     1.000000,    -0.264265,     1.264787,    -0.000831,     0.000310), \
F(x, 0.000000,     0.258824,     0.000502,     1.740189,     0.372609,    -1.140405), \
F(x, 0.258824,     0.521569,     0.294585,     0.168953,     2.224051,    -2.019746), \
F(x, 0.521569,     1.000000,     0.427054,     0.465070,     0.118419,    -0.025793), \
F(x, 0.000000,     0.407843,    -0.000319,     0.297254,     0.186624,    -0.857099), \
F(x, 0.407843,     0.886275,    -0.121618,     1.027933,    -1.611968,     0.962846), \
F(x, 0.886275,     1.000000,    -0.107974,     0.338385,     0.000000,    -0.000000); exit

#plot [:1] [0:1] $moreland \
   u ($1*dx):2 ev nev w p pt 7 ps 2 lc 'red'  , \
'' u ($1*dx):3 ev nev w p pt 7 ps 2 lc 'green', \
'' u ($1*dx):4 ev nev w p pt 7 ps 2 lc 'blue' , \
F(x, 0.000000,     0.616813,     0.337271,     0.983489,     0.747265,    -1.195076), \
F(x, 0.616813,     1.000000,    -0.632934,     5.066503,    -4.609700,     0.870305), \
F(x, 0.000000,     0.522972,     0.283621,     1.787108,    -0.273675,    -1.969690), \
F(x, 0.522972,     0.953079,     0.509720,     2.032810,    -2.808791,     0.356746), \
F(x, 0.953079,     1.000000,  1127.813666, -3504.201247,  3632.816785, -1256.425801), \
F(x, 0.000000,     0.577713,     0.753360,     1.710547,    -2.941316,    -0.067188), \
F(x, 0.577713,     1.000000,     1.240366,     0.386625,    -3.038613,     1.566723); exit

#plot [] [0:1] $thermometer \
   u ($1*dx):2 ev nev w p pt 7 ps 2 lc 'red'  , \
'' u ($1*dx):3 ev nev w p pt 7 ps 2 lc 'green', \
'' u ($1*dx):4 ev nev w p pt 7 ps 2 lc 'blue' , \
F(x, 0.000000,     0.125490,     0.162174,     1.201690,    -6.399066,    51.243518), \
F(x, 0.125490,     0.431373,     0.148005,     0.893940,     4.075881,    -6.121149), \
F(x, 0.431373,     0.549020,    -4.033245,    27.676482,   -53.245771,    34.957323), \
F(x, 0.549020,     0.682353,     6.764309,   -30.526849,    52.649881,   -30.087859), \
F(x, 0.682353,     0.823529,    -9.803570,    42.258774,   -54.801709,    23.209690), \
F(x, 0.823529,     0.964706,    11.780482,   -36.707047,    42.042428,   -16.601556), \
F(x, 0.964706,     1.000000,     2.266182,    -1.864046,     0.199584,    -0.067639), \
F(x, 0.000000,     0.149020,     0.120337,     2.121807,    -1.069806,    14.995654), \
F(x, 0.149020,     0.274510,    -0.223638,     7.153696,   -21.865513,    31.179040), \
F(x, 0.274510,     0.372549,     1.546500,    -9.998720,    36.221021,   -38.417792), \
F(x, 0.372549,     0.466667,     4.898614,   -30.601930,    76.109977,   -61.896098), \
F(x, 0.466667,     0.576471,    10.130934,   -54.795573,   108.692312,   -72.115484), \
F(x, 0.576471,     0.670588,    15.937621,   -74.007429,   122.359430,   -68.317410), \
F(x, 0.670588,     0.776471,     4.301756,   -14.376467,    21.324263,   -11.667634), \
F(x, 0.776471,     0.960784,   -11.106730,    45.205518,   -56.022964,    22.036255), \
F(x, 0.960784,     1.000000,     1.807516,    -1.709471,    -0.019174,     0.006442), \
F(x, 0.000000,     0.133333,     0.794128,     0.966570,     4.388059,   -37.453981), \
F(x, 0.133333,     0.286275,     0.720815,     2.417972,    -8.750889,    10.817765), \
F(x, 0.286275,     0.419608,     1.358732,    -3.883291,    12.716874,   -14.600037), \
F(x, 0.419608,     1.000000,     0.636442,     2.380296,    -5.193566,     2.345535); exit

#plot [] [0:1] $redbluetones \
   u ($1*dx):2 ev nev w p pt 7 ps 2 lc 'red'  , \
'' u ($1*dx):3 ev nev w p pt 7 ps 2 lc 'green', \
'' u ($1*dx):4 ev nev w p pt 7 ps 2 lc 'blue' , \
F(x, 0.000000,     0.274510,     0.449268,     1.574639,    -0.822878,    -1.491296), \
F(x, 0.274510,     0.403922,    -0.288228,     8.474660,   -22.200837,    20.561929), \
F(x, 0.403922,     0.509804,     3.707488,   -19.016191,    42.456321,   -31.670207), \
F(x, 0.509804,     0.635294,     9.680930,   -46.704589,    83.339750,   -50.422752), \
F(x, 0.635294,     0.878431,    -2.137951,    13.824286,   -20.185904,     8.656232), \
F(x, 0.878431,     1.000000,    34.681746,  -106.281568,   110.295455,   -38.557173), \
F(x, 0.000000,     0.113725,     0.156994,     1.203165,    -4.497128,    33.013271), \
F(x, 0.113725,     0.321569,     0.143463,     0.889417,     4.005165,    -5.745487), \
F(x, 0.321569,     0.435294,     2.876165,   -21.447213,    64.655291,   -60.582508), \
F(x, 0.435294,     0.529412,    10.474115,   -64.301217,   140.750865,  -101.324183), \
F(x, 0.529412,     0.635294,    11.899055,   -59.479146,   106.631277,   -63.671925), \
F(x, 0.635294,     0.752941,     3.869591,   -13.953571,    22.184265,   -12.226311), \
F(x, 0.752941,     0.905882,    -7.963197,    32.488971,   -39.126611,    15.003588), \
F(x, 0.905882,     1.000000,     2.422365,    -2.104851,    -0.009003,     0.003156), \
F(x, 0.000000,     0.184314,     0.219042,     0.690306,     0.608545,     2.946494), \
F(x, 0.184314,     0.529412,     0.260646,     0.314043,     2.383175,    -2.250743), \
F(x, 0.529412,     0.639216,     8.263630,   -42.062465,    77.213158,   -46.359035), \
F(x, 0.639216,     0.741176,    11.232975,   -48.568134,    75.137634,   -38.554181), \
F(x, 0.741176,     0.850980,     4.206315,   -14.675529,    21.410637,   -10.502424), \
F(x, 0.850980,     0.929412,    72.226097,  -246.169093,   283.448543,  -109.133180), \
F(x, 0.929412,     1.000000,     2.153883,    -1.589228,    -0.021327,     0.007324); exit

#plot [] [0:1] $tmap \
   u ($1*dx):2 ev nev w p pt 7 ps 2 lc 'red'  , \
'' u ($1*dx):3 ev nev w p pt 7 ps 2 lc 'green', \
'' u ($1*dx):4 ev nev w p pt 7 ps 2 lc 'blue' , \
F(x, 0.000000,     0.184314,     0.178983,     1.552415,    -0.207975,     5.223602), \
F(x, 0.184314,     0.341176,     0.066012,     2.264477,     0.673474,    -2.095856), \
F(x, 0.341176,     0.447059,     4.749953,   -33.884109,    94.216232,   -83.714029), \
F(x, 0.447059,     0.670588,    -0.868992,     8.918165,   -14.178095,     7.483969), \
F(x, 0.670588,     1.000000,    -0.109994,     4.838318,    -6.513456,     2.604411), \
F(x, 0.000000,     0.184314,     0.303925,     1.772992,    -2.122434,     6.139617), \
F(x, 0.184314,     0.349020,     0.277610,     1.481387,     2.270264,    -4.510611), \
F(x, 0.349020,     0.450980,     4.078043,   -27.598191,    76.434953,   -67.732938), \
F(x, 0.450980,     0.666667,    -0.655808,     8.086928,   -13.137827,     7.055502), \
F(x, 0.666667,     0.756863,    17.528872,   -67.906478,    94.487329,   -44.770083), \
F(x, 0.756863,     0.921569,     4.773772,   -10.130974,     9.493433,    -3.915789), \
F(x, 0.921569,     1.000000,     3.920833,    -3.780528,    -0.009403,     0.003226), \
F(x, 0.000000,     0.419608,     0.935393,     0.005125,     0.678124,    -0.845400), \
F(x, 0.419608,     0.505882,     6.416300,   -34.137703,    72.988688,   -53.484807), \
F(x, 0.505882,     0.588235,    11.099459,   -52.804586,    93.846613,   -57.963352), \
F(x, 0.588235,     0.666667,     2.844359,    -3.625769,    -0.006251,     0.003290), \
F(x, 0.666667,     0.752941,    -6.371378,    32.213485,   -48.309479,    22.921856), \
F(x, 0.752941,     1.000000,     1.358718,    -2.895860,     2.697215,    -0.996030); exit

#plot [] [0:1] $ltmap \
   u ($1*dx):2 ev nev w p pt 7 ps 2 lc 'red'  , \
'' u ($1*dx):3 ev nev w p pt 7 ps 2 lc 'green', \
'' u ($1*dx):4 ev nev w p pt 7 ps 2 lc 'blue' , \
F(x, 0.000000,     0.533333,     0.163818,     1.791632,    -0.273066,    -0.474513), \
F(x, 0.533333,     0.568627,   132.161025,  -731.752508,  1358.253910,  -838.937696), \
F(x, 0.568627,     1.000000,     0.374538,     2.632874,    -3.441293,     1.279874), \
F(x, 0.000000,     0.121569,     0.283064,     2.509624,     3.097213,   -20.814460), \
F(x, 0.121569,     0.262745,     0.526602,    -1.130885,    18.376538,   -36.628242), \
F(x, 0.262745,     0.368627,     1.975347,   -13.981096,    52.456257,   -59.971000), \
F(x, 0.368627,     0.482353,     4.081605,   -24.751935,    63.659714,   -53.109146), \
F(x, 0.482353,     0.607843,     3.013910,   -12.490787,    25.671804,   -17.528325), \
F(x, 0.607843,     0.698039,    22.796447,  -103.139979,   163.020009,   -86.219083), \
F(x, 0.698039,     0.815686,    21.113658,   -80.795116,   108.914800,   -49.613536), \
F(x, 0.815686,     0.933333,    21.071994,   -69.419374,    80.770813,   -32.128173), \
F(x, 0.933333,     1.000000,     2.643322,    -2.266569,    -0.011028,     0.003803), \
F(x, 0.000000,     0.450980,     0.936797,     0.075986,     0.411047,    -0.704823), \
F(x, 0.450980,     0.909804,     1.411502,    -0.307129,    -1.841263,     0.895902), \
F(x, 0.909804,     1.000000,     1.269224,    -1.085101,     0.010742,    -0.003716); exit

#plot [0.:] [0.:1.1] $cividis \
   u ($1*dx):2 ev nev w p pt 7 ps 2 lc 'red'  , \
'' u ($1*dx):3 ev nev w p pt 7 ps 2 lc 'green', \
'' u ($1*dx):4 ev nev w p pt 7 ps 2 lc 'blue' , \
F(x, 0.000000,    0.121569,   -0.000426,    0.060312,    -1.553906,   10.160156), \
F(x, 0.121569,    0.196078,   -1.982815,   32.605251,  -172.417296,  319.227164), \
F(x, 0.196078,    0.666667,   -0.133935,    2.107203,    -2.734498,    2.026416), \
F(x, 0.666667,    0.988235,    0.321580,   -0.218041,     1.387263,   -0.470313), \
F(x, 0.988235,    1.000000,    1.000001,   -0.000002,     0.000002,   -0.000001), \
F(x, 0.000000,    1.000000,    0.125938,    0.709731,    -0.053306,    0.135030), \
F(x, 0.000000,    0.078431,    0.302000,    1.452622,     8.254845,  -63.440633), \
F(x, 0.078431,    0.509804,    0.463205,   -0.353857,     0.722687,    0.040738), \
F(x, 0.509804,    0.976471,    0.546672,   -0.581206,     1.466601,   -1.184898), \
F(x, 0.976471,    0.992157,-6481.511854,19864.890017,-20291.986582, 6908.908418), \
F(x, 0.992157,    1.000000,   15.165100,  -31.148250,    16.256250,    0.000000); exit

# 拟合代码

array Beg[N]; array End[N];
array A0[N];  array A1[N]; array A2[N]; array A3[N];
n=1; idx=1
while(n<N) {
	beg=(n-1)*dt
	i=N; FIT_STDFIT=1E99

	if(N-n<Nord) {
		a0=1; a1=0; a2=0; a3=0;
		if(N-n==0){             fit [beg:1] f(x) data u ($1*dx):col via a0 }
		if(N-n==1){ a1=1;       fit [beg:1] f(x) data u ($1*dx):col via a0,a1 }
		if(N-n==2){ a1=1; a2=1; fit [beg:1] f(x) data u ($1*dx):col via a0,a1,a2 }
		Beg[idx]=beg; End[idx]=1; FIT_STDFIT=0
		A0[idx]=a0; A1[idx]=a1; A2[idx]=a2; A3[idx]=a3
	}

	while(i>n && FIT_STDFIT>maxSTD) {
		end=(i-1)*dt
		a0=1; a1=1; a2=1; a3=1;
		fit [beg:end] f(x) data u ($1*dx):col via a0,a1,a2,a3

		print ">>>> ", n, i, i-n, beg, end, FIT_STDFIT

		if(FIT_STDFIT<maxSTD || i-n==Nord) {
			FIT_STDFIT=0
			Beg[idx]=beg; End[idx]=end
			A0[idx]=a0; A1[idx]=a1; A2[idx]=a2; A3[idx]=a3
		} else {
			i=i-1
		}
	}
	n=i
	idx=idx+1
}

# 输出拟合结果
do for [i=1:idx] {
	print sprintf("F(x, %13.6f,%13.6f,%13.6f,%13.6f,%13.6f,%13.6f)", \
		Beg[i], End[i], A0[i], A1[i], A2[i], A3[i])
}

示例

python常用的viridis颜色映射来说吧, 原始数据和拟合结果对比如下

常见映射方案的颜色标尺如下

解析表达式的单函数形式

如果所用的编程语言支持 ?: 三元条件运算, 可以将得到的拟合结果写成单个函数的形式, 例如viridis对应的函数三个函数如下:

r=  x<=0.455523 ? cubic(x,   0.282113,     0.585552,   -3.199668,    3.273305) : \
    x<=0.700880 ? cubic(x,  -2.797136,    18.565436,  -37.531934,   24.521874) : \
    x<=1.000000 ? cubic(x,   4.563142,   -18.410284,   23.559766,   -8.759172) : 0

g=  x<=0.074291 ? cubic(x,  -0.000438,    -0.351490,   48.354922, -360.945985) : \
    x<=0.702835 ? cubic(x,  -0.014997,     1.564580,   -1.154700,    0.700928) : \
    x<=1.000000 ? cubic(x,  -0.960756,     4.427797,   -3.482215,    0.926468) : 0

b=  x<=0.476051 ? cubic(x,   0.327782,     1.734561,   -4.300750,    3.500006) : \
    x<=0.932551 ? cubic(x,   0.218208,     1.243982,   -0.736753,   -0.849211) : \
    x<=1.000000 ? cubic(x, 348.076269, -1047.336499, 1048.282355, -348.910600) : 0

python不支持 ?: 运算符, 或者老老实实地if elif, 或者使用一些奇技淫巧.

随意赞赏

微信

支付宝
◆本文地址: , 转载请注明◆
◆评论问题: https://jerkwin.herokuapp.com/category/3/博客, 欢迎留言◆


前一篇: 谷歌的改版
后一篇: GROMACS的蛋白构象聚类分析

访问人次(2015年7月 9日起): | 最后更新: 2019-08-05 22:56:18 UTC | 版权所有 © 2008 - 2019 Jerkwin