Skip to content
Snippets Groups Projects
Unverified Commit 17742e3d authored by Recolic Keghart's avatar Recolic Keghart
Browse files

Support smooth curly line fitting without any poly-fitting limits!

parent 3f0eb66b
No related branches found
No related tags found
No related merge requests found
**/__pycache__/**
**/.vscode/**
test.dat
......@@ -3,3 +3,8 @@
Pick what you want.
You can also get data processor at https://recolic.net/phy and https://recolic.net/phy2.
''' sh
Tip: folder name with leading lowercase letter indicates: this experiment was made at first semester.
with leading uppercase letter : second
'''
#!/usr/bin/env python3
minutes = [4, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30, 32, 34, 36, 38, 40, 42, 44, 46, 48, 50, 52, 54, 56, 58, 60, 62, 64, 66, 68, 70, 72, 74, 76, 78, 80, 82, 84, 86, 88, 90, 92]
minutes = [0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30, 32, 34, 36, 38, 40, 42, 44, 46, 48, 50, 52, 54, 56, 58, 60, 62, 64, 66, 68, 70, 72, 74, 76, 78, 80, 82, 84, 86, 88, 90, 92]
rawTa = [18.7, 36.0, 47.8, 50.2, 50.3, 50.2, 50.3, 50.4, 50.4, 50.5, 50.5, 50.5, 50.5, 50.6, 50.5, 50.6, 50.5, 50.7, 50.8, 50.8, 50.8, 50.8, 50.8, 50.8, 50.9, 50.9, 50.9, 50.9, 50.9, 50.9, 50.9, 50.9, 50.9, 51.0, 51.0, 51.0, 51.1, 51.1, 51.1, 51.1, 51.1, 51.1, 51.1, 51.1, 51.1, 51.2, 51.2]
rawTa = [Ti - 0.3 for Ti in rawTa]
rawTc = [20.0, 21.6, 23.3, 25.0, 26.8, 28.3, 29.4, 30.7, 31.8, 32.9, 33.5, 34.3, 34.9, 35.6, 36.1, 36.7, 37.0, 37.5, 37.8, 38.2, 38.5, 38.8, 39.0, 39.3, 39.5, 39.6, 39.8, 39.9, 40.1, 40.2, 40.3, 40.4, 40.6, 40.6, 40.7, 40.8, 40.8, 40.9, 41.0, 41.0, 41.0, 41.0, 41.0, 41.0, 41.0]
rawTc = [19.4, 19.4, 20.0, 21.6, 23.3, 25.0, 26.8, 28.3, 29.4, 30.7, 31.8, 32.9, 33.5, 34.3, 34.9, 35.6, 36.1, 36.7, 37.0, 37.5, 37.8, 38.2, 38.5, 38.8, 39.0, 39.3, 39.5, 39.6, 39.8, 39.9, 40.1, 40.2, 40.3, 40.4, 40.6, 40.6, 40.7, 40.8, 40.8, 40.9, 41.0, 41.0, 41.0, 41.0, 41.0, 41.0, 41.0]
import sys
sys.path.append("..")
import quickmap
# Junk old version poly fitting......
##for powerTest in range(16):
## quickmap.GetMap(minutes, rawTa, line=True, maxXPower=powerTest)#, inverseK=True)
##for powerTest in range(16):
......@@ -18,5 +19,5 @@ import quickmap
# quickmap.GetMap(minutes, rawTc, line=True, maxXPower=powerTest)#, inverseK=True)
##for powerTest in range(16):
## quickmap.GetMap(minutes, rawTc, line=True, maxXPower=powerTest, inverseK=True)
quickmap.GetMap(minutes, rawTa)
quickmap.GetMap(minutes, rawTc, line=True, maxXPower=9)
quickmap.GetMap(minutes, rawTa, smoothLine=True)
quickmap.GetMap(minutes, rawTc, smoothLine=True)#line=True, maxXPower=9)
#!/usr/bin/env python3
import sys
sys.path.append('..')
import quickmap
x,y = quickmap.DataFileToXYArray('without_damping.dat')
quickmap.GetMap(x,y,smoothLine=True)
x,y = quickmap.DataFileToXYArray('with_damping.dat')
quickmap.GetMap(x,y,smoothLine=True)
0 5*2 10*2 15*2 20*2 25*2
261.417 227.469 203.045 185.083 170.032 159.016
1.496 1.493 1.346 1.320 1.224 0.962
14.633 19.327 24.256 29.192 34.589 39.547
259.176 0.069
259.576 0.081
259.976 0.101
260.376 0.140
260.876 0.328
261.076 0.777
261.096 0.879
261.116 0.992
261.136 1.105
261.156 1.193
261.166 1.213
261.176 1.222 #
261.186 1.208
261.196 1.183
261.216 1.094
261.236 0.965
261.266 0.806
261.306 0.631
261.346 0.510
261.386 0.424
261.466 0.306
261.546 0.237
261.646 0.183
261.846 0.121
262.146 0.077
262.546 0.048
262.946 0.033
263.176 0.027
259.417 0.076
259.617 0.082
259.817 0.089
260.217 0.112
260.617 0.156
261.017 0.286
261.117 0.371
261.167 0.436
261.207 0.505
261.247 0.610
261.287 0.752
261.317 0.911
261.337 1.043
261.357 1.194
261.377 1.349
261.397 1.466
261.407 1.496
261.417 1.498 #
261.427 1.473
261.437 1.419
261.457 1.290
261.477 1.144
261.507 0.936
261.547 0.729
261.587 0.586
261.657 0.427
261.777 0.280
262.037 0.152
262.337 0.095
262.737 0.059
263.137 0.040
263.417 0.032
......@@ -7,6 +7,7 @@ import numpy
from scipy.optimize import leastsq
import matplotlib.pyplot as plt
from matplotlib import rcParams
from scipy.interpolate import spline
def dotMultiply(vctA, vctB):
if len(vctA) != len(vctB):
......@@ -17,42 +18,49 @@ def dotMultiply(vctA, vctB):
ans += a * b
return ans
def GetMap(parrX, parrY, windowSizeX=12, windowSizeY=8, extendXRate=1, extendYRate=1, line=False, passO=False, maxXPower=1, inverseK=False):
def GetMap(arrX, arrY, windowSizeX=12, windowSizeY=8, extendXRate=1, extendYRate=1, polyLine=False, poly_passO=False, poly_maxXPower=1, poly_inverseK=False, smoothLine=False):
'''
Arguments:
parrX and parrY: array of coordinates of points. Ex: GetMap([1,2,3,4,5], [1,2,3,4,5]) -> y=x
line: Should I draw a fitting line?
arrX and arrY: array of coordinates of points. Ex: GetMap([1,2,3,4,5], [1,2,3,4,5]) -> y=x
polyLine(conflict with polyLine): Should I draw a fitting line by polynomial fitting?
smoothLine(conflict with smoothLine): Should I draw a fitting line by connecting all points and smooth them?
passO: Should the fitting line pass (0,0)? That's saying, should k0 be zero?
maxXPower: If I should draw a fitting line, what polynomial function should I use? Ex: GetMap([1,2,3,4,5], [1,4,9,16,25], line=True, maxXPower) -> y=x^2
inverseK: Usually, I'm drawing the curl `y=KX`, while K=[k0,k1,k2,...], X=[x^0,x^1,x^2,...]. If this switch is set, I'll drawing the curl `KY=x`. Don't worry, this switch is transparent to you.
Ex: GetMap([0,1,1,4,4,9,9], [0,1,-1,2,-2,3,-3], maxXPower=2, line=True, inverseK=True) -> y^2=x
Ex: GetMap([0,1,1,4,4,9,9], [0,1,-1,2,-2,3,-3], poly_maxXPower=2, polyLine=True, poly_inverseK=True) -> y^2=x
ReturnValue:
void
'''
if polyLine and smoothLine:
raise RuntimeError("bad argument: polyLine and smoothLine can't been set simultaneously.")
if len(arrX) != len(arrY):
raise ValueError("arrX and arrY must in same size, but {} != {}.".format(len(arrX), len(arrY)))
print('Generic-GetMap by Recolic.')
arrX, arrY = parrX, parrY
maxX, maxY = max(arrX) * extendXRate, max(arrY) * extendYRate
minX, minY = min(arrX) * extendXRate, min(arrY) * extendYRate
# Do calculate
# y = [k0 k1 k2 ...] dot [x^0 x^1 x^2 ...]
print('Your input: ', arrX, '|', arrY)
print('Data collection done. Generating result...')
X, Y = numpy.array(arrX), numpy.array(arrY)
print('Your input: ', arrX, '|', arrY)
if polyLine:
print('You want a polynomial fitting, with maxXPower = {}, passO = {}, inverseK = {}.'.format(poly_maxXPower, poly_passO, poly_inverseK))
if smoothLine:
print('You want a smooth fitting.')
print('Generating result...')
# y = [k0 k1 k2 ...] dot [x^0 x^1 x^2 ...]
def lineFunc(k, x):
vctX = [x ** power for power in range(maxXPower + 1)]
if passO:
vctX = [x ** power for power in range(poly_maxXPower + 1)]
if poly_passO:
vctX[0] = 0
return dotMultiply(k, vctX)
def lossFunc(k, x, y): return abs(lineFunc(k, x) - y)
# Fire!
if line:
kInit = [1 for _ in range(maxXPower + 1)]
if polyLine:
kInit = [1 for _ in range(poly_maxXPower + 1)]
kInit[0] = 0 # guarantee passO.
if inverseK:
if poly_inverseK:
kFinal, _ = leastsq(lossFunc, kInit, args=(Y, X))
print('Fitting line done. k^-1=', kFinal)
else:
......@@ -67,13 +75,17 @@ def GetMap(parrX, parrY, windowSizeX=12, windowSizeY=8, extendXRate=1, extendYRa
rcParams['grid.linewidth'] = 0.2
plt.figure(figsize=(windowSizeX, windowSizeY))
plt.scatter(X, Y, color="red", label="Sample Point", linewidth=3)
if line:
if inverseK:
if polyLine:
if poly_inverseK:
py = numpy.linspace(minY, maxY, 1000)
px = dotMultiply(kFinal, [py ** power for power in range(maxXPower + 1)])
px = dotMultiply(kFinal, [py ** power for power in range(poly_maxXPower + 1)])
else:
px = numpy.linspace(minX, maxX, 1000)
py = dotMultiply(kFinal, [px ** power for power in range(maxXPower + 1)])
py = dotMultiply(kFinal, [px ** power for power in range(poly_maxXPower + 1)])
plt.plot(px, py, color="orange", label="Fitting Line", linewidth=2)
if smoothLine:
px = numpy.linspace(minX, maxX, 1000)
py = spline(X, Y, px)
plt.plot(px, py, color="orange", label="Fitting Line", linewidth=2)
plt.legend()
plt.grid()
......@@ -84,3 +96,56 @@ def GetMap(parrX, parrY, windowSizeX=12, windowSizeY=8, extendXRate=1, extendYRa
return 0.0
else:
return float(sstr)
def _str_remove_extra_space(s):
begin, end = 0, len(s)
for char in s:
if char == ' ':
begin += 1
else:
break
for char in s[::-1]:
if char == ' ':
end -= 1
else:
break
return s[begin:end]
def _str_remove_comments(s, commentSym = '#'):
pos = s.find(commentSym)
if pos == -1:
return s
else:
return s[0:pos]
def DataFileToXYArray(fname, lineDelimiter = '\n', wordDelimiter = ' ', commentSym = '#', _DataType = float):
'''
give file name, return tuple: ([_DataType], [_DataType])
(xArray, yArray)
Default _DataType is python float, and I think that's enough.
Example:
x, y = DataFileToXYArray('test.dat')
GenMap(x, y)
'''
with open(fname, 'r') as fd:
s = fd.read()
xArray, yArray = [], []
for ori_line in s.split(lineDelimiter):
line = _str_remove_comments(ori_line, commentSym)
line = _str_remove_extra_space(line)
if len(line) == '':
continue
ar = line.split(wordDelimiter)
if len(ar) < 2:
print('Warning: Bad data line "{}" skipped.'.format(ori_line))
continue
if len(ar) > 2:
print('Warning: Invalid data line "{}" parsed as "{} {}"'.format(ori_line, ar[0], ar[1]))
try:
xArray.append(_DataType(ar[0]))
yArray.append(_DataType(ar[1]))
except:
print('At data line "{}":'.format(ori_line))
return xArray, yArray
\ No newline at end of file
......@@ -5,7 +5,18 @@
#
#quickmap.GetMap(xrc, yrc, line=True, maxXPower=16, inverseK=True)
#
import virtualtype
#import virtualtype
#virtualtype.virtual_type_array([99.78,8.26,99.82,9.86,99.72,8.4,99.7,9.9,99.72,8.42,99.72,9.88,99.68,8.34,99.72,9.74,99.7,8.2,99.86,9.82])
virtualtype.virtual_type_array([42.5,42.2,41.9,41.5,41.2,40.9,40.5,40.3,40.0,39.7])
\ No newline at end of file
#virtualtype.virtual_type_array([42.5,42.2,41.9,41.5,41.2,40.9,40.5,40.3,40.0,39.7])
import quickmap
x, y = quickmap.DataFileToXYArray('test.dat')
quickmap.GetMap(x,y,polyLine=True)
x=[50,80,110,140,170,200]
y=[63.7,75.0,93.7,112.0,126.1,143.0]
quickmap.GetMap(x,y,smoothLine=True)
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment