diff --git a/.github/workflows/ci_PR.yml b/.github/workflows/ci_PR.yml index 44f9b528..63de2eec 100644 --- a/.github/workflows/ci_PR.yml +++ b/.github/workflows/ci_PR.yml @@ -6,12 +6,13 @@ on: jobs: tests: + name: python test runs-on: ${{ matrix.os }} strategy: fail-fast: false matrix: os: [ubuntu-latest, macos-latest, windows-latest] - python-version: [3.9] + python-version: ['3.10','3.11'] steps: - uses: actions/checkout@v3 diff --git a/.github/workflows/ci_scheduled.yml b/.github/workflows/ci_scheduled.yml index 5c94e264..0659c65f 100644 --- a/.github/workflows/ci_scheduled.yml +++ b/.github/workflows/ci_scheduled.yml @@ -12,7 +12,7 @@ jobs: fail-fast: false matrix: os: [ubuntu-latest, macos-latest, windows-latest] - python-version: [3.8, 3.9, "3.10", "3.11"] + python-version: [3.8, 3.9, "3.10", "3.11", "3.12"] steps: - uses: actions/checkout@v3 diff --git a/.github/workflows/deploy_ghpages.yml b/.github/workflows/deploy_ghpages.yml index 7b89f15e..0e28a743 100644 --- a/.github/workflows/deploy_ghpages.yml +++ b/.github/workflows/deploy_ghpages.yml @@ -12,10 +12,10 @@ jobs: steps: - name: Checkout uses: actions/checkout@v1 - - name: Set up Python 3.9 + - name: Set up Python 3.10 uses: actions/setup-python@v1 with: - python-version: 3.9 + python-version: '3.10' - uses: actions/cache@v2 with: path: | diff --git a/.github/workflows/docs_PR.yml b/.github/workflows/docs_PR.yml index dbeaf7e6..2d5faea1 100644 --- a/.github/workflows/docs_PR.yml +++ b/.github/workflows/docs_PR.yml @@ -16,10 +16,10 @@ jobs: steps: - name: Checkout uses: actions/checkout@v1 - - name: Set up Python 3.9 + - name: Set up Python 3.10 uses: actions/setup-python@v1 with: - python-version: 3.9 + python-version: '3.10' - uses: actions/cache@v2 with: path: | diff --git a/.github/workflows/examples_PR.yml b/.github/workflows/examples_PR.yml index 34896559..e38c8609 100644 --- a/.github/workflows/examples_PR.yml +++ b/.github/workflows/examples_PR.yml @@ -16,10 +16,10 @@ jobs: steps: - name: Checkout uses: actions/checkout@v1 - - name: Set up Python 3.9 + - name: Set up Python 3.10 uses: actions/setup-python@v1 with: - python-version: 3.9 + python-version: '3.10' - uses: actions/cache@v2 with: path: | diff --git a/README.md b/README.md index 99d2d626..3ade50f1 100644 --- a/README.md +++ b/README.md @@ -16,7 +16,7 @@ The early versions of DeerLab (up to version 0.9.2) are written in MATLAB. The o ## Requirements -DeerLab is available for Windows, Mac and Linux systems and requires **Python 3.8**, **3.9**, or **3.10**. +DeerLab is available for Windows, Mac and Linux systems and requires **Python 3.8**, **3.9**, **3.10**, or **3.11**. All additional dependencies are automatically downloaded and installed during the setup. diff --git a/VERSION b/VERSION index 0ec25f75..570c7965 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v1.0.0 +v1.0.2 diff --git a/deerlab/model.py b/deerlab/model.py index e977fda9..70aa9f70 100644 --- a/deerlab/model.py +++ b/deerlab/model.py @@ -1569,11 +1569,11 @@ def fit(model_, y, *constants, par0=None, penalties=None, bootstrap=0, noiselvl= nnlsSolver : string, optional Solver used to solve a non-negative least-squares problem (if applicable): - * ``'qp'`` - Optimization of the NNLS problem using the ``quadprog`` package. + * ``'qp'`` - Optimization of the NNLS problem using the ``quadprog`` package. Only Python <= 3.10. * ``'cvx'`` - Optimization of the NNLS problem using the ``cvxopt`` package. * ``'fnnls'`` - Optimization using the fast NNLS algorithm. - The default is ``'qp'``. + The default is ``'cvx'``. verbose : scalar integer, optional Level of verbosity during the analysis: @@ -1754,7 +1754,7 @@ def bootstrap_fcn(ysim): # Dictionary of parameter names and fit uncertainties FitResult_paramuq = {f'{key}Uncert': model._getparamuq(fitresults.paramUncert,idx) for key,idx in zip(keys,param_idx)} # Dictionary of other fit quantities of interest - FitResult_dict = {key: getattr(fitresults,key) for key in ['param','paramUncert','model','cost','plot','residuals','stats','regparam']} + FitResult_dict = {key: getattr(fitresults,key) for key in ['param','paramUncert','model','cost','plot','residuals','stats','regparam','regparam_stats']} _paramlist = model._parameter_list('vector') # Prepare the propagate() and evaluate() methods diff --git a/deerlab/selregparam.py b/deerlab/selregparam.py index ac281d97..27a800cf 100644 --- a/deerlab/selregparam.py +++ b/deerlab/selregparam.py @@ -99,7 +99,7 @@ def selregparam(y, A, solver, method='aic', algorithm='brent', noiselvl=None, residuals : ndarray Values of the residual norms evaluated during the search. Returned if full_output is True. - residuals : ndarray + penalties : ndarray Values of the penalty norms evaluated during the search. Returned if full_output is True. """ diff --git a/deerlab/solvers.py b/deerlab/solvers.py index 845c56ef..0b8360a4 100644 --- a/deerlab/solvers.py +++ b/deerlab/solvers.py @@ -14,7 +14,6 @@ import time from functools import partial from copy import deepcopy -from quadprog import solve_qp def timestamp(): @@ -136,7 +135,6 @@ def _plot(ys,yfits,yuqs,noiselvls,axis=None,xlabel=None,gof=False,fontsize=13): # Adjust fontsize for ax in axs: for label in (ax.get_xticklabels() + ax.get_yticklabels()): - label.set_fontname('Calibri') label.set_fontsize(fontsize) return fig @@ -388,7 +386,7 @@ def _insertfrozen(parfit,parfrozen,frozen): # =========================================================================================== -def snlls(y, Amodel, par0=None, lb=None, ub=None, lbl=None, ubl=None, nnlsSolver='qp', reg='auto', weights=None, verbose=0, +def snlls(y, Amodel, par0=None, lb=None, ub=None, lbl=None, ubl=None, nnlsSolver='cvx', reg='auto', weights=None, verbose=0, regparam='aic', regparamrange=None, multistart=1, regop=None, alphareopt=1e-3, extrapenalty=None, subsets=None, ftol=1e-8, xtol=1e-8, max_nfev=1e8, lin_tol=1e-15, lin_maxiter=1e4, noiselvl=None, lin_frozen=None, mask=None, nonlin_frozen=None, uq=True, modeluq=False): @@ -482,11 +480,11 @@ def snlls(y, Amodel, par0=None, lb=None, ub=None, lbl=None, ubl=None, nnlsSolver nnlsSolver : string, optional Solver used to solve a non-negative least-squares problem (if applicable): - * ``'qp'`` - Optimization of the NNLS problem using the ``quadprog`` package. + * ``'qp'`` - Optimization of the NNLS problem using the ``quadprog`` package. Only Python <= 3.10. * ``'cvx'`` - Optimization of the NNLS problem using the ``cvxopt`` package. * ``'fnnls'`` - Optimization using the fast NNLS algorithm. - The default is ``'qp'``. + The default is ``'cvx'``. noiselvl : array_like, optional Noise standard deviation of the input signal(s), if not specified it is estimated automatically. @@ -714,10 +712,18 @@ def linear_problem(y,A,optimize_alpha,alpha): yfrozen = (A[:,lin_frozen]@lin_parfrozen[lin_frozen]).astype(float) # Optimiza the regularization parameter only if needed + + if optimize_alpha: - alpha = dl.selregparam((y-yfrozen)[mask], Ared[mask,:], linSolver, regparam, + output = dl.selregparam((y-yfrozen)[mask], Ared[mask,:], linSolver, regparam, weights=weights[mask], regop=L, candidates=regparamrange, - noiselvl=noiselvl,searchrange=regparamrange) + noiselvl=noiselvl,searchrange=regparamrange,full_output=True) + alpha = output[0] + alpha_stats['alphas_evaled'] = output[1] + alpha_stats['functional'] = output[2] + alpha_stats['residuals'] = output[3] + alpha_stats['penalties'] = output[4] + # Components for linear least-squares AtA, Aty = _lsqcomponents((y-yfrozen)[mask], Ared[mask,:], L, alpha, weights=weights[mask]) @@ -732,8 +738,8 @@ def linear_problem(y,A,optimize_alpha,alpha): # Insert back the frozen linear parameters xfit = _insertfrozen(xfit,lin_parfrozen,lin_frozen) - return xfit, alpha, Ndof + #=========================================================================== def ResidualsFcn(p): @@ -745,8 +751,10 @@ def ResidualsFcn(p): non-linear least-squares solver. """ - nonlocal par_prev, check, regparam_prev, xfit, alpha, Ndof, Ndof_lin - + nonlocal par_prev, check, regparam_prev, xfit, alpha, alpha_stats, Ndof, Ndof_lin + + + # Non-linear model evaluation A = Amodel(p) @@ -803,6 +811,8 @@ def ResidualsFcn(p): return res #=========================================================================== + alpha_stats = {'alphas_evaled':[],'functional':[],'residuals':[],'penalties':[]} + # ------------------------------------------------------------------- # Only linear parameters # ------------------------------------------------------------------- @@ -981,7 +991,7 @@ def ymodel(n): return FitResult(nonlin=nonlinfit, lin=linfit, param=parfit, model=modelfit, nonlinUncert=paramuq_nonlin, linUncert=paramuq_lin, paramUncert=paramuq, modelUncert=modelfituq, regparam=alpha, plot=plotfcn, - stats=stats, cost=fvals, residuals=res, noiselvl=noiselvl) + stats=stats, cost=fvals, residuals=res, noiselvl=noiselvl,regparam_stats=alpha_stats) # =========================================================================================== @@ -1232,6 +1242,15 @@ def qpnnls(AtA, Atb): Mathematical Programming, 27, 1-33. """ + + try: + from quadprog import solve_qp + except ModuleNotFoundError: + raise ModuleNotFoundError( + 'The quadprog package is required for this function.'+ + 'Install it with "pip install quadprog".') + + N = np.shape(AtA)[1] I = np.eye(N) lb = np.zeros(N) diff --git a/docsrc/source/changelog.rst b/docsrc/source/changelog.rst index 78d1fab2..006fe0fc 100644 --- a/docsrc/source/changelog.rst +++ b/docsrc/source/changelog.rst @@ -14,6 +14,35 @@ Release Notes - |fix| : Something which was not working as expected or leading to errors has been fixed. - |api| : This will require changes in your scripts or code. +Release ``v1.0.2`` - July 2023 +------------------------------------------ +- |fix| : Fixes errors in documentation (:pr:`429`). + + * Changes the file name of figures 'modelling*.png` to `modeling*.png`. To keep all spelling consistent with american english. + * Adds a missing `)` in `fitting_guide` + * Corrects the time axis in `echo_crossing` example. + +- |fix| : Fixes an errors in tests (:pr:`429`). + + * The test `test_algorithms` had an incorrect search range. + +- |fix| : Removes the default font from the `fit` function due to conflicts on some systems (:pr:`429`). + +Release ``v1.0.1`` - March 2023 +------------------------------------------ +- |fix| : Fixes some minor bugs in the documentation. + + * The file modelling_guide.rst is renamed to modeling_guide.rst to keep spelling consistency. + * The "Simulating a two-pathway 5-pulse DEER signal" and "Simulating a three-pathway 4-pulse DEER signal" examples now run correctly. + * +- |fix| : Fixes issues with CVXOPT in tests. + + * The testing will now use quadprog as the default solver. To account for the change tested values are now generated using the grid method + +- |fix| : Updates GitHub actions to use latest packages. +- |api| : Hardcodes out Python 3.11 support. This will remain until quadprog is fixed. +- |api| : Removes hard-wired RNG seeding + Release ``v1.0.0`` - December 2022 ------------------------------------------ diff --git a/docsrc/source/fitting_guide.rst b/docsrc/source/fitting_guide.rst index fc81f13e..483a5348 100644 --- a/docsrc/source/fitting_guide.rst +++ b/docsrc/source/fitting_guide.rst @@ -66,7 +66,7 @@ Fitting multi-dataset models Models merged using the ``merge`` function (see :ref:`here ` for details) can describe multiple datasets with a single model object and a common parameter set. To fit such a merged model to multiple datasets, the ``fit`` function can be used as above by passing a list of datasets ``[y1,y2,...,yN]`` instead of a single dataset :: # Fit the model to multiple datasets - result = dl.fit(model, [y1,y2,y3] + result = dl.fit(model, [y1,y2,y3]) The number of datasets must match the number of responses returned by the model. Additionally, the ordering in the list of datasets must match the order of responses from the model, i.e. ``response1`` of ``model`` will be fitted to ``y1``, and so on. diff --git a/docsrc/source/images/modelling_guide_graphics.svg b/docsrc/source/images/modeling_guide_graphics.svg similarity index 100% rename from docsrc/source/images/modelling_guide_graphics.svg rename to docsrc/source/images/modeling_guide_graphics.svg diff --git a/docsrc/source/images/modelling_guide_lincombine.png b/docsrc/source/images/modeling_guide_lincombine.png similarity index 100% rename from docsrc/source/images/modelling_guide_lincombine.png rename to docsrc/source/images/modeling_guide_lincombine.png diff --git a/docsrc/source/images/modelling_guide_lincombine1.png b/docsrc/source/images/modeling_guide_lincombine1.png similarity index 100% rename from docsrc/source/images/modelling_guide_lincombine1.png rename to docsrc/source/images/modeling_guide_lincombine1.png diff --git a/docsrc/source/images/modelling_guide_lincombine2.png b/docsrc/source/images/modeling_guide_lincombine2.png similarity index 100% rename from docsrc/source/images/modelling_guide_lincombine2.png rename to docsrc/source/images/modeling_guide_lincombine2.png diff --git a/docsrc/source/images/modelling_guide_link.png b/docsrc/source/images/modeling_guide_link.png similarity index 100% rename from docsrc/source/images/modelling_guide_link.png rename to docsrc/source/images/modeling_guide_link.png diff --git a/docsrc/source/images/modelling_guide_merge.png b/docsrc/source/images/modeling_guide_merge.png similarity index 100% rename from docsrc/source/images/modelling_guide_merge.png rename to docsrc/source/images/modeling_guide_merge.png diff --git a/docsrc/source/images/modelling_guide_relate.png b/docsrc/source/images/modeling_guide_relate.png similarity index 100% rename from docsrc/source/images/modelling_guide_relate.png rename to docsrc/source/images/modeling_guide_relate.png diff --git a/docsrc/source/installation.rst b/docsrc/source/installation.rst index 13832f9f..dd974a46 100644 --- a/docsrc/source/installation.rst +++ b/docsrc/source/installation.rst @@ -8,7 +8,7 @@ Requirements To install DeerLab, first install Python on your computer. Python can be downloaded from the `official Python distribution `_. There are many online tutorials to guide you through the installation and setup (see `here `_ for example). Make sure you install -one of the Python versions compatible with DeerLab, either **Python 3.8**, **3.9**, or **3.10**. +one of the Python versions compatible with DeerLab, either **Python 3.8**, **3.9**, **3.10**, or **3.11**. .. rubric:: Windows systems diff --git a/examples/intermediate/ex_crossing_echoes_masking.py b/examples/intermediate/ex_crossing_echoes_masking.py index de7796a9..db2cba66 100644 --- a/examples/intermediate/ex_crossing_echoes_masking.py +++ b/examples/intermediate/ex_crossing_echoes_masking.py @@ -27,6 +27,8 @@ # Load the experimental data t,Vexp = dl.deerload('../data/example_4pdeer_5.DTA') +t *= 1e3 # convert from ms to us + # Experimental parameters tau1 = 0.5 # First inter-pulse time delay, μs tau2 = 4.5 # Second inter-pulse time delay, μs diff --git a/examples/intermediate/ex_selregparam.py b/examples/intermediate/ex_selregparam.py new file mode 100644 index 00000000..7bd1b0e1 --- /dev/null +++ b/examples/intermediate/ex_selregparam.py @@ -0,0 +1,162 @@ +# %% [markdown] +""" +Analysing the selection of regularisation parameter +------------------------------------------------------------------------- + +This example demonstates how to generate plots from the regularisation parameter selection, +and how to use this infomation. + +""" + +import numpy as np +import matplotlib.pyplot as plt +import deerlab as dl + +# %% + +# File location +path = '../data/' +file = 'example_4pdeer_1.DTA' + +# Experimental parameters +tau1 = 0.3 # First inter-pulse delay, μs +tau2 = 4.0 # Second inter-pulse delay, μs +deadtime = 0.1 # Acquisition deadtime, μs + +# Load the experimental data +t,Vexp = dl.deerload(path + file) + +# Pre-processing +Vexp = dl.correctphase(Vexp) # Phase correction +Vexp = Vexp/np.max(Vexp) # Rescaling (aesthetic) +t = t + deadtime # Account for deadtime + +# Distance vector +r = np.linspace(1.5,7,50) # nm + +# Construct the model +Vmodel = dl.dipolarmodel(t,r, experiment=dl.ex_4pdeer(tau1,tau2, pathways=[1])) + +# Fit the model to the data with compactness criterion +results= dl.fit(Vmodel,Vexp,regparam='bic') +print(results) + +""" +The regularisation parameter in DeerLab can be selected using a variety of criteria. +The default criterion is the Akaike complexity criterion (aic) however other +criterion exists and can be selected. + +Each criterion has its own functional, which is minimised. These functionals +are often based on the residuals of the fit vs the raw data, such that a minimal functional value +will occur at the location of the best fit. Some methods such as the L-Curve-based methods do not follow this approach. + +Traditionally the L-Curve has been used to investigate and select the regularisation parameter. +The L-Curve is a plot of the Residual Norm against the Penalty Norm. Each point represents a +different regularisation parameter. Normally the optimal regularisation parameter can be found at the kink +of the curve, i.e. the place that has both a low Residual Norm and a low Pentalty Norm. +Recently, this approach has taken a back foot as the existence of an L-shape or kink is not guaranteed. +Nonetheless, it can be useful to diagnose problems in the selection of the regularisation parameter. + +""" +# %% + +fig, axs =plt.subplots(1,3, figsize=(9,4),width_ratios=(1,1,0.1)) +fig.tight_layout() +alphas = results.regparam_stats['alphas_evaled'][1:] +funcs = results.regparam_stats['functional'][1:] + + +idx = np.argsort(alphas) + +axs[0].semilogx(alphas[idx], funcs[idx],marker='.',ls='-') +axs[0].set_title(r"$\alpha$ selection functional"); +axs[0].set_xlabel("Regularisation Parameter") +axs[0].set_ylabel("Functional Value ") + +# Just the final L-Curve + +cmap = plt.get_cmap('plasma') +import matplotlib as mpl + +x = results.regparam_stats['residuals'] +y = results.regparam_stats['penalties'] +idx = np.argsort(x) + + +axs[1].loglog(x[idx],y[idx]) + +n_points = results.regparam_stats['alphas_evaled'].shape[-1] +lams = results.regparam_stats['alphas_evaled'] +norm = mpl.colors.LogNorm(vmin=lams[1:].min(), vmax=lams.max()) +for i in range(n_points): + axs[1].plot(x[i], y[i],marker = '.', ms=8, color=cmap(norm(lams[i]))) + +i_optimal = np.argmin(np.abs(lams - results.regparam)) +axs[1].annotate(fr"$\lambda =$ {results.regparam:.2g}", xy = (x[i_optimal],y[i_optimal]),arrowprops=dict(facecolor='black', shrink=0.05, width=5), xytext=(20, 20),textcoords='offset pixels') +fig.colorbar(mpl.cm.ScalarMappable(norm=norm, cmap=cmap),cax=axs[2]) +axs[1].set_ylabel("Penalties") +axs[2].set_ylabel("Regularisation Parameter") +axs[1].set_xlabel("Residuals") +axs[1].set_title("L-Curve"); + +# %% +""" +Over and Under selection of the regularisation parameter +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + +Here we will demonstrate the effect of selecting a regularisation parameter +that is either too small or too large. + + + +""" + +result_high= dl.fit(Vmodel,Vexp,regparam=1.0) + +result_low= dl.fit(Vmodel,Vexp,regparam=1e-4) + +green = '#3cb4c6' +red = '#f84862' +fig, axs =plt.subplots(1,2, figsize=(9,4),width_ratios=(1,1)) +fig.tight_layout() + +axs[0].set_xlabel("Time $t$ (μs)") +axs[0].set_ylabel('$V(t)$ (arb.u.)') +axs[0].plot(t,Vexp,'.',color='grey',label='Data') +axs[0].plot(t,result_high.model,linewidth=3,color=green,label='High regparam') +axs[0].plot(t,result_low.model,linewidth=3,color=red,label='Low regparam') +axs[0].legend(frameon=False,loc='best') + +Pfit_h = result_high.P +Pci95_h = result_high.PUncert.ci(95) + +Pfit_l = result_low.P +Pci95_l = result_low.PUncert.ci(95) + + + +axs[1].plot(r,Pfit_h,linewidth=3,color=green,label='High regparam') +axs[1].fill_between(r,Pci95_h[:,0],Pci95_h[:,1],alpha=0.3,color=green,linewidth=0) +axs[1].plot(r,Pfit_l,linewidth=3,color=red,label='Low regparam') +axs[1].fill_between(r,Pci95_l[:,0],Pci95_l[:,1],alpha=0.3,color=red,linewidth=0) +axs[1].set_ylim(0,max([Pfit_h.max(),Pfit_l.max()])) +axs[1].legend(frameon=False,loc='best') +axs[1].set_xlabel('Distance $r$ (nm)') +axs[1].set_ylabel('$P(r)$ (nm$^{-1}$)') + + +# %% +""" +As we can see when the regularisation parameter is too small we still get a high +quality fit in the time domain, however, our distance domain data is now way too +spikey and non-physical. + +In contrast when the regularisation parameter is too large we struggle to get +a good fit, however, we get a much smoother distance distribution. + +This could have been seen from the selection functional above. The effect of +lower regularisation parameter had a smaller effect on the functional than the +effect of going to a larger one. + + +""" diff --git a/setup.py b/setup.py index f7934a11..6a10e370 100644 --- a/setup.py +++ b/setup.py @@ -11,7 +11,7 @@ 'Documentation': 'https://jeschkelab.github.io/DeerLab/', 'Source': 'https://github.com/JeschkeLab/DeerLab', }, - python_requires='>=3.8,<3.11', + python_requires='>=3.8', license='LICENSE.txt', include_package_data = True, keywords='data analysis modeling least-squares EPR spectroscopy DEER PELDOR'.split(), @@ -30,10 +30,10 @@ 'pytest>=6.2.2', 'setuptools>=53.0.0', 'numexpr>=2.7.3', - 'quadprog>=0.1.11', + 'quadprog>=0.1.11; python_version <= "3.10"', ], classifiers=[ - 'Development Status :: 4 - Beta', + 'Development Status :: 5 - Production/Stable', 'Intended Audience :: Science/Research', 'License :: OSI Approved :: MIT License', 'Operating System :: Microsoft :: Windows', @@ -42,6 +42,7 @@ 'Programming Language :: Python :: 3.8', 'Programming Language :: Python :: 3.9', 'Programming Language :: Python :: 3.10', - 'Topic :: Scientific/Engineering', + 'Programming Language :: Python :: 3.11', + 'Topic :: Scientific/Engineering::Chemistry', ] ) diff --git a/test/test_selregparam.py b/test/test_selregparam.py index e50f70c1..30a64ee9 100644 --- a/test/test_selregparam.py +++ b/test/test_selregparam.py @@ -4,7 +4,7 @@ from deerlab import dipolarkernel, regoperator, selregparam, whitegaussnoise from deerlab.dd_models import dd_gauss,dd_gauss2 from deerlab.utils import assert_docstring -from deerlab.solvers import qpnnls +from deerlab.solvers import cvxnnls import pytest # Fixtures @@ -46,15 +46,15 @@ def dataset(design_matrix): ]) def test_selection_criteria_values(dataset, design_matrix, regularization_matrix, criterion,log_reference): "Check that the value returned by the selection criteria are correct" - alpha = selregparam(dataset, design_matrix, qpnnls, method=criterion, searchrange=[1e-10,1e-6], noiselvl=0, regop=regularization_matrix) - assert abs(1-np.log10(alpha)/log_reference) < 0.2 + alpha = selregparam(dataset, design_matrix, cvxnnls, method=criterion, searchrange=[1e-10,1e-6], noiselvl=0, regop=regularization_matrix) + assert abs(1-np.log10(alpha)/log_reference) < 0.25 #======================================================================= #======================================================================= def test_algorithms(dataset, design_matrix, regularization_matrix): "Check that the value returned by the the grid and Brent algorithms coincide" - alpha_grid = selregparam(dataset,design_matrix,qpnnls,method='aic',algorithm='grid',regop=regularization_matrix) - alpha_brent = selregparam(dataset,design_matrix,qpnnls,method='aic',algorithm='brent',regop=regularization_matrix) + alpha_grid = selregparam(dataset,design_matrix,cvxnnls,method='aic',algorithm='grid',regop=regularization_matrix, searchrange=[1e-9,1e-5]) + alpha_brent = selregparam(dataset,design_matrix,cvxnnls,method='aic',algorithm='brent',regop=regularization_matrix, searchrange=[1e-9,1e-5]) assert abs(1-alpha_grid/alpha_brent) < 0.2 #======================================================================= @@ -66,7 +66,7 @@ def test_nonuniform_r(dataset): K = dipolarkernel(t,r) L = regoperator(r,2) - logalpha = np.log10(selregparam(dataset,K,qpnnls,method='aic',regop=L)) + logalpha = np.log10(selregparam(dataset,K,cvxnnls,method='aic',regop=L)) logalpharef = -4.3325 assert abs(1 - logalpha/logalpharef) < 0.2 @@ -76,7 +76,7 @@ def test_nonuniform_r(dataset): @pytest.mark.parametrize('algorithm',['brent','grid']) def test_full_output_behavior_with_algorithms(dataset,design_matrix,regularization_matrix,algorithm): "Check that the full output argument works using all algorithms" - alpha,alphas_evaled,functional,residuals,penalties = selregparam(dataset,design_matrix,qpnnls,method='aic',algorithm=algorithm,full_output=True,regop=regularization_matrix) + alpha,alphas_evaled,functional,residuals,penalties = selregparam(dataset,design_matrix,cvxnnls,method='aic',algorithm=algorithm,full_output=True,regop=regularization_matrix) errors = [] if np.size(alpha)!=1: @@ -103,15 +103,15 @@ def test_unconstrained(dataset,design_matrix,regularization_matrix): def test_manual_candidates(dataset,design_matrix,regularization_matrix): "Check that the alpha-search range can be manually passed" alphas = np.linspace(-8,2,60) - alpha_manual = np.log10(selregparam(dataset,design_matrix,qpnnls,method='aic',candidates=alphas,regop=regularization_matrix)) - alpha_auto = np.log10(selregparam(dataset,design_matrix,qpnnls,method='aic',regop=regularization_matrix)) + alpha_manual = np.log10(selregparam(dataset,design_matrix,cvxnnls,method='aic',candidates=alphas,regop=regularization_matrix)) + alpha_auto = np.log10(selregparam(dataset,design_matrix,cvxnnls,method='aic',regop=regularization_matrix)) assert abs(alpha_manual-alpha_auto)<1e-4 #======================================================================= #======================================================================= def test_tikh_value(dataset,design_matrix,regularization_matrix): "Check that the value returned by Tikhonov regularization" - alpha = selregparam(dataset + whitegaussnoise(t,0.01,seed=1),design_matrix,qpnnls,method='aic',regop=regularization_matrix) + alpha = selregparam(dataset + whitegaussnoise(t,0.01,seed=1),design_matrix,cvxnnls,method='aic',regop=regularization_matrix) loga = np.log10(alpha) logaref = -3.667 # Computed with DeerLab-Matlab (1.0.0) assert abs(1-loga/logaref) < 0.02 # less than 2% error diff --git a/test/test_snlls_rlls.py b/test/test_snlls_rlls.py index 2c494e69..66c077a7 100644 --- a/test/test_snlls_rlls.py +++ b/test/test_snlls_rlls.py @@ -31,7 +31,14 @@ def mock_data(design_matrix, reference): # ============================================================ -@pytest.mark.parametrize('solver',['qp','fnnls','cvx']) + +try : + import quadprog + solvers = ['qp','fnnls','cvx'] +except ImportError: + solvers = ['fnnls','cvx'] + +@pytest.mark.parametrize('solver',solvers) def test_RLLS_fit_solvers(mock_data, design_matrix, reference, solver): "Check that the RLLS problem is correctly solved by all numerical solvers" fit = snlls(mock_data,design_matrix,lbl=lbl,nnlsSolver=solver,uq=False)