[demos] avoid running demos on import

parent 64f2c543
Pipeline #68578 canceled with stages
in 1 minute and 49 seconds
......@@ -14,116 +14,121 @@ import matplotlib.pyplot as plt
from pymor.basic import *
# # Model
def run_demo():
# # Model
# In[ ]:
# In[ ]:
def H(s, mu):
tau = mu['tau'][0]
return np.array([[np.exp(-s) / (tau * s + 1)]])
def H(s, mu):
tau = mu['tau'][0]
return np.array([[np.exp(-s) / (tau * s + 1)]])
def dH(s, mu):
tau = mu['tau'][0]
return np.array([[-(tau * s + tau + 1) * np.exp(-s) / (tau * s + 1) ** 2]])
def dH(s, mu):
tau = mu['tau'][0]
return np.array([[-(tau * s + tau + 1) * np.exp(-s) / (tau * s + 1) ** 2]])
# In[ ]:
# In[ ]:
fom = TransferFunction(1, 1,
H, dH,
parameters={'tau': 1})
fom = TransferFunction(1, 1,
H, dH,
parameters={'tau': 1})
# # Magnitude plot
# # Magnitude plot
# In[ ]:
# In[ ]:
mu_list_short = [0.01, 0.1, 1]
mu_list_short = [0.01, 0.1, 1]
# In[ ]:
# In[ ]:
w = np.logspace(-2, 4, 100)
w = np.logspace(-2, 4, 100)
fig, ax = plt.subplots()
for mu in mu_list_short:
fom.mag_plot(w, ax=ax, mu=mu, label=fr'$\tau = {mu}$')
ax.legend()
plt.show()
fig, ax = plt.subplots()
for mu in mu_list_short:
fom.mag_plot(w, ax=ax, mu=mu, label=fr'$\tau = {mu}$')
ax.legend()
plt.show()
# In[ ]:
# In[ ]:
w_list = np.logspace(-2, 4, 100)
mu_list = np.logspace(-2, 0, 50)
w_list = np.logspace(-2, 4, 100)
mu_list = np.logspace(-2, 0, 50)
fom_w_mu = np.zeros((len(w_list), len(mu_list)))
for i, mu in enumerate(mu_list):
fom_w_mu[:, i] = spla.norm(fom.freq_resp(w_list, mu=mu), axis=(1, 2))
fom_w_mu = np.zeros((len(w_list), len(mu_list)))
for i, mu in enumerate(mu_list):
fom_w_mu[:, i] = spla.norm(fom.freq_resp(w_list, mu=mu), axis=(1, 2))
# In[ ]:
# In[ ]:
fig, ax = plt.subplots()
out = ax.contourf(w_list, mu_list, fom_w_mu.T,
norm=mpl.colors.LogNorm(),
levels=np.logspace(np.log10(fom_w_mu.min()), np.log10(fom_w_mu.max()), 100))
ax.set_xlabel(r'Frequency $\omega$')
ax.set_ylabel(r'Parameter $\mu$')
ax.set_xscale('log')
ax.set_yscale('log')
fig.colorbar(out, ticks=np.logspace(-4, 1, 6))
plt.show()
fig, ax = plt.subplots()
out = ax.contourf(w_list, mu_list, fom_w_mu.T,
norm=mpl.colors.LogNorm(),
levels=np.logspace(np.log10(fom_w_mu.min()), np.log10(fom_w_mu.max()), 100))
ax.set_xlabel(r'Frequency $\omega$')
ax.set_ylabel(r'Parameter $\mu$')
ax.set_xscale('log')
ax.set_yscale('log')
fig.colorbar(out, ticks=np.logspace(-4, 1, 6))
plt.show()
# # TF-IRKA
# # TF-IRKA
# In[ ]:
# In[ ]:
r = 10
roms_tf_irka = []
for mu in mu_list_short:
tf_irka = TFIRKAReductor(fom, mu=mu)
rom = tf_irka.reduce(r, conv_crit='h2', maxit=1000, num_prev=5)
roms_tf_irka.append(rom)
r = 10
roms_tf_irka = []
for mu in mu_list_short:
tf_irka = TFIRKAReductor(fom, mu=mu)
rom = tf_irka.reduce(r, conv_crit='h2', maxit=1000, num_prev=5)
roms_tf_irka.append(rom)
# In[ ]:
# In[ ]:
fig, ax = plt.subplots()
for mu, rom in zip(mu_list_short, roms_tf_irka):
poles = rom.poles()
ax.plot(poles.real, poles.imag, '.', label=fr'$\tau = {mu}$')
ax.set_title("Poles of TF-IRKA's ROMs")
ax.legend()
plt.show()
fig, ax = plt.subplots()
for mu, rom in zip(mu_list_short, roms_tf_irka):
poles = rom.poles()
ax.plot(poles.real, poles.imag, '.', label=fr'$\tau = {mu}$')
ax.set_title("Poles of TF-IRKA's ROMs")
ax.legend()
plt.show()
# In[ ]:
# In[ ]:
fig, ax = plt.subplots()
for mu, rom in zip(mu_list_short, roms_tf_irka):
rom.mag_plot(w, ax=ax, label=fr'$\tau = {mu}$')
ax.set_title("Magnitude plot of TF-IRKA's ROMs")
ax.legend()
plt.show()
fig, ax = plt.subplots()
for mu, rom in zip(mu_list_short, roms_tf_irka):
rom.mag_plot(w, ax=ax, label=fr'$\tau = {mu}$')
ax.set_title("Magnitude plot of TF-IRKA's ROMs")
ax.legend()
plt.show()
# In[ ]:
# In[ ]:
fig, ax = plt.subplots()
for mu, rom in zip(mu_list_short, roms_tf_irka):
(fom - rom).mag_plot(w, ax=ax, mu=mu, label=fr'$\tau = {mu}$')
ax.set_title("Magnitude plot of error systems")
ax.legend()
plt.show()
fig, ax = plt.subplots()
for mu, rom in zip(mu_list_short, roms_tf_irka):
(fom - rom).mag_plot(w, ax=ax, mu=mu, label=fr'$\tau = {mu}$')
ax.set_title("Magnitude plot of error systems")
ax.legend()
plt.show()
if __name__ == "__main__":
run_demo()
This diff is collapsed.
This diff is collapsed.
......@@ -20,211 +20,215 @@ from pymor.reductors.bt import BTReductor
from pymor.reductors.h2 import IRKAReductor
# # Model
#
# https://morwiki.mpi-magdeburg.mpg.de/morwiki/index.php/Synthetic_parametric_model
def run_demo():
# # Model
#
# https://morwiki.mpi-magdeburg.mpg.de/morwiki/index.php/Synthetic_parametric_model
# In[ ]:
# In[ ]:
n = 100 # order of the resulting system
n = 100 # order of the resulting system
# set coefficients
a = -np.linspace(1e1, 1e3, n // 2)
b = np.linspace(1e1, 1e3, n // 2)
c = np.ones(n // 2)
d = np.zeros(n // 2)
# set coefficients
a = -np.linspace(1e1, 1e3, n // 2)
b = np.linspace(1e1, 1e3, n // 2)
c = np.ones(n // 2)
d = np.zeros(n // 2)
# build 2x2 submatrices
aa = np.empty(n)
aa[::2] = a
aa[1::2] = a
bb = np.zeros(n)
bb[::2] = b
# build 2x2 submatrices
aa = np.empty(n)
aa[::2] = a
aa[1::2] = a
bb = np.zeros(n)
bb[::2] = b
# set up system matrices
Amu = sps.diags(aa, format='csc')
A0 = sps.diags([bb, -bb], [1, -1], shape=(n, n), format='csc')
B = np.zeros((n, 1))
B[::2, 0] = 2
C = np.empty((1, n))
C[0, ::2] = c
C[0, 1::2] = d
# set up system matrices
Amu = sps.diags(aa, format='csc')
A0 = sps.diags([bb, -bb], [1, -1], shape=(n, n), format='csc')
B = np.zeros((n, 1))
B[::2, 0] = 2
C = np.empty((1, n))
C[0, ::2] = c
C[0, 1::2] = d
# In[ ]:
# In[ ]:
A0 = NumpyMatrixOperator(A0)
Amu = NumpyMatrixOperator(Amu)
B = NumpyMatrixOperator(B)
C = NumpyMatrixOperator(C)
A0 = NumpyMatrixOperator(A0)
Amu = NumpyMatrixOperator(Amu)
B = NumpyMatrixOperator(B)
C = NumpyMatrixOperator(C)
# In[ ]:
# In[ ]:
A = A0 + Amu * ProjectionParameterFunctional('mu')
A = A0 + Amu * ProjectionParameterFunctional('mu')
# In[ ]:
# In[ ]:
lti = LTIModel(A, B, C)
lti = LTIModel(A, B, C)
# # Magnitude plot
# # Magnitude plot
# In[ ]:
# In[ ]:
mu_list_short = [1/50, 1/20, 1/10, 1/5, 1/2, 1]
mu_list_short = [1/50, 1/20, 1/10, 1/5, 1/2, 1]
# In[ ]:
# In[ ]:
w = np.logspace(0.5, 3.5, 200)
w = np.logspace(0.5, 3.5, 200)
fig, ax = plt.subplots()
for mu in mu_list_short:
lti.mag_plot(w, ax=ax, mu=mu, label=fr'$\mu = {mu}$')
ax.legend()
plt.show()
fig, ax = plt.subplots()
for mu in mu_list_short:
lti.mag_plot(w, ax=ax, mu=mu, label=fr'$\mu = {mu}$')
ax.legend()
plt.show()
# In[ ]:
# In[ ]:
w_list = np.logspace(0.5, 3.5, 200)
mu_list = np.linspace(1/50, 1, 50)
w_list = np.logspace(0.5, 3.5, 200)
mu_list = np.linspace(1/50, 1, 50)
lti_w_mu = np.zeros((len(w_list), len(mu_list)))
for i, mu in enumerate(mu_list):
lti_w_mu[:, i] = spla.norm(lti.freq_resp(w_list, mu=mu), axis=(1, 2))
lti_w_mu = np.zeros((len(w_list), len(mu_list)))
for i, mu in enumerate(mu_list):
lti_w_mu[:, i] = spla.norm(lti.freq_resp(w_list, mu=mu), axis=(1, 2))
# In[ ]:
# In[ ]:
fig, ax = plt.subplots()
out = ax.contourf(w_list, mu_list, lti_w_mu.T,
norm=mpl.colors.LogNorm(),
levels=np.logspace(np.log10(lti_w_mu.min()), np.log10(lti_w_mu.max()), 100))
ax.set_xlabel(r'Frequency $\omega$')
ax.set_ylabel(r'Parameter $\mu$')
ax.set_xscale('log')
#ax.set_yscale('log')
fig.colorbar(out, ticks=np.logspace(-2, 1, 7))
plt.show()
fig, ax = plt.subplots()
out = ax.contourf(w_list, mu_list, lti_w_mu.T,
norm=mpl.colors.LogNorm(),
levels=np.logspace(np.log10(lti_w_mu.min()), np.log10(lti_w_mu.max()), 100))
ax.set_xlabel(r'Frequency $\omega$')
ax.set_ylabel(r'Parameter $\mu$')
ax.set_xscale('log')
#ax.set_yscale('log')
fig.colorbar(out, ticks=np.logspace(-2, 1, 7))
plt.show()
# # Hankel singular values
# # Hankel singular values
# In[ ]:
# In[ ]:
fig, ax = plt.subplots()
for mu in mu_list_short:
hsv = lti.hsv(mu=mu)
ax.semilogy(range(1, len(hsv) + 1), hsv, '.-', label=fr'$\mu = {mu}$')
ax.set_title('Hankel singular values')
ax.legend()
plt.show()
fig, ax = plt.subplots()
for mu in mu_list_short:
hsv = lti.hsv(mu=mu)
ax.semilogy(range(1, len(hsv) + 1), hsv, '.-', label=fr'$\mu = {mu}$')
ax.set_title('Hankel singular values')
ax.legend()
plt.show()
# # System norms
# # System norms
# In[ ]:
# In[ ]:
fig, ax = plt.subplots()
mu_fine = np.linspace(1/50, 1, 20)
h2_norm_mu = [lti.h2_norm(mu=mu) for mu in mu_fine]
ax.plot(mu_fine, h2_norm_mu, '.-', label=r'$\mathcal{H}_2$-norm')
fig, ax = plt.subplots()
mu_fine = np.linspace(1/50, 1, 20)
h2_norm_mu = [lti.h2_norm(mu=mu) for mu in mu_fine]
ax.plot(mu_fine, h2_norm_mu, '.-', label=r'$\mathcal{H}_2$-norm')
if config.HAVE_SLYCOT:
hinf_norm_mu = [lti.hinf_norm(mu=mu) for mu in mu_fine]
ax.plot(mu_fine, hinf_norm_mu, '.-', label=r'$\mathcal{H}_\infty$-norm')
if config.HAVE_SLYCOT:
hinf_norm_mu = [lti.hinf_norm(mu=mu) for mu in mu_fine]
ax.plot(mu_fine, hinf_norm_mu, '.-', label=r'$\mathcal{H}_\infty$-norm')
hankel_norm_mu = [lti.hankel_norm(mu=mu) for mu in mu_fine]
ax.plot(mu_fine, hankel_norm_mu, '.-', label='Hankel norm')
hankel_norm_mu = [lti.hankel_norm(mu=mu) for mu in mu_fine]
ax.plot(mu_fine, hankel_norm_mu, '.-', label='Hankel norm')
ax.set_xlabel(r'$\mu$')
ax.set_title('System norms')
ax.legend()
plt.show()
ax.set_xlabel(r'$\mu$')
ax.set_title('System norms')
ax.legend()
plt.show()
# # Balanced truncation
# # Balanced truncation
# In[ ]:
# In[ ]:
def reduction_errors(lti, r, mu_fine, method):
h2_err_mu = []
hinf_err_mu = []
hankel_err_mu = []
for mu in mu_fine:
rom_mu = method(lti, r, mu=mu)
h2_err_mu.append((lti - rom_mu).h2_norm(mu=mu) / lti.h2_norm(mu=mu))
if config.HAVE_SLYCOT:
hinf_err_mu.append((lti - rom_mu).hinf_norm(mu=mu) / lti.hinf_norm(mu=mu))
hankel_err_mu.append((lti - rom_mu).hankel_norm(mu=mu) / lti.hankel_norm(mu=mu))
return h2_err_mu, hinf_err_mu, hankel_err_mu
def reduction_errors(lti, r, mu_fine, method):
h2_err_mu = []
hinf_err_mu = []
hankel_err_mu = []
for mu in mu_fine:
rom_mu = method(lti, r, mu=mu)
h2_err_mu.append((lti - rom_mu).h2_norm(mu=mu) / lti.h2_norm(mu=mu))
if config.HAVE_SLYCOT:
hinf_err_mu.append((lti - rom_mu).hinf_norm(mu=mu) / lti.hinf_norm(mu=mu))
hankel_err_mu.append((lti - rom_mu).hankel_norm(mu=mu) / lti.hankel_norm(mu=mu))
return h2_err_mu, hinf_err_mu, hankel_err_mu
# In[ ]:
# In[ ]:
r = 20
mu_fine = np.linspace(1/50, 1, 10)
(
h2_bt_err_mu,
hinf_bt_err_mu,
hankel_bt_err_mu
) = reduction_errors(lti, r, mu_fine,
lambda lti, r, mu=None: BTReductor(lti, mu=mu).reduce(r))
r = 20
mu_fine = np.linspace(1/50, 1, 10)
(
h2_bt_err_mu,
hinf_bt_err_mu,
hankel_bt_err_mu
) = reduction_errors(lti, r, mu_fine,
lambda lti, r, mu=None: BTReductor(lti, mu=mu).reduce(r))
# In[ ]:
# In[ ]:
fig, ax = plt.subplots()
ax.semilogy(mu_fine, h2_bt_err_mu, '.-', label=r'$\mathcal{H}_2$')
if config.HAVE_SLYCOT:
ax.semilogy(mu_fine, hinf_bt_err_mu, '.-', label=r'$\mathcal{H}_\infty$')
ax.semilogy(mu_fine, hankel_bt_err_mu, '.-', label='Hankel')
fig, ax = plt.subplots()
ax.semilogy(mu_fine, h2_bt_err_mu, '.-', label=r'$\mathcal{H}_2$')
if config.HAVE_SLYCOT:
ax.semilogy(mu_fine, hinf_bt_err_mu, '.-', label=r'$\mathcal{H}_\infty$')
ax.semilogy(mu_fine, hankel_bt_err_mu, '.-', label='Hankel')
ax.set_xlabel(r'$\mu$')
ax.set_title('Balanced truncation errors')
ax.legend()
plt.show()
ax.set_xlabel(r'$\mu$')
ax.set_title('Balanced truncation errors')
ax.legend()
plt.show()
# # Iterative Rational Krylov Algorithm (IRKA)
# # Iterative Rational Krylov Algorithm (IRKA)
# In[ ]:
# In[ ]:
(
h2_irka_err_mu,
hinf_irka_err_mu,
hankel_irka_err_mu
) = reduction_errors(lti, r, mu_fine,
lambda lti, r, mu=mu: IRKAReductor(lti, mu=mu).reduce(r, conv_crit='h2'))
(
h2_irka_err_mu,
hinf_irka_err_mu,
hankel_irka_err_mu
) = reduction_errors(lti, r, mu_fine,
lambda lti, r, mu=mu: IRKAReductor(lti, mu=mu).reduce(r, conv_crit='h2'))
# In[ ]:
# In[ ]:
fig, ax = plt.subplots()
ax.semilogy(mu_fine, h2_irka_err_mu, '.-', label=r'$\mathcal{H}_2$')
if config.HAVE_SLYCOT:
ax.semilogy(mu_fine, hinf_irka_err_mu, '.-', label=r'$\mathcal{H}_\infty$')
ax.semilogy(mu_fine, hankel_irka_err_mu, '.-', label='Hankel')
fig, ax = plt.subplots()
ax.semilogy(mu_fine, h2_irka_err_mu, '.-', label=r'$\mathcal{H}_2$')
if config.HAVE_SLYCOT:
ax.semilogy(mu_fine, hinf_irka_err_mu, '.-', label=r'$\mathcal{H}_\infty$')
ax.semilogy(mu_fine, hankel_irka_err_mu, '.-', label='Hankel')
ax.set_xlabel(r'$\mu$')
ax.set_title('IRKA errors')
ax.legend()
plt.show()
ax.set_xlabel(r'$\mu$')
ax.set_title('IRKA errors')
ax.legend()
plt.show()
if __name__ == "__main__":
run_demo()
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment