Skip to content
Closed
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension


Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion .github/workflows/test-notebooks.yml
Original file line number Diff line number Diff line change
Expand Up @@ -42,7 +42,7 @@ jobs:
done

- name: Run notebook tests
run: pytest --nbval-lax --current-env tests/
run: pytest --nbval-lax --nbval-current-env tests/

- name: Upload executed notebooks on failure
if: failure()
Expand Down
3 changes: 3 additions & 0 deletions .gitignore
Original file line number Diff line number Diff line change
Expand Up @@ -30,5 +30,8 @@ __pycache__
*.Identifier
!requirements.txt
.jupyter_cache

# Ignore log files
tests/**/*.log
tests/multiprocess_log.txt
.vscode
5 changes: 2 additions & 3 deletions pyAMARES/__init__.py
Original file line number Diff line number Diff line change
@@ -1,8 +1,7 @@
__author__ = "Jia Xu, MR Research Facility, University of Iowa"
__version__ = "0.3.33"
__version__ = "0.3.34dev"

# print("Current pyAMARES version is %s" % __version__)
# print("Author: %s" % __author__)
# print(f"Author: {__author__)}"

from .fileio import * # noqa: F403
from .kernel import * # noqa: F403
Expand Down
10 changes: 3 additions & 7 deletions pyAMARES/fileio/readfidall.py
Original file line number Diff line number Diff line change
Expand Up @@ -2,12 +2,9 @@

import mat73
import numpy as np
from loguru import logger
from scipy import io

from ..libs.logger import get_logger

logger = get_logger(__name__)


def is_mat_file_v7_3(filename):
with open(filename, "rb") as f:
Expand Down Expand Up @@ -132,10 +129,9 @@ def read_fidall(filename):

if len(data.shape) != 1:
logger.warning(
"Note pyAMARES.fitAMARES only fits 1D MRS data, however, your data shape is {data.shape}. Is it MRSI or raw MRS data that needs to be coil-combined?"
f"Note pyAMARES.fitAMARES only fits 1D MRS data, however, your data shape is {data.shape}. Is it MRSI or raw MRS data that needs to be coil-combined?"
)

# print("data.shape=", data.shape)
logger.debug("data.shape=%s", data.shape)
logger.debug(f"data.shape={data.shape}")

return header, data
12 changes: 2 additions & 10 deletions pyAMARES/fileio/readmat.py
Original file line number Diff line number Diff line change
@@ -1,12 +1,10 @@
import mat73
import numpy as np
from loguru import logger
from scipy import io

from ..libs.logger import get_logger
from .readfidall import is_mat_file_v7_3

logger = get_logger(__name__)


def readmrs(filename):
"""
Expand Down Expand Up @@ -43,28 +41,23 @@ def readmrs(filename):
- For MATLAB files, both traditional (.mat) and V7.3 (.mat) files are supported, but the variable must be named ``fid`` or ``data``.
"""
if filename.endswith("csv"):
# print("Try to load 2-column CSV")
logger.debug("Try to load 2-column CSV")
data = np.loadtxt(filename, delimiter=",")
data = data[:, 0] + 1j * data[:, 1]
elif filename.endswith("txt"):
# print("Try to load 2-column ASCII data")
logger.debug("Try to load 2-column ASCII data")
data = np.loadtxt(filename, delimiter=" ")
data = data[:, 0] + 1j * data[:, 1]
elif filename.endswith("npy"):
# print("Try to load python NPY file")
logger.debug("Try to load python NPY file")
data = np.load(filename)
elif filename.endswith("mat"):
if is_mat_file_v7_3(filename):
# print("Try to load Matlab V7.3 mat file with the var saved as fid or data")
logger.debug(
"Try to load Matlab V7.3 mat file with the var saved as fid or data"
)
matdic = mat73.loadmat(filename)
else:
# print("Try to load Matlab mat file with the var saved as fid or data")
logger.debug(
"Try to load Matlab mat file with the var saved as fid or data"
)
Expand All @@ -87,6 +80,5 @@ def readmrs(filename):
"Note pyAMARES.fitAMARES only fits 1D MRS data, however, your data shape is {data.shape}. Is it MRSI or raw MRS data that needs to be coil-combined?"
)

# print("data.shape=", data.shape)
logger.debug("data.shape=%s", data.shape)
logger.debug(f"data.shape={data.shape}")
return data
5 changes: 1 addition & 4 deletions pyAMARES/fileio/readnifti.py
Original file line number Diff line number Diff line change
@@ -1,10 +1,7 @@
import argparse

import numpy as np

from ..libs.logger import get_logger

logger = get_logger(__name__)
from loguru import logger


def read_nifti(filename):
Expand Down
98 changes: 41 additions & 57 deletions pyAMARES/kernel/PriorKnowledge.py
Original file line number Diff line number Diff line change
Expand Up @@ -5,13 +5,16 @@
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

# Disable automatic string dtype inference for pandas 3.0 and above
if hasattr(pd.options, "future") and hasattr(pd.options.future, "infer_string"):
pd.options.future.infer_string = False

from lmfit import Parameters
from loguru import logger

from ..libs.logger import get_logger
from .fid import fft_params

logger = get_logger(__name__)


def safe_convert_to_numeric(x):
try:
Expand Down Expand Up @@ -181,6 +184,9 @@ def unitconverter(df_ini, MHz=120.0):
pandas.DataFrame: A DataFrame with converted unit values in specified rows.
"""
df = deepcopy(df_ini)
df = df.apply(
pd.to_numeric, errors="raise", downcast="float"
) # By this point the values should only be numeric
if "chemicalshift" in df.index:
df.loc["chemicalshift", df.notna().loc["chemicalshift"]] *= MHz

Expand All @@ -189,7 +195,7 @@ def unitconverter(df_ini, MHz=120.0):

if "phase" in df.index:
df.loc["phase", df.notna().loc["phase"]] = np.deg2rad(
df.loc["phase"][df.notna().loc["phase"]].astype(float)
df.loc["phase"][df.notna().loc["phase"]]
)

return df
Expand Down Expand Up @@ -262,30 +268,27 @@ def assert_peak_format(input_str):
if re.search(r"\.\d+$", input_str):
logger.error(msg)
raise ValueError(
"The peak name %s cannot end with a floating-point number!" % input_str
f"The peak name {input_str} cannot end with a floating-point number!"
)
if re.search(r"\d+[\D]", input_str) or re.search(r"^\d+", input_str):
logger.error(msg)
raise ValueError(
"The peak name %s cannot contain numbers at the beginning or in the middle!"
% input_str
f"The peak name {input_str} cannot contain numbers at the beginning or in the middle!"
)


def find_header_row(filename, comment_char="#"):
"""Determine the index of the first non-commented line."""
with open(filename, "r") as file:
logger.info("Checking comment lines in the prior knowledge file")
logger.debug("Checking comment lines in the prior knowledge file")
for i, line in enumerate(file):
if "#" in line:
logger.info("Comment: in line %d: %s", i, line)
logger.debug(f"Comment: in line {i}: {line}")
with open(filename, "r") as file:
for i, line in enumerate(file):
processedline = line.replace('"', "").replace("'", "").strip()
if not processedline.startswith(comment_char):
return i
# else:
# print("Comment:", processedline)
return None # Return None if all lines are comments or file is empty


Expand Down Expand Up @@ -326,41 +329,36 @@ def generateparameter(
else:
raise NotImplementedError("file format must be Excel (xlsx) or CSV!")

# Compatible with pandas both older and newer than 2.1.0
pk = (
pk.map(safe_convert_to_numeric)
if hasattr(pk, "map")
else pk.applymap(safe_convert_to_numeric)
) # To be compatible with CSV
def backward_compatible_map(df, func):
# Check if the newer 'map' method exists (pandas >= 2.1.0)
if hasattr(df, "map"):
return df.map(func)
# Fallback to the older 'applymap' (pandas < 2.1.0)
elif hasattr(df, "applymap"):
return df.applymap(func) # type: ignore
else:
raise AttributeError(
"Pandas DataFrame has neither 'map' nor 'applymap' method."
)

pk = backward_compatible_map(pk, safe_convert_to_numeric)
peaklist = pk.columns.to_list() # generate a peak list directly from the
[assert_peak_format(x) for x in peaklist]
dfini = extractini(pk, MHz=MHz) # Parse initial values
dfini2 = unitconverter(
dfini, MHz=MHz
) # Convert ppm to Hz, convert FWHM to dk, convert degree to radians.
# print(f"{dfini2=}")
df_lb, df_ub = parse_bounds(pk) # Parse bounds
df_expr = extract_expr(pk, MHz=MHz) # Parse expression
df_lb2 = unitconverter(df_lb, MHz=MHz)
df_ub2 = unitconverter(df_ub, MHz=MHz)
# Make sure the bounds are numeric
# Compatible with pandas both older and newer than 2.1.0
df_lb2 = (
df_lb2.map(safe_convert_to_numeric)
if hasattr(df_lb2, "map")
else df_lb2.applymap(safe_convert_to_numeric)
)
df_ub2 = (
df_ub2.map(safe_convert_to_numeric)
if hasattr(df_ub2, "map")
else df_ub2.applymap(safe_convert_to_numeric)
)
df_lb2 = backward_compatible_map(df_lb2, safe_convert_to_numeric)
df_ub2 = backward_compatible_map(df_ub2, safe_convert_to_numeric)
if g_global is False:
logger.debug(
"Parameter g will be fit with the initial value set in the file %s" % fname
f"Parameter g will be fit with the initial value set in the file {fname}"
)
# print(f"Parameter g will be fit with the initial value set in the file {fname}")
allpara = Parameters()
for peak in dfini2.columns:
for i, para in enumerate(paramter_prefix):
Expand Down Expand Up @@ -492,27 +490,21 @@ def initialize_FID(
deadtime = float(deadtime)
dwelltime = 1.0 / sw
if truncate_initial_points > 0:
logger.info(
"Truncating %i points from the beginning of the FID signal"
% truncate_initial_points
logger.debug(
f"Truncating {truncate_initial_points} points from the beginning of the FID signal"
)
deadtime_old = deadtime * 1.0
deadtime = deadtime + truncate_initial_points * dwelltime
fid = fid[truncate_initial_points:]
logger.info(
"The deadtime is changing from %f seconds to %f seconds"
% (deadtime_old, deadtime)
logger.debug(
f"The deadtime is changing from {deadtime} seconds to {deadtime_old} seconds"
)
fidpt = len(fid)
# TD = fidpt * 2
# at = TD / (2 * sw)

ppm = np.linspace(-sw / 2, sw / 2, fidpt) / np.abs(MHz)
Hz = np.linspace(-sw / 2, sw / 2, fidpt)
# print(f"{sw=}")
# print(f"{np.max(ppm)=} {np.min(ppm)=}")
# print(f"{np.max(Hz)=} {np.min(Hz)=}")
# print(f"{-sw/2=}")

opts = argparse.Namespace()
opts.deadtime = deadtime
Expand All @@ -523,7 +515,7 @@ def initialize_FID(
# This must be done before the shifting FID for carrier.
fid = np.conj(fid)
if carrier != 0:
logger.info("Shift FID so that center frequency is at %s ppm!" % carrier)
logger.debug(f"Shift FID so that center frequency is at {carrier} ppm!")
fid = fid * np.exp(1j * 2 * np.pi * carrier * MHz * opts.timeaxis)
# ppm = ppm + carrier
# Hz = Hz + carrier / np.abs(MHz)
Expand Down Expand Up @@ -584,7 +576,7 @@ def initialize_FID(
timeaxis=opts.timeaxis, params=opts.initialParams, fid=True
)
if ppm_offset != 0:
logger.info("Shifting the ppm by ppm_offset=%2.2f ppm" % ppm_offset)
logger.debug(f"Shifting the ppm by ppm_offset={ppm_offset:.2f} ppm")
for p in opts.initialParams:
if p.startswith("freq"):
hz_offset = opts.ppm_offset * opts.MHz
Expand All @@ -597,25 +589,17 @@ def initialize_FID(
): # Check if there's an upper bound set
opts.initialParams[p].max += hz_offset
logger.debug(
"before opts.initialParams[%s].value=%s"
% (p, opts.initialParams[p].value)
f"before opts.initialParams[{p}].value={opts.initialParams[p].value}"
)
logger.debug(
"new value should be opts.initialParams[%s].value + opts.ppm_offset * opts.MHz=%s"
% (p, opts.initialParams[p].value + opts.ppm_offset * opts.MHz)
f"new value should be opts.initialParams[{p}].value + opts.ppm_offset * opts.MHz="
f"{opts.initialParams[p].value + opts.ppm_offset * opts.MHz}"
)

# print(f"before {opts.initialParams[p].value=}")
# print(
# f"new value should be {opts.initialParams[p].value + opts.ppm_offset * opts.MHz=}"
# )
opts.initialParams[p].value = (
opts.initialParams[p].value + hz_offset
)
# print(f"after {opts.initialParams[p].value=}")
logger.debug(
"after opts.initialParams[%s].value=%s"
% (p, opts.initialParams[p].value)
f"after opts.initialParams[{p}].value={opts.initialParams[p].value}"
)

opts.allpara = opts.initialParams # obsolete API, will be removed
Expand All @@ -635,12 +619,12 @@ def initialize_FID(
plt.xlabel("ppm")
plt.show()
if priorknowledgefile is not None:
logger.info("Printing the Prior Knowledge File %s" % priorknowledgefile)
logger.debug(f"Printing the Prior Knowledge File {priorknowledgefile}")
try:
from IPython.display import display

display(opts.PK_table) # display table
except ImportError:
logger.info(opts.PK_table)
logger.debug(opts.PK_table)

return opts
Loading