from doctest import master
from gettext import ngettext
import from doctest import master
from gettext import ngettext
import locale
import csdmpy as cp
import matplotlib.pyplot as plt
from lmfit import Minimizer
import nmrglue as ng
from mrsimulator import Site, SpinSystem, Simulator
from mrsimulator.method.lib import BlochDecaySpectrum
from mrsimulator import signal_processor as sp
from mrsimulator.utils import spectral_fitting as sf
from mrsimulator.utils import get_spectral_dimensions

# initialize nmrglue converter object 

ax = plt.subplot(projection="csdm")
# read in the JEOL dataset file

exp_spectrum = cp.load ("/Users/pdellis/Desktop/mrsim_fit/SnO_Transformed.csdf") 
fig, ax = plt.subplots(1, 2, figsize=(9, 3.5), subplot_kw={"projection": "csdm"})

ax[0].plot(exp_spectrum.real)
ax[0].set_title("Full Spectrum")
ax[0].grid()
ax[1].plot(exp_spectrum.real,color= "red", label="real")
ax[1].set_title("Zoomed Spectrum")
ax[1].set_xlim(1.5e6, 1.72e6)
ax[1].grid()
plt.tight_layout()
plt.legend()
plt.show()

spectral_dims = get_spectral_dimensions(exp_spectrum)
_ = [item.to("Hz", "nmr_frequency_ratio") for item in exp_spectrum.dimensions]

Tin_119 = Site(
    isotope="119Sn",
    isotropic_chemical_shift = 1.6e6,
    shielding_symmetric={"zeta":1.8e5,"eta": 0.0}, #zeta in Hz
)



Dict = {'origin_offset': 'Hz', 'reference_offset': 'Hz', 'spectral_width': 'Hz', 'reference_offset': 'Hz', 'origin_offset': 'Hz'},
spectral_dims = get_spectral_dimensions(exp_spectrum)
bloch = BlochDecaySpectrum(
    channels = ["119Sn"],
    magnetic_flux_density = 11.75,  # in T
    rotor_angle = 0.9553234193716164, # in rad
    rotor_frequency = 6600,  # in Hz  
    spectral_dimensions=spectral_dims,
    experiment = exp_spectrum,  # add the measurement to the method
    )

# Optimize the script by pre-setting the transition pathways for each spin system from
# the method.

# Use only the real part of the spectrum
exp_spectrum = exp_spectrum.real

sys = [SpinSystem(sites=[Tin_119]), ]
sim = Simulator(spin_systems = sys, methods=[bloch])
sim.config.decompose_spectrum = "spin_system"
sim.run()

relative_intensity_factor = exp_spectrum.max() / sim.methods[0].simulation.max()
# Post Simulation Processing

processor = sp.SignalProcessor(
    operations=[
        sp.IFFT(),
        sp.apodization.Gaussian(FWHM="50 Hz"),
        sp.FFT(),        
        sp.Scale(factor=relative_intensity_factor)
    ]
)
_ = [item.to("Hz", "nmr_frequency_ratio") for item in exp_spectrum.dimensions]
dataset = exp_spectrum
processed_dataset = processor.apply_operations(dataset = sim.methods[0].simulation).real

# Plot the spectrum
plt.figure(figsize=(4.25, 3.0))
ax = plt.subplot(projection="csdm")
ax.plot(exp_spectrum, linewidth=1, label="Experiment")
ax.plot(processed_dataset, linewidth=1, label="Simulation")
ax.set_xlim(1.5e6, 1.72e6)
plt.grid()
plt.legend()
plt.tight_layout()
plt.show()


params = sf.make_LMFIT_params(sim, processor, include={"rotor_frequency"})

fit_parameters = sf.make_LMFIT_params(sim, processor)
print(fit_parameters.pretty_print(columns=["value", "min", "max", "vary", "expr"]))

# standard deviation of noise from the dataset
sigma = 1.523217
from lmfit import Minimizer
minner = Minimizer(sf.LMFIT_min_function, fit_parameters, fcn_args=(sim, processor, sigma))
result = minner.minimize()
result

best_fit = sf.bestfit(sim, processor)[0].real
residuals = sf.residuals(sim, processor)[0].real

# Plot the spectrum
plt.figure(figsize=(4.25, 3.0))
ax = plt.subplot(projection="csdm")
ax.plot(exp_spectrum, linewidth=1, label="Experiment")
ax.plot(best_fit, linewidth=1, label="Best Fit")
ax.plot(residuals, linewidth=1, label="Residuals")
ax.set_xlim(1.5e6, 1.72e6)
plt.grid()
plt.legend()
plt.tight_layout()
plt.show()

import csdmpy as cp
import matplotlib.pyplot as plt
from lmfit import Minimizer
import nmrglue as ng
from mrsimulator import Site, SpinSystem, Simulator
from mrsimulator.method.lib import BlochDecaySpectrum
from mrsimulator import signal_processor as sp
from mrsimulator.utils import spectral_fitting as sf
from mrsimulator.utils import get_spectral_dimensions

# initialize nmrglue converter object 

ax = plt.subplot(projection="csdm")
# read in the JEOL dataset file

exp_spectrum = cp.load ("/Users/pdellis/Desktop/mrsim_fit/SnO_Transformed.csdf") 
fig, ax = plt.subplots(1, 2, figsize=(9, 3.5), subplot_kw={"projection": "csdm"})

ax[0].plot(exp_spectrum.real)
ax[0].set_title("Full Spectrum")
ax[0].grid()
ax[1].plot(exp_spectrum.real,color= "red", label="real")
ax[1].set_title("Zoomed Spectrum")
ax[1].set_xlim(1.5e6, 1.72e6)
ax[1].grid()
plt.tight_layout()
plt.legend()
plt.show()

spectral_dims = get_spectral_dimensions(exp_spectrum)
_ = [item.to("Hz", "nmr_frequency_ratio") for item in exp_spectrum.dimensions]

Tin_119 = Site(
    isotope="119Sn",
    isotropic_chemical_shift = 1.6e6,
    shielding_symmetric={"zeta":1.8e5,"eta": 0.0}, #zeta in Hz
)



Dict = {'origin_offset': 'Hz', 'reference_offset': 'Hz', 'spectral_width': 'Hz', 'reference_offset': 'Hz', 'origin_offset': 'Hz'},
spectral_dims = get_spectral_dimensions(exp_spectrum)
bloch = BlochDecaySpectrum(
    channels = ["119Sn"],
    magnetic_flux_density = 11.75,  # in T
    rotor_angle = 0.9553234193716164, # in rad
    rotor_frequency = 6600,  # in Hz  
    spectral_dimensions=spectral_dims,
    experiment = exp_spectrum,  # add the measurement to the method
    )

# Optimize the script by pre-setting the transition pathways for each spin system from
# the method.

# Use only the real part of the spectrum
exp_spectrum = exp_spectrum.real

sys = [SpinSystem(sites=[Tin_119]), ]
sim = Simulator(spin_systems = sys, methods=[bloch])
sim.config.decompose_spectrum = "spin_system"
sim.run()

relative_intensity_factor = exp_spectrum.max() / sim.methods[0].simulation.max()
# Post Simulation Processing

processor = sp.SignalProcessor(
    operations=[
        sp.IFFT(),
        sp.apodization.Gaussian(FWHM="50 Hz"),
        sp.FFT(),        
        sp.Scale(factor=relative_intensity_factor)
    ]
)
_ = [item.to("Hz", "nmr_frequency_ratio") for item in exp_spectrum.dimensions]
dataset = exp_spectrum
processed_dataset = processor.apply_operations(dataset = sim.methods[0].simulation).real

# Plot the spectrum
plt.figure(figsize=(4.25, 3.0))
ax = plt.subplot(projection="csdm")
ax.plot(exp_spectrum, linewidth=1, label="Experiment")
ax.plot(processed_dataset, linewidth=1, label="Simulation")
ax.set_xlim(1.5e6, 1.72e6)
plt.grid()
plt.legend()
plt.tight_layout()
plt.show()


params = sf.make_LMFIT_params(sim, processor, include={"rotor_frequency"})

fit_parameters = sf.make_LMFIT_params(sim, processor)
print(fit_parameters.pretty_print(columns=["value", "min", "max", "vary", "expr"]))

# standard deviation of noise from the dataset
sigma = 1.523217
from lmfit import Minimizer
minner = Minimizer(sf.LMFIT_min_function, fit_parameters, fcn_args=(sim, processor, sigma))
result = minner.minimize()
result

best_fit = sf.bestfit(sim, processor)[0].real
residuals = sf.residuals(sim, processor)[0].real

# Plot the spectrum
plt.figure(figsize=(4.25, 3.0))
ax = plt.subplot(projection="csdm")
ax.plot(exp_spectrum, linewidth=1, label="Experiment")
ax.plot(best_fit, linewidth=1, label="Best Fit")
ax.plot(residuals, linewidth=1, label="Residuals")
ax.set_xlim(1.5e6, 1.72e6)
plt.grid()
plt.legend()
plt.tight_layout()
plt.show()

/usr/local/bin/python3 /Users/pdellis/Desktop/mrsim_fit/SnO_Fitting
pdellis@dotynmr ~ % /usr/local/bin/python3 /Users/pdellis/Desktop/mrsim_fit/SnO_Fitting
/Library/Frameworks/Python.framework/Versions/3.11/lib/python3.11/site-packages/astropy/units/quantity.py:673: RuntimeWarning: divide by zero encountered in divide
  result = super().__array_ufunc__(function, method, *arrays, **kwargs)
/Library/Frameworks/Python.framework/Versions/3.11/lib/python3.11/site-packages/astropy/units/quantity.py:673: RuntimeWarning: invalid value encountered in divide
  result = super().__array_ufunc__(function, method, *arrays, **kwargs)
/Library/Frameworks/Python.framework/Versions/3.11/lib/python3.11/site-packages/astropy/units/quantity.py:1327: ComplexWarning: Casting complex values to real discards the imaginary part
  return float(self.to_value(dimensionless_unscaled))
/Library/Frameworks/Python.framework/Versions/3.11/lib/python3.11/site-packages/csdmpy/csdm.py:411: RuntimeWarning: invalid value encountered in multiply
  item.components *= other
Name                                      Value      Min      Max     Vary     Expr
SP_0_operation_1_Gaussian_FWHM               50     -inf      inf     True     None
SP_0_operation_3_Scale_factor               inf     -inf      inf     True     None
sys_0_abundance                             100        0      100    False      100
sys_0_site_0_isotropic_chemical_shift   1.6e+06     -inf      inf     True     None
sys_0_site_0_shielding_symmetric_eta          0        0        1     True     None
sys_0_site_0_shielding_symmetric_zeta   1.8e+05     -inf      inf     True     None
None
Traceback (most recent call last):
  File "/Users/pdellis/Desktop/mrsim_fit/SnO_Fitting", line 102, in <module>
    result = minner.minimize()
             ^^^^^^^^^^^^^^^^^
  File "/Library/Frameworks/Python.framework/Versions/3.11/lib/python3.11/site-packages/lmfit/minimizer.py", line 2379, in minimize
    return function(**kwargs)
           ^^^^^^^^^^^^^^^^^^
  File "/Library/Frameworks/Python.framework/Versions/3.11/lib/python3.11/site-packages/lmfit/minimizer.py", line 1703, in leastsq
    lsout = scipy_leastsq(self.__residual, variables, **lskws)
            ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/Library/Frameworks/Python.framework/Versions/3.11/lib/python3.11/site-packages/scipy/optimize/_minpack_py.py", line 413, in leastsq
    shape, dtype = _check_func('leastsq', 'func', func, x0, args, n)
                   ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/Library/Frameworks/Python.framework/Versions/3.11/lib/python3.11/site-packages/scipy/optimize/_minpack_py.py", line 26, in _check_func
    res = atleast_1d(thefunc(*((x0[:numinputs],) + args)))
                     ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/Library/Frameworks/Python.framework/Versions/3.11/lib/python3.11/site-packages/lmfit/minimizer.py", line 604, in __residual
    return _nan_policy(np.asfarray(out).ravel(),
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/Library/Frameworks/Python.framework/Versions/3.11/lib/python3.11/site-packages/lmfit/minimizer.py", line 2463, in _nan_policy
    raise ValueError(msg)
ValueError: NaN values detected in your input data or the output of your objective/model function - fitting algorithms cannot handle this! Please read https://lmfit.github.io/lmfit-py/faq.html#i-get-errors-from-nan-in-my-fit-what-can-i-do for more information.
pdellis@dotynmr ~ %