From 2bf919554f0e36849954b89459fcd4a127be7409 Mon Sep 17 00:00:00 2001 From: James Tocknell Date: Fri, 10 Jun 2022 13:07:42 +1000 Subject: [PATCH 1/8] Add compute_jacobian --- src/disc_solver/analyse/compute_jacobian.py | 78 +++++++++++++++++++++ tests/test_run.py | 4 ++ 2 files changed, 82 insertions(+) create mode 100644 src/disc_solver/analyse/compute_jacobian.py diff --git a/src/disc_solver/analyse/compute_jacobian.py b/src/disc_solver/analyse/compute_jacobian.py new file mode 100644 index 00000000..b934e17a --- /dev/null +++ b/src/disc_solver/analyse/compute_jacobian.py @@ -0,0 +1,78 @@ +# -*- coding: utf-8 -*- +""" +""" + +from numpy import zeros + +from ..float_handling import float_type +from ..solve.solution import ode_system +from ..utils import get_solutions +from .utils import ( + analyse_main_wrapper, analysis_func_wrapper, +) + + +def compute_jacobian( + *, γ, a_0, norm_kepler_sq, init_con, θ_scale, η_derivs, use_E_r, θ, params, + eps, +): + rhs_eq, _ = ode_system( + γ=γ, a_0=a_0, norm_kepler_sq=norm_kepler_sq, init_con=init_con, + θ_scale=θ_scale, with_taylor=False, η_derivs=η_derivs, + store_internal=False, use_E_r=use_E_r, v_θ_sonic_crit=None, + after_sonic=None, deriv_v_θ_func=None, check_params=False + ) + + solution_length = params.shape[0] + ode_size = params.shape[1] + + J = zeros([solution_length, ode_size, ode_size], dtype=float_type) + + # compute column for each param + for i in range(ode_size): + derivs_h = zeros([solution_length, ode_size], dtype=float_type) + derivs_l = zeros([solution_length, ode_size], dtype=float_type) + params_h = params.copy() + params_l = params.copy() + + # offset only the value associated with this column + params_h[:, i] += eps + params_l[:, i] -= eps + + # we don't check the validity of inputs as we have these from the + # solution + rhs_eq(θ, params_h, derivs_h) + rhs_eq(θ, params_l, derivs_l) + + J[:, :, i] = (derivs_h - derivs_l) / (2 * eps) + return J + + +def compute_jacobian_from_solution(soln, *, eps, θ_scale=float_type(1)): + solution = soln.solution + angles = soln.angles + cons = soln.initial_conditions + soln_input = soln.solution_input + + init_con = cons.init_con + γ = cons.γ + a_0 = cons.a_0 + norm_kepler_sq = cons.norm_kepler_sq + + η_derivs = soln_input.η_derivs + use_E_r = soln_input.use_E_r + + return compute_jacobian( + γ=γ, a_0=a_0, norm_kepler_sq=norm_kepler_sq, init_con=init_con, + θ_scale=θ_scale, η_derivs=η_derivs, use_E_r=use_E_r, θ=angles, + params=solution, eps=eps, + ) + + +@analysis_func_wrapper +def compute_jacobian_from_file( + soln_file, *, soln_range=None, eps, θ_scale=float_type(1), **kwargs +): + + soln = get_solutions(soln_file, soln_range) + return compute_jacobian_from_solution(soln, eps=eps, θ_scale=θ_scale) diff --git a/tests/test_run.py b/tests/test_run.py index 344ef68f..7da68ef8 100644 --- a/tests/test_run.py +++ b/tests/test_run.py @@ -7,6 +7,7 @@ from disc_solver.analyse.combine_plot import combine_plot from disc_solver.analyse.compare_plot import compare_plot from disc_solver.analyse.component_plot import plot as component_plot +from disc_solver.analyse.compute_jacobian import compute_jacobian_from_file from disc_solver.analyse.conserve_plot import conserve_main from disc_solver.analyse.derivs_plot import derivs_plot from disc_solver.analyse.diverge_plot import diverge_main @@ -519,6 +520,9 @@ def test_vert_plot_file(self, solution, plot_file): def test_stats(self, solution, tmpdir): return stats_main([solution, '--file', str(tmpdir / 'stats.csv')]) + def test_compute_jacobian_from_solution(self, solution): + compute_jacobian_from_file(solution, eps=1e-10) + def test_taylor_space(mpl_interactive): taylor_space_main() From ff3d4d44385f3cf65abd82b9a47d15179a621af0 Mon Sep 17 00:00:00 2001 From: James Tocknell Date: Fri, 24 Jun 2022 12:28:37 +1000 Subject: [PATCH 2/8] WIP --- setup.py | 2 +- tests/test_run.py | 5 +++-- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/setup.py b/setup.py index 84c854f3..e4ff8d2b 100644 --- a/setup.py +++ b/setup.py @@ -12,7 +12,7 @@ packages = setuptools.find_packages('src'), package_dir = {'': 'src'}, install_requires = [ - "numpy", + "numpy<1.23", "matplotlib>=2.2,<3.4a0", "scikits.odes>=2.3.0dev0", "scipy", diff --git a/tests/test_run.py b/tests/test_run.py index 7da68ef8..8427f635 100644 --- a/tests/test_run.py +++ b/tests/test_run.py @@ -15,6 +15,7 @@ from disc_solver.analyse.dump_config import dump_cfg from disc_solver.analyse.hydro_check_plot import hydro_check_plot from disc_solver.analyse.info import info +from disc_solver.analyse.jacobian_plot import jacobian_plot from disc_solver.analyse.j_e_plot import j_e_plot from disc_solver.analyse.params_plot import params_plot from disc_solver.analyse.phys_ratio_plot import plot as phys_ratio_plot @@ -485,10 +486,10 @@ def test_hydro_check_file_no_internal( hydro_check_plot(solution_no_internal, plot_filename=plot_file) def test_jacobian_show(self, solution, mpl_interactive): - plot(solution, show=True) + jacobian_plot(solution, show=True) def test_jacobian_file(self, solution, plot_file): - return plot(solution, plot_filename=plot_file) + return jacobian_plot(solution, plot_filename=plot_file) def test_diverge_show(self, solution, mpl_interactive): diverge_main([solution, '--show']) From 0d15fd4a7f6276455f1b39eddbde81844495f38e Mon Sep 17 00:00:00 2001 From: James Tocknell Date: Fri, 24 Jun 2022 15:53:51 +1000 Subject: [PATCH 3/8] WIP --- setup.py | 2 + src/disc_solver/analyse/compute_jacobian.py | 125 +++++++++++++++++++- tests/test_run.py | 12 +- 3 files changed, 137 insertions(+), 2 deletions(-) diff --git a/setup.py b/setup.py index e4ff8d2b..515c2564 100644 --- a/setup.py +++ b/setup.py @@ -57,6 +57,8 @@ "ds-params-plot = disc_solver.analyse.params_plot:params_main", "ds-taylor-plot = disc_solver.analyse.taylor_plot:taylor_main", "ds-combine-plot = disc_solver.analyse.combine_plot:combine_main", + "ds-jacobian-plot = disc_solver.analyse.jacobian_plot:jacobian_main", + "ds-jacobian-eigenvalue-plot = disc_solver.analyse.compute_jacobian:jacobian_eigenvalues_main", "ds-validate-plot = disc_solver.analyse.validate_plot:validate_plot_main", "ds-hydro-check-plot = disc_solver.analyse.hydro_check_plot:hydro_check_plot_main", "ds-acc-plot = disc_solver.analyse.acc_plot:acc_main", diff --git a/src/disc_solver/analyse/compute_jacobian.py b/src/disc_solver/analyse/compute_jacobian.py index b934e17a..a19ed0cb 100644 --- a/src/disc_solver/analyse/compute_jacobian.py +++ b/src/disc_solver/analyse/compute_jacobian.py @@ -1,17 +1,61 @@ # -*- coding: utf-8 -*- """ """ +from functools import wraps +import numpy as np +from numpy import degrees from numpy import zeros +from scipy.linalg import eigvals from ..float_handling import float_type from ..solve.solution import ode_system from ..utils import get_solutions +from ..utils import ODEIndex from .utils import ( - analyse_main_wrapper, analysis_func_wrapper, + single_solution_plotter, common_plotting_options, analyse_main_wrapper, + get_common_plot_args, analysis_func_wrapper, plot_output_wrapper, + DEFAULT_MPL_STYLE, ) +def plot_parser(parser): + """ + Add arguments for plot command to parser + """ + common_plotting_options(parser) + parser.add_argument( + "--eps", action='store', default=1e-10, type=float + ) + return parser + + +def get_plot_args(args): + """ + Parse plot args + """ + return { + "eps": float(args.get("eps", 1e-10)), + } + + +@analyse_main_wrapper( + "Plot eigenvalues of jacobians for DiscSolver", + plot_parser, + cmd_parser_splitters={ + "common_plot_args": get_common_plot_args, + "plot_args": get_plot_args, + } +) +def jacobian_eigenvalues_main(soln, *, soln_range, common_plot_args, plot_args): + """ + Entry point for ds-derivs-plot + """ + return jacobian_eigenvalues_plot( + soln, soln_range=soln_range, **common_plot_args, **plot_args + ) + + def compute_jacobian( *, γ, a_0, norm_kepler_sq, init_con, θ_scale, η_derivs, use_E_r, θ, params, eps, @@ -76,3 +120,82 @@ def compute_jacobian_from_file( soln = get_solutions(soln_file, soln_range) return compute_jacobian_from_solution(soln, eps=eps, θ_scale=θ_scale) + + +def jacobian_single_solution_plot_wrapper(func): + @wraps(func) + def jacobian_wrapper( + fig, soln, *args, eps, θ_scale=float_type(1), **kwargs + ): + jacobians = compute_jacobian_from_solution( + soln, eps=eps, θ_scale=θ_scale + ) + return func( + fig, *args, jacobians=jacobians, angles=soln.angles, **kwargs + ) + + return jacobian_wrapper + + +@analysis_func_wrapper +def jacobian_eigenvalues_plot( + soln, *, soln_range=None, plot_filename=None, show=False, start=0, stop=90, + figargs=None, linestyle='.', title=None, close=True, filename, + mpl_style=DEFAULT_MPL_STYLE, with_version=True, eps +): + """ + Show derivatives + """ + # pylint: disable=too-many-function-args,unexpected-keyword-arg + fig = plot_jacobian_eigenvalues( + soln, soln_range, linestyle=linestyle, start=start, stop=stop, + figargs=figargs, title=title, filename=filename, + mpl_style=mpl_style, with_version=with_version, eps=eps + ) + + return plot_output_wrapper( + fig, file=plot_filename, show=show, close=close + ) + + +def compute_eigenvalues(jacobian): + if not np.isfinite(jacobian).all(): + return np.full(jacobian.shape[0], np.nan) + return eigvals(jacobian) + + +@single_solution_plotter +@jacobian_single_solution_plot_wrapper +def plot_jacobian_eigenvalues( + fig, *, jacobians, angles, use_E_r, linestyle='.', figargs=None, start=0, + stop=90 +): + """ + Generate plot of jacobians + """ + if figargs is None: + figargs = {} + + npnot = np.logical_not + npand = np.logical_and + indexes = degrees(angles) <= stop + + data = np.array([compute_eigenvalues(j) for j in jacobians]).real + + axes = fig.subplots(nrows=2, ncols=6, sharex=True, **figargs) + axes.shape = 12 + for param in ODEIndex: + ax = axes[param] + pos_data = data[:, param] >= 0 + ax.plot( + degrees(angles[npand(pos_data, indexes)]), + data[npand(pos_data, indexes), param], linestyle + "b", + ) + ax.plot( + degrees(angles[npand(npnot(pos_data), indexes)]), + - data[npand(npnot(pos_data), indexes), param], linestyle + "g", + ) + ax.set_xlabel("angle from plane (°)") + ax.set_yscale("log") + ax.set_title(param.name) + return fig diff --git a/tests/test_run.py b/tests/test_run.py index 8427f635..417ae1d6 100644 --- a/tests/test_run.py +++ b/tests/test_run.py @@ -7,7 +7,9 @@ from disc_solver.analyse.combine_plot import combine_plot from disc_solver.analyse.compare_plot import compare_plot from disc_solver.analyse.component_plot import plot as component_plot -from disc_solver.analyse.compute_jacobian import compute_jacobian_from_file +from disc_solver.analyse.compute_jacobian import ( + compute_jacobian_from_file, jacobian_eigenvalues_plot +) from disc_solver.analyse.conserve_plot import conserve_main from disc_solver.analyse.derivs_plot import derivs_plot from disc_solver.analyse.diverge_plot import diverge_main @@ -524,6 +526,14 @@ def test_stats(self, solution, tmpdir): def test_compute_jacobian_from_solution(self, solution): compute_jacobian_from_file(solution, eps=1e-10) + def test_jacobian_eigenvalues_show(self, solution, mpl_interactive): + jacobian_eigenvalues_plot(solution, show=True, eps=1e-10) + + def test_jacobian_eigenvalues_file(self, solution, plot_file): + return jacobian_eigenvalues_plot( + solution, plot_filename=plot_file, eps=1e-10 + ) + def test_taylor_space(mpl_interactive): taylor_space_main() From 6fe27e1f9a5772dea3c9a376ec2a93c123f63e5d Mon Sep 17 00:00:00 2001 From: James Tocknell Date: Fri, 15 Jul 2022 14:21:15 +1000 Subject: [PATCH 4/8] WIP --- src/disc_solver/analyse/compute_jacobian.py | 17 +++++++---------- 1 file changed, 7 insertions(+), 10 deletions(-) diff --git a/src/disc_solver/analyse/compute_jacobian.py b/src/disc_solver/analyse/compute_jacobian.py index a19ed0cb..59f79ac1 100644 --- a/src/disc_solver/analyse/compute_jacobian.py +++ b/src/disc_solver/analyse/compute_jacobian.py @@ -182,20 +182,17 @@ def plot_jacobian_eigenvalues( data = np.array([compute_eigenvalues(j) for j in jacobians]).real - axes = fig.subplots(nrows=2, ncols=6, sharex=True, **figargs) - axes.shape = 12 - for param in ODEIndex: - ax = axes[param] - pos_data = data[:, param] >= 0 + ax = fig.subplots(**figargs) + ax.set_xlabel("angle from plane (°)") + ax.set_yscale("log") + for i in range(data.shape[1]): + pos_data = data[:, i] >= 0 ax.plot( degrees(angles[npand(pos_data, indexes)]), - data[npand(pos_data, indexes), param], linestyle + "b", + data[npand(pos_data, indexes), i], linestyle + "b", ) ax.plot( degrees(angles[npand(npnot(pos_data), indexes)]), - - data[npand(npnot(pos_data), indexes), param], linestyle + "g", + - data[npand(npnot(pos_data), indexes), i], linestyle + "g", ) - ax.set_xlabel("angle from plane (°)") - ax.set_yscale("log") - ax.set_title(param.name) return fig From f30625ee74017ba703641f2049b4b11f3d38c10c Mon Sep 17 00:00:00 2001 From: James Tocknell Date: Sat, 23 Jul 2022 00:23:47 +1000 Subject: [PATCH 5/8] wip --- src/disc_solver/analyse/compute_jacobian.py | 64 +++++++++++++++------ 1 file changed, 46 insertions(+), 18 deletions(-) diff --git a/src/disc_solver/analyse/compute_jacobian.py b/src/disc_solver/analyse/compute_jacobian.py index 59f79ac1..ddee6bfb 100644 --- a/src/disc_solver/analyse/compute_jacobian.py +++ b/src/disc_solver/analyse/compute_jacobian.py @@ -3,10 +3,12 @@ """ from functools import wraps -import numpy as np -from numpy import degrees -from numpy import zeros +from numpy import ( + degrees, zeros, logical_not as npnot, logical_and as npand, isfinite, nan, + full, array as nparray +) from scipy.linalg import eigvals +import matplotlib.pyplot as plt from ..float_handling import float_type from ..solve.solution import ode_system @@ -159,11 +161,40 @@ def jacobian_eigenvalues_plot( def compute_eigenvalues(jacobian): - if not np.isfinite(jacobian).all(): - return np.full(jacobian.shape[0], np.nan) + if not isfinite(jacobian).all(): + return full(jacobian.shape[0], nan) return eigvals(jacobian) +def plot_log_complex_data(ax, *, x_data, y_data, **kwargs): + pos_real_data_slice = y_data.real >= 0 + pos_imag_data_slice = y_data.imag >= 0 + ax.plot( + x_data[pos_real_data_slice], + y_data[pos_real_data_slice].real, + marker="triangle_up", + **kwargs + ) + ax.plot( + x_data[npnot(pos_real_data_slice)], + y_data[npnot(pos_real_data_slice)].real, + marker="triangle_up", + **kwargs + ) + ax.plot( + x_data[pos_imag_data_slice], + y_data[pos_imag_data_slice].imag, + marker="tri_up", + **kwargs + ) + ax.plot( + x_data[npnot(pos_imag_data_slice)], + y_data[npnot(pos_imag_data_slice)].imag, + marker="tri_up", + **kwargs + ) + + @single_solution_plotter @jacobian_single_solution_plot_wrapper def plot_jacobian_eigenvalues( @@ -176,23 +207,20 @@ def plot_jacobian_eigenvalues( if figargs is None: figargs = {} - npnot = np.logical_not - npand = np.logical_and indexes = degrees(angles) <= stop - data = np.array([compute_eigenvalues(j) for j in jacobians]).real + data = nparray([compute_eigenvalues(j) for j in jacobians]) ax = fig.subplots(**figargs) ax.set_xlabel("angle from plane (°)") ax.set_yscale("log") - for i in range(data.shape[1]): - pos_data = data[:, i] >= 0 - ax.plot( - degrees(angles[npand(pos_data, indexes)]), - data[npand(pos_data, indexes), i], linestyle + "b", - ) - ax.plot( - degrees(angles[npand(npnot(pos_data), indexes)]), - - data[npand(npnot(pos_data), indexes), i], linestyle + "g", - ) + with plt.style.context({ + "axes.prop_cycle": cycler.cycler(color=plt.get_cmap("tab20").colors) + }): + for i in range(data.shape[1]): + plot_log_complex_data( + ax, x_data=degrees(angles[indexes]), + y_data=data[indexes, i], + color="C" + str(i), + ) return fig From 84ca1e78f8189bbde2a39604c6becf49f046762c Mon Sep 17 00:00:00 2001 From: James Tocknell Date: Fri, 12 Aug 2022 18:41:33 +1000 Subject: [PATCH 6/8] wip --- setup.py | 1 + src/disc_solver/analyse/compute_jacobian.py | 140 +++++++++++++++++--- 2 files changed, 123 insertions(+), 18 deletions(-) diff --git a/setup.py b/setup.py index 515c2564..50909ef0 100644 --- a/setup.py +++ b/setup.py @@ -59,6 +59,7 @@ "ds-combine-plot = disc_solver.analyse.combine_plot:combine_main", "ds-jacobian-plot = disc_solver.analyse.jacobian_plot:jacobian_main", "ds-jacobian-eigenvalue-plot = disc_solver.analyse.compute_jacobian:jacobian_eigenvalues_main", + "ds-jacobian-eigenvector-plot = disc_solver.analyse.compute_jacobian:jacobian_eigenvectors_main", "ds-validate-plot = disc_solver.analyse.validate_plot:validate_plot_main", "ds-hydro-check-plot = disc_solver.analyse.hydro_check_plot:hydro_check_plot_main", "ds-acc-plot = disc_solver.analyse.acc_plot:acc_main", diff --git a/src/disc_solver/analyse/compute_jacobian.py b/src/disc_solver/analyse/compute_jacobian.py index ddee6bfb..a90fb36a 100644 --- a/src/disc_solver/analyse/compute_jacobian.py +++ b/src/disc_solver/analyse/compute_jacobian.py @@ -3,23 +3,27 @@ """ from functools import wraps +import cycler from numpy import ( - degrees, zeros, logical_not as npnot, logical_and as npand, isfinite, nan, - full, array as nparray + degrees, zeros, logical_not as npnot, isfinite, nan, full, + array as nparray, ) -from scipy.linalg import eigvals +from scipy.linalg import eigvals, eig +import matplotlib as mpl import matplotlib.pyplot as plt from ..float_handling import float_type from ..solve.solution import ode_system from ..utils import get_solutions -from ..utils import ODEIndex from .utils import ( single_solution_plotter, common_plotting_options, analyse_main_wrapper, get_common_plot_args, analysis_func_wrapper, plot_output_wrapper, DEFAULT_MPL_STYLE, ) +mpl.rcParams['agg.path.chunksize'] = 1000000 +COLOURS = plt.get_cmap("tab20").colors + def plot_parser(parser): """ @@ -49,7 +53,9 @@ def get_plot_args(args): "plot_args": get_plot_args, } ) -def jacobian_eigenvalues_main(soln, *, soln_range, common_plot_args, plot_args): +def jacobian_eigenvalues_main( + soln, *, soln_range, common_plot_args, plot_args +): """ Entry point for ds-derivs-plot """ @@ -58,6 +64,25 @@ def jacobian_eigenvalues_main(soln, *, soln_range, common_plot_args, plot_args): ) +@analyse_main_wrapper( + "Plot eigenvectors of jacobians for DiscSolver", + plot_parser, + cmd_parser_splitters={ + "common_plot_args": get_common_plot_args, + "plot_args": get_plot_args, + } +) +def jacobian_eigenvectors_main( + soln, *, soln_range, common_plot_args, plot_args +): + """ + Entry point for ds-derivs-plot + """ + return jacobian_eigenvectors_plot( + soln, soln_range=soln_range, **common_plot_args, **plot_args + ) + + def compute_jacobian( *, γ, a_0, norm_kepler_sq, init_con, θ_scale, η_derivs, use_E_r, θ, params, eps, @@ -160,37 +185,69 @@ def jacobian_eigenvalues_plot( ) +@analysis_func_wrapper +def jacobian_eigenvectors_plot( + soln, *, soln_range=None, plot_filename=None, show=False, start=0, stop=90, + figargs=None, linestyle='.', title=None, close=True, filename, + mpl_style=DEFAULT_MPL_STYLE, with_version=True, eps +): + """ + Show derivatives + """ + # pylint: disable=too-many-function-args,unexpected-keyword-arg + fig = plot_jacobian_eigenvectors( + soln, soln_range, linestyle=linestyle, start=start, stop=stop, + figargs=figargs, title=title, filename=filename, + mpl_style=mpl_style, with_version=with_version, eps=eps + ) + + return plot_output_wrapper( + fig, file=plot_filename, show=show, close=close + ) + + def compute_eigenvalues(jacobian): if not isfinite(jacobian).all(): return full(jacobian.shape[0], nan) return eigvals(jacobian) +def compute_eigenvalues_and_eigenvectors(jacobian): + if not isfinite(jacobian).all(): + M = jacobian.shape[0] + return full(M, nan), full([M, M], nan) + return eig(jacobian) + + def plot_log_complex_data(ax, *, x_data, y_data, **kwargs): pos_real_data_slice = y_data.real >= 0 pos_imag_data_slice = y_data.imag >= 0 ax.plot( x_data[pos_real_data_slice], y_data[pos_real_data_slice].real, - marker="triangle_up", + marker="^", + linestyle='None', **kwargs ) ax.plot( x_data[npnot(pos_real_data_slice)], - y_data[npnot(pos_real_data_slice)].real, - marker="triangle_up", + - y_data[npnot(pos_real_data_slice)].real, + marker="v", + linestyle='None', **kwargs ) ax.plot( x_data[pos_imag_data_slice], y_data[pos_imag_data_slice].imag, - marker="tri_up", + marker="2", + linestyle='None', **kwargs ) ax.plot( x_data[npnot(pos_imag_data_slice)], - y_data[npnot(pos_imag_data_slice)].imag, - marker="tri_up", + - y_data[npnot(pos_imag_data_slice)].imag, + marker="1", + linestyle='None', **kwargs ) @@ -207,20 +264,67 @@ def plot_jacobian_eigenvalues( if figargs is None: figargs = {} - indexes = degrees(angles) <= stop + indexes = (start <= degrees(angles)) & (degrees(angles) <= stop) data = nparray([compute_eigenvalues(j) for j in jacobians]) ax = fig.subplots(**figargs) ax.set_xlabel("angle from plane (°)") ax.set_yscale("log") - with plt.style.context({ - "axes.prop_cycle": cycler.cycler(color=plt.get_cmap("tab20").colors) - }): - for i in range(data.shape[1]): + for i in range(data.shape[1]): + plot_log_complex_data( + ax, x_data=degrees(angles[indexes]), + y_data=data[indexes, i], + color=COLOURS[i], + alpha=.5, + label=str(i) + ) + ax.legend() + return fig + + +@single_solution_plotter +@jacobian_single_solution_plot_wrapper +def plot_jacobian_eigenvectors( + fig, *, jacobians, angles, use_E_r, linestyle='.', start=0, + stop=90 +): + """ + Generate plot of jacobians + """ + indexes = (start <= degrees(angles)) & (degrees(angles) <= stop) + + eigvalues, eigvecs = zip(*[ + compute_eigenvalues_and_eigenvectors(j) for j in jacobians + ]) + eigvalues = nparray(eigvalues) + eigvecs = nparray(eigvecs) + + axes = fig.subplots( + nrows=2, ncols=6, sharex=True, gridspec_kw=dict(hspace=0), + ) + + # only add label to bottom plots + for ax in axes[1]: + ax.set_xlabel("angle from plane (°)") + + axes.shape = 12 + for j in range(eigvalues.shape[1]): + ax = axes[j] + ax.set_yscale("log") + plot_log_complex_data( + ax, x_data=degrees(angles[indexes]), + y_data=eigvalues[indexes, j], + color="black", + alpha=.5, + ) + for i in range(eigvecs.shape[2]): plot_log_complex_data( ax, x_data=degrees(angles[indexes]), - y_data=data[indexes, i], - color="C" + str(i), + y_data=eigvecs[indexes, i, j], + color=COLOURS[i], + alpha=.5, + label=str(i) ) + ax.legend() return fig From a7c8107e969ff5bc9f58d1977df41c884b5f90d9 Mon Sep 17 00:00:00 2001 From: James Tocknell Date: Sat, 13 Aug 2022 22:18:08 +1000 Subject: [PATCH 7/8] wip --- src/disc_solver/analyse/compute_jacobian.py | 52 ++++++++++++--------- tests/conftest.py | 12 +++++ tests/test_run.py | 24 +++------- 3 files changed, 49 insertions(+), 39 deletions(-) diff --git a/src/disc_solver/analyse/compute_jacobian.py b/src/disc_solver/analyse/compute_jacobian.py index a90fb36a..b9b270db 100644 --- a/src/disc_solver/analyse/compute_jacobian.py +++ b/src/disc_solver/analyse/compute_jacobian.py @@ -1,9 +1,9 @@ # -*- coding: utf-8 -*- """ +Analysis functions related to the jacobians of existing solutions """ from functools import wraps -import cycler from numpy import ( degrees, zeros, logical_not as npnot, isfinite, nan, full, array as nparray, @@ -14,7 +14,6 @@ from ..float_handling import float_type from ..solve.solution import ode_system -from ..utils import get_solutions from .utils import ( single_solution_plotter, common_plotting_options, analyse_main_wrapper, get_common_plot_args, analysis_func_wrapper, plot_output_wrapper, @@ -87,11 +86,14 @@ def compute_jacobian( *, γ, a_0, norm_kepler_sq, init_con, θ_scale, η_derivs, use_E_r, θ, params, eps, ): + """ + Compute jacobian from solution + """ rhs_eq, _ = ode_system( γ=γ, a_0=a_0, norm_kepler_sq=norm_kepler_sq, init_con=init_con, θ_scale=θ_scale, with_taylor=False, η_derivs=η_derivs, store_internal=False, use_E_r=use_E_r, v_θ_sonic_crit=None, - after_sonic=None, deriv_v_θ_func=None, check_params=False + after_sonic=None, deriv_v_θ_func=None, derivs_post_solution=True, ) solution_length = params.shape[0] @@ -120,6 +122,9 @@ def compute_jacobian( def compute_jacobian_from_solution(soln, *, eps, θ_scale=float_type(1)): + """ + Compute jacobian from solution + """ solution = soln.solution angles = soln.angles cons = soln.initial_conditions @@ -140,20 +145,15 @@ def compute_jacobian_from_solution(soln, *, eps, θ_scale=float_type(1)): ) -@analysis_func_wrapper -def compute_jacobian_from_file( - soln_file, *, soln_range=None, eps, θ_scale=float_type(1), **kwargs -): - - soln = get_solutions(soln_file, soln_range) - return compute_jacobian_from_solution(soln, eps=eps, θ_scale=θ_scale) - - def jacobian_single_solution_plot_wrapper(func): + """ + Helper wrapper for solutions working with jacobians + """ @wraps(func) def jacobian_wrapper( - fig, soln, *args, eps, θ_scale=float_type(1), **kwargs + fig, soln, *args, eps, use_E_r, θ_scale=float_type(1), **kwargs ): + # pylint: disable=unused-argument jacobians = compute_jacobian_from_solution( soln, eps=eps, θ_scale=θ_scale ) @@ -171,9 +171,10 @@ def jacobian_eigenvalues_plot( mpl_style=DEFAULT_MPL_STYLE, with_version=True, eps ): """ - Show derivatives + Plot eigenvalues of jacobians of a solution """ # pylint: disable=too-many-function-args,unexpected-keyword-arg + # pylint: disable=missing-kwoa fig = plot_jacobian_eigenvalues( soln, soln_range, linestyle=linestyle, start=start, stop=stop, figargs=figargs, title=title, filename=filename, @@ -195,6 +196,7 @@ def jacobian_eigenvectors_plot( Show derivatives """ # pylint: disable=too-many-function-args,unexpected-keyword-arg + # pylint: disable=missing-kwoa fig = plot_jacobian_eigenvectors( soln, soln_range, linestyle=linestyle, start=start, stop=stop, figargs=figargs, title=title, filename=filename, @@ -207,12 +209,18 @@ def jacobian_eigenvectors_plot( def compute_eigenvalues(jacobian): + """ + Compute eigenvalues of a jacobian matrix + """ if not isfinite(jacobian).all(): return full(jacobian.shape[0], nan) return eigvals(jacobian) def compute_eigenvalues_and_eigenvectors(jacobian): + """ + Compute eigenvalues and eigenvectors of a jacobian matrix + """ if not isfinite(jacobian).all(): M = jacobian.shape[0] return full(M, nan), full([M, M], nan) @@ -220,6 +228,9 @@ def compute_eigenvalues_and_eigenvectors(jacobian): def plot_log_complex_data(ax, *, x_data, y_data, **kwargs): + """ + Plot complex data logarithmically + """ pos_real_data_slice = y_data.real >= 0 pos_imag_data_slice = y_data.imag >= 0 ax.plot( @@ -255,20 +266,17 @@ def plot_log_complex_data(ax, *, x_data, y_data, **kwargs): @single_solution_plotter @jacobian_single_solution_plot_wrapper def plot_jacobian_eigenvalues( - fig, *, jacobians, angles, use_E_r, linestyle='.', figargs=None, start=0, - stop=90 + fig, *, jacobians, angles, start=0, stop=90, linestyle=',', ): """ Generate plot of jacobians """ - if figargs is None: - figargs = {} - + # pylint: disable=unused-argument indexes = (start <= degrees(angles)) & (degrees(angles) <= stop) data = nparray([compute_eigenvalues(j) for j in jacobians]) - ax = fig.subplots(**figargs) + ax = fig.subplots() ax.set_xlabel("angle from plane (°)") ax.set_yscale("log") for i in range(data.shape[1]): @@ -286,12 +294,12 @@ def plot_jacobian_eigenvalues( @single_solution_plotter @jacobian_single_solution_plot_wrapper def plot_jacobian_eigenvectors( - fig, *, jacobians, angles, use_E_r, linestyle='.', start=0, - stop=90 + fig, *, jacobians, angles, start=0, stop=90, linestyle=',', ): """ Generate plot of jacobians """ + # pylint: disable=unused-argument indexes = (start <= degrees(angles)) & (degrees(angles) <= stop) eigvalues, eigvecs = zip(*[ diff --git a/tests/conftest.py b/tests/conftest.py index ab3b2484..11938c8b 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -49,6 +49,13 @@ USE_E_R_SOLUTIONS = [ "mod_hydro_solution_use_E_r", ] +NON_APPROX_SOLUTIONS = [ + "single_solution_default", + "sonic_root_solution_default", + "step_solution_default", + "single_solution_no_internal", + "sonic_root_solution_no_internal", +] @pytest.fixture(scope="session") @@ -213,6 +220,11 @@ def solution_use_E_r(request): return request.getfixturevalue(request.param) +@pytest.fixture(scope="session", params=NON_APPROX_SOLUTIONS) +def solution_not_approx(request): + return request.getfixturevalue(request.param) + + @pytest.fixture() def mpl_interactive(request): import matplotlib as mpl diff --git a/tests/test_run.py b/tests/test_run.py index 417ae1d6..c1d31bdf 100644 --- a/tests/test_run.py +++ b/tests/test_run.py @@ -7,9 +7,7 @@ from disc_solver.analyse.combine_plot import combine_plot from disc_solver.analyse.compare_plot import compare_plot from disc_solver.analyse.component_plot import plot as component_plot -from disc_solver.analyse.compute_jacobian import ( - compute_jacobian_from_file, jacobian_eigenvalues_plot -) +from disc_solver.analyse.compute_jacobian import jacobian_eigenvalues_plot from disc_solver.analyse.conserve_plot import conserve_main from disc_solver.analyse.derivs_plot import derivs_plot from disc_solver.analyse.diverge_plot import diverge_main @@ -17,7 +15,6 @@ from disc_solver.analyse.dump_config import dump_cfg from disc_solver.analyse.hydro_check_plot import hydro_check_plot from disc_solver.analyse.info import info -from disc_solver.analyse.jacobian_plot import jacobian_plot from disc_solver.analyse.j_e_plot import j_e_plot from disc_solver.analyse.params_plot import params_plot from disc_solver.analyse.phys_ratio_plot import plot as phys_ratio_plot @@ -487,12 +484,6 @@ def test_hydro_check_file_no_internal( with pytest.raises(AnalysisError): hydro_check_plot(solution_no_internal, plot_filename=plot_file) - def test_jacobian_show(self, solution, mpl_interactive): - jacobian_plot(solution, show=True) - - def test_jacobian_file(self, solution, plot_file): - return jacobian_plot(solution, plot_filename=plot_file) - def test_diverge_show(self, solution, mpl_interactive): diverge_main([solution, '--show']) @@ -523,15 +514,14 @@ def test_vert_plot_file(self, solution, plot_file): def test_stats(self, solution, tmpdir): return stats_main([solution, '--file', str(tmpdir / 'stats.csv')]) - def test_compute_jacobian_from_solution(self, solution): - compute_jacobian_from_file(solution, eps=1e-10) - - def test_jacobian_eigenvalues_show(self, solution, mpl_interactive): - jacobian_eigenvalues_plot(solution, show=True, eps=1e-10) + def test_jacobian_eigenvalues_show( + self, solution_not_approx, mpl_interactive + ): + jacobian_eigenvalues_plot(solution_not_approx, show=True, eps=1e-10) - def test_jacobian_eigenvalues_file(self, solution, plot_file): + def test_jacobian_eigenvalues_file(self, solution_not_approx, plot_file): return jacobian_eigenvalues_plot( - solution, plot_filename=plot_file, eps=1e-10 + solution_not_approx, plot_filename=plot_file, eps=1e-10 ) From d86712fb40493c9985eed43239a346edfbc78d6b Mon Sep 17 00:00:00 2001 From: James Tocknell Date: Sun, 14 Aug 2022 11:14:55 +1000 Subject: [PATCH 8/8] wip --- tests/test_run.py | 14 +++++++++++++- 1 file changed, 13 insertions(+), 1 deletion(-) diff --git a/tests/test_run.py b/tests/test_run.py index c1d31bdf..39e5d39e 100644 --- a/tests/test_run.py +++ b/tests/test_run.py @@ -7,7 +7,9 @@ from disc_solver.analyse.combine_plot import combine_plot from disc_solver.analyse.compare_plot import compare_plot from disc_solver.analyse.component_plot import plot as component_plot -from disc_solver.analyse.compute_jacobian import jacobian_eigenvalues_plot +from disc_solver.analyse.compute_jacobian import ( + jacobian_eigenvalues_plot, jacobian_eigenvectors_plot, +) from disc_solver.analyse.conserve_plot import conserve_main from disc_solver.analyse.derivs_plot import derivs_plot from disc_solver.analyse.diverge_plot import diverge_main @@ -524,6 +526,16 @@ def test_jacobian_eigenvalues_file(self, solution_not_approx, plot_file): solution_not_approx, plot_filename=plot_file, eps=1e-10 ) + def test_jacobian_eigenvectors_show( + self, solution_not_approx, mpl_interactive + ): + jacobian_eigenvectors_plot(solution_not_approx, show=True, eps=1e-10) + + def test_jacobian_eigenvectors_file(self, solution_not_approx, plot_file): + return jacobian_eigenvectors_plot( + solution_not_approx, plot_filename=plot_file, eps=1e-10 + ) + def test_taylor_space(mpl_interactive): taylor_space_main()