我对此有几个想法，并尝试了整体曲率、弧长等以及混合组合。结果证明一个简单的chi**2 效果最好。所以它就这么简单：
- 通过样条线获得一些结以适应每个块的给定精度
- 加入一切
- 减少结以避免接触组中非常接近的结，这会导致较大的曲率。
- 在整个集合上使用最小平方拟合，并在连接和减少的结集合上使用样条来查找变化。
理论上可以玩/修改：
- 样条顺序
- 最小结密度
- 最大结密度
- 如何处理相邻的集合
- 为大间隙添加结
- 等
（注意：在一些随机数据中，splrev 会产生错误信息。由于这些信息大多不是很有帮助，我只能说这段代码不是 100% 健壮的。）
代码如下
import matplotlib.pyplot as plt
import numpy as np
from scipy.interpolate import interp1d, splrep, splev
from scipy.optimize import fmin, leastsq
def reduce_knots( inList, dist ):
outList=[]
addList=[]
for i in inList:
try:
if abs( i - addList[ -1 ] ) < dist:
addList += [ i ]
else:
outList += [ addList ]
addList = [ i ]
except IndexError:### basically the first
addList = [ i]
outList += [ addList ]
return [ sum( x ) / len( x ) for x in outList ]
def adaptive_knots( inX, inY, thresh=.005 ):
ll = len( inX )
sup = ll - 4
assert sup > 3
nN = 3
test = True
while test:
testknots = np.linspace( 1, len( inX ) - 2, nN, dtype=np.int )
testknots = [ inX[ x ] for x in testknots ]
myTCK= splrep( inX , inY, t=testknots )
newY = splev( inX , myTCK )
chi2 = np.sum( ( newY - inY )**2 ) / ll
if chi2 > thresh:
nN += 1
if nN > sup:
test = False
else:
test = False
return testknots
def global_residuals( shiftList, xBlocks, yBlocks, allTheKnots ):# everything shifted (1 is redundant by global offset) Blocks must be ordered an np.arrays
localYBlocks = [ s + yList for s, yList in zip( shiftList, yBlocks ) ]
allTheX = np.concatenate( xBlocks )
allTheY = np.concatenate( localYBlocks )
tck = splrep( allTheX, allTheY, t=allTheKnots )
yList = splev( allTheX, tck )
diff = yList - allTheY
return diff
#~ np.random.seed( 28561 )
np.random.seed( 5561 )
#~ np.random.seed( 733437 )
### python way for test data
knots = np.random.rand( 8, 2 )
knots = np.array( sorted( [ [ 0, np.random.rand() ] ] + list( knots ) + [ [ 1, np.random.rand() ] ], key=lambda x: x[ 0 ] ) )
dummyX = np.linspace( 0, 1, 3e4 )
f = interp1d( knots[ :, 0 ], knots[ :, 1 ], 'cubic' )
dummyY = np.fromiter( ( f( x ) for x in dummyX ), np.float )
chunk = np.append( [ 0 ], np.append( np.sort( np.random.randint( 7, high=len( dummyX ) - 10 , size= 10, dtype=np.int ) ), len( dummyX ) ) )
xDataDict = dict()
yDataDict = dict()
allX = np.array( [] )
allY = np.array( [] )
allK = np.array( [] )
allS = []
for i, val in enumerate(chunk[ : -1 ] ):
if np.random.rand() < .75: ## 25% of not appearing
xDataDict[ i ] = dummyX[ val:chunk[ i + 1 ] ]
realShift = 1.5 * ( 1 - 2 * np.random.rand() )
allS += [ realShift ]
yDataDict[ i ] = dummyY[ val:chunk[ i + 1 ] ] + realShift
yDataDict[ i ] = np.fromiter( ( np.random.normal( scale=.05, loc=y ) for y in yDataDict[ i ] ), np.float )
allX = np.append( allX, xDataDict[ i ] )
allY = np.append( allY, yDataDict[ i ] )
### Plotting
fig = plt.figure()
ax = fig.add_subplot( 3, 1, 1 )
ax.plot( knots[ :, 0 ],knots[ :, 1 ], ls='', c='r', marker='o')
ax.plot( dummyX , dummyY, '--' )
for key in xDataDict.keys():
ax.plot(xDataDict[ key ], yDataDict[ key ] )
myKnots = adaptive_knots( xDataDict[ key ], yDataDict[ key ] )
allK = np.append( allK, myKnots )
myTCK = splrep( xDataDict[ key ], yDataDict[ key ], t=myKnots )
ax.plot( xDataDict[ key ], splev( xDataDict[ key ] , myTCK ) )
myTCK = splrep( allX, allY, t=allK )
ax.plot( allX, splev( allX, myTCK ) )
for x in allK:
ax.axvline( x=x, linestyle=':', color='#AAAAAA', linewidth=1 )
### now fitting
myXBlockList = []
myYBlockList = []
for key in sorted( xDataDict.keys() ):
myXBlockList += [ xDataDict[ key ] ]
myYBlockList += [ yDataDict[ key ] ]
#start values
s = [ 0 ]
for i,y in enumerate( myYBlockList[ :-1 ] ):
ds = myYBlockList[ i + 1 ][ 0 ] - y[ -1 ]
s += [ -ds ]
startShift = np.cumsum( s )
allK = reduce_knots( allK, .01 )
sol, ierr = leastsq( global_residuals, x0=startShift, args=( myXBlockList, myYBlockList, allK ), maxfev=10000 )
sol = np.array(sol) - sol[ 0 ]
print "solution: ", -sol
print "real: ", np.array( allS ) - allS[ 0 ]
### Plotting solutions
bx = fig.add_subplot( 3, 1, 3, sharex=ax )
for x, y, s in zip( myXBlockList, myYBlockList, sol ):
bx.plot( x, y + s )
localYBlocks = [ s + yList for s,yList in zip( sol, myYBlockList ) ]
allTheX = np.concatenate( myXBlockList )
allTheY = np.concatenate( localYBlocks )
tck = splrep( allTheX, allTheY, t=allK )
dx = allTheX[ 1 ] - allTheX[ 0 ]
testX = np.arange( allTheX[ 0 ], allTheX[ -1 ], dx )
finalyList = splev( testX, tck)
bx.plot( testX, finalyList , 'k--' )
mean = sum( dummyY ) / len( dummyY ) - sum( finalyList ) / len( finalyList )
bx.plot( dummyX, dummyY - mean, '--' )
for x in allK:
bx.axvline( x=x, linestyle=':', color='#AAAAAA', linewidth=1 )
cx = fig.add_subplot( 3, 1, 2, sharex=ax )
for x, y, s in zip( myXBlockList, myYBlockList, startShift ):
cx.plot( x, y + s )
plt.show()
对于小的差距，这在测试数据上效果很好
上图将原始样条曲线及其节点显示为红点。这产生了数据。此外，它显示了嘈杂的移位块、初始拟合结为垂直线和相应的样条拟合。
中间图显示了由预先计算的起始值移动的块 - 对齐的末端。
下图显示了原始样条、拟合样条、减少的结位置以及根据拟合解移动的块。
当然，间隙越大，解决方案越偏离原始解决方案
...不过还是不错的。