diff --git a/.github/workflows/docker-image.yml b/.github/workflows/docker-image.yml index f754f6f..904209d 100644 --- a/.github/workflows/docker-image.yml +++ b/.github/workflows/docker-image.yml @@ -15,8 +15,10 @@ jobs: registry: ghcr.io username: ${{github.actor}} password: ${{secrets.GITHUB_TOKEN}} - - name: Build and push - run: | - docker build . --tag ghcr.io/LemurPwned/cmtj:latest - docker push ghcr.io/LemurPwned/cmtj:latest - + - name: Build and push image + uses: docker/build-push-action@v2 + with: + push: true + file: docker/Dockerfile + tags: | + ghcr.io/lemurpwned/cmtj:latest \ No newline at end of file diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml index 6b6c037..32ee4a8 100644 --- a/.github/workflows/main.yml +++ b/.github/workflows/main.yml @@ -35,7 +35,7 @@ jobs: # Steps represent a sequence of tasks that will be executed as part of the job steps: # Checks-out your repository under $GITHUB_WORKSPACE, so your job can access it - - uses: actions/checkout@v2 + - uses: actions/checkout@v3 - name: Python wheels manylinux stable build uses: RalfG/python-wheels-manylinux-build@v0.5.0 with: diff --git a/.github/workflows/python-test.yml b/.github/workflows/python-test.yml index f9a776d..82bcaa9 100644 --- a/.github/workflows/python-test.yml +++ b/.github/workflows/python-test.yml @@ -8,10 +8,10 @@ on: branches: [ "master" ] pull_request: branches: [ "master" ] + workflow_dispatch: jobs: - build: - + build-linux: runs-on: ubuntu-latest strategy: fail-fast: true @@ -33,3 +33,26 @@ jobs: - name: Test with pytest run: | pytest + + build-windows: + runs-on: windows-latest + strategy: + fail-fast: true + matrix: + python-version: ["3.9", "3.10", "3.11"] + + steps: + - uses: actions/checkout@v3 + - name: Set up Python ${{ matrix.python-version }} + uses: actions/setup-python@v3 + with: + python-version: ${{ matrix.python-version }} + - name: Install dependencies + run: | + python -m pip install --upgrade pip + python -m pip install flake8 pytest + python -m pip install -e .[utils,test] + + - name: Test with pytest + run: | + pytest \ No newline at end of file diff --git a/CHANGELOG.md b/CHANGELOG.md index 59fcaed..52ce9cd 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,16 @@ # Changelog +# 1.4.0 + +- Adding new, dynamic symbolic model compatible with `Solver` class. It is now possible to use the `Solver` class with `LayerDynamic` to solve the LLG equation. +- Added tests for procedures and operators. +- Added missing operators for `CVector` class in Python. +- `CVector` is now subscriptable. +- Added new `CVector` and `AxialDriver` initialisations. +- `VSD` and `PIMM` procedures accept additional new parameters. +- Added some optimization utilities like `coordinate_descent`. +- Added a `streamlit` service for an example PIMM simulation. + # 1.3.2 - Added new `ScalarDrivers` -- Gaussian impulse and Gaussian step. diff --git a/README.md b/README.md index 8703e1d..594b751 100644 --- a/README.md +++ b/README.md @@ -6,6 +6,7 @@ [![pages-build-deployment](https://github.com/LemurPwned/cmtj/actions/workflows/pages/pages-build-deployment/badge.svg?branch=gh-pages)](https://github.com/LemurPwned/cmtj/actions/workflows/pages/pages-build-deployment) [![Version](https://img.shields.io/pypi/v/cmtj)](https://pypi.org/project/cmtj/) [![License](https://img.shields.io/pypi/l/cmtj.svg)](https://github.com/LemurPwned/cmtj/blob/master/LICENSE) +[![Streamlit](https://static.streamlit.io/badges/streamlit_badge_black_white.svg)](http://cmtj-simulations.streamlit.app/) ![Downloads](https://img.shields.io/pypi/dm/cmtj.svg) ## Short description @@ -14,6 +15,10 @@ A name may be misleading -- the MTJ (Magnetic Tunnel Junctions) are not the only The library allows for macromagnetic simulation of various multilayer spintronic structures. The package uses C++ implementation of (s)LLGS (stochastic Landau-Lifschitz-Gilbert-Slonczewski) equation with various field contributions included for instance: anisotropy, interlayer exchange coupling, demagnetisation, dipole fields etc. It is also possible to connect devices in parallel or in series to have electrically coupled arrays. +## Demo + +Check out the [streamlit hosted demo here](http://cmtj-simulations.streamlit.app/). + ## Quickstart #### Installation :rocket: @@ -97,6 +102,11 @@ Many thanks to professor Jack Sankey for his help with the development of therma All contributions are welcome, please leave an issue if you've encountered any trouble with setup or running the library. +## Docker + +In the `docker` directory there's a `Dockerfile` that can be used to build a docker image with the library installed. +`Dockerfile.app` is used for streamlit development. + ## Precommit There's a `.pre-commit-config.yaml` that does some basic python and cpp lints and checks. More static analysis to come in the future. diff --git a/cmtj/models/__init__.py b/cmtj/models/__init__.py index e69de29..70d0afb 100644 --- a/cmtj/models/__init__.py +++ b/cmtj/models/__init__.py @@ -0,0 +1,8 @@ +from .general_sb import LayerDynamic, LayerSB, Solver, VectorObj + +__all__ = [ + "LayerSB", + "LayerDynamic", + "Solver", + "VectorObj", +] diff --git a/cmtj/models/ensemble.py b/cmtj/models/ensemble.py index 3bc3412..963c7d6 100644 --- a/cmtj/models/ensemble.py +++ b/cmtj/models/ensemble.py @@ -47,5 +47,5 @@ def mixed_lorentz(H, dH, Hr, Va, Vas): :param Vs: amplitude of symmetric Lorentzian :param Vas: amplitude of antisymmetric Lorentzian """ - return symmetric_lorentz(H, dH, Hr, - Va) + Vb * antisymmetric_lorentz(H, dH, Hr, Vas) + return symmetric_lorentz(H, dH, Hr, Va) + antisymmetric_lorentz( + H, dH, Hr, Vas) diff --git a/cmtj/models/general_sb.py b/cmtj/models/general_sb.py index c4b2839..c7c28ae 100644 --- a/cmtj/models/general_sb.py +++ b/cmtj/models/general_sb.py @@ -13,6 +13,8 @@ from ..utils import VectorObj, gamma, gamma_rad, mu0, perturb_position from ..utils.solvers import RootFinder +EPS = np.finfo('float64').resolution + def real_deocrator(fn): """Using numpy real cast is way faster than sympy.""" @@ -24,9 +26,11 @@ def wrap_fn(*args): @njit -def fast_norm(x): +def fast_norm(x): # sourcery skip: for-index-underscore, sum-comprehension """Fast norm function for 1D arrays.""" - sum_ = sum(x_**2 for x_ in x) + sum_ = 0 + for x_ in x: + sum_ += x_**2 return math.sqrt(sum_) @@ -41,8 +45,7 @@ def general_hessian_functional(N: int): # indx_i = str(i + 1) # for display purposes indx_i = str(i) all_symbols.extend( - (sym.Symbol(r"\theta_" + indx_i), sym.Symbol(r"\phi_" + indx_i)) - ) + (sym.Symbol(r"\theta_" + indx_i), sym.Symbol(r"\phi_" + indx_i))) energy_functional_expr = sym.Function("E")(*all_symbols) return get_hessian_from_energy_expr( N, energy_functional_expr), energy_functional_expr @@ -50,6 +53,14 @@ def general_hessian_functional(N: int): @lru_cache def get_hessian_from_energy_expr(N: int, energy_functional_expr: sym.Expr): + """ + Computes the Hessian matrix of the energy functional expression with respect to the spin angles and phases. + + :param N (int): The number of spins. + :param energy_functional_expr (sympy.Expr): The energy functional expression. + + returns: sympy.Matrix: The Hessian matrix of the energy functional expression. + """ hessian = [[0 for _ in range(2 * N)] for _ in range(2 * N)] for i in range(N): # indx_i = str(i + 1) # for display purposes @@ -155,7 +166,7 @@ def __post_init__(self): raise ValueError("Only up to 10 layers supported.") self.theta = sym.Symbol(r"\theta_" + str(self._id)) self.phi = sym.Symbol(r"\phi_" + str(self._id)) - self.m = sym.Matrix([ + self.m = sym.ImmutableMatrix([ sym.sin(self.theta) * sym.cos(self.phi), sym.sin(self.theta) * sym.sin(self.phi), sym.cos(self.theta) @@ -169,7 +180,8 @@ def get_m_sym(self): """Returns the magnetisation vector.""" return self.m - def symbolic_layer_energy(self, H: sym.Matrix, J1top: float, + @lru_cache(3) + def symbolic_layer_energy(self, H: sym.ImmutableMatrix, J1top: float, J1bottom: float, J2top: float, J2bottom: float, top_layer: "LayerSB", down_layer: "LayerSB"): """Returns the symbolic expression for the energy of the layer. @@ -191,12 +203,14 @@ def symbolic_layer_energy(self, H: sym.Matrix, J1top: float, other_m) - (J2bottom / self.thickness) * m.dot(other_m)**2 return eng_non_interaction + top_iec_energy + bottom_iec_energy - def no_iec_symbolic_layer_energy(self, H: sym.Matrix): + def no_iec_symbolic_layer_energy(self, H: sym.ImmutableMatrix): """Returns the symbolic expression for the energy of the layer. Coupling contribution comes only from the bottom layer (top-down crawl)""" m = self.get_m_sym() - alpha = sym.Matrix([sym.cos(self.Kv.phi), sym.sin(self.Kv.phi), 0]) + alpha = sym.ImmutableMatrix( + [sym.cos(self.Kv.phi), + sym.sin(self.Kv.phi), 0]) field_energy = -mu0 * self.Ms * m.dot(H) surface_anistropy = (-self.Ks + @@ -208,10 +222,49 @@ def sb_correction(self): omega = sym.Symbol(r'\omega') return (omega / gamma) * self.Ms * sym.sin(self.theta) * self.thickness + def __hash__(self) -> int: + return hash(str(self)) + + def __eq__(self, __value: "LayerSB") -> bool: + return self._id == __value._id and self.thickness == __value.thickness and self.Kv == __value.Kv and self.Ks == __value.Ks and self.Ms == __value.Ms + @dataclass -class SolverSB: - layers: List[LayerSB] +class LayerDynamic(LayerSB): + alpha: float + + def rhs_llg(self, H: sym.Matrix, J1top: float, J1bottom: float, + J2top: float, J2bottom: float, top_layer: "LayerSB", + down_layer: "LayerSB"): + """Returns the symbolic expression for the RHS of the spherical LLG equation. + Coupling contribution comes only from the bottom layer (top-down crawl)""" + U = self.symbolic_layer_energy(H, + J1top=J1top, + J1bottom=J1bottom, + J2top=J2top, + J2bottom=J2bottom, + top_layer=top_layer, + down_layer=down_layer) + # sum all components + prefac = gamma_rad / (1. + self.alpha)**2 + inv_sin = 1. / (sym.sin(self.theta) + EPS) + dUdtheta = sym.diff(U, self.theta) + dUdphi = sym.diff(U, self.phi) + + dtheta = (-inv_sin * dUdphi - self.alpha * dUdtheta) + dphi = (inv_sin * dUdtheta - self.alpha * dUdphi * (inv_sin)**2) + return prefac * sym.ImmutableMatrix([dtheta, dphi]) / self.Ms + + def __eq__(self, __value: "LayerDynamic") -> bool: + return super().__eq__(__value) and self.alpha == __value.alpha + + def __hash__(self) -> int: + return super().__hash__() + + +@dataclass +class Solver: + layers: List[Union[LayerSB, LayerDynamic]] J1: List[float] J2: List[float] H: VectorObj = None @@ -225,7 +278,8 @@ def __post_init__(self): id_sets = {layer._id for layer in self.layers} ideal_set = set(range(len(self.layers))) if id_sets != ideal_set: - raise ValueError("Layer ids must be 0, 1, 2, ... and unique") + raise ValueError("Layer ids must be 0, 1, 2, ... and unique." + "Ids must start from 0.") def get_layer_references(self, layer_indx, interaction_constant): """Returns the references to the layers above and below the layer @@ -243,8 +297,26 @@ def get_layer_references(self, layer_indx, interaction_constant): 1], interaction_constant[layer_indx - 1], interaction_constant[layer_indx] + def compose_llg_jacobian(self, H: VectorObj): + """Create a symbolic jacobian of the LLG equation in spherical coordinates.""" + # has order theta0, phi0, theta1, phi1, ... + if isinstance(H, VectorObj): + H = sym.ImmutableMatrix(H.get_cartesian()) + + symbols, fns = [], [] + for i, layer in enumerate(self.layers): + symbols.extend((layer.theta, layer.phi)) + top_layer, bottom_layer, Jtop, Jbottom = self.get_layer_references( + i, self.J1) + _, _, J2top, J2bottom = self.get_layer_references(i, self.J2) + fns.append( + layer.rhs_llg(H, Jtop, Jbottom, J2top, J2bottom, top_layer, + bottom_layer)) + jac = sym.ImmutableMatrix(fns).jacobian(symbols) + return jac, symbols + def create_energy(self, - H: Union[VectorObj, sym.Matrix] = None, + H: Union[VectorObj, sym.ImmutableMatrix] = None, volumetric: bool = False): """Creates the symbolic energy expression. @@ -257,7 +329,7 @@ def create_energy(self, """ if H is None: h = self.H.get_cartesian() - H = sym.Matrix(h) + H = sym.ImmutableMatrix(h) energy = 0 if volumetric: # volumetric energy -- DO NOT USE IN GENERAL @@ -327,7 +399,7 @@ def create_energy_hessian(self, equilibrium_position: List[float]): hessian[2 * i + 1][2 * j] = expr hessian[2 * j][2 * i + 1] = expr - hes = sym.Matrix(hessian) + hes = sym.ImmutableMatrix(hessian) _, U, _ = hes.LUdecomposition() return U.det().subs(subs) @@ -338,7 +410,8 @@ def get_gradient_expr(self, accel="math"): symbols = [] for layer in self.layers: (theta, phi) = layer.get_coord_sym() - grad_vector.extend((sym.diff(energy, theta), sym.diff(energy, phi))) + grad_vector.extend((sym.diff(energy, theta), sym.diff(energy, + phi))) symbols.extend((theta, phi)) return sym.lambdify(symbols, grad_vector, accel) @@ -417,7 +490,8 @@ def solve(self, perturbation: float = 1e-3, ftol: float = 0.01e9, max_freq: float = 80e9, - force_single_layer: bool = False): + force_single_layer: bool = False, + force_sb: bool = False): """Solves the system. 1. Computes the energy functional. 2. Computes the gradient of the energy functional. @@ -435,6 +509,11 @@ def solve(self, :param pertubarion: the perturbation to use for the numerical gradient computation. :param ftol: tolerance for the frequency search. [numerical only] :param max_freq: maximum frequency to search for. [numerical only] + :param force_single_layer: whether to force the computation of the frequencies + for each layer individually. + :param force_sb: whether to force the computation of the frequencies. + Takes effect only if the layers are LayerDynamic, not LayerSB. + :return: equilibrium position and frequencies in [GHz] (and eigenvectors if LayerDynamic instead of LayerSB). """ if self.H is None: raise ValueError( @@ -447,9 +526,12 @@ def solve(self, first_momentum_decay=first_momentum_decay, second_momentum_decay=second_momentum_decay, perturbation=perturbation) + if not force_sb and isinstance(self.layers[0], LayerDynamic): + eigenvalues, eigenvectors = self.dynamic_layer_solve(eq) + return eq, eigenvalues / 1e9, eigenvectors N = len(self.layers) if N == 1: - return eq, self.single_layer_resonance(0, eq) / 1e9 + return eq, [self.single_layer_resonance(0, eq) / 1e9] if force_single_layer: frequencies = [] for indx in range(N): @@ -458,6 +540,20 @@ def solve(self, return eq, frequencies return self.num_solve(eq, ftol=ftol, max_freq=max_freq) + def dynamic_layer_solve(self, eq: List[float]): + """Return the FMR frequencies and modes for N layers using the + dynamic RHS model + :param eq: the equilibrium position of the system. + :return: frequencies and eigenmode vectors.""" + jac, symbols = self.compose_llg_jacobian(self.H) + subs = {symbols[i]: eq[i] for i in range(len(eq))} + jac = jac.subs(subs) + jac = np.asarray(jac, dtype=np.float32) + eigvals, eigvecs = np.linalg.eig(jac) + eigvals_im = np.imag(eigvals) / (2 * np.pi) + indx = np.argwhere(eigvals_im > 0).ravel() + return eigvals_im[indx], eigvecs[indx] + def num_solve(self, eq: List[float], ftol: float = 0.01e9, diff --git a/cmtj/utils/filters.py b/cmtj/utils/filters.py index 72e2ca4..355c01a 100644 --- a/cmtj/utils/filters.py +++ b/cmtj/utils/filters.py @@ -28,10 +28,10 @@ def butter_bandpass_filter(data: np.ndarray, btype="bandpass", analog=False, ) - except ValueError: + except ValueError as e: print(fs, pass_freq, nyq, 0.9 * pass_freq / nyq, pass_freq / nyq) - raise ValueError("Error in filtering") + raise ValueError("Error in filtering") from e elif isinstance(pass_freq, tuple): b, a = butter(order, [pass_freq[0], pass_freq[1]], btype="bandpass", diff --git a/cmtj/utils/general.py b/cmtj/utils/general.py index a8b7c4b..3844f1e 100644 --- a/cmtj/utils/general.py +++ b/cmtj/utils/general.py @@ -37,6 +37,12 @@ def __rmul__(self, other: Union["VectorObj", float]): def __repr__(self) -> str: return f"VectorObj(theta={self.theta}, phi={self.phi}, mag={self.mag})" + def __hash__(self) -> int: + return hash(str(self)) + + def __eq__(self, __value: "VectorObj") -> bool: + return self.theta == __value.theta and self.phi == __value.phi and self.mag == __value.mag + def _componentwise_mul(self, other): coors = self.get_cartesian() other_coords = other.get_cartesian() diff --git a/cmtj/utils/optimization.py b/cmtj/utils/optimization.py new file mode 100644 index 0000000..222a1db --- /dev/null +++ b/cmtj/utils/optimization.py @@ -0,0 +1,32 @@ +from typing import Callable, Dict + +import numpy as np +from tqdm import tqdm + + +def coordinate_descent(operating_point: Dict[str, float], + fn: Callable, + best_mse: float = float("-inf"), + granularity: int = 10, + percentage: float = 0.05): + """Performs coordinate descent on the operating point. + :param operating_point: operating point to be optimised. Order of that dict matters. + :param fn: function to be optimised + :param best_mse: best mse so far + :param granularity: granularity of the search + :param percentage: percentage of the search + :returns: best operating point, best mse + """ + opt_params = operating_point + for k, org_v in tqdm(operating_point.items(), desc="Coordinate descent"): + new_params = operating_point.copy() + for v in tqdm(np.linspace((1 - percentage) * org_v, + (1 + percentage) * org_v, granularity), + desc=f"Optimising {k}", + leave=False): + new_params[k] = v + mse = fn(**new_params) + if mse > best_mse: + opt_params = new_params + best_mse = mse + return opt_params, best_mse diff --git a/cmtj/utils/procedures.py b/cmtj/utils/procedures.py index 32610c2..093f502 100644 --- a/cmtj/utils/procedures.py +++ b/cmtj/utils/procedures.py @@ -1,25 +1,23 @@ import math from collections import defaultdict from dataclasses import dataclass -from typing import Any, Dict, List, Tuple +from typing import Any, Callable, Dict, List, Tuple import numpy as np from scipy.fft import fft, fftfreq from tqdm import tqdm from cmtj import AxialDriver, Axis, CVector, Junction, NullDriver, ScalarDriver -from cmtj.utils.filters import Filters from .resistance import calculate_resistance_series, compute_sd @dataclass class ResistanceParameters: - """A data holder for resistance parameters. Not all have to be filled""" + """A data holder for resistance parameters. Not all have to be filled in.""" + Rxx0: float = 0 Rxy0: float = 0 - Rp: float = 100 - Rap: float = 200 Rahe: float = 0 Rsmr: float = 0 Ramr: float = 0 @@ -29,6 +27,7 @@ class ResistanceParameters: def compute_spectrum_strip(input_m: np.ndarray, int_step: float, max_frequency: float): + """Compute the spectrum of a given magnetization trajectory.""" yf = np.abs(fft(input_m)) freqs = fftfreq(len(yf), int_step) freqs = freqs[:len(freqs) // 2] @@ -42,7 +41,7 @@ def compute_spectrum_strip(input_m: np.ndarray, int_step: float, def PIMM_procedure( - junction: 'Junction', + junction: "Junction", Hvecs: np.ndarray, int_step: float, resistance_params: List[ResistanceParameters], @@ -52,10 +51,16 @@ def PIMM_procedure( simulation_duration: float = 5e-9, wait_time: float = 0e-9, max_frequency: float = 80e9, + resistance_fn: Callable = calculate_resistance_series, disturbance: float = 1e-3, - full_output: bool = False + take_last_n: int = 100, + full_output: bool = False, + disable_tqdm: bool = False, + static_only: bool = False, ) -> Tuple[np.ndarray, np.ndarray, Dict[str, Any]]: - """Procedure for computing Pulse Induced Microwave Magnetometry + """Procedure for computing Pulse Induced Microwave Magnetometry. + It computes both PIMM and Resistance (for instance AHE loops). + Set `static_only` to True to only compute the static resistance. :param junction: junction to be simulated. :param Hvecs: list of cartesian vectors. (use FieldScan.amplitude_scan or alike) :param int_step: integration step [s]. @@ -65,8 +70,13 @@ def PIMM_procedure( :param wait_time: time to wait before taking vector for the fft [s]. :param Hoe_duration: duration of Hoe excitation in multiples of in step :param max_frequency: maximum frequency -- larger will be dropped [Hz]. + :param resistance_fn: function to be used to compute the resistance + (either calculate_resistance_series or calculate_resistance_parallel). :param disturbance: disturbance to be applied to the magnetization (std of normal distribution). + :param take_last_n: number of last time steps to be taken for the compuation. :param full_output: if True, return the full trajectories and per layer spectra. + :param disable_tqdm: if True, disable tqdm progress bar. + :param static_only: if True, only compute the static resistance. :return: (spectrum, frequencies, other_data) other_data is a dictionary with the following keys: - 'H': Hext field [A/m] @@ -80,82 +90,103 @@ def PIMM_procedure( spectrum = [] extraction_m_component = None if Hoe_direction == Axis.zaxis: - extraction_m_component = 'z' + extraction_m_component = "z" oedriver = AxialDriver( - NullDriver(), NullDriver(), + NullDriver(), + NullDriver(), ScalarDriver.getStepDriver(0, Hoe_excitation, 0, - int_step * Hoe_duration)) + int_step * Hoe_duration), + ) elif Hoe_direction == Axis.yaxis: - extraction_m_component = 'y' + extraction_m_component = "y" oedriver = AxialDriver( NullDriver(), ScalarDriver.getStepDriver(0, Hoe_excitation, 0, - int_step * Hoe_duration), NullDriver()) + int_step * Hoe_duration), + NullDriver(), + ) else: - extraction_m_component = 'x' + extraction_m_component = "x" oedriver = AxialDriver( ScalarDriver.getStepDriver(0, Hoe_excitation, 0, - int_step * Hoe_duration), NullDriver(), - NullDriver()) + int_step * Hoe_duration), + NullDriver(), + NullDriver(), + ) # get layer strings layer_ids = junction.getLayerIds() + if len(layer_ids) != len(resistance_params): + raise ValueError( + "The number of layers in the junction must match the number of resistance parameters!" + ) output = defaultdict(list) normalising_factor = np.sum( [layer.thickness * layer.Ms for layer in junction.layers]) - for H in tqdm(Hvecs, desc="Computing PIMM"): + freqs = None # in case of static_only + for H in tqdm(Hvecs, desc="Computing PIMM", disable=disable_tqdm): junction.clearLog() junction.setLayerExternalFieldDriver( "all", - AxialDriver(ScalarDriver.getConstantDriver(H[0]), - ScalarDriver.getConstantDriver(H[1]), - ScalarDriver.getConstantDriver(H[2]))) + AxialDriver( + ScalarDriver.getConstantDriver(H[0]), + ScalarDriver.getConstantDriver(H[1]), + ScalarDriver.getConstantDriver(H[2]), + ), + ) junction.setLayerOerstedFieldDriver("all", oedriver) if disturbance: for layer_id in layer_ids: old_mag = junction.getLayerMagnetisation(layer_id) - new_mag = CVector(old_mag.x + np.random.normal(0, disturbance), - old_mag.y + np.random.normal(0, disturbance), - old_mag.z + np.random.normal(0, disturbance)) + new_mag = CVector( + old_mag.x + np.random.normal(0, disturbance), + old_mag.y + np.random.normal(0, disturbance), + old_mag.z + np.random.normal(0, disturbance), + ) new_mag.normalize() junction.setLayerMagnetisation(layer_id, new_mag) junction.runSimulation(simulation_duration, int_step, int_step) log = junction.getLog() - indx = np.argwhere(np.asarray(log['time']) >= wait_time).ravel() + indx = np.argwhere(np.asarray(log["time"]) >= wait_time).ravel() m_traj = np.asarray([ np.asarray([ - log[f'{layer.id}_mx'], log[f'{layer.id}_my'], - log[f'{layer.id}_mz'] + log[f"{layer.id}_mx"], + log[f"{layer.id}_my"], + log[f"{layer.id}_mz"], ]) * layer.thickness * layer.Ms / normalising_factor for layer in junction.layers ]) - m = m_traj[:, :, -100:] # all layers, all x, y, z, last 100 steps - mixed = np.asarray([ - np.asarray(log[f"{layer.id}_m{extraction_m_component}"])[indx] * - layer.thickness * layer.Ms / normalising_factor - for layer in junction.layers - ]) - mixed = np.squeeze(mixed) - mixed_sum = mixed.sum(axis=0) - yf, freqs = compute_spectrum_strip(mixed_sum, int_step, max_frequency) - - spectrum.append(yf) - Rx, Ry = calculate_resistance_series( - [r.Rxx0 - for r in resistance_params], [r.Rxy0 for r in resistance_params], - [r.Ramr - for r in resistance_params], [r.Rahe for r in resistance_params], + m = m_traj[:, :, + -take_last_n:] # all layers, all x, y, z, last 100 steps + Rx, Ry = resistance_fn( + [r.Rxx0 for r in resistance_params], + [r.Rxy0 for r in resistance_params], + [r.Ramr for r in resistance_params], + [r.Rahe for r in resistance_params], [r.Rsmr for r in resistance_params], m, l=[r.l for r in resistance_params], - w=[r.w for r in resistance_params]) + w=[r.w for r in resistance_params], + ) + if not static_only: + mixed = np.asarray([ + np.asarray(log[f"{layer.id}_m{extraction_m_component}"])[indx] + * layer.thickness * layer.Ms / normalising_factor + for layer in junction.layers + ]) + mixed_sum = mixed.sum(axis=0) + yf, freqs = compute_spectrum_strip(mixed_sum, int_step, + max_frequency) + + spectrum.append(yf) + # fill the output dict - output['H'].append(H) - output['Rx'].append(Rx) - output['Ry'].append(Ry) - output['m_avg'].append(m_traj[:, :, -1].sum(0)) - if full_output: - output['m_traj'].append(m_traj) + output["H"].append(H) + output["Rx"].append(Rx) + output["Ry"].append(Ry) + output["m_avg"].append(m_traj[:, :, -1].sum(0)) + if full_output and not static_only: + output["m_traj"].append(m_traj) for li, layer_id in enumerate(layer_ids): y, _ = compute_spectrum_strip(mixed[li], int_step, max_frequency) @@ -167,16 +198,20 @@ def PIMM_procedure( return spectrum, freqs, output -def VSD_procedure(junction: Junction, - Hvecs: np.ndarray, - frequencies: np.ndarray, - int_step: float, - resistance_params, - Hoe_direction: Axis = Axis.yaxis, - Hoe_excitation: float = 50, - simulation_duration: float = 30e-9, - disturbance: float = 1e-3, - Rtype: str = 'Rz'): +def VSD_procedure( + junction: Junction, + Hvecs: np.ndarray, + frequencies: np.ndarray, + int_step: float, + resistance_params: List[ResistanceParameters] = [], + Hoe_direction: Axis = Axis.yaxis, + Hoe_excitation: float = 50, + simulation_duration: float = 30e-9, + disturbance: float = 1e-3, + Rtype: str = "Rz", + resistance_fn: Callable = calculate_resistance_series, + disable_tqdm: bool = False, +): """Procedure for computing Voltage-Spin Diode. We use the Oersted field sine exctitation to excite the system. :param junction: junction to be simulated. @@ -188,54 +223,85 @@ def VSD_procedure(junction: Junction, :param Hoe_excitation: excitation amplitude of Hoe [A/m]. :param simulation_duration: duration of simulation [s]. :param disturbance: disturbance to be applied to the magnetization (std of normal distribution). + :param resistance_fn: function to be used to compute the resistance + (either calculate_resistance_series or calculate_resistance_parallel). Rz forces standard magnetores. :param Rtype: type of resistance to be used. (Rx Ry or Rz) + :param disable_tqdm: if True, disable tqdm progress bar. """ layer_ids = junction.getLayerIds() + if Rtype == "Rz" and len(layer_ids) > 2: + raise ValueError( + "Rz can only be used for 2 layer junctions. Use Rx or Ry instead.") + elif len(resistance_params) != len(layer_ids): + raise ValueError( + "The number of layers in the junction must match the number of resistance parameters!" + ) - def simulate_VSD(H: np.ndarray, frequency: float, resistance_params): - + def simulate_VSD(H: np.ndarray, frequency: float, + resistance_params: ResistanceParameters): if Hoe_direction == Axis.zaxis: oedriver = AxialDriver( - NullDriver(), NullDriver(), - ScalarDriver.getSineDriver(0, Hoe_excitation, frequency, 0)) + NullDriver(), + NullDriver(), + ScalarDriver.getSineDriver(0, Hoe_excitation, frequency, 0), + ) elif Hoe_direction == Axis.yaxis: oedriver = AxialDriver( NullDriver(), ScalarDriver.getSineDriver(0, Hoe_excitation, frequency, 0), - NullDriver()) + NullDriver(), + ) else: oedriver = AxialDriver( ScalarDriver.getSineDriver(0, Hoe_excitation, frequency, 0), - NullDriver(), NullDriver()) + NullDriver(), + NullDriver(), + ) junction.clearLog() junction.setLayerExternalFieldDriver( "all", - AxialDriver(ScalarDriver.getConstantDriver(H[0]), - ScalarDriver.getConstantDriver(H[1]), - ScalarDriver.getConstantDriver(H[2]))) + AxialDriver( + ScalarDriver.getConstantDriver(H[0]), + ScalarDriver.getConstantDriver(H[1]), + ScalarDriver.getConstantDriver(H[2]), + ), + ) junction.setLayerOerstedFieldDriver("all", oedriver) if disturbance: for layer_id in layer_ids: old_mag = junction.getLayerMagnetisation(layer_id) - new_mag = CVector(old_mag.x + np.random.normal(0, disturbance), - old_mag.y + np.random.normal(0, disturbance), - old_mag.z + np.random.normal(0, disturbance)) + new_mag = CVector( + old_mag.x + np.random.normal(0, disturbance), + old_mag.y + np.random.normal(0, disturbance), + old_mag.z + np.random.normal(0, disturbance), + ) new_mag.normalize() junction.setLayerMagnetisation(layer_id, new_mag) junction.runSimulation(simulation_duration, int_step, int_step) log = junction.getLog() m_traj = np.asarray([[ - log[f'{layer_ids[i]}_mx'], log[f'{layer_ids[i]}_my'], - log[f'{layer_ids[i]}_mz'] + log[f"{layer_ids[i]}_mx"], + log[f"{layer_ids[i]}_my"], + log[f"{layer_ids[i]}_mz"], ] for i in range(len(layer_ids))]) - if Rtype == 'Rz': - if isinstance(resistance_params, - list) and len(resistance_params) > 1: - resistance_params = resistance_params[0] - R = log[f'R_{layer_ids[0]}_{layer_ids[1]}'] + if Rtype == "Rz": + if len(layer_ids) > 2: + raise ValueError( + "Rz can only be used for 2 layer junctions. One layer can be fictisious." + ) + elif len(layer_ids) == 2: + R = log[f"R_{layer_ids[0]}_{layer_ids[1]}"] + elif len(layer_ids) == 1: + R = log["Resistance"] + else: + raise ValueError( + "Resistance definition ambiguous!" + "If you want to use Rz, you must provide" + "a single resistance parameter set or set Rp Rap" + " at junction creation.") else: - Rx, Ry = calculate_resistance_series( + Rx, Ry = resistance_fn( [r.Rxx0 for r in resistance_params], [r.Rxy0 for r in resistance_params], [r.Ramr for r in resistance_params], @@ -243,19 +309,21 @@ def simulate_VSD(H: np.ndarray, frequency: float, resistance_params): [r.Rsmr for r in resistance_params], m_traj, l=[r.l for r in resistance_params], - w=[r.w for r in resistance_params]) - if Rtype == 'Rx': + w=[r.w for r in resistance_params], + ) + if Rtype == "Rx": R = Rx - elif Rtype == 'Ry': + elif Rtype == "Ry": R = Ry else: raise ValueError("Rtype must be either Rx or Ry or Rz") - dynamicI = np.sin(2 * math.pi * frequency * np.asarray(log['time'])) + dynamicI = np.sin(2 * math.pi * frequency * np.asarray(log["time"])) vmix = compute_sd(R, dynamicI, int_step) return vmix spectrum = np.zeros((len(Hvecs), len(frequencies))) - for hindx, H in enumerate(tqdm(Hvecs, "Computing VSD")): + for hindx, H in enumerate( + tqdm(Hvecs, "Computing VSD", disable=disable_tqdm)): for findx, f in enumerate(frequencies): spectrum[hindx, findx] = simulate_VSD(H, f, resistance_params) return spectrum diff --git a/core/cvector.hpp b/core/cvector.hpp index 81f8902..c576b6a 100644 --- a/core/cvector.hpp +++ b/core/cvector.hpp @@ -11,7 +11,6 @@ class CVector { public: - // friend std::ostream& operator<<(std::ostream &o, const CVector &obj); T x, y, z; CVector() { @@ -19,6 +18,7 @@ class CVector this->y = 0.0; this->z = 0.0; } + explicit CVector(std::vector vec) { if (vec.size() != 3) @@ -36,6 +36,7 @@ class CVector this->y = y; this->z = z; } + CVector(const CVector& v) { this->x = v.x; @@ -66,6 +67,7 @@ class CVector this->z -= v.z; return *this; } + CVector operator+(CVector v) { CVector res( @@ -120,6 +122,7 @@ class CVector y = v.y; z = v.z; } + bool operator==(const CVector& v) { if ( diff --git a/core/drivers.hpp b/core/drivers.hpp index 4e8c5f5..657dbca 100644 --- a/core/drivers.hpp +++ b/core/drivers.hpp @@ -422,6 +422,15 @@ class AxialDriver : public Driver { } + explicit AxialDriver( + const T x, const T y, const T z + ) : AxialDriver( + ScalarDriver::getConstantDriver(x), + ScalarDriver::getConstantDriver(y), + ScalarDriver::getConstantDriver(z)) + { + } + explicit AxialDriver(std::vector> axialDrivers) { if (axialDrivers.size() != 3) @@ -460,6 +469,7 @@ class AxialDriver : public Driver this->drivers[1].constantValue, this->drivers[2].constantValue); } + /** * Returns the mask for the Axial Driver. * For instance: a vector (1213, 123, 0) returns (1, 1, 0) diff --git a/core/junction.hpp b/core/junction.hpp index 1f2dd3d..4643c5b 100644 --- a/core/junction.hpp +++ b/core/junction.hpp @@ -1020,6 +1020,8 @@ class Junction // std::string fileSave; unsigned int logLength = 0; unsigned int layerNo; + std::string Rtag = "R"; + Junction() {} /** @@ -1057,6 +1059,8 @@ class Junction this->Rp = Rp; this->Rap = Rap; this->MR_mode = CLASSIC; + // A string representing the tag for the junction's resistance value. + this->Rtag = "R_" + this->layers[0].id + "_" + this->layers[1].id; } /** @@ -1432,7 +1436,7 @@ class Junction { const auto magnetoresistance = calculateMagnetoresistance(c_dot(this->layers[0].mag, this->layers[1].mag)); - this->log["R_free_bottom"].emplace_back(magnetoresistance); + this->log[this->Rtag].emplace_back(magnetoresistance); } else if (MR_mode == STRIP) { diff --git a/Dockerfile b/docker/Dockerfile similarity index 64% rename from Dockerfile rename to docker/Dockerfile index 3729885..a832e16 100644 --- a/Dockerfile +++ b/docker/Dockerfile @@ -1,5 +1,4 @@ -FROM ubuntu:20.04 - +FROM python:3.9 # let's copy all the necessary files # we install our only 2 dependencies :) and vim for nice workflow @@ -7,12 +6,7 @@ RUN apt-get update && \ apt-get install -y build-essential libfftw3-dev git python3 python3-pip vim RUN python3 -m pip install --upgrade pip WORKDIR /scratch -# COPY cmtj /scratch/cmtj -# COPY core /scratch/core -# COPY python /scratch/python -# COPY setup.py setup.cfg ./ -# RUN python3 -m pip install . -RUN git clone -b logging-fixes https://github.com/LemurPwned/cmtj.git && \ +RUN git clone https://github.com/LemurPwned/cmtj.git && \ cd cmtj && \ python3 -m pip install numpy scipy && \ python3 -m pip install . diff --git a/docker/Dockerfile.app b/docker/Dockerfile.app new file mode 100644 index 0000000..69771ad --- /dev/null +++ b/docker/Dockerfile.app @@ -0,0 +1,12 @@ +FROM python:3.11 +EXPOSE 8501 + +WORKDIR /app +COPY view/ /app/ +RUN ls +RUN python3 -m pip install streamlit && \ + git clone https://github.com/LemurPwned/cmtj.git && \ + cd cmtj && \ + python3 -m pip install .[utils] + +CMD streamlit run streamlit_app.py diff --git a/python/cmtj.cpp b/python/cmtj.cpp index f7897b2..5bc7df1 100644 --- a/python/cmtj.cpp +++ b/python/cmtj.cpp @@ -44,9 +44,15 @@ PYBIND11_MODULE(cmtj, m) // operators .def(py::self + py::self) .def(py::self += py::self) + .def(py::self - py::self) + .def(py::self -= py::self) .def(py::self *= double()) + .def(py::self == py::self) + .def(py::self != py::self) .def(double() * py::self) .def(py::self * double()) + .def("__getitem__", [](const DVector& v, const int key) { return v[key]; }) + .def("__len__", [](const DVector& v) { return 3; }) .def("__str__", &DVector::toString) .def("__repr__", &DVector::toString); @@ -127,6 +133,8 @@ PYBIND11_MODULE(cmtj, m) py::class_(m, "AxialDriver") .def(py::init()) .def(py::init>>()) + .def(py::init()) + .def(py::init()) .def("getVectorAxialDriver", &DAxialDriver::getVectorAxialDriver) .def("getCurrentAxialDrivers", &DAxialDriver::getCurrentAxialDrivers) diff --git a/setup.cfg b/setup.cfg index f4ff2ef..046e0f0 100644 --- a/setup.cfg +++ b/setup.cfg @@ -1,7 +1,7 @@ [metadata] -description-file = CMTJ.md -long_description = file: CMTJ.md +long_description = file: README.md long_description_content_type = text/markdown +license_file = LICENSE classifiers = Programming Language :: Python :: 3 Programming Language :: C++ @@ -10,12 +10,17 @@ classifiers = Topic :: Scientific/Engineering Topic :: Scientific/Engineering :: Physics Typing :: Stubs Only + Typing :: Typed [options.extras_require] - utils = matplotlib - numpy - scipy - pandas - tqdm - multiprocess - numba +utils = + matplotlib + numpy + scipy + pandas + tqdm + multiprocess + numba + sympy +test = + pytest diff --git a/setup.py b/setup.py index 2cf1fee..51e8cd3 100644 --- a/setup.py +++ b/setup.py @@ -1,11 +1,11 @@ +import contextlib import os -import sys import setuptools from setuptools import Extension, find_namespace_packages, setup from setuptools.command.build_ext import build_ext -__version__ = "1.3.2" +__version__ = "1.4.0" """ As per https://github.com/pybind/python_example @@ -24,7 +24,6 @@ class get_pybind_include(object): until it is actually installed, so that the ``get_include()`` method can be invoked. """ - def __str__(self): import pybind11 @@ -68,10 +67,8 @@ def has_flag(compiler, flagname): except setuptools.distutils.errors.CompileError: return False finally: - try: + with contextlib.suppress(OSError): os.remove(fname) - except OSError: - pass return True @@ -122,7 +119,8 @@ def build_extensions(self): opts.append("-fvisibility=hidden") for ext in self.extensions: - ext.define_macros = [("VERSION_INFO", f'"{self.distribution.get_version()}"')] + ext.define_macros = [("VERSION_INFO", + f'"{self.distribution.get_version()}"')] ext.extra_compile_args = opts ext.extra_link_args = link_opts build_ext.build_extensions(self) @@ -136,7 +134,6 @@ def build_extensions(self): author_email="mojsieju@agh.edu.pl", url="https://github.com/LemurPwned/cmtj", description="CMTJ - C Magnetic Tunnel Junctions.", - long_description="Fast library for simulating magnetic multilayers.", ext_modules=ext_modules, include_package_data=True, namespace_packages=["cmtj"], diff --git a/tests/conftest.py b/tests/conftest.py new file mode 100644 index 0000000..0fc65e4 --- /dev/null +++ b/tests/conftest.py @@ -0,0 +1,228 @@ +import pytest +from typing import Tuple +from cmtj.utils import mu0 +from cmtj import Junction, CVector, Layer, ScalarDriver +from cmtj.utils.procedures import ResistanceParameters +from cmtj.models import LayerDynamic, VectorObj + +rp = ResistanceParameters(Rxx0=100, + Rxy0=1, + Rsmr=-0.46, + Rahe=-2.7, + Ramr=-0.24, + l=30, + w=20) + + +@pytest.fixture +def arg(request): + return request.getfixturevalue(request.param) + + +@pytest.fixture +def two_layer_symbolic_dyn() -> Tuple[LayerDynamic]: + layerA = LayerDynamic( + 0, + thickness=2.3e-9, + Kv=VectorObj(4e3, 0), + Ks=100, + Ms=1.8 / mu0, + alpha=1e-3, + ) + layerB = LayerDynamic( + 1, + thickness=3.2e-9, + Kv=VectorObj(340e3, 0), + Ks=100, + Ms=1.2 / mu0, + alpha=1e-3, + ) + return (layerA, layerB) + + +@pytest.fixture +def single_layer_mtj() -> Tuple[Junction, ResistanceParameters]: + # approximate demagnetisation tensor + demag = [CVector(0, 0, 0), CVector(0, 0, 0), CVector(0, 0, 1)] + alpha = 0.005 + Kdir = CVector(1, 0, 0) + l1 = Layer("free", + mag=CVector(0, 0, 1), + anis=Kdir, + Ms=1.65, + thickness=3e-9, + cellSurface=0, + demagTensor=demag, + damping=alpha) + K1 = 1.05e3 + junction = Junction([l1]) + junction.setLayerAnisotropyDriver("free", + ScalarDriver.getConstantDriver(K1)) + + return (junction, [rp]) + + +@pytest.fixture +def single_layer_mtj_fictious() -> Tuple[Junction, ResistanceParameters]: + # approximate demagnetisation tensor + demag = [CVector(0, 0, 0), CVector(0, 0, 0), CVector(0, 0, 1)] + alpha = 0.005 + Kdir = CVector(1, 0, 0) + l1 = Layer("free", + mag=CVector(0, 0, 1), + anis=Kdir, + Ms=1.65, + thickness=3e-9, + cellSurface=0, + demagTensor=demag, + damping=alpha) + K1 = 1.05e3 + l1.setReferenceLayer(CVector(1, 0, 0)) + junction = Junction([l1], 100, 200) + junction.setLayerAnisotropyDriver("free", + ScalarDriver.getConstantDriver(K1)) + + return (junction, None) + + +@pytest.fixture +def two_layer_mtj() -> Tuple[Junction, ResistanceParameters]: + # approximate demagnetisation tensor + demag = [CVector(0, 0, 0), CVector(0, 0, 0), CVector(0, 0, 1)] + alpha = 0.005 + Kdir = CVector(1, 0, 0) + l1 = Layer("free", + mag=CVector(0, 0, 1), + anis=Kdir, + Ms=1.65, + thickness=3e-9, + cellSurface=0, + demagTensor=demag, + damping=alpha) + + l2 = Layer("bottom", + mag=CVector(0, 0, 1), + anis=Kdir, + Ms=1.45, + thickness=4e-9, + cellSurface=0, + demagTensor=demag, + damping=alpha) + J1 = -1.78e-3 + J2 = -1.69e-4 + + K1 = 1.05e3 + K2 = 50e3 + + junction = Junction([l1, l2]) + + junction.setLayerAnisotropyDriver("free", + ScalarDriver.getConstantDriver(K1)) + + junction.setLayerAnisotropyDriver("bottom", + ScalarDriver.getConstantDriver(K2)) + junction.setIECDriver("free", "bottom", ScalarDriver.getConstantDriver(J1)) + junction.setQuadIECDriver("free", "bottom", + ScalarDriver.getConstantDriver(J2)) + return (junction, [rp, rp]) + + +@pytest.fixture +def two_layer_mtj_rz() -> Tuple[Junction, ResistanceParameters]: + # approximate demagnetisation tensor + demag = [CVector(0, 0, 0), CVector(0, 0, 0), CVector(0, 0, 1)] + alpha = 0.005 + Kdir = CVector(1, 0, 0) + l1 = Layer("free", + mag=CVector(0, 0, 1), + anis=Kdir, + Ms=1.65, + thickness=3e-9, + cellSurface=0, + demagTensor=demag, + damping=alpha) + + l2 = Layer("bottom", + mag=CVector(0, 0, 1), + anis=Kdir, + Ms=1.45, + thickness=4e-9, + cellSurface=0, + demagTensor=demag, + damping=alpha) + J1 = -1.78e-3 + J2 = -1.69e-4 + + K1 = 1.05e3 + K2 = 50e3 + + junction = Junction([l1, l2], 100, 200) + + junction.setLayerAnisotropyDriver("free", + ScalarDriver.getConstantDriver(K1)) + + junction.setLayerAnisotropyDriver("bottom", + ScalarDriver.getConstantDriver(K2)) + junction.setIECDriver("free", "bottom", ScalarDriver.getConstantDriver(J1)) + junction.setQuadIECDriver("free", "bottom", + ScalarDriver.getConstantDriver(J2)) + return (junction, None) + + +@pytest.fixture +def tri_layer_mtj() -> Tuple[Junction, ResistanceParameters]: + # approximate demagnetisation tensor + demag = [CVector(0, 0, 0), CVector(0, 0, 0), CVector(0, 0, 1)] + alpha = 0.005 + Kdir = CVector(1, 0, 0) + l1 = Layer("free", + mag=CVector(0, 0, 1), + anis=Kdir, + Ms=1.65, + thickness=3e-9, + cellSurface=0, + demagTensor=demag, + damping=alpha) + + l2 = Layer("middle", + mag=CVector(0, 0, 1), + anis=Kdir, + Ms=1.45, + thickness=4e-9, + cellSurface=0, + demagTensor=demag, + damping=alpha) + + l3 = Layer("bottom", + mag=CVector(0., 0, 1.), + anis=CVector(0, 0, 1.), + Ms=0.7, + thickness=1e-9, + cellSurface=0, + demagTensor=demag, + damping=alpha) + J11 = -1.78e-3 + J12 = -1.69e-4 + + J21 = 3e-4 + J22 = 2e-6 + + K1 = 1.05e3 + K2 = 50e3 + K3 = 1e4 + junction = Junction([l1, l2, l3]) + + for l, anis_val in zip(["free", "middle", "bottom"], [K1, K2, K3]): + junction.setLayerAnisotropyDriver( + l, ScalarDriver.getConstantDriver(anis_val)) + + junction.setIECDriver("free", "middle", + ScalarDriver.getConstantDriver(J11)) + junction.setQuadIECDriver("free", "middle", + ScalarDriver.getConstantDriver(J12)) + + junction.setIECDriver("middle", "bottom", + ScalarDriver.getConstantDriver(J21)) + junction.setQuadIECDriver("middle", "bottom", + ScalarDriver.getConstantDriver(J22)) + return (junction, [rp, rp, rp]) \ No newline at end of file diff --git a/tests/test_drivers.py b/tests/test_drivers.py new file mode 100644 index 0000000..2ee5c0a --- /dev/null +++ b/tests/test_drivers.py @@ -0,0 +1,39 @@ +from cmtj import AxialDriver, CVector + + +def test_cvector_operators(): + vec1 = (1.0, 2.0, 3.0) + vec2 = (1.0, 2.0, 3.0) + vec3 = (2, 4, 6) + + cvec1 = CVector(*vec1) + cvec2 = CVector(*vec2) + cvec3 = CVector(*vec3) + + # test assignment operators + assert vec1[0] == cvec1[0] + assert vec1[1] == cvec1[1] + assert vec1[2] == cvec1[2] + + # test ops + assert cvec1 == cvec2 + assert cvec1 != cvec3 + assert (cvec1 + cvec2) == cvec3 + assert (cvec3 - cvec1) == cvec1 + + # test self-assignment operators + assert 2 * cvec1 == cvec3 + assert cvec3 * 0.5 == cvec1 + cvec1 += CVector(1.0, 2.0, 3.0) + assert cvec1 == cvec3 + cvec3 -= cvec2 + assert cvec3 == cvec2 + + +def test_axial_definitions(): + vec = (1.0, 2.0, 3.0) + cvec = CVector(*vec) + d1 = AxialDriver(*vec) + assert d1.getCurrentAxialDrivers(0.0) == cvec + d1 = AxialDriver(cvec) + assert d1.getCurrentAxialDrivers(0.0) == cvec diff --git a/tests/test_models.py b/tests/test_models.py new file mode 100644 index 0000000..45d4ee9 --- /dev/null +++ b/tests/test_models.py @@ -0,0 +1,33 @@ +import numpy as np +from cmtj.models import Solver, LayerDynamic, VectorObj +from typing import Tuple + + +def test_sb_dynamic(two_layer_symbolic_dyn: Tuple[LayerDynamic]): + J1 = -1e-3 + J2 = 1e-4 + + for Hmag in np.linspace(-300e3, 300e3, 8): + H = VectorObj(np.deg2rad(88), np.deg2rad(0.1), Hmag) + solver_dyn = Solver(layers=two_layer_symbolic_dyn, + J1=[J1], + J2=[J2], + H=H) + pos = np.asarray( + [np.deg2rad(88), + np.deg2rad(0.1), + np.deg2rad(88), + np.deg2rad(0.1)]) + # set perturbation to 0 to avoid numerical errors + eq_sb, f_sb = solver_dyn.solve(init_position=pos, + perturbation=0, + max_steps=1e8, + force_sb=True) + eq_dyn, f_dyn, _ = solver_dyn.solve(init_position=pos, + max_steps=1e8, + perturbation=0) + f_sb.sort() + f_dyn.sort() + assert np.allclose(eq_sb, eq_dyn) + assert np.allclose(f_sb, f_dyn, atol=0.2) + pos = eq_dyn \ No newline at end of file diff --git a/tests/test_procedures.py b/tests/test_procedures.py new file mode 100644 index 0000000..3a607d4 --- /dev/null +++ b/tests/test_procedures.py @@ -0,0 +1,85 @@ +from itertools import product +from typing import Tuple + +import numpy as np +import pytest + +from cmtj import Axis, Junction +from cmtj.utils import FieldScan +from cmtj.utils.procedures import (PIMM_procedure, ResistanceParameters, + VSD_procedure) +from cmtj.utils.resistance import (calculate_resistance_parallel, + calculate_resistance_series) + + +def product_dict(inp): + return (dict(zip(inp.keys(), values)) for values in product(*inp.values())) + + +def test_base_pimm(single_layer_mtj: Tuple[Junction, ResistanceParameters]): + junction, res = single_layer_mtj + _, Hvecs = FieldScan.amplitude_scan(-500e3, + 500e3, + steps=50, + theta=90, + phi=0) + parameter_variants = { + "static_only": [True, False], + "full_output": [True, False], + "resistance_fn": + [calculate_resistance_series, calculate_resistance_parallel], + "Hoe_direction": [Axis.xaxis, Axis.yaxis] + } + for param_inp in product_dict(parameter_variants): + spec, f, out = PIMM_procedure(junction, + Hvecs=Hvecs, + resistance_params=res, + int_step=1e-12, + simulation_duration=10e-9, + disable_tqdm=True, + **param_inp) + if param_inp["static_only"]: + assert f is None + assert len(spec) == 0 + elif param_inp['full_output']: + assert 'm_traj' in out + + +@pytest.mark.parametrize( + 'arg', ['single_layer_mtj', 'two_layer_mtj', 'tri_layer_mtj'], + indirect=True) +def test_pimm(arg: Tuple[Junction, ResistanceParameters]): + junction, res = arg + _, Hvecs = FieldScan.amplitude_scan(-500e3, + 500e3, + steps=50, + theta=90, + phi=0) + PIMM_procedure(junction, + Hvecs=Hvecs, + int_step=1e-12, + resistance_params=res, + simulation_duration=10e-9, + disable_tqdm=True) + + +@pytest.mark.parametrize( + 'arg', ['single_layer_mtj', 'two_layer_mtj', 'tri_layer_mtj'], + indirect=True) +def test_vsd(arg: Tuple[Junction, ResistanceParameters]): + junction, res = arg + _, Hvecs = FieldScan.amplitude_scan(-500e3, + 500e3, + steps=50, + theta=90, + phi=0) + # TODO: add more tests for Rz + for rtype in ['Rx']: + VSD_procedure(junction, + Hvecs=Hvecs, + resistance_params=res, + frequencies=np.arange(1e9, 15e9, 1e9), + int_step=1e-12, + simulation_duration=10e-9, + Rtype=rtype, + disable_tqdm=True) diff --git a/tests/test_symbolic.py b/tests/test_symbolic.py new file mode 100644 index 0000000..4652725 --- /dev/null +++ b/tests/test_symbolic.py @@ -0,0 +1,54 @@ +import numpy as np +import sympy as sym +from cmtj.models.general_sb import LayerSB, Solver, VectorObj +import sympy as sym + + +def test_layer_energy(): + # create a test layer + layer = LayerSB(_id=1, + thickness=1.0, + Kv=VectorObj(theta=0, phi=0.0, mag=1.0), + Ks=1.0, + Ms=1.0) + + # create test values for the input parameters + H = sym.ImmutableMatrix([0, 0, 1]) + J1top = 0.0 + J1bottom = 0.0 + J2top = 0.0 + J2bottom = 0.0 + top_layer = None + down_layer = None + + # calculate the energy of the layer + energy = layer.symbolic_layer_energy(H, J1top, J1bottom, J2top, J2bottom, + top_layer, down_layer) + + # check that the energy is a sympy expression + assert isinstance(energy, sym.Expr) + + +def test_solver_init(): + # create test layers + layer1 = LayerSB(_id=0, + thickness=1.0, + Kv=VectorObj(theta=00, phi=0.0, mag=1.0), + Ks=1.0, + Ms=1.0) + layer2 = LayerSB(_id=1, + thickness=1.0, + Kv=VectorObj(theta=0, phi=0.0, mag=1.0), + Ks=1.0, + Ms=1.0) + layers = [layer1, layer2] + + # create test values for J1 and J2 + J1 = [0.0] + J2 = [0.0] + + # create a solver object + solver = Solver(layers=layers, J1=J1, J2=J2) + + # check that the solver object was created successfully + assert isinstance(solver, Solver) diff --git a/view/autofit.py b/view/autofit.py new file mode 100644 index 0000000..5a83b00 --- /dev/null +++ b/view/autofit.py @@ -0,0 +1,94 @@ +from typing import List + +import numpy as np +import streamlit as st +from hebo.design_space.design_space import DesignSpace +from hebo.optimizers.hebo import HEBO +from helpers import read_data +from simulation_fns import compute_sb_mse, simulate_sb + + +def hebo_fit( + hvals: List[float], + y: List[float], + design_space: DesignSpace, + N: int, + n_suggestions: int = 1, +): + def target_fn(yhat, y): + return np.asarray([compute_sb_mse(target=y, data=yhat)]) + + opt = HEBO(design_space) + prog_bar = st.progress(0, text="Optimisation") + try: + for i in range(1, N + 1): + rec = opt.suggest(n_suggestions=n_suggestions) + param_values = rec.to_dict("records") + for k, v in param_values[0].items(): + st.session_state[k] = v + result_dictionary = simulate_sb(hvals=hvals) + opt.observe(rec, target_fn(yhat=result_dictionary, y=y)) + try: + val = opt.y.min() + progres_text = f"({i}/{N}) MSE: {val:.2f}" + prog = int((i * 100) / N) + prog_bar.progress(prog, text=progres_text) + except Exception as e: + print(f"Error printing progress: {e}") + except KeyboardInterrupt: + print("Optimisation interrupted") + return opt + return opt + + +def autofit(placeholder=None): + N = st.session_state.N + cfg = [] + # keep the same units as in the GUI + bounds = { + "Ms": (0.2, 2.0), + "K": (10, 10e3), + } + for param_name in ("Ms", "K"): + cfg.extend( + { + "name": f"{param_name}{i}", + "type": "num", + "lb": bounds[param_name][0], + "ub": bounds[param_name][1], + } + for i in range(N) + ) + cfg.extend( + { + "name": f"J{i}", + "type": "num", + "lb": -1e3, + "ub": 1e3, + } + for i in range(N - 1) + ) + try: + target_data = read_data() + h, f = target_data + f = np.asarray(f).ravel().squeeze().tolist() + target_data = (h, f) + # TODO: sort + # target_data = sorted(zip(*target_data)) + except AttributeError: + st.write("Upload the data to start optimization!") + else: + result = hebo_fit( + hvals=h, + y=f, + design_space=DesignSpace().parse(cfg), + N=st.session_state.n_iters, + ) + placeholder.markdown( + f"""## OPTIMISATION COMPLETE\n + Best MSE: {result.y.min():.2f}\n + Best parameters: {result.best_x.iloc[0].to_dict()} + """ + ) + for k, v in result.best_x.iloc[0].to_dict().items(): + st.session_state[k] = v diff --git a/view/helpers.py b/view/helpers.py new file mode 100644 index 0000000..7e79466 --- /dev/null +++ b/view/helpers.py @@ -0,0 +1,106 @@ +import matplotlib.pyplot as plt +import numpy as np +import streamlit as st +from simulation_fns import get_pimm_data, get_vsd_data + + +def read_data(): + filedata = st.session_state.upload.read().decode("utf-8") + lines = filedata.split("\n") + fields, freqs = [], [] + for line in lines[1:]: + if line.startswith("#"): + continue + vals = line.split() + if len(vals): + fields.append(float(vals[0])) + freqs.append([float(val) for val in vals[1:]]) + return fields, freqs + + +def plot_data(Hscan, freqs, spec, mag=None, title="Resonance spectrum"): + with plt.style.context(["dark_background"]): + if mag is not None: + fig = plt.figure(dpi=300) + gs = plt.GridSpec(3, 6) + indx = 4 + ax1 = fig.add_subplot(gs[:, :indx]) + prev_ax = [] + for i, lab in enumerate("xyz"): + _ax = fig.add_subplot(gs[i, indx:]) + prev_ax.append(_ax) + _ax.plot(Hscan / 1e3, mag[:, i], color="white") + # move ticks to the right side + _ax.yaxis.tick_right() + _ax.yaxis.set_label_position("right") + _ax.set_ylim(-1.1, 1.1) + if lab != "z": + # remove x ticks + _ax.set_xticks([]) + _ax.set_ylabel(rf"$m_{lab}$", rotation=270, labelpad=15) + prev_ax[-1].set_xlabel("H (kA/m)") + fig.subplots_adjust(wspace=0.05, hspace=0.05) + fig.align_ylabels() + else: + fig, ax1 = plt.subplots(dpi=300) + ax1.pcolormesh( + Hscan / 1e3, + freqs / 1e9, + 10 * np.log10(np.abs(spec.T)), + shading="auto", + cmap="inferno", + rasterized=True, + ) + ax1.set_xlabel("H (kA/m)") + ax1.set_ylabel("Frequency (GHz)") + fig.suptitle(title) + try: + fields, freqs = read_data() + freqs = np.asarray(freqs) + fields = np.asarray(fields) + cmap = plt.colormaps["Pastel2"] + for f_indx in range(freqs.shape[1]): + ax1.plot( + fields / 1e3, + freqs[:, f_indx] / 1e9, + "o", + markeredgecolor="white", + color=cmap(f_indx), + label=f"Data {f_indx}", + ) + except (ValueError, AttributeError) as e: + print(f"Error plotting data: {e}") + st.pyplot(fig) + + +def simulate_vsd(): + with st.spinner("Simulating VSD..."): + spec, freqs, Hscan = get_vsd_data( + int_step=st.session_state.int_step, + fmin=st.session_state.fmin * 1e9, + fmax=st.session_state.fmax * 1e9, + H_axis=st.session_state.H_axis, + Hmin=st.session_state.Hmin * 1e3, + Hmax=st.session_state.Hmax * 1e3, + Hsteps=st.session_state.Hsteps, + Hoex=st.session_state.Hoex, + fstep=st.session_state.fstep * 1e9, + Rtype=st.session_state.res_type, + sim_time=st.session_state.sim_time * 1e-9, + Hoex_mag=st.session_state.Hoex_mag * 1e3, + ) + plot_data(Hscan, freqs, spec, title="VSD spectrum") + + +def simulate_pimm(): + with st.spinner("Simulating PIMM..."): + spec, freqs, output, Hscan = get_pimm_data( + H_axis=st.session_state.H_axis, + Hmin=st.session_state.Hmin * 1e3, + Hmax=st.session_state.Hmax * 1e3, + Hsteps=st.session_state.Hsteps, + int_step=st.session_state.int_step, + sim_time=st.session_state.sim_time * 1e-9, + ) + mag = np.asarray(output["m_avg"]) + plot_data(Hscan, freqs, spec, mag=mag, title="PIMM spectrum") diff --git a/view/requirements.txt b/view/requirements.txt new file mode 100644 index 0000000..09cd4d5 --- /dev/null +++ b/view/requirements.txt @@ -0,0 +1,9 @@ +git+https://github.com/huawei-noah/HEBO.git#subdirectory=HEBO +git+https://github.com/lemurpwned/cmtj.git@cmtj-updates +matplotlib +networkx +numba +numpy +scipy +streamlit +tqdm diff --git a/view/simulation_fns.py b/view/simulation_fns.py new file mode 100644 index 0000000..606dc69 --- /dev/null +++ b/view/simulation_fns.py @@ -0,0 +1,209 @@ +from collections import defaultdict +from itertools import groupby +from typing import List + +import numpy as np +import streamlit as st + +from cmtj import * +from cmtj.models import LayerSB, Solver +from cmtj.utils import FieldScan, VectorObj, mu0 +from cmtj.utils.procedures import (PIMM_procedure, ResistanceParameters, + VSD_procedure) + + +def create_single_layer(id_: str) -> tuple: + """Do not forget to rescale the units!""" + demag = [CVector(0, 0, 0), CVector(0, 0, 0), CVector(0, 0, 1)] + Kdir1 = get_axis_cvector(st.session_state[f"anisotropy_axis{id_}"]) + layer = Layer( + id=f"layer_{id_}", + mag=Kdir1, + anis=Kdir1, + Ms=st.session_state[f"Ms{id_}"], + thickness=st.session_state[f"thickness{id_}"] * 1e-9, + cellSurface=1e-16, + demagTensor=demag, + damping=st.session_state[f"alpha{id_}"], + ) + layer.setAnisotropyDriver( + ScalarDriver.getConstantDriver(st.session_state[f"K{id_}"] * 1e3) + ) + rp = ResistanceParameters( + Rxx0=100, + Rxy0=1, + Rsmr=-0.46, + Rahe=-2.7, + Ramr=-0.24, + l=st.session_state[f"length{id_}"] * 1e-6, + w=st.session_state[f"width{id_}"] * 1e-6, + ) + return layer, rp + + +def get_axis_cvector(axis: str): + if axis == "x": + return CVector(1, 0, 0) + elif axis == "y": + return CVector(0, 1, 0) + elif axis == "z": + return CVector(0, 0, 1) + else: + raise ValueError(f"Invalid axis {axis}") + + +def get_axis(axis: str) -> Axis: + if axis == "x": + return Axis.xaxis + elif axis == "y": + return Axis.yaxis + elif axis == "z": + return Axis.zaxis + else: + raise ValueError(f"Invalid axis {axis}") + + +def get_axis_angles(axis: str): + if axis == "x": + return 0, 0 + elif axis == "y": + return 0, 90 + elif axis == "z": + return 90, 0 + else: + raise ValueError(f"Invalid axis {axis}") + + +def prepare_simulation(): + layers = [] + rparams = [] + N = st.session_state["N"] + for i in range(N): + layer, rp = create_single_layer(i) + layers.append(layer) + rparams.append(rp) + j = Junction(layers=layers) + for jvals in range(N - 1): + J = st.session_state[f"J{jvals}"] * 1e-3 # rescale GUI units + l1_name = layers[jvals].id + l2_name = layers[jvals + 1].id + j.setIECDriver(l1_name, l2_name, ScalarDriver.getConstantDriver(J)) + return j, rparams + + +# @st.cache_data +def get_pimm_data( + H_axis, + Hmin, + Hmax, + Hsteps, + int_step=1e-12, + sim_time=16e-9, +): + htheta, hphi = get_axis_angles(H_axis) + Hscan, Hvecs = FieldScan.amplitude_scan(Hmin, Hmax, Hsteps, htheta, hphi) + j, rparams = prepare_simulation() + spec, freqs, out = PIMM_procedure( + j, + Hvecs=Hvecs, + int_step=int_step, + resistance_params=rparams, + max_frequency=60e9, + simulation_duration=sim_time, + disturbance=1e-6, + wait_time=4e-9, + ) + return spec, freqs, out, Hscan + + +def get_vsd_data( + H_axis, + Hmin, + Hmax, + Hsteps, + Hoex, + fmin=0, + fmax=30e9, + fstep=0.5e9, + int_step=1e-12, + sim_time=16e-9, + Rtype="Ry", + Hoex_mag=500, +): + htheta, hphi = get_axis_angles(H_axis) + Hscan, Hvecs = FieldScan.amplitude_scan(Hmin, Hmax, Hsteps, htheta, hphi) + j, rparams = prepare_simulation() + frequencies = np.arange(fmin, fmax, step=fstep) + spec = VSD_procedure( + j, + Hvecs=Hvecs, + int_step=int_step, + frequencies=frequencies, + resistance_params=rparams, + simulation_duration=sim_time, + Rtype=Rtype, + Hoe_excitation=Hoex_mag, + Hoe_direction=get_axis(Hoex), + disturbance=1e-6, + ) + return spec, frequencies, Hscan + + +def compute_sb_mse(target, data): + # if there are multiple layers, and target is a single list + # we need to take the closest resonance + if len(data["frequency"]) < len(target): + return float("inf") + + mse = 0 + for i, (_, f) in enumerate( + groupby(zip(data["Hmag"], data["frequency"]), key=lambda x: x[0]) + ): + # f is a list of tuples (Hmag, frequency) + contrib = min(((target[i] / 1e9 - f_[1]) ** 2).sum() for f_ in f) + mse += contrib + + return mse + + +def simulate_sb(hvals: List[float]): + layers = [] + N = st.session_state.N + init_pos = [] + for i in range(N): + ktheta, kphi = get_axis_angles(st.session_state[f"anisotropy_axis{i}"]) + Kval = st.session_state[f"K{i}"] * 1e3 # rescale GUI units + if ktheta == 0: + Ks = Kval + Kv = 10 + else: + Ks = 10 + Kv = Kval + layer = LayerSB( + _id=i, + thickness=st.session_state[f"thickness{i}"] * 1e-9, # rescale GUI units + Kv=VectorObj(np.deg2rad(0.0), np.deg2rad(0), Kv), + Ks=Ks, + Ms=st.session_state[f"Ms{i}"] / mu0, + # alpha=st.session_state[f"alpha{i}"], + ) + init_pos.extend([np.deg2rad(ktheta), np.deg2rad(kphi)]) + layers.append(layer) + Js = [st.session_state[f"J{i}"] * 1e-3 for i in range(N - 1)] + result_dictionary = defaultdict(list) + # we perform a sweep over the field magnitude + htheta, hphi = get_axis_angles(st.session_state.H_axis) + for i in range(len(hvals)): + solver = Solver( + layers=layers, + J1=Js, + J2=[0 for _ in Js], + H=VectorObj(theta=htheta, phi=hphi, mag=hvals[i]), + ) + eq, frequencies = solver.solve(init_position=init_pos, perturbation=1e-4) + for freq in frequencies: + result_dictionary["Hmag"].append(hvals[i] / 1e3) + result_dictionary["frequency"].append(freq) + init_pos = eq + + return result_dictionary diff --git a/view/streamlit_app.py b/view/streamlit_app.py new file mode 100644 index 0000000..e332f8c --- /dev/null +++ b/view/streamlit_app.py @@ -0,0 +1,200 @@ +# Authors: LemurPwned +from functools import partial + +import streamlit as st +from autofit import autofit +from helpers import simulate_pimm, simulate_vsd + +apptitle = "CMTJ simulator" + +st.set_page_config(page_title=apptitle, page_icon=":rocket:") +st.title(apptitle) +container = st.container() +N = container.number_input( + "Number of layers", min_value=1, max_value=10, value=1, key="N", format="%d" +) +container.markdown( + """ + ## Data Upload + If you want to upload data, to overlay it on the plot, please upload + a file with \t separated columns: H, f1, f2, ... Put H in (A/m) and f in (Hz). + They will be rescaled to (kA/m) and (GHz) automatically. + """ +) +container.file_uploader( + "Upload your data here", + help="Upload your data here. Must be `\t` separated values and have H and f headers.", + type=["txt", "dat"], + accept_multiple_files=False, + key="upload", +) + +with st.sidebar: + st.markdown("## Simulation parameters") + st.markdown("### Layer parameters") + for i in range(N): + st.markdown(f"#### Layer {i+1}") + st.slider( + f"Ms ({i+1}) (T)", + min_value=0.2, + max_value=2.0, + value=0.52, + step=0.01, + key=f"Ms{i}", + ) + st.number_input( + f"K ({i+1}) (kJ/m^3)", + min_value=0.1, + max_value=10e3, + value=150.0, + step=10.0, + key=f"K{i}", + ) + st.number_input( + f"alpha ({i+1})", + min_value=1e-3, + max_value=0.1, + value=1e-3, + key=f"alpha{i}", + format="%.3f", + ) + st.number_input( + f"thickness ({i+1}) (nm)", + min_value=1.0, + max_value=10.0, + value=1.0, + key=f"thickness{i}", + ) + st.number_input( + f"width ({i+1}) (um)", + min_value=1.0, + max_value=500.0, + value=10.0, + key=f"width{i}", + ) + st.number_input( + f"length ({i+1}) (um)", + min_value=1.0, + max_value=500.0, + value=10.0, + key=f"length{i}", + ) + st.radio( + f"anisotropy axis ({i+1})", + options=["x", "y", "z"], + key=f"anisotropy_axis{i}", + index=2, + ) + st.markdown("-----\n") + + st.markdown("### Interlayer parameters") + for j in range(N - 1): + st.number_input( + f"J ({j+1}<-->{j+2}) (mJ/m^2)", + min_value=-1.0, + max_value=1.0, + value=0.0, + key=f"J{j}", + format="%.2f", + ) + st.markdown("-----\n") + st.markdown("## Control parameters") + st.markdown("### External field") + st.selectbox("H axis", options=["x", "y", "z"], key="H_axis", index=2) + st.number_input( + "Hmin (kA/m)", min_value=-1000.0, max_value=1000.0, value=-400.0, key="Hmin" + ) + st.number_input( + "Hmax (kA/m)", min_value=0.0, max_value=1000.0, value=400.0, key="Hmax" + ) + st.number_input( + "H steps", min_value=1, max_value=1000, value=50, key="Hsteps", format="%d" + ) + st.number_input( + "int_step", + min_value=1e-14, + max_value=1e-12, + value=1e-12, + key="int_step", + format="%.1e", + ) + st.number_input( + "sim_time (ns)", + min_value=1, + max_value=500, + value=16, + key="sim_time", + format="%d", + ) + + +pimm_tab, vsd_tab, opt_tab = st.tabs(["PIMM", "VSD", "Optimization"]) +with opt_tab: + st.markdown( + """ + ## Optimization + + Run Bayesian optimisation -- fitting data is source from file upload. + Select the number of iterations. + The only optimised values are: Ms, and K. + All other parameters are treated as constants. + """ + ) + st.number_input( + "Number of iterations", + min_value=1, + max_value=100, + value=10, + key="n_iters", + format="%d", + ) + placeholder = st.empty() + st.button( + "Run optimization", + on_click=partial(autofit, placeholder=placeholder), + key="opt_btn", + ) + +with vsd_tab: + st.number_input( + "Frequency min (GHz)", min_value=0.0, max_value=50.0, value=0.0, key="fmin" + ) + st.number_input( + "Frequency max (GHz)", min_value=0.0, max_value=50.0, value=20.0, key="fmax" + ) + st.number_input( + "Frequency step (GHz)", + min_value=0.1, + max_value=10.0, + value=0.5, + key="fstep", + step=0.1, + ) + st.number_input( + "Excitation (kA/m)", min_value=0.5, max_value=50.0, value=5.0, key="Hoex_mag" + ) + st.selectbox( + "Resistance type", + options=["Rx", "Ry"], + key="res_type", + index=1, + ) + st.selectbox("excitation axis", options=["x", "y", "z"], key="Hoex", index=1) + st.markdown( + """### Simulation info + + This is Voltage Spin Diode experiment (WIP). + """ + ) + + st.button("Simulate VSD", on_click=simulate_vsd, key="VSD_btn") + +with pimm_tab: + fn = simulate_pimm + st.markdown( + """### Simulation info + + This app simulates PIMM experiment. + """ + ) + st.button("Simulate PIMM", on_click=simulate_pimm, key="PIMM_btn") diff --git a/view/utils.py b/view/utils.py new file mode 100644 index 0000000..e9fa515 --- /dev/null +++ b/view/utils.py @@ -0,0 +1,34 @@ +from typing import List + +import numpy as np + + +def extract_max_resonance_lines( + spectrum, + h_vals: List[float], + frequencies: List[float], + N_layers: int, + mask_size_x=0, + mask_size_y=0, +): + masked_spectrum = spectrum[mask_size_x:, mask_size_y:] + frequencies = frequencies[mask_size_y:].squeeze() + h_vals = h_vals[mask_size_x:] + linespectra = {} + for i in range(len(h_vals)): + grad = np.gradient(masked_spectrum[i, :]) + sign = np.sign(grad) + diff = np.diff(sign) + indices = np.where(diff != 0)[0] + grad = grad / np.max(np.abs(grad)) + + # take two largest values at the indices + amp_vals = masked_spectrum[i, indices] + amp_max_indx = np.argsort(amp_vals) + indices = indices[amp_max_indx] + + indices = indices[-N_layers:] + n_diff = N_layers - len(indices) + freqs = frequencies[indices].ravel() + [None] * n_diff + linespectra[h_vals[i]] = freqs + return linespectra