# -*-Python-*-
# Created by grierson at 19 May 2018 07:01
"""
This script demonstrates the use of the GASpline class for fitting profiles.
defaultVars parameters
----------------------
:param demo: demo number, from 1-14
1) single autoknot spline to core
2) single manual knot spline to core
3) single autoknot spline to core + edge
4) single manual knot spline to core + edge
5) monte-carlo autoknot spline to core
6) monte-carlo manual knot spline to core
7) monte-carlo autoknot spline to core + edge
8) monte-carlo manual knot spline to core + edge
9) monte-carlo manual knot spline to core + edge with fixed y(1)=1.7
10) single auto-knot spline to H-mode profile core
11) single auto-knot spline to H-mode profile core+edge
12) single biased auto-knot spline to H-mode profile core
13) single biased auto-knot spline to H-mode profile core+edge
14) single biased auto-knot spline to H-mode profile core+edge w/continuous y'
15) monte-carlo baised auto-knot spline to H-mode profile core+edge
"""
defaultVars(demo=1)
def mtanh(x, A, B, xsym, hwid, core, edge):
"""
Groebner APS 2004
y = A * MTANH(alpha,z) + B
MTANH(alpha,z) = [(1+alpha*z)*exp(z) - exp(-z)] / [exp(z) + exp(-z)]
generalized to have a arbitrary order poly on numerator exp(z) and exp(-z)
to allow for more variation in the core and edge. Re-formulated as
MTANH(alpha,z) = [nc*exp(z) - ne*exp(-z)] / [exp(z) + exp(-z)]
where nc and ne are the core and edge polynomials and z = (xsym-x)/hwid
Pedestal height is A+B
Pedestal offset is A-B
Location of the knee is xsym-hwid
Location of the foot is xsym+hwid
:param x: x-axis values
:param A: Amplitude
:param B: Offset
:param xsym: Symmetry location
:param hwid: half-width
:param core: core polynomial coefficients [1.0, ...]
:param edge: edge polynomial coefficients [1.0, ...]
:return: Function y(x)
"""
# Z shifted function
z = (xsym - x) / hwid
# Core polynomial
nc = core[0]
for i in range(1, len(core)):
nc += core[i] * z ** (i)
# Edge polynomial
ne = edge[0]
for i in range(1, len(edge)):
ne += edge[i] * z ** (i)
# mtanh in all its glory
mt = (nc * np.exp(z) - ne * np.exp(-z)) / (np.exp(z) + np.exp(-z))
# Final function
y = A * mt + B
return y
# ----------
# Begin demo
# ----------
if demo in [1, 2, 3, 4, 5, 6, 7, 8, 9]:
nxdata = 51
xdata = linspace(0, 1.2, nxdata)
coremask = xdata > 1.0
ydata = np.sin(2.0 * np.pi * xdata) + 2.0
ydata[coremask] = ydata[(xdata < 1.0)][-1] * exp((1 - xdata[coremask]) / 0.05)
yerr = np.repeat(0.1, nxdata)
# Keyword options
if demo == 1:
kw = {'numknot': 4, 'autoknot': True, 'sol': False, 'doPlot': False, 'verbose': True}
if demo == 2:
kw = {'numknot': 4, 'autoknot': False, 'sol': False, 'doPlot': False, 'verbose': True}
if demo == 3:
kw = {'numknot': 4, 'autoknot': True, 'sol': True, 'doPlot': False, 'verbose': True}
if demo == 4:
kw = {'numknot': 4, 'autoknot': False, 'sol': True, 'doPlot': False, 'verbose': True}
if demo == 5:
kw = {'numknot': 4, 'autoknot': True, 'sol': False, 'monte': 10, 'doPlot': False, 'verbose': True}
if demo == 6:
kw = {'numknot': 4, 'autoknot': False, 'sol': False, 'monte': 10, 'doPlot': False, 'verbose': True}
if demo == 7:
kw = {'numknot': 4, 'autoknot': True, 'sol': True, 'monte': 10, 'doPlot': False, 'verbose': True}
if demo == 8:
kw = {'numknot': 4, 'autoknot': False, 'sol': True, 'monte': 100, 'doPlot': False, 'verbose': True}
if demo == 9:
kw = {'numknot': 4, 'autoknot': True, 'sol': True, 'monte': 10, 'y0p': 0.0, 'y1': 1.7, 'doPlot': False, 'verbose': True}
# MTANH data with a fit using an edge bias
if demo in [10, 11, 12, 13, 14]:
xsym = 0.91
hwid = 0.02
core = np.array([1.0, 2e-2])
edge = np.array([1.0, -1e-1, -3e-4])
pedestal = 5.0
offset = 1.0
B = (pedestal + offset) / 2.0
A = B - offset
xdata = linspace(0.0, 1.2, 101)
coremask = xdata > 1.0
ydata = np.abs(mtanh(xdata, A, B, xsym, hwid, core, edge) + 0.5 * random(len(xdata)))
ydata[coremask] = ydata[(xdata < 1.0)][-1] * exp((1 - xdata[coremask]) / 0.05)
yerr = 0.1 * ydata + 0.1
numknot = 5
if demo == 10:
kw = {'numknot': numknot, 'knotbias': 0, 'autoknot': True, 'sol': False}
if demo == 11:
kw = {'numknot': numknot, 'knotbias': 0, 'autoknot': True, 'sol': True}
if demo == 12:
kw = {'numknot': numknot, 'knotbias': 2, 'autoknot': True, 'sol': False}
if demo == 13:
kw = {'numknot': numknot, 'knotbias': 2, 'autoknot': True, 'sol': True}
if demo == 14:
kw = {'numknot': numknot, 'knotbias': 2, 'autoknot': True, 'sol': True, 'solc': True}
if demo == 15:
kw = {'numknot': numknot, 'knotbias': 2, 'autoknot': True, 'monte': 10}
xfit = linspace(0, 1.2, 221)
fit_obj = GASpline(xdata, ydata, yerr, xfit, **kw)
fit_obj.plot()
