From 17742e3d61bc5f81a81a12e223d047d41512eba1 Mon Sep 17 00:00:00 2001 From: Recolic Keghart <root@recolic.net> Date: Mon, 11 Dec 2017 15:56:33 +0800 Subject: [PATCH] Support smooth curly line fitting without any poly-fitting limits! --- .gitignore | 1 + README.md | 5 ++ ThermalConductivity/draw1.py | 9 +-- TuningFork/draw.py | 13 +++++ TuningFork/var_m.dat | 4 ++ TuningFork/with_damping.dat | 28 +++++++++ TuningFork/without_damping.dat | 32 +++++++++++ quickmap.py | 101 +++++++++++++++++++++++++++------ test.py | 15 ++++- 9 files changed, 184 insertions(+), 24 deletions(-) create mode 100755 TuningFork/draw.py create mode 100644 TuningFork/var_m.dat create mode 100644 TuningFork/with_damping.dat create mode 100644 TuningFork/without_damping.dat diff --git a/.gitignore b/.gitignore index c380da6..d5bf0a6 100644 --- a/.gitignore +++ b/.gitignore @@ -1,2 +1,3 @@ **/__pycache__/** **/.vscode/** +test.dat diff --git a/README.md b/README.md index 85dc3b9..621dc1b 100644 --- a/README.md +++ b/README.md @@ -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 +''' diff --git a/ThermalConductivity/draw1.py b/ThermalConductivity/draw1.py index 3103ff6..dfb4a7b 100755 --- a/ThermalConductivity/draw1.py +++ b/ThermalConductivity/draw1.py @@ -1,14 +1,15 @@ #!/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) diff --git a/TuningFork/draw.py b/TuningFork/draw.py new file mode 100755 index 0000000..a3aafac --- /dev/null +++ b/TuningFork/draw.py @@ -0,0 +1,13 @@ +#!/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) + + diff --git a/TuningFork/var_m.dat b/TuningFork/var_m.dat new file mode 100644 index 0000000..c28a390 --- /dev/null +++ b/TuningFork/var_m.dat @@ -0,0 +1,4 @@ +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 diff --git a/TuningFork/with_damping.dat b/TuningFork/with_damping.dat new file mode 100644 index 0000000..a86c072 --- /dev/null +++ b/TuningFork/with_damping.dat @@ -0,0 +1,28 @@ +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 diff --git a/TuningFork/without_damping.dat b/TuningFork/without_damping.dat new file mode 100644 index 0000000..a1d3031 --- /dev/null +++ b/TuningFork/without_damping.dat @@ -0,0 +1,32 @@ +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 diff --git a/quickmap.py b/quickmap.py index 6c455bc..b036571 100644 --- a/quickmap.py +++ b/quickmap.py @@ -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 diff --git a/test.py b/test.py index 2054255..6f6ccf3 100644 --- a/test.py +++ b/test.py @@ -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) -- GitLab