diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index d656046..e6d5cc4 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -26,7 +26,7 @@ jobs: run: | python -m pip install --upgrade pip python -m pip install flake8 pytest coverage - if [ -f requirements.txt ]; then pip install -r requirements.txt; fi + python -m pip install . - name: Lint with flake8 run: | # stop the build if there are Python syntax errors or undefined names diff --git a/.github/workflows/version.yml b/.github/workflows/version.yml new file mode 100644 index 0000000..f3c2eae --- /dev/null +++ b/.github/workflows/version.yml @@ -0,0 +1,18 @@ +name: Version_Checks + +on: + pull_request: + branches: [master] + +jobs: + version-checks: + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v3 + - uses: actions/setup-python@v4 + - name: Upgrade pip and install linters + run: | + python -m pip install --upgrade pip + python -m pip install packaging + - name: Check version number + run: python ./bin/check_version.py \ No newline at end of file diff --git a/Basin-hopping_Nelder-Mead/BHNM/msf_fit.py b/Basin-hopping_Nelder-Mead/BHNM/msf_fit.py deleted file mode 100644 index 967239c..0000000 --- a/Basin-hopping_Nelder-Mead/BHNM/msf_fit.py +++ /dev/null @@ -1,255 +0,0 @@ -import numpy as np -from scipy.optimize import minimize, basinhopping -import sys - - -class max_fit_BH(object): - def __init__(self, x, y, N, mid_point, **kwargs): - - self.step_size = kwargs.pop('step_size', 0.5) - self.temp = kwargs.pop('temp', 1) - self.interval = kwargs.pop('interval', 50) - self.normalisation = kwargs.pop('data_norm', True) - - if self.normalisation is True: - self.x = x/x.max() - self.true_x = x - self.y = y/y.std() + y.mean() / y.std() - self.true_y = y - else: - self.x = x - self.y = y - - self.N = N - self.mid_point = mid_point - - self.fit_result = self.fit() - - def fit(self): - - def func1(params): - y_sum = self.y[self.mid_point] * np.sum([ - params[i] * (self.x / self.x[self.mid_point])**i - for i in range(self.N)], axis=0) - message = 'test' - const = constraint(params, message) - if const == 1: - return np.sum((self.y - y_sum)**2) - if const == -1: - return np.sum((self.y - y_sum)**2)+1e80 - - def constraint(params, message): - - m = np.arange(0, self.N, 1) - deriv = [] - for j in range(len(m)): - if m[j] >= 2: - dif = [] - for i in range(self.N - m[j]): - if i <= (self.N - m[j]): - dif_m_bit = ( - self.y[self.mid_point] / - self.x[self.mid_point]) * \ - np.math.factorial(m[j]+i) / \ - np.math.factorial(i) * \ - params[m[j]+i]*(self.x)**i / \ - (self.x[self.mid_point])**(i+1) - dif.append(dif_m_bit) - dif = np.array(dif) - derivative = dif.sum(axis=0) - deriv.append([m[j], derivative]) - deriv = np.array(deriv) - derive = [] - for i in range(deriv.shape[0]): - derive.append(deriv[i, 1]) - derive = np.array(derive) - pass_fail = [] # 1==pass, 0==fail - for i in range(derive.shape[0]): - if np.any(derive[i, :] == 1e-7): - pass_fail.append(0) - elif np.any(derive[i, :] > -1e-7) and \ - np.any(derive[i, :] < 1e-7): - pass_fail.append(0) - else: - pass_fail.append(1) - pass_fail = np.array(pass_fail) - - if np.any(pass_fail == 0): - const = -1 # failed - else: - const = 1 # satisfied - - return const - - params0 = [(self.y[-1]-self.y[0])/2]*(self.N) - res = basinhopping( - func1, params0, niter=10000, niter_success=100*self.N, - stepsize=self.step_size, T=self.temp, interval=self.interval, - seed=1) - print('msf fit params', res) - - parameters = res.x.copy() - if self.normalisation is True: - for i in range(len(parameters)): - if i == 0: - parameters[i] = parameters[i] * \ - (1+self.true_y.mean()/self.true_y[self.mid_point]) \ - - self.true_y.mean()/(self.true_y[self.mid_point]) - else: - parameters[i] = parameters[i] * \ - (1+self.true_y.mean()/self.true_y[self.mid_point]) - - message = 'summary' - const = constraint(parameters, message) - if const != 1: - print('Error: condition violated') - sys.exit(1) - - def fitting(params): - if self.normalisation is True: - y_sum = self.true_y[self.mid_point]*np.sum([ - params[i]*(self.true_x/self.true_x[self.mid_point])**i - for i in range(self.N)], axis=0) - else: - y_sum = self.y[self.mid_point]*np.sum([ - params[i]*(self.x/self.x[self.mid_point])**i - for i in range(self.N)], axis=0) - return y_sum - - fitted_y = fitting(parameters) - if self.normalisation is True: - print('chi BH', np.sum((self.true_y-fitted_y)**2)) - else: - print('chi BH', np.sum((self.y-fitted_y)**2)) - - return res.x - - -class max_fit_NM(object): - def __init__(self, x, y, N, BH_params, mid_point, **kwargs): - - self.normalisation = kwargs.pop('data_norm', True) - - if self.normalisation is True: - self.x = x/x.max() - self.true_x = x - self.y = y/y.std() + y.mean()/y.std() - self.true_y = y - else: - self.x = x - self.y = y - - self.N = N - self.BH_params = BH_params - self.mid_point = mid_point - self.fit_result, self.fit_h, self.chi = self.fit() - - def fit(self): - print('-------------------------------------------------------------') - - def func1(params): - y_sum = self.y[self.mid_point]*np.sum([ - params[i]*(self.x/self.x[self.mid_point])**i - for i in range(self.N)], axis=0) - message = 'test' - const, h = constraint(params, message) - if const == 1: - return np.sum((self.y-y_sum)**2) - if const == -1: - return np.sum((self.y-y_sum)**2)+1e80 - - def constraint(params, message): - - m = np.arange(0, self.N, 1) - deriv = [] - for j in range(len(m)): - if m[j] >= 2: - dif = [] - for i in range(self.N-m[j]): - if i <= (self.N-m[j]): - dif_m_bit = ( - self.y[self.mid_point] / - self.x[self.mid_point]) * \ - np.math.factorial(m[j]+i) / \ - np.math.factorial(i) * \ - params[m[j]+i]*(self.x)**i / \ - (self.x[self.mid_point])**(i+1) - dif.append(dif_m_bit) - dif = np.array(dif) - derivative = dif.sum(axis=0) - deriv.append([m[j], derivative]) - deriv = np.array(deriv) - derive = [] - for i in range(deriv.shape[0]): - derive.append(deriv[i, 1]) - derive = np.array(derive) - if message == 'summary': - print('min derivative', derive.min()) - h = np.abs(derive).min()/2 - if message == 'test': - h = None - pass - pass_fail = [] # 1==pass, 0==fail - for i in range(derive.shape[0]): - if np.any(derive[i, :] == 1e-7): - pass_fail.append(0) - elif np.any(derive[i, :] > -1e-7) and \ - np.any(derive[i, :] < 1e-7): - pass_fail.append(0) - else: - pass_fail.append(1) - pass_fail = np.array(pass_fail) - - if np.any(pass_fail == 0): - const = -1 # failed - else: - const = 1 # satisfied - - return const, h - - def fitting(params): - if self.normalisation is True: - y_sum = self.true_y[self.mid_point]*np.sum([ - params[i]*(self.true_x/self.true_x[self.mid_point])**i - for i in range(self.N)], axis=0) - else: - y_sum = self.y[self.mid_point]*np.sum([ - params[i]*(self.x/self.x[self.mid_point])**i - for i in range(self.N)], axis=0) - return y_sum - - params0 = self.BH_params - - res = minimize( - func1, params0, - options={ - 'maxiter': 100000, 'adaptive': True, 'fatol': 1e-7, - 'xatol': 1e-7}, method='Nelder-Mead') - print(res) - - parameters = res.x - if self.normalisation is True: - for i in range(len(parameters)): - if i == 0: - parameters[i] = parameters[i] * \ - (1+self.true_y.mean()/self.true_y[self.mid_point]) \ - - self.true_y.mean()/(self.true_y[self.mid_point]) - else: - parameters[i] = parameters[i] * \ - (1+self.true_y.mean()/self.true_y[self.mid_point]) - - fitted_y = fitting(parameters) - if self.normalisation is True: - print('chi NM', np.sum((self.true_y-fitted_y)**2)) - chi = np.sum((self.true_y-fitted_y)**2) - else: - print('chi NM', np.sum((self.y-fitted_y)**2)) - chi = np.sum((self.y-fitted_y)**2) - - message = 'summary' - const, h = constraint(parameters, message) - if const != 1: - print('Error: condition violated') - sys.exit(1) - - return parameters, h, chi diff --git a/Basin-hopping_Nelder-Mead/times_chis.py b/Basin-hopping_Nelder-Mead/times_chis.py deleted file mode 100644 index bc8f232..0000000 --- a/Basin-hopping_Nelder-Mead/times_chis.py +++ /dev/null @@ -1,120 +0,0 @@ -import numpy as np -import matplotlib.pyplot as plt -from maxsmooth.DCF import smooth -from BHNM.msf_fit import max_fit_BH, max_fit_NM -import time - - -def fit(params, mid_point, x, y, N): - y_sum = y[mid_point]*np.sum([ - params[i]*(x/x[mid_point])**i - for i in range(N)], axis=0) - return y_sum - - -"""x = np.linspace(1, 5, 100) -noise = np.random.normal(0,0.5,len(x)) -y = 0.6+2*x+4*x**(3)+9*x**(4)+noise""" - -x = np.loadtxt('x_data_for_comp.txt') -y = np.loadtxt('y_data_for_comp.txt') - -N = [3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13] - -# Testing all signs -times_qp = [] -chi_qp = [] -for i in range(len(N)): - s = time.time() - result = smooth( - x, y, N[i], model_type='normalised_polynomial', fit_type='qp') - e = time.time() - times_qp.append(e-s) - chi_qp.append(result.optimum_chi) - -# Using the maxsmooth sign navigating algorithm -times_qpsf = [] -chi_qpsf = [] -for i in range(len(N)): - s = time.time() - result = smooth(x, y, N[i], model_type='normalised_polynomial') - e = time.time() - times_qpsf.append(e-s) - chi_qpsf.append(result.optimum_chi) - -np.savetxt('times_qp_for_comp.txt', times_qp) -np.savetxt('times_qpsf_for_comp.txt', times_qpsf) -np.savetxt('chi_qp_for_comp.txt', chi_qp) -np.savetxt('chi_qpsf_for_comp.txt', chi_qpsf) - -# Basin-hopping followed by a Nelder-Mead routine -N_BHNM = np.arange(3, 8, 1) -mid_point = len(x)//2 -BHNM_chis_defualt = [] -BHNM_times_defualt = [] -h = [] -for i in range(len(N_BHNM)): - s = time.time() - BH = max_fit_BH( - x, y, N_BHNM[i], mid_point, step_size=10*N_BHNM[i], - interval=50, temp=5 * N_BHNM[i]) - NM = max_fit_NM(x, y, N_BHNM[i], BH.fit_result, mid_point) - fitted_y = fit(NM.fit_result, mid_point, x, y, N_BHNM[i]) - e = time.time() - print( - 'N:', N_BHNM[i], 'TIME:', e-s, 'CHI_RATIO:', - np.sum((y-fitted_y)**2)/chi_qpsf[i]) - print('CHI_BHNM:', np.sum((y-fitted_y)**2), 'CHI_QP:', chi_qpsf[i]) - chi = np.sum((y-fitted_y)**2) - BHNM_chis_defualt.append(chi) - BHNM_times_defualt.append(e-s) - -np.save('BHNM_times.npy', BHNM_times_defualt) -np.save('BHNM_chis.npy', BHNM_chis_defualt) - -chi_qp_plotting = [] -for i in range(len(N)): - if N[i] <= N_BHNM.max(): - chi_qp_plotting.append(chi_qpsf[i]) - -ratio = [] -for i in range(len(BHNM_chis_defualt)): - ratio.append(BHNM_chis_defualt[i]/chi_qp_plotting[i]) - -plt.subplot(111) -plt.plot(N_BHNM, ratio, ls='-', c='k') -plt.xlabel('N', fontsize=12) -plt.ylabel( - r'$\frac{X^2_{Basinhopping + Nelder-Mead}}{X^2_{QP Sampling Sign Space}}$', - fontsize=12) -plt.tight_layout() -plt.savefig('chi_ratio.pdf') -plt.close() - -plt.subplot(111) -plt.plot(N, times_qp, ls='-', label='QP Testing All Sign Combinations', c='g') -plt.plot(N, times_qpsf, ls='-', label='QP Sampling Sign Space', c='k') -plt.plot( - N_BHNM, BHNM_times_defualt, - ls='-', label='Basinhopping + Nelder-Mead\n', c='r') -plt.xlabel('N', fontsize=12) -plt.ylabel(r'$t$ [s]', fontsize=12) -plt.yscale('log') -plt.legend(loc=0) -plt.tight_layout() -plt.savefig('Times.pdf') -plt.close() - -plt.subplot(111) -plt.plot(N, chi_qp, ls='-', label='QP Testing All Sign Combinations', c='g') -plt.plot(N, chi_qpsf, ls='-', label='QP Sampling Sign Space', c='k') -plt.plot( - N_BHNM, BHNM_chis_defualt, - ls='-', label='Basinhopping + Nelder-Mead', c='r') -plt.xlabel('N', fontsize=12) -plt.ylabel(r'$X^2$', fontsize=12) -plt.yscale('log') -plt.legend(loc=0) -plt.tight_layout() -plt.savefig('chi.pdf') -plt.close() diff --git a/Basin-hopping_Nelder-Mead/x_data_for_comp.txt b/Basin-hopping_Nelder-Mead/x_data_for_comp.txt deleted file mode 100644 index 9f718a4..0000000 --- a/Basin-hopping_Nelder-Mead/x_data_for_comp.txt +++ /dev/null @@ -1,100 +0,0 @@ -1.000000000000000000e+00 -1.040404040404040442e+00 -1.080808080808080884e+00 -1.121212121212121104e+00 -1.161616161616161547e+00 -1.202020202020201989e+00 -1.242424242424242431e+00 -1.282828282828282873e+00 -1.323232323232323315e+00 -1.363636363636363757e+00 -1.404040404040403978e+00 -1.444444444444444420e+00 -1.484848484848484862e+00 -1.525252525252525304e+00 -1.565656565656565746e+00 -1.606060606060605966e+00 -1.646464646464646631e+00 -1.686868686868686851e+00 -1.727272727272727293e+00 -1.767676767676767735e+00 -1.808080808080808177e+00 -1.848484848484848619e+00 -1.888888888888888840e+00 -1.929292929292929504e+00 -1.969696969696969724e+00 -2.010101010101010388e+00 -2.050505050505050608e+00 -2.090909090909090828e+00 -2.131313131313131493e+00 -2.171717171717171713e+00 -2.212121212121211933e+00 -2.252525252525252597e+00 -2.292929292929293261e+00 -2.333333333333333481e+00 -2.373737373737373701e+00 -2.414141414141414366e+00 -2.454545454545454586e+00 -2.494949494949494806e+00 -2.535353535353535470e+00 -2.575757575757576134e+00 -2.616161616161616355e+00 -2.656565656565656575e+00 -2.696969696969697239e+00 -2.737373737373737459e+00 -2.777777777777777679e+00 -2.818181818181818343e+00 -2.858585858585859008e+00 -2.898989898989899228e+00 -2.939393939393939448e+00 -2.979797979797980112e+00 -3.020202020202020332e+00 -3.060606060606060996e+00 -3.101010101010101216e+00 -3.141414141414141437e+00 -3.181818181818182101e+00 -3.222222222222222321e+00 -3.262626262626262985e+00 -3.303030303030303205e+00 -3.343434343434343425e+00 -3.383838383838384090e+00 -3.424242424242424310e+00 -3.464646464646464974e+00 -3.505050505050505194e+00 -3.545454545454545858e+00 -3.585858585858586078e+00 -3.626262626262626299e+00 -3.666666666666666963e+00 -3.707070707070707183e+00 -3.747474747474747847e+00 -3.787878787878788067e+00 -3.828282828282828731e+00 -3.868686868686868952e+00 -3.909090909090909172e+00 -3.949494949494949836e+00 -3.989898989898990056e+00 -4.030303030303031164e+00 -4.070707070707070940e+00 -4.111111111111110716e+00 -4.151515151515152269e+00 -4.191919191919192045e+00 -4.232323232323232709e+00 -4.272727272727273373e+00 -4.313131313131313149e+00 -4.353535353535353813e+00 -4.393939393939394478e+00 -4.434343434343434254e+00 -4.474747474747474918e+00 -4.515151515151515582e+00 -4.555555555555555358e+00 -4.595959595959596911e+00 -4.636363636363636687e+00 -4.676767676767676463e+00 -4.717171717171718015e+00 -4.757575757575757791e+00 -4.797979797979798455e+00 -4.838383838383839120e+00 -4.878787878787878896e+00 -4.919191919191919560e+00 -4.959595959595960224e+00 -5.000000000000000000e+00 diff --git a/Basin-hopping_Nelder-Mead/y_data_for_comp.txt b/Basin-hopping_Nelder-Mead/y_data_for_comp.txt deleted file mode 100644 index 0b4c8c2..0000000 --- a/Basin-hopping_Nelder-Mead/y_data_for_comp.txt +++ /dev/null @@ -1,100 +0,0 @@ -1.538480190059443231e+01 -1.803994844410931719e+01 -2.074192622407534614e+01 -2.340100984140073592e+01 -2.621510063558880432e+01 -2.887939230577280370e+01 -3.267760037624015723e+01 -3.624293586100588982e+01 -4.183230771978232099e+01 -4.472484203235710254e+01 -4.901901109503933185e+01 -5.458384801381982498e+01 -6.044395282904604727e+01 -6.710163163497058747e+01 -7.331305934512305100e+01 -8.038417701150292771e+01 -8.828937593277244389e+01 -9.658444786378277058e+01 -1.042045848601506322e+02 -1.146416307471176879e+02 -1.235112892305587451e+02 -1.345268586166006060e+02 -1.450430823504431430e+02 -1.581643332442427266e+02 -1.717913514750528066e+02 -1.835712719659944128e+02 -1.979150031939569203e+02 -2.139361084172272456e+02 -2.295680700846963020e+02 -2.465373692074055327e+02 -2.641143352719346922e+02 -2.831127401702224233e+02 -3.022997542458672910e+02 -3.228565015249954513e+02 -3.445193178888176249e+02 -3.676949607788102412e+02 -3.913593215276684418e+02 -4.168242323983324695e+02 -4.426368502745826845e+02 -4.702519040672153210e+02 -4.999275768833090865e+02 -5.289023404506694988e+02 -5.603229923493009892e+02 -5.931106443932800403e+02 -6.271303491261907084e+02 -6.633693290327342993e+02 -7.017142035218076899e+02 -7.397376073564449825e+02 -7.806386443959225971e+02 -8.228136223978054886e+02 -8.656723069434885929e+02 -9.102137208256813210e+02 -9.580326031355238001e+02 -1.007291992119634074e+03 -1.058029537451610622e+03 -1.111047519335722654e+03 -1.166443910326044715e+03 -1.222809821319707908e+03 -1.281110152112206606e+03 -1.342530443746999026e+03 -1.405175264221236603e+03 -1.471403152526458143e+03 -1.536878719515193097e+03 -1.608101102592653660e+03 -1.679812150650684544e+03 -1.754709199343822320e+03 -1.832464343650941828e+03 -1.911709046912803842e+03 -1.994382101504000275e+03 -2.078654055364936994e+03 -2.166598274438367753e+03 -2.255622982675686217e+03 -2.348654835621773600e+03 -2.444463611754800695e+03 -2.543845966417291038e+03 -2.645305246365588118e+03 -2.749656581282604748e+03 -2.858029640524120168e+03 -2.968665196804036896e+03 -3.083208241631706187e+03 -3.198777704047548468e+03 -3.320655724687552720e+03 -3.444593388278271505e+03 -3.572310906936728770e+03 -3.704479259388221635e+03 -3.838077949702294973e+03 -3.975746853596534947e+03 -4.118423306260904610e+03 -4.263847070864699162e+03 -4.413647032393576410e+03 -4.567425210209468787e+03 -4.724884307810286373e+03 -4.886477595104389366e+03 -5.051962838680952700e+03 -5.220691820972367168e+03 -5.395668600079392490e+03 -5.574189572375433272e+03 -5.756885988219393766e+03 -5.943798781023513584e+03 -6.134840236402888877e+03 diff --git a/README.rst b/README.rst index 4e07d20..06bd150 100644 --- a/README.rst +++ b/README.rst @@ -10,7 +10,7 @@ Introduction :maxsmooth: Derivative Constrained Function Fitting :Author: Harry Thomas Jones Bevins -:Version: 1.2.3 +:Version: 1.2.4 :Homepage: https://github.com/htjb/maxsmooth :Documentation: https://maxsmooth.readthedocs.io/ diff --git a/bin/check_version.py b/bin/check_version.py new file mode 100644 index 0000000..deeed0b --- /dev/null +++ b/bin/check_version.py @@ -0,0 +1,118 @@ +#!/usr/bin/env python +"""Check Version. + +Verify version has been incremented. + +Based on check_version.py from the anesthetic GitHub repository: +https://github.com/handley-lab/anesthetic +""" + + +import sys +import subprocess +from packaging import version + + +# Filestructure +readme_file = "README.rst" + + +# Utility functions +def run_on_commandline(*args): + """Run the given arguments as a command on the command line.""" + return subprocess.run(args, text=True, stdout=subprocess.PIPE, + stderr=subprocess.PIPE).stdout + + +def unit_incremented(version_a: str, version_b: str) -> bool: + """Check if version_a is one version larger than version_b. + + Parameters + ---------- + version_a: str + New version, version number + version_b: str + Old version, version number + + Returns + ------- + unit_incremented: bool + Whether, version number has been unit incremented + """ + # Convert to version objects + version_a = version.parse(version_a) + version_b = version.parse(version_b) + + # Check increment for pre-release versions + if version_a.pre is not None and version_b.pre is not None: + # Matching pre-release levels + if version_a.pre[0] == version_b.pre[0]: + return ((version_a.pre[1] == (int(version_b.pre[1]) + 1)) and + (version_a.base_version == version_b.base_version)) + # Differing pre-release levels + else: + return (version_a.pre[1] == 0 and + version_a.pre[0] > version_b.pre[0] and + version_a.base_version == version_b.base_version) + + # New pre-release level + elif version_a.pre is not None: + return (version_a.base_version > version_b.base_version and + version_a.pre[1] == 0) + + # Full release + elif version_b.pre is not None: + return version_a.base_version == version_b.base_version + + # Standard version major, minor and micro increments + else: + return (version_a.micro == version_b.micro + 1 and + version_a.minor == version_b.minor and + version_a.major == version_b.major or + version_a.micro == 0 and + version_a.minor == version_b.minor + 1 and + version_a.major == version_b.major or + version_a.micro == 0 and + version_a.minor == 0 and + version_a.major == version_b.major+1) + +def get_current_version() -> str: + """Get current version of package from README.rst""" + current_version = run_on_commandline("grep", ":Version:", readme_file) + current_version = current_version.split(":")[-1].strip() + return current_version + +def main(): + """Check version is consistent and incremented correctly.""" + # Get current version from readme + current_version = get_current_version() + + # Get previous version from main branch of code + run_on_commandline("git", "fetch", "origin", "master") + readme_contents = run_on_commandline("git", "show", + "remotes/origin/master:" + + readme_file) + + previous_version = None + for line in readme_contents.splitlines(): + if ":Version:" in line: + previous_version = line.split(":")[-1].strip() + break + + if previous_version is None: + raise ValueError("Could not find version in README.rst on master branch") + + # Check versions have been incremented + if not unit_incremented(current_version, previous_version): + sys.stderr.write(("Version must be incremented by one:\n" + "HEAD: {},\n" + "master: {}.\n").format(current_version, + previous_version)) + sys.exit(1) + + # No issues found, exit happily :) + sys.exit(0) + + +if __name__ == "__main__": + main() \ No newline at end of file diff --git a/maxsmooth/DCF.py b/maxsmooth/DCF.py index c38a1d4..97379b4 100755 --- a/maxsmooth/DCF.py +++ b/maxsmooth/DCF.py @@ -8,7 +8,7 @@ """ -from maxsmooth.qp import qp_class +from maxsmooth.qp import qp_class, precompute_matrices from maxsmooth.Models import Models_class from maxsmooth.derivatives import derivative_class from maxsmooth.Data_save import save, save_optimum @@ -369,12 +369,16 @@ def qp(x, y, pivot_point): # Testing all signs append_params, append_chi, append_zc_dict, append_passed_signs = \ params.append, chi_squared.append, zc_dict.append, \ passed_signs.append + precomputed = precompute_matrices( + x, y, self.N, pivot_point, self.model_type, + self.zero_crossings, self.constraints, self.new_basis) for j in range(signs.shape[0]): fit = qp_class( x, y, self.N, signs[j, :], pivot_point, self.model_type, self.cvxopt_maxiter, self.zero_crossings, - self.initial_params, self.constraints, self.new_basis) + self.initial_params, self.constraints, self.new_basis, + precomputed=precomputed) if self.print_output == 2: print('-'*50) @@ -499,11 +503,15 @@ def qp_sign_flipping(x, y, pivot_point): # Sign Sampling for i in range(len(array_signs)): if i == r: tested_indices.append(i) + precomputed = precompute_matrices( + x, y, self.N, pivot_point, self.model_type, + self.zero_crossings, self.constraints, self.new_basis) fit = qp_class( x, y, self.N, signs, pivot_point, self.model_type, self.cvxopt_maxiter, self.zero_crossings, - self.initial_params, self.constraints, self.new_basis) + self.initial_params, self.constraints, self.new_basis, + precomputed=precomputed) chi_squared.append(fit.chi_squared) tested_signs.append(signs) parameters.append(fit.parameters) @@ -570,7 +578,8 @@ def qp_sign_flipping(x, y, pivot_point): # Sign Sampling x, y, self.N, signs, pivot_point, self.model_type, self.cvxopt_maxiter, self.zero_crossings, self.initial_params, - self.constraints, self.new_basis) + self.constraints, self.new_basis, + precomputed=precomputed) if fit.chi_squared < chi_squared_old: chi_squared_new = fit.chi_squared previous_signs = signs @@ -659,7 +668,8 @@ def qp_sign_flipping(x, y, pivot_point): # Sign Sampling x, y, self.N, signs, pivot_point, self.model_type, self.cvxopt_maxiter, self.zero_crossings, self.initial_params, - self.constraints, self.new_basis) + self.constraints, self.new_basis, + precomputed=precomputed) chi_down = fit.chi_squared chi_squared.append(fit.chi_squared) tested_signs.append(signs) @@ -721,7 +731,8 @@ def qp_sign_flipping(x, y, pivot_point): # Sign Sampling x, y, self.N, signs, pivot_point, self.model_type, self.cvxopt_maxiter, self.zero_crossings, self.initial_params, - self.constraints, self.new_basis) + self.constraints, self.new_basis, + precomputed=precomputed) chi_up = fit.chi_squared chi_squared.append(fit.chi_squared) tested_signs.append(signs) diff --git a/maxsmooth/_version.py b/maxsmooth/_version.py index cd7e2da..fecd8c0 100644 --- a/maxsmooth/_version.py +++ b/maxsmooth/_version.py @@ -1 +1 @@ -__version__ = "1.2.3" \ No newline at end of file +__version__ = "1.2.4" \ No newline at end of file diff --git a/maxsmooth/derivatives.py b/maxsmooth/derivatives.py index 00767ee..3888a4e 100755 --- a/maxsmooth/derivatives.py +++ b/maxsmooth/derivatives.py @@ -27,7 +27,7 @@ def derivatives_func(self): def mth_order_derivatives(m): if self.derivatives_function is None: - if np.any(self.model_type != ['legendre', 'exponential']): + if self.model_type not in ['legendre', 'exponential']: mth_order_derivative = [] for i in range(self.N-m): if self.model_type == 'normalised_polynomial': diff --git a/maxsmooth/qp.py b/maxsmooth/qp.py index 00618c9..7da900f 100755 --- a/maxsmooth/qp.py +++ b/maxsmooth/qp.py @@ -9,11 +9,177 @@ warnings.simplefilter('always', UserWarning) +def _constraint_prefactors(m, x, y, N, pivot_point, model_type, + derivative_pres, args): + """Compute the constraint matrix prefactors for derivative order m.""" + if derivative_pres is None: + if model_type not in ['legendre', 'exponential']: + derivatives = [] + for i in range(N): + if i <= m - 1: + derivatives.append([0]*len(x)) + for i in range(N-m): + if model_type == 'normalised_polynomial': + mth_order_derivative_term = ( + y[pivot_point] / + x[pivot_point]) \ + * math.factorial(m + i) \ + / math.factorial(i) * \ + x**i / x[pivot_point]**(i + 1) + derivatives.append(mth_order_derivative_term) + if model_type == 'polynomial': + mth_order_derivative_term = math.factorial(m+i)\ + / math.factorial(i) * x**i + derivatives.append(mth_order_derivative_term) + if model_type == 'log_polynomial': + mth_order_derivative_term = math.factorial(m+i)\ + / math.factorial(i) * \ + np.log10(x/x[pivot_point])**i + derivatives.append(mth_order_derivative_term) + if model_type == 'loglog_polynomial': + mth_order_derivative_term = math.factorial(m+i)\ + / math.factorial(i) * np.log10(x)**i + derivatives.append(mth_order_derivative_term) + if model_type == 'difference_polynomial': + mth_order_derivative_term = math.factorial(m+i)\ + / math.factorial(i) * ( + x - x[pivot_point])**i + derivatives.append(mth_order_derivative_term) + + if derivative_pres is not None: + if args is None: + derivatives = derivative_pres(m, x, y, N, pivot_point) + if args is not None: + derivatives = derivative_pres(m, x, y, N, pivot_point, *args) + + if model_type == 'legendre': + interval = np.linspace(-0.999, 0.999, len(x)) + alps = [] + for i in range(N): + alps.append(lpmv(m, i, interval)) + alps = np.array(alps) + derivatives = [] + for row in range(len(alps)): + derivatives.append( + ((alps[row, :]*(-1)**(m))/(1-interval**2)**(m/2))) + derivatives = np.array(derivatives) + if model_type == 'exponential': + derivatives = np.empty([N, len(x)]) + for i in range(N): + for col in range(len(x)): + derivatives[i, col] = \ + y[pivot_point] * \ + (np.exp(-i*x[col]/x[pivot_point])) * \ + (-i/x[pivot_point])**m + derivatives = np.array(derivatives) + + derivatives = np.array(derivatives).astype(np.double) + derivatives = matrix(derivatives) + if derivatives.size == (len(x), N): + pass + else: + derivatives = derivatives.T + return derivatives + + +def precompute_matrices(x, y, N, pivot_point, model_type, zero_crossings, + constraints, new_basis): + """Precompute QP matrices that do not depend on the sign vector. + + Returns a dict with base_derivatives, phi, Q, q, h. Pass as + precomputed=... to qp_class when calling it in a loop over sign + combinations to avoid redundant computation. + """ + derivative_pres = new_basis['der_pres'] + basis_functions = new_basis['basis_function'] + args = new_basis['args'] + + m_vals = np.arange(0, N, 1) + base_derivatives = [] + for i in range(len(m_vals)): + if m_vals[i] >= constraints: + if zero_crossings is not None: + if m_vals[i] not in set(zero_crossings): + dp = _constraint_prefactors( + m_vals[i], x, y, N, pivot_point, model_type, + derivative_pres, args) + if dp != []: + base_derivatives.append(dp) + else: + dp = _constraint_prefactors( + m_vals[i], x, y, N, pivot_point, model_type, + derivative_pres, args) + if dp != []: + base_derivatives.append(dp) + + if basis_functions is None: + phi = np.empty([len(x), N]) + if model_type != 'legendre': + for row in range(len(x)): + for i in range(N): + if model_type == 'normalised_polynomial': + phi[row, i] = y[pivot_point] * ( + x[row] / x[pivot_point])**i + if model_type == 'polynomial': + phi[row, i] = (x[row])**i + if model_type == 'log_polynomial': + phi[row, i] = \ + np.log10(x[row]/x[pivot_point])**i + if model_type == 'loglog_polynomial': + phi[row, i] = np.log10(x[row])**i + if model_type == 'difference_polynomial': + phi[row, i] = (x[row]-x[pivot_point])**i + if model_type == 'exponential': + phi[row, i] = y[pivot_point] * \ + np.exp(-i*x[row]/x[pivot_point]) + if model_type == 'legendre': + interval = np.linspace(-0.999, 0.999, len(x)) + phi = [] + for i in range(N): + P = legendre(i) + phi.append(P(interval)) + phi = np.array(phi).T + phi = matrix(phi) + else: + if args is None: + phi = basis_functions(x, y, pivot_point, N) + phi = matrix(phi) + if args is not None: + phi = basis_functions(x, y, pivot_point, N, *args) + phi = matrix(phi) + + if model_type == 'loglog_polynomial': + data_matrix = matrix( + np.log10(y).astype(np.double), (len(y), 1), 'd') + else: + data_matrix = matrix( + y.astype(np.double), (len(y), 1), 'd') + + Q = phi.T * phi + q = -phi.T * data_matrix + + if zero_crossings is None: + h = matrix(0.0, ((N - constraints) * len(x), 1), 'd') + else: + h = matrix( + 0.0, ( + (N - constraints - len(zero_crossings)) + * len(x), 1), 'd') + + return { + 'base_derivatives': base_derivatives, + 'phi': phi, + 'Q': Q, + 'q': q, + 'h': h, + } + + class qp_class(object): def __init__( self, x, y, N, signs, pivot_point, model_type, cvxopt_maxiter, zero_crossings, initial_params, - constraints, new_basis): + constraints, new_basis, precomputed=None): self.model_type = model_type self.pivot_point = pivot_point self.y = y @@ -30,6 +196,7 @@ def __init__( self.args = new_basis['args'] self.new_basis = new_basis self.constraints = constraints + self.precomputed = precomputed self.parameters, self.chi_squared, self.zc_dict = self.fit() def fit(self): @@ -37,158 +204,23 @@ def fit(self): solvers.options['maxiters'] = self.cvxopt_maxiter solvers.options['show_progress'] = False - def constraint_prefactors(m): - # Derivative prefactors on parameters - if self.derivative_pres is None: - if np.any(self.model_type != ['legendre', 'exponential']): - derivatives = [] - for i in range(self.N): - if i <= m - 1: - derivatives.append([0]*len(self.x)) - for i in range(self.N-m): - if self.model_type == 'normalised_polynomial': - mth_order_derivative_term = ( - self.y[self.pivot_point] / - self.x[self.pivot_point]) \ - * math.factorial(m + i) \ - / math.factorial(i) * \ - (self.x)**i/(self.x[self.pivot_point])**(i + 1) - derivatives.append(mth_order_derivative_term) - if self.model_type == 'polynomial': - mth_order_derivative_term = math.factorial(m+i)\ - / math.factorial(i) * (self.x)**i - derivatives.append(mth_order_derivative_term) - if self.model_type == 'log_polynomial': - mth_order_derivative_term = math.factorial(m+i)\ - / math.factorial(i) * \ - np.log10(self.x/self.x[self.pivot_point])**i - derivatives.append(mth_order_derivative_term) - if self.model_type == 'loglog_polynomial': - mth_order_derivative_term = math.factorial(m+i)\ - / math.factorial(i) * np.log10(self.x)**i - derivatives.append(mth_order_derivative_term) - if self.model_type == 'difference_polynomial': - mth_order_derivative_term = math.factorial(m+i)\ - / math.factorial(i) * ( - self.x - self.x[self.pivot_point])**i - derivatives.append(mth_order_derivative_term) - - if self.derivative_pres is not None: - if self.args is None: - derivatives = self.derivative_pres( - m, self.x, self.y, self.N, self.pivot_point) - if self.args is not None: - derivatives = self.derivative_pres( - m, self.x, self.y, self.N, self.pivot_point, - *self.args) - - if self.model_type == 'legendre': - interval = np.linspace(-0.999, 0.999, len(self.x)) - alps = [] - for i in range(self.N): - alps.append(lpmv(m, i, interval)) - alps = np.array(alps) - derivatives = [] - for h in range(len(alps)): - derivatives.append( - ((alps[h, :]*(-1)**(m))/(1-interval**2)**(m/2))) - derivatives = np.array(derivatives) - if self.model_type == 'exponential': - derivatives = np.empty([self.N, len(self.x)]) - for i in range(self.N): - for h in range(len(self.x)): - derivatives[i, h] = \ - self.y[self.pivot_point] * \ - (np.exp(-i*self.x[h]/self.x[self.pivot_point])) * \ - (-i/self.x[self.pivot_point])**m - derivatives = np.array(derivatives) - - derivatives = np.array(derivatives).astype(np.double) - derivatives = matrix(derivatives) - if derivatives.size == (len(self.x), self.N): - pass - else: - derivatives = derivatives.T - return derivatives - - m = np.arange(0, self.N, 1) - derivatives = [] - signs = matrix(self.signs) - for i in range(len(m)): - if m[i] >= self.constraints: - if self.zero_crossings is not None: - if m[i] not in set(self.zero_crossings): - derivative_prefactors = constraint_prefactors(m[i]) - if derivative_prefactors != []: - derivatives.append(derivative_prefactors) - else: - derivative_prefactors = constraint_prefactors(m[i]) - if derivative_prefactors != []: - derivatives.append(derivative_prefactors) - - for i in range(len(derivatives)): - derivatives[i] *= signs[i] - - G = matrix(derivatives) - - if self.basis_functions is None: - phi = np.empty([len(self.x), self.N]) - if self.model_type != 'legendre': - for h in range(len(self.x)): - for i in range(self.N): - if self.model_type == 'normalised_polynomial': - phi[h, i] = self.y[self.pivot_point] * ( - self.x[h] / self.x[self.pivot_point])**i - if self.model_type == 'polynomial': - phi[h, i] = (self.x[h])**i - if self.model_type == 'log_polynomial': - phi[h, i] = \ - np.log10(self.x[h]/self.x[self.pivot_point])**i - if self.model_type == 'loglog_polynomial': - phi[h, i] = np.log10(self.x[h])**i - if self.model_type == 'difference_polynomial': - phi[h, i] = (self.x[h]-self.x[self.pivot_point])**i - if self.model_type == 'exponential': - phi[h, i] = self.y[self.pivot_point] * \ - np.exp(-i*self.x[h]/self.x[self.pivot_point]) - if self.model_type == 'legendre': - interval = np.linspace(-0.999, 0.999, len(self.x)) - phi = [] - for i in range(self.N): - P = legendre(i) - phi.append(P(interval)) - phi = np.array(phi).T - phi = matrix(phi) - if self.basis_functions is not None: - if self.args is None: - phi = self.basis_functions( - self.x, self.y, self.pivot_point, self.N) - phi = matrix(phi) - if self.args is not None: - phi = self.basis_functions( - self.x, self.y, self.pivot_point, self.N, *self.args) - phi = matrix(phi) - - if self.model_type == 'loglog_polynomial': - data_matrix = matrix( - np.log10(self.y).astype(np.double), (len(self.y), 1), - 'd') - else: - data_matrix = matrix( - self.y.astype(np.double), (len(self.y), 1), - 'd') - - if self.zero_crossings is None: - h = matrix(0.0, ((self.N-self.constraints)*len(self.x), 1), 'd') + if self.precomputed is not None: + base_derivatives = self.precomputed['base_derivatives'] + Q = self.precomputed['Q'] + q = self.precomputed['q'] + h = self.precomputed['h'] else: - h = matrix( - 0.0, ( - (self.N-self.constraints-len(self.zero_crossings)) - * len(self.x), 1), 'd') - - Q = phi.T*phi + pc = precompute_matrices( + self.x, self.y, self.N, self.pivot_point, self.model_type, + self.zero_crossings, self.constraints, self.new_basis) + base_derivatives = pc['base_derivatives'] + Q = pc['Q'] + q = pc['q'] + h = pc['h'] - q = -phi.T*data_matrix + signs = matrix(self.signs) + G = matrix([base_derivatives[i] * signs[i] + for i in range(len(base_derivatives))]) if self.initial_params is None: qpfit = solvers.qp(Q, q, G, h) diff --git a/pyproject.toml b/pyproject.toml new file mode 100644 index 0000000..8aab0bc --- /dev/null +++ b/pyproject.toml @@ -0,0 +1,34 @@ +[build-system] +requires = ["setuptools>=77.0", "wheel"] +build-backend = "setuptools.build_meta" + +[project] +name = "maxsmooth" +version = "1.2.4" +description = "Derivative-constrained function fitting." +readme = "README.rst" +requires-python = ">=3.6" +dependencies = [ + 'matplotlib', + 'numpy', + 'CVXOPT', + 'scipy', + 'progressbar', +] + +[tool.setuptools.packages.find] +where = ["."] +include = ["maxsmooth*"] + +[project.optional-dependencies] +dev = [ + "pytest", + "ruff", + "pre-commit", + "packaging" +] +docs = [ + 'sphinx', + 'sphinx_rtd_theme', + 'numpydoc' +] \ No newline at end of file diff --git a/requirements.txt b/requirements.txt deleted file mode 100644 index 0552818..0000000 --- a/requirements.txt +++ /dev/null @@ -1,5 +0,0 @@ -matplotlib -numpy -CVXOPT -scipy -progressbar diff --git a/setup.py b/setup.py deleted file mode 100644 index 5420251..0000000 --- a/setup.py +++ /dev/null @@ -1,37 +0,0 @@ -from setuptools import setup, find_packages - - -def readme(short=False): - with open('README.rst') as f: - if short: - return f.readlines()[1].strip() - else: - return f.read() - -setup( - name='maxsmooth', - version='1.2.2', - description='maxsmooth:Derivative Constrained Function Fitting', - long_description=readme(), - author='Harry T. J. Bevins', - author_email='htjb2@cam.ac.uk', - url='https://github.com/htjb/maxsmooth', - packages=find_packages(), - install_requires=['numpy', 'scipy', 'cvxopt', 'matplotlib', 'progressbar'], - license='MIT', - extras_require={ - 'docs': ['sphinx', 'sphinx_rtd_theme', 'numpydoc'], - }, - tests_require=['pytest'], - classifiers=[ - 'Intended Audience :: Science/Research', - 'Natural Language :: English', - 'License :: OSI Approved :: MIT License', - 'Programming Language :: Python :: 3.6', - 'Programming Language :: Python :: 3.7', - 'Programming Language :: Python :: 3.8', - 'Topic :: Scientific/Engineering', - 'Topic :: Scientific/Engineering :: Astronomy', - 'Topic :: Scientific/Engineering :: Physics', - ], -)