diff --git a/.github/workflows/python-app.yml b/.github/workflows/python-app.yml new file mode 100644 index 0000000..47d0ba7 --- /dev/null +++ b/.github/workflows/python-app.yml @@ -0,0 +1,39 @@ +# This workflow will install Python dependencies, run tests and lint with a single version of Python +# For more information see: https://docs.github.com/en/actions/automating-builds-and-tests/building-and-testing-python + +name: KinOpt application + +on: + push: + branches: [ "master" ] + pull_request: + branches: [ "master" ] + +permissions: + contents: read + +jobs: + build: + + runs-on: ubuntu-latest + + steps: + - uses: actions/checkout@v4 + - name: Set up Python 3.10 + uses: actions/setup-python@v3 + with: + python-version: "3.10" + - name: Install dependencies + run: | + python -m pip install --upgrade pip + pip install flake8 pytest + if [ -f requirements.txt ]; then pip install -r requirements.txt; fi + - name: Lint with flake8 + run: | + # stop the build if there are Python syntax errors or undefined names + flake8 . --count --select=E9,F63,F7,F82 --show-source --statistics + # exit-zero treats all errors as warnings. The GitHub editor is 127 chars wide + flake8 . --count --exit-zero --max-complexity=10 --max-line-length=127 --statistics + - name: Test with pytest + run: | + pytest diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md new file mode 100644 index 0000000..9ae218f --- /dev/null +++ b/CONTRIBUTING.md @@ -0,0 +1,62 @@ +# Community Guidelines for KinOpt + +Welcome to the KinOpt community! This project is open to contributors from all backgrounds, and we aim to foster a collaborative and respectful environment focused on chemical kinetics and kinetic modeling. + +This document provides clear guidance for third parties who wish to: + +--- + +## 1. Contribute to the Software + +We welcome contributions that improve the functionality, performance, or usability of KinOpt. If you would like to contribute: + +- For most users with a chemical kinetics background, your primary interest will likely be in adding new kinetic models. These should be added to the file: + `kinopt/src/kinetic_models.py` + +- When adding a model: + - Create a function named in the format: + `rate_for_[name_of_the_model]` + - The first two arguments **must** be: + - `extent`: extent of reaction + - `T`: temperature of the reaction + - Any additional arguments should correspond to **user input parameters** + +- Please ensure your code is: + - Well-documented + - Properly tested (include tests if possible) + - Consistent with the existing codebase and style + +- Submit contributions via a **Pull Request** on GitHub: + [https://github.com/alan-tabore/KinOpt](https://github.com/alan-tabore/KinOpt) + +--- + +## 2. Report Issues or Problems + +We encourage users to report bugs, unexpected behavior, or any other problems with the software. + +- Please check the [Issues page](https://github.com/alan-tabore/KinOpt/issues) to see if your concern has already been raised. +- When reporting a new issue: + - Be clear and concise + - Include your operating system, Python version, and KinOpt version + - If possible, provide a minimal reproducible example + +--- + +## 3. Seek Support + +If you need help using KinOpt or understanding how to contribute: + +- Feel free to open a **discussion or issue** on the GitHub page: + [https://github.com/alan-tabore/KinOpt](https://github.com/alan-tabore/KinOpt) + +- When asking for help: + - Describe your problem clearly + - Include relevant code or context + - Be respectful of the maintainers’ time + +--- + +## General Conduct + +While our community is primarily technical, we expect all interactions to remain professional and respectful. Disruptive or abusive behavior will not be tolerated. \ No newline at end of file diff --git a/README.md b/README.md index 8eeec43..a85f41e 100644 --- a/README.md +++ b/README.md @@ -41,13 +41,10 @@ Once you’ve downloaded the project, you can install the required python module Open a command prompt in the KinOpt folder and execute: ``` bash -python -m pip install requirements.txt -``` -or with conda: -``` bash -conda install requirements.txt +python -m pip install -r requirements.txt ``` + ### Tutorials Tutorials are available on Youtube to show you how to install and use the software: diff --git a/docs/source/launch_data_extraction.rst b/docs/source/launch_data_extraction.rst index 41cc685..29a91af 100644 --- a/docs/source/launch_data_extraction.rst +++ b/docs/source/launch_data_extraction.rst @@ -8,5 +8,9 @@ Store them in the same folder and: 1. Click on the "Add file(s)" button of the data extraction dock 2. Navigate to the folder containing your data -3. Selected all your files -4. Make sure that the file list becomes green (indicating the extraction was successful) +3. Select all your files +4. Click the "open" button. (your files should now be loaded) +5. Check the box "has header" if the first line of your file is a header. (if you have multiple lines to skip at the beginning of your file, indicate the number of lines to skip in the entry box "Number of lines to skip". +6. Indicate the delimiter of your data in the "Delimiter" entry box. If no delimiter is indicated, the default delimiter (a comma ",") will be used. The delimiter used in the example data files is a comma ",". The tab delimiter can be specified as "\t" or "tab" +7. Click the extract button. +8. Make sure that the file list becomes green (indicating the extraction was successful) diff --git a/kinopt/src/__init__.py b/kinopt/src/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/kinopt/src/kinetic_models.py b/kinopt/src/kinetic_models.py index 95b77e5..8614d3c 100644 --- a/kinopt/src/kinetic_models.py +++ b/kinopt/src/kinetic_models.py @@ -21,11 +21,11 @@ def arrhenius_rate_constant(T,A,Ea): Parameters ---------- T : Float - Temperature of reaction. + Temperature of reaction in Kelvin. A : Float Pre-exponential factor. Ea : Float - Activation energy. + Activation energy in J/mol. Returns ------- @@ -53,11 +53,11 @@ def rate_for_nth_order(extent,T,A1,E1,n): extent : 1-D array Extent of reaction. T : 1-D array - Temperature of reaction. + Temperature of reaction in Kelvin. A : Float Pre-exponential factor of the reaction. Ea : Float - Activation energy of the reaction. + Activation energy of the reaction in J/mol. n : Float Order of reaction. @@ -67,9 +67,48 @@ def rate_for_nth_order(extent,T,A1,E1,n): Rate of reaction for a nth order reaction. """ - return arrhenius_rate_constant(T,A1,E1)*(1-extent)**n + return arrhenius_rate_constant(T,A1,E1)*((1-extent)**n) +def rate_for_autocatalytic(extent,T,A,Ea,m,n): + r""" + Compute the rate of reaction for an autocatalytic reaction. + + Parameters + ---------- + extent : ndarray + Extent of reaction. + T : ndarray + Temperature of reaction in Kelvin. + A : float + Pre-exponential factor of the reaction. + Ea : float + Activation energy of the reaction in J/mol. + m : float + Order of reaction for the autocatalyzed reaction. + n : float + Order of reaction for the regular reaction. + + Returns + ------- + rate : ndarray + Rate of reaction for an autocatalytic equation. + + Notes + ----- + The rate for an autocatalytic reaction is given by: + + .. math:: \frac{d\alpha}{dt} = A e^{ \left( \frac{-E_a}{RT} \right)} \alpha^m (1-\alpha)^n + + References + ---------- + [1] M. R. Keenan, « Autocatalytic cure kinetics from DSC measurements: Zero initial cure rate », J. Appl. Polym. Sci., vol. 33, nᵒ 5, p. 1725‑1734, avr. 1987, doi: 10.1002/app.1987.070330525. + + """ + if np.any(extent) <= 0: + raise ValueError("Be careful, for an autocatalytic model the extent can't be inferior or equal to 0 !!! \n It would result in a rate of reaction equal to 0 and no reaction would occur.") + return arrhenius_rate_constant(T,A,Ea) * extent**m * (1 - extent)**n + def rate_for_kamal(extent,T,A1,E1,A2,E2,m,n): r""" @@ -80,15 +119,15 @@ def rate_for_kamal(extent,T,A1,E1,A2,E2,m,n): extent : ndarray Extent of reaction. T : ndarray - Temperature of reaction. + Temperature of reaction in Kelvin. A1 : float Pre-exponential factor of the regular reaction. E1 : float - Activation energy of the regular reaction. + Activation energy of the regular reaction in J/mol. A2 : float Pre-exponential factor of the autocatalyzed reaction. E2 : float - Activation energy of the autocatalyzed reaction. + Activation energy of the autocatalyzed reaction in J/mol. m : float Order of reaction for the autocatalyzed reaction. n : float @@ -153,45 +192,6 @@ def rate_for_kamal(extent,T,A1,E1,A2,E2,m,n): """ return (arrhenius_rate_constant(T,A1,E1) + (arrhenius_rate_constant(T,A2,E2) * extent**m)) * (1 - extent)**n -def rate_for_autocatalytic(extent,T,A,Ea,m,n): - r""" - Compute the rate of reaction for an autocatalytic reaction. - - Parameters - ---------- - extent : ndarray - Extent of reaction. - T : ndarray - Temperature of reaction. - A : float - Pre-exponential factor of the reaction. - Ea : float - Activation energy of the reaction. - m : float - Order of reaction for the autocatalyzed reaction. - n : float - Order of reaction for the regular reaction. - - Returns - ------- - rate : ndarray - Rate of reaction for an autocatalytic equation. - - Notes - ----- - The rate for an autocatalytic reaction is given by: - - .. math:: \frac{d\alpha}{dt} = A e^{ \left( \frac{-E_a}{RT} \right)} \alpha^m (1-\alpha)^n - - References - ---------- - [1] M. R. Keenan, « Autocatalytic cure kinetics from DSC measurements: Zero initial cure rate », J. Appl. Polym. Sci., vol. 33, nᵒ 5, p. 1725‑1734, avr. 1987, doi: 10.1002/app.1987.070330525. - - """ - if np.any(extent) <= 0: - raise ValueError("Be careful, for an autocatalytic model the extent can't be inferior or equal to 0 !!! \n It would result in a rate of reaction equal to 0 and no reaction would occur.") - return arrhenius_rate_constant(T,A,Ea) * extent**m * (1 - extent)**n - def vitrification_WLF_rate(T,Tg,Ad,C1,C2): r""" @@ -240,6 +240,7 @@ def vitrification_WLF_rate(T,Tg,Ad,C1,C2): #If so, the vitrification term is computed, else the rate is equal to 0 return Ad*np.exp( (C1*(T-Tg)) / (C2+abs(T-Tg))) + def vitrification_WLF_rate_no_reaction_below_Tg(T,Tg,Ad,C1,C2): r""" Compute the vitrification term for a WLF-like model when the reaction temperature is above Tg. @@ -288,16 +289,16 @@ def tg_diBennedetto(extent,Tg_0,Tg_inf,coeff): alpha : array-like conversion or extent of reaction Tg_0 : Float - glass transition temperature of unreacted material + glass transition temperature of unreacted material in Kelvin Tg_inf : Float - glass transition temperature of fully reacted material + glass transition temperature of fully reacted material in Kelvin coeff : Float ratio of the changes in isobaric heat capacities at Tg of the fully reacted material and of the initial unreacted material Returns ------- Tg : array-like - Glass transition temperature + Glass transition temperature in Kelvin Notes ----- @@ -314,6 +315,7 @@ def tg_diBennedetto(extent,Tg_0,Tg_inf,coeff): """ return Tg_0 + (Tg_inf - Tg_0)*((coeff*extent)/(1-(1-coeff)*extent)) + def coupling_harmonic_mean(kc,kv,experimental_parameters=None): r""" Return the harmonic mean of a chemical rate and vitrification rate. @@ -347,6 +349,7 @@ def coupling_harmonic_mean(kc,kv,experimental_parameters=None): return rate + def coupling_product(kc,kv,experimental_parameters=None): r""" Return the product of a chemical rate and vitrification rate. @@ -369,10 +372,7 @@ def coupling_product(kc,kv,experimental_parameters=None): Rate : 1-D array Rate of reaction """ - #If kc or kv are equal to 0, then it returns 0 - rate = 1 / ((1/kc)+(1/kv)) - - return rate + return kc*kv def jac_for_rate_for_kamal(extent, T, A1, E1, A2, E2, m, n): @@ -540,7 +540,7 @@ def compute_extent_and_rate(time, temperature, rate_law=None, rate_law_args=None time : array-like List or array containing the times during the reaction. temperature : array-like - List or array containing the temperatures during the reaction. + List or array containing the temperatures (in Kelvin) during the reaction. rate_law : function The rate law function that calculates the rate of reaction. rate_law_args : tuple @@ -571,7 +571,7 @@ def compute_extent_and_rate(time, temperature, rate_law=None, rate_law_args=None vitrification_term : ndarray, optional Evolution of the vitrification rate of reaction (if vitrification parameters are provided). tg : ndarray, optional - Evolution of the Tg (if Tg parameters are provided). + Evolution of the Tg in Kelvin (if Tg parameters are provided). """ # Importation of a module to display a progress bar for the finite difference from tqdm import tqdm @@ -661,93 +661,6 @@ def compute_extent_and_rate(time, temperature, rate_law=None, rate_law_args=None - -def compute_extent_and_rate_using_scipy(t0, tf, temperature_program, rate_law, rate_law_args, vitrification_law=None, vitrification_args=None, tg_law=None, tg_law_args=None, coupling_law=None, coupling_law_args=None, initial_extent=0): - """ - Compute extent and rate using scipy's solve_ivp. - - Parameters - ---------- - t0 : float - Initial time. - tf : float - Final time. - temperature_program : callable - A function that takes time as input and returns temperature. - rate_law : callable - Rate law function. - rate_law_args : tuple - Arguments for the rate law function. - vitrification_law : callable - Vitrification law function. - vitrification_args : tuple - Arguments for the vitrification law function. - tg_law : callable - Tg law function. - tg_law_args : tuple - Arguments for the Tg law function. - coupling_law : callable - Coupling law function. - coupling_law_args : tuple - Arguments for the coupling law function. - initial_extent : float, optional - Initial extent. Default is 0. - - Returns - ------- - scipy.integrate.OdeResult - Solution object from solve_ivp. - """ - if (vitrification_law is not None) and (tg_law is None): - raise ValueError("Please make sure to indicate a tg law since you have selected a vitrification law.") - if (vitrification_law is not None) and (coupling_law is None): - raise ValueError("Please make sure to indicate a tg law since you have selected a vitrification law.") - if (coupling_law is not None) and (vitrification_law is None): - raise ValueError("Please make sure to indicate a tg law since you have selected a vitrification law.") - if (coupling_law is not None) and (tg_law is None): - raise ValueError("Please make sure to indicate a tg law since you have selected a vitrification law.") - if (tg_law is not None) and (vitrification_law is None): - raise ValueError("Please make sure to indicate a tg law since you have selected a vitrification law.") - if (tg_law is not None) and (coupling_law is None): - raise ValueError("Please make sure to indicate a tg law since you have selected a vitrification law.") - - def rate(t, extent, temperature, *args): - """ - Compute the rate at a given time. - - Parameters - ---------- - t : float - Time. - extent : float - Current extent. - temperature : float - Current temperature. - *args : tuple - Additional arguments. - - Returns - ------- - float - Rate. - """ - rate = rate_law(extent, temperature_program(t), *rate_law_args) - - if vitrification_law is not None: - tg = tg_law(extent, *tg_law_args) - vitrification_term = vitrification_law(temperature_program(t), tg, *vitrification_args) - rate = coupling_law(rate, vitrification_term, *coupling_law_args) - - return rate - - initial_extent = np.atleast_1d(initial_extent) - return scipy.integrate.solve_ivp(rate, [t0, tf], initial_extent, - args=(temperature_program, *rate_law_args, - vitrification_law, *vitrification_args, - tg_law, *tg_law_args, - coupling_law, *coupling_law_args)) - - #%% Example 1 - Kamal if __name__=="__main__": diff --git a/kinopt/src/kinopt_interface.py b/kinopt/src/kinopt_interface.py index 0763452..b48b617 100644 --- a/kinopt/src/kinopt_interface.py +++ b/kinopt/src/kinopt_interface.py @@ -285,7 +285,7 @@ def setupUi(self, MainWindow): MainWindow.setStatusBar(self.statusbar) self.dockWidgetData_extraction = QDockWidget(MainWindow) self.dockWidgetData_extraction.setObjectName(u"dockWidgetData_extraction") - self.dockWidgetData_extraction.setMinimumSize(QSize(352, 281)) + self.dockWidgetData_extraction.setMinimumSize(QSize(356, 281)) self.dockWidgetData_extraction.setFeatures(QDockWidget.AllDockWidgetFeatures) self.dockWidgetOpenFileContents = QWidget() self.dockWidgetOpenFileContents.setObjectName(u"dockWidgetOpenFileContents") @@ -477,7 +477,7 @@ def setupUi(self, MainWindow): self.toolBox_models_and_optimization_parameters.setFrameShape(QFrame.NoFrame) self.page_models_and_parameters = QWidget() self.page_models_and_parameters.setObjectName(u"page_models_and_parameters") - self.page_models_and_parameters.setGeometry(QRect(0, 0, 271, 747)) + self.page_models_and_parameters.setGeometry(QRect(0, 0, 214, 747)) self.verticalLayout_2 = QVBoxLayout(self.page_models_and_parameters) self.verticalLayout_2.setObjectName(u"verticalLayout_2") self.scrollArea_models_and_parameters = QScrollArea(self.page_models_and_parameters) @@ -487,7 +487,7 @@ def setupUi(self, MainWindow): self.scrollAreaWidgetContents_models_and_parameters = QWidget() self.scrollAreaWidgetContents_models_and_parameters.setObjectName(u"scrollAreaWidgetContents_models_and_parameters") self.scrollAreaWidgetContents_models_and_parameters.setEnabled(True) - self.scrollAreaWidgetContents_models_and_parameters.setGeometry(QRect(0, 0, 251, 680)) + self.scrollAreaWidgetContents_models_and_parameters.setGeometry(QRect(0, 0, 194, 680)) self.verticalLayout_3 = QVBoxLayout(self.scrollAreaWidgetContents_models_and_parameters) self.verticalLayout_3.setObjectName(u"verticalLayout_3") self.label_rate_model = QLabel(self.scrollAreaWidgetContents_models_and_parameters) @@ -578,7 +578,7 @@ def setupUi(self, MainWindow): self.toolBox_models_and_optimization_parameters.addItem(self.page_models_and_parameters, u"Models and parameters") self.page_optimization_methods_and_parameters = QWidget() self.page_optimization_methods_and_parameters.setObjectName(u"page_optimization_methods_and_parameters") - self.page_optimization_methods_and_parameters.setGeometry(QRect(0, 0, 271, 747)) + self.page_optimization_methods_and_parameters.setGeometry(QRect(0, 0, 214, 747)) self.verticalLayout_5 = QVBoxLayout(self.page_optimization_methods_and_parameters) self.verticalLayout_5.setObjectName(u"verticalLayout_5") self.scrollArea_optimization_methods = QScrollArea(self.page_optimization_methods_and_parameters) @@ -588,7 +588,7 @@ def setupUi(self, MainWindow): self.scrollArea_optimization_methods.setWidgetResizable(True) self.scrollAreaWidgetContents_optimization_methods = QWidget() self.scrollAreaWidgetContents_optimization_methods.setObjectName(u"scrollAreaWidgetContents_optimization_methods") - self.scrollAreaWidgetContents_optimization_methods.setGeometry(QRect(0, 0, 251, 680)) + self.scrollAreaWidgetContents_optimization_methods.setGeometry(QRect(0, 0, 194, 680)) self.verticalLayout_4 = QVBoxLayout(self.scrollAreaWidgetContents_optimization_methods) self.verticalLayout_4.setObjectName(u"verticalLayout_4") self.label_global_optimization_methods = QLabel(self.scrollAreaWidgetContents_optimization_methods) diff --git a/kinopt/src/main.py b/kinopt/src/main.py index 857eab2..2c4575e 100644 --- a/kinopt/src/main.py +++ b/kinopt/src/main.py @@ -1482,7 +1482,9 @@ def launch_optimization(self): self.experimental_args_for_cost_function, self.max_iter ) - self.ui.pushButton_launch_optimization.setEnabled(False) + self.ui.pushButton_launch_optimization.setText("Cancel optimization") + self.ui.pushButton_launch_optimization.clicked.disconnect() + self.ui.pushButton_launch_optimization.clicked.connect(self.cancel_optimization) self.optimization_thread.start() self.optimization_thread.update_progress_bar_signal.connect(self.update_progress_bar) self.optimization_thread.update_graph_signal.connect(self.update_graph) @@ -1494,10 +1496,42 @@ def launch_optimization(self): QMessageBox.critical(self,"Error", f"An error occurred: {str(e)}") self.optimization_thread.terminate() self.ui.pushButton_launch_optimization.setEnabled(True) - self.ui.progressBar.setValue(100) + self.ui.progressBar.setValue(0) self.ui.label_remaing_time.setText("Remaining time: (No optimization runnning)") return - + + def cancel_optimization(self): + """ + Cancel the ongoing optimization process. + + This function stops the optimization thread and resets the GUI elements related to the optimization. + + Parameters + ---------- + self : object + The object instance. + + Returns + ------- + None + """ + try: + if hasattr(self, 'optimization_thread') and self.optimization_thread.isRunning(): + self.ui.pushButton_launch_optimization.setText("Start optimization") + self.ui.pushButton_launch_optimization.clicked.disconnect() + self.ui.pushButton_launch_optimization.clicked.connect(self.launch_optimization) + self.optimization_thread.terminate() + self.ui.pushButton_launch_optimization.setEnabled(True) + self.ui.progressBar.setValue(0) + self.ui.label_remaing_time.setText("Remaining time: (No optimization runnning)") + self.ui.textEdit_output_of_optimization.insertPlainText("Optimization cancelled.\n") + self.ui.textEdit_output_of_optimization.moveCursor(QTextCursor.End) + self.ax_optimization.clear() + self.canvas_optimization.draw_idle() + except Exception as e: + # Handle other exceptions with a generic error message + QMessageBox.critical(self,"Error", f"An error occurred: {str(e)}") + def update_progress_bar(self,progress): """ Update the progress bar with the given progress value. diff --git a/kinopt/tests/test_interpolation.py b/kinopt/tests/test_interpolation.py new file mode 100644 index 0000000..cd9c249 --- /dev/null +++ b/kinopt/tests/test_interpolation.py @@ -0,0 +1,86 @@ +import pytest +import numpy as np +from kinopt.src import interpolation as interp + +def test_linear_interpolation(): + """ + Test the linear_interpolation function to ensure it correctly interpolates + conversion, time, temperature, and rate arrays to a specified number of points. + + The test sets up a simple case with known input arrays for conversions, times, + temperatures, and rates, and specifies 9 desired interpolation points. It then + verifies that the output arrays match the expected interpolated values using + np.allclose for numerical comparison. + """ + conversions = [np.array([0.0, 0.25, 0.5, 0.75, 1.0])] + times = [np.array([0.0, 2.5, 5.0, 7.5, 10.0])] + temperatures = [np.array([300, 325, 350, 350, 300])] + rates = [np.array([0.0, 0.25, 0.5, 0.25, 0.0])] + + num_points = 9 + + expected_conversions = np.linspace(0.0, 1, num_points) + expected_times = np.linspace(0.0, 10.0, num_points) + expected_temperatures = [300, 312.5, 325, 337.5, 350, 350, 350, 325, 300] + expected_rates = [0.0, 0.125, 0.25, 0.375, 0.5, 0.375, 0.25, 0.125, 0.0] + + new_conversions, new_times, new_temperatures, new_rates = interp.linear_interpolation( + conversions, times, temperatures, rates, num_points + ) + + assert np.allclose(new_conversions[0], expected_conversions) + assert np.allclose(new_times[0], expected_times) + assert np.allclose(new_temperatures[0], expected_temperatures) + assert np.allclose(new_rates[0], expected_rates) + +def test_linear_interpolation_multiple_limits(): + """ + Test the linear_interpolation_multiple_limits function to ensure it correctly interpolates + conversion, time, temperature, and rate arrays to a specified number of points. + + + """ + conversions = [np.array([0.0, 0.2, 0.4, 0.6, 0.7, 0.9]), + np.array([0, 0.2, 0.4, 0.6, 0.9]), + np.array([0.1, 0.4, 0.6, 0.8, 1.0])] + + times = [np.array([0.0, 2.0, 4.0, 6.0, 8.0, 10.0]), + np.array([0, 5, 15, 30, 45]), + np.array([1.0, 2.0, 4.0, 7.0, 20])] + + temperatures = [np.array([300, 340, 380, 400, 340, 300]), + np.array([300, 340, 360, 360, 300]), + np.array([300, 360, 360, 340, 300])] + + rates = [np.array([0.0, 0.2, 0.6, 0.4, 0.2, 0.0]), + np.array([0.0, 0.4, 0.6, 0.3, 0.0]), + np.array([0.0, 0.3, 0.4, 0.6, 0.1]),] + + num_points = 10 + + # conversion is an increasing function so np.linspace is used + expected_conversions = [np.linspace(0.0, 0.9, num_points), + np.linspace(0.0, 0.9, num_points), + np.linspace(0.1, 1.0, num_points)] + + # time is an increasing function so np.linspace is used + expected_times = [[0, 1, 2, 3, 4, 5, 6, 8, 9, 10], + [0, 2.5, 5, 10, 15, 22.5, 30, 35, 40, 45], + [1, 4/3, 5/3, 2, 3, 4, 5.5, 7, 13.5, 20]] + # trivial solution to interpolation + expected_temperatures = [[300, 320, 340, 360, 380, 390, 400, 340, 320, 300], + [300, 320, 340, 350, 360, 360, 360, 340, 320, 300], + [300, 320, 340, 360, 360, 360, 350, 340, 320, 300]] + # trivial solution to interpolation + expected_rates = [[0, 0.1, 0.2, 0.4, 0.6, 0.5, 0.4, 0.2, 0.1 , 0], + [0, 0.2 , 0.4, 0.5, 0.6, 0.45, 0.3, 0.2, 0.1, 0], + [0,0.1,0.2, 0.3, 0.35, 0.4, 0.5, 0.6, 0.35, 0.1]] + + new_conversions, new_times, new_temperatures, new_rates = interp.linear_interpolation_multiple_limits( + conversions, times, temperatures, rates, num_points + ) + + assert np.allclose(new_conversions, expected_conversions) + assert np.allclose(new_times, expected_times) + assert np.allclose(new_temperatures, expected_temperatures) + assert np.allclose(new_rates, expected_rates) \ No newline at end of file diff --git a/kinopt/tests/test_kinetic_models.py b/kinopt/tests/test_kinetic_models.py index 090cc95..5fe34c8 100644 --- a/kinopt/tests/test_kinetic_models.py +++ b/kinopt/tests/test_kinetic_models.py @@ -6,43 +6,173 @@ """ import pytest import numpy as np -import src.kinetic_models as km +import inspect +from kinopt.src import kinetic_models as km + +def test_arrhenius_rate_constant(): + A = 1.0 # Pre-exponential factor + Ea = 831.446261815324 # Activation energy in J/mol + T = 100.0 # Temperature in Kelvin -def test_vitrification_WLF_rate_no_reaction_below_Tg(): + expected_result = np.exp(-1) + result = km.arrhenius_rate_constant(T, A, Ea) + assert np.allclose(result, expected_result), "Arrhenius rate constant calculation failed." + +def test_rate_for_nth_order(): + extent = np.array([0, 0.1, 0.3, 0.5, 0.6, 1]) + T = np.array([100.0, 200.0, 300.0, 400.0, 500.0, 600.0]) + A = 1.0 # Pre-exponential factor + Ea = 831.446261815324 # Activation energy in J/mol + n = 2 # Order of the reaction + + expected_results = [np.exp(-1), ((0.9)**2)*np.exp(-1/2), ((0.7)**2)*np.exp(-1/3), ((0.5)**2)*np.exp(-1/4), ((0.4)**2)*np.exp(-1/5), 0] + results = km.rate_for_nth_order(extent, T, A, Ea, n) + assert np.allclose(results, expected_results), "Rate for nth order reaction calculation failed." + +def test_rate_for_autocatalytic(): + extent = np.array([0, 0.1, 0.3, 0.5, 0.6, 1]) + T = np.array([100.0, 200.0, 300.0, 400.0, 500.0, 600.0]) + A = 1.0 # Pre-exponential factor + Ea = 831.446261815324 # Activation energy in J/mol + n = 2 # Order of the reaction + m = 0.3 # Autocatalytic order + + expected_results = [0, (0.1**0.3)*((0.9)**2)*np.exp(-1/2), (0.3**0.3)*((0.7)**2)*np.exp(-1/3), (0.5**0.3)*((0.5)**2)*np.exp(-1/4), (0.6**0.3)*((0.4)**2)*np.exp(-1/5), 0] + results = km.rate_for_autocatalytic(extent, T, A, Ea, m, n) + assert np.allclose(results, expected_results), "Rate for nth order reaction calculation failed." + +def test_rate_for_kamal(): + extent = np.array([0, 0.1, 0.3, 0.5, 0.6, 1]) + T = np.array([100.0, 200.0, 300.0, 400.0, 500.0, 600.0]) + A1 = 1.0 # Pre-exponential factor + E1 = 831.446261815324 # Activation energy in J/mol + A2 = 2 # Pre-exponential factor for second reaction + E2 = 415.723130907662 # Activation energy for second reaction in J/mol + n = 2 # Order of the reaction + m = 0.3 # Autocatalytic order + + expected_results = [np.exp(-1), + (np.exp(-1/2)+2*np.exp(-1/4)*0.1**0.3)*(0.9**2), + (np.exp(-1/3)+2*np.exp(-1/6)*0.3**0.3)*(0.7**2), + (np.exp(-1/4)+2*np.exp(-1/8)*0.5**0.3)*(0.5**2), + (np.exp(-1/5)+2*np.exp(-1/10)*0.6**0.3)*(0.4**2), + 0] + results = km.rate_for_kamal(extent,T,A1,E1,A2,E2,m,n) + assert np.allclose(results, expected_results), "Rate for kamal reaction calculation failed." + + +def test_vitrification_WLF_rate(): # Test case 1: All temperatures above Tg temperature = np.array([350, 400, 450]) Ad = 1.0 - C1 = 0.5 - C2 = 0.2 + C1 = 1.0 + C2 = 10 Tg = 320 - expected_result = np.array([1.64327096, 1.64666679, 1.64745546]) - result = km.vitrification_WLF_rate_no_reaction_below_Tg(temperature, Ad, C1, C2, Tg) - assert np.allclose(result, expected_result), "Test case 1 failed" + expected_result = np.array([np.exp(3/4), np.exp(8/9), np.exp(13/14)]) + result = km.vitrification_WLF_rate(temperature, Tg, Ad, C1, C2) + assert np.allclose(result, expected_result), "Test case 1 failed. When the reaction temperature is above the glass transition temperature, this vitrification rate should return a non-zero value." # Test case 2: All temperatures below Tg temperature = np.array([280, 300, 310]) Ad = 1.0 - C1 = 0.5 - C2 = 0.2 + C1 = 1.0 + C2 = 10 Tg = 320 - expected_result = np.array([0, 0, 0]) - result = km.vitrification_WLF_rate_no_reaction_below_Tg(temperature, Ad, C1, C2, Tg) - assert np.allclose(result, expected_result), "Test case 2 failed" + expected_result = np.array([np.exp(-4/5), np.exp(-2/3), np.exp(-1/2)]) + result = km.vitrification_WLF_rate(temperature, Tg, Ad, C1, C2) + assert np.allclose(result, expected_result), "Test case 2 failed. When the reaction temperature is below the glass transition temperature, this vitrification rate should return a rate equal to 0." - # Test case 3: Large values for Ad, C1, and C2 + +def test_vitrification_WLF_rate_no_reaction_below_Tg(): + # Test case 1: All temperatures above Tg temperature = np.array([350, 400, 450]) - Ad = 1e6 - C1 = 1e6 - C2 = 1e6 + Ad = 1.0 + C1 = 1.0 + C2 = 10 Tg = 320 - expected_result = np.array([1.06768614e+19, 5.50527845e+40, 2.82451591e+62]) - result = km.vitrification_WLF_rate_no_reaction_below_Tg(temperature, Ad, C1, C2, Tg) - assert np.allclose(result, expected_result), "Test case 4 failed" + expected_result = np.array([np.exp(3/4), np.exp(8/9), np.exp(13/14)]) + result = km.vitrification_WLF_rate_no_reaction_below_Tg(temperature, Tg, Ad, C1, C2) + assert np.allclose(result, expected_result), "Test case 1 failed. When the reaction temperature is above the glass transition temperature, this vitrification rate should return a non-zero value." + + # Test case 2: All temperatures below Tg + temperature = np.array([280, 300, 310]) + Ad = 1.0 + C1 = 1.0 + C2 = 10 + Tg = 320 + expected_result = np.array([0, 0, 0]) + result = km.vitrification_WLF_rate_no_reaction_below_Tg(temperature, Tg, Ad, C1, C2) + assert np.allclose(result, expected_result), "Test case 2 failed. When the reaction temperature is below the glass transition temperature, this vitrification rate should return a rate equal to 0." + + +def test_tg_diBennedetto(): + Tg_0 = -100 + Tg_inf = 100 + coeff = 0.5 + extent = np.array([0, 0.1, 0.3, 0.5, 0.6, 1]) + + expected_result = np.array([-100, -100+200*1/19, -100+200*3/17, -100+200*1/3, -100+200*3/7, 100]) + result = km.tg_diBennedetto(extent, Tg_0, Tg_inf, coeff) + assert np.allclose(result, expected_result), "DiBenedetto Tg calculation failed." + + +def test_coupling_harmonic_mean(): + kc = 2 + kv = 3 + + expected_result = 1.2 + result = km.coupling_harmonic_mean(kc, kv) + assert np.allclose(result, expected_result), "Harmonic mean coupling calculation failed." + + +def test_coupling_product(): + kc = 2 + kv = 3 + + expected_result = 6 + result = km.coupling_product(kc, kv) + assert np.allclose(result, expected_result), "Product coupling calculation failed." + + +def test_rate_functions_signature(): + for name, func in inspect.getmembers(km, inspect.isfunction): + if name.startswith("rate_"): + sig = inspect.signature(func) + params = list(sig.parameters.values()) + assert len(params) >= 2, f"{name} should have at least 2 parameters" + assert params[0].name == "extent", f"{name}: first argument should be 'extent'" + assert params[1].name == "T", f"{name}: second argument should be 'T'" + +def test_vitrification_functions_signature(): + for name, func in inspect.getmembers(km, inspect.isfunction): + if name.startswith("vitrification_"): + sig = inspect.signature(func) + params = list(sig.parameters.values()) + assert len(params) >= 2, f"{name} should have at least 2 parameters" + assert params[0].name == "T", f"{name}: first argument should be 'T'" + assert params[1].name == "Tg", f"{name}: second argument should be 'Tg'" + +def test_tg_functions_signature(): + for name, func in inspect.getmembers(km, inspect.isfunction): + if name.startswith("tg_"): + sig = inspect.signature(func) + params = list(sig.parameters.values()) + assert len(params) >= 1, f"{name} should have at least 1 parameter" + assert params[0].name == "extent", f"{name}: first argument should be 'extent'" - print("All test cases passed!") +print("All test cases passed!") +def test_coupling_functions_signature(): + for name, func in inspect.getmembers(km, inspect.isfunction): + if name.startswith("coupling_"): + sig = inspect.signature(func) + params = list(sig.parameters.values()) + assert len(params) >= 3, f"{name} should have at least 3 parameters: kc, kv, and experimental_parameters (that can be set to 'None' if no experimental parameter is required for the coupling) " + assert params[0].name == "kc", f"{name}: first argument should be 'kc' (purely chemical rate)" + assert params[1].name == "kv", f"{name}: second argument should be 'kv' (vitrification rate)" + if __name__=="__main__": diff --git a/kinopt/tests/test_optimization.py b/kinopt/tests/test_optimization.py index ad4067a..b453292 100644 --- a/kinopt/tests/test_optimization.py +++ b/kinopt/tests/test_optimization.py @@ -8,42 +8,21 @@ import pytest import inspect import numpy as np -import src.optimization as opt - -# Define some sample data for testing -experimental_rate = np.array([0.5, 0.7, 1.0, 1.5]) -x = np.array([0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, ]) -fraction_to_amplify = 0.5 -amplification_factor = 2.0 - -# Define some mock rate laws for testing -def mock_rate_law(*args): - return np.array([0.3, 0.6, 0.9, 1.2]) - -def mock_vitrification_law(*args): - return np.array([0.1, 0.2, 0.3, 0.4]) - -def mock_tg_law(*args): - return 100.0 - -def mock_coupling_law(*args): - return np.array([0.2, 0.3, 0.4, 0.5]) - -# Test rss_increase_of_small_rates_impact_with_zones function -def test_rss_increase_of_small_rates_impact_with_zones(): - rss = opt.rss_increase_of_small_rates_impact_with_zones(x, experimental_rate, mock_rate_law, (), 0, mock_vitrification_law, (), 0, mock_coupling_law, (), mock_tg_law, (), (), fraction_to_amplify, amplification_factor) - assert np.isclose(rss, 0.025, atol=1e-3) - -# Test rss_increase_of_small_extents_impact function -def test_rss_increase_of_small_extents_impact(): - extent = np.array([0.05, 0.1, 0.15, 0.2]) - extent_limit = 0.1 - rss = opt.rss_increase_of_small_extents_impact(x, experimental_rate, mock_rate_law, (), 0, mock_vitrification_law, (), 0, mock_coupling_law, (), mock_tg_law, (), (), extent, extent_limit, amplification_factor) - assert np.isclose(rss, 0.0325, atol=1e-3) +from kinopt.src import optimization as opt def test_rss_functions_arguments(): - rss_functions = [func for func in locals().values() if inspect.isfunction(func) and func.__name__.startswith('rss')] + """ + Test that all RSS functions have the correct argument names and order. + """ + # Get all RSS functions from the optimization module + rss_functions = [func for func in dict(inspect.getmembers(opt, inspect.isfunction)).values() + if func.__name__.startswith('rss')] + + # Check that we found some functions + assert len(rss_functions) > 0, "No RSS functions found in optimization module" + + # Define the expected argument names in the correct order expected_arguments = [ 'x', 'experimental_rate', @@ -59,6 +38,145 @@ def test_rss_functions_arguments(): 'experimental_args_for_tg', 'tg_args' ] + for func in rss_functions: - actual_arguments = inspect.signature(func).parameters.keys() - assert list(actual_arguments) == expected_arguments, f"Function {func.__name__} has incorrect arguments: {list(actual_arguments)}" \ No newline at end of file + # Get actual parameter names while preserving their order + params = list(inspect.signature(func).parameters.items()) + actual_arguments = [param[0] for param in params[:13]] + + # Compare each parameter name and position + for expected, actual in zip(expected_arguments, actual_arguments): + assert expected == actual, f"In function {func.__name__}, expected parameter '{expected}' but got '{actual}'" + +def test_rss_standard_with_known_function(): + """ + Test RSS standard calculation using a simple quadratic function with carefully chosen values: + - Using just 3 time points: t = [0, 1, 2] + - True function: f(t) = t² + t + 1 (a=1, b=1, c=1) + - Test function: f(t) = 2t² + t + 1 (a=2, b=1, c=1) + + This gives us: + t=0: diff = 0 → (0)² = 0 + t=1: diff = 1 → (1)² = 1 + t=2: diff = 4 → (4)² = 16 + + Total RSS = 0 + 1 + 16 = 17 + """ + def simple_rate_law(t, temp, a, b, c, *args): + t = np.asarray(t) + return a * t**2 + b * t + c + + t_data = np.array([0, 1, 2]) + temp_data = 298 * np.ones_like(t_data) + true_params = [1.0, 1.0, 1.0] + experimental_rate = simple_rate_law(t_data, temp_data, *true_params) + x_test = [2.0, 1.0, 1.0] + + rss = opt.rss_standard(x_test, experimental_rate, simple_rate_law, + (t_data, temp_data), 3, + None, None, 0, None, None, None, None, None) + + assert np.isclose(rss, 17.0, rtol=1e-10) + +def test_rss_mean_with_known_function(): + """ + Test RSS mean calculation with same function but divided by n points: + Total RSS = (0 + 1 + 16) / 3 = 5.6667 + """ + def simple_rate_law(t, temp, a, b, c, *args): + t = np.asarray(t) + return a * t**2 + b * t + c + + t_data = np.array([0, 1, 2]) + temp_data = 298 * np.ones_like(t_data) + true_params = [1.0, 1.0, 1.0] + experimental_rate = simple_rate_law(t_data, temp_data, *true_params) + x_test = [2.0, 1.0, 1.0] + + rss = opt.rss_mean(x_test, experimental_rate, simple_rate_law, + (t_data, temp_data), 3, + None, None, 0, None, None, None, None, None) + + assert np.isclose(rss, 17.0/3, rtol=1e-10) + +def test_rss_relative_with_known_function(): + """ + Test RSS relative calculation with simple values: + - True values: [1, 2, 4] + - Test values: [0.9, 1.8, 3.6] + + Relative differences: (0.9-1)/1 = -0.1, (1.8-2)/2 = -0.1, (3.6-4)/4 = -0.1 + RSS = (-0.1)² + (-0.1)² + (-0.1)² = 0.03 + """ + def simple_rate_law(t, temp, k, *args): + return k * np.array([0.9, 1.8, 3.6]) + + t_data = np.array([1.0, 2.0, 3.0]) + temp_data = 298 * np.ones_like(t_data) + experimental_rate = np.array([1.0, 2.0, 4.0]) + + rss = opt.rss_relative([1.0], experimental_rate, simple_rate_law, + (t_data, temp_data), 1, + None, None, 0, None, None, None, None, None) + + assert np.isclose(rss, 0.03, rtol=1e-10) + +def test_rss_small_extents_impact(): + """ + Test RSS with small extents impact: + - Experimental rate: [1, 1, 1] + - Model rate: [0.9, 0.9, 0.9] + - Extents: [0.1, 0.5, 0.9] + - Extent limit: 0.3 + - Amplification factor: 2 + + Differences: all are 0.1 + First point amplified (extent < 0.3): (0.1 * 2)² = 0.04 + Other points normal: (0.1)² + (0.1)² = 0.02 + Mean RSS = (0.04 + 0.02) / 3 = 0.02 + """ + def simple_rate_law(t, temp, k, *args): + return k * np.array([0.9, 0.9, 0.9]) + + t_data = np.array([1.0, 2.0, 3.0]) + temp_data = 298 * np.ones_like(t_data) + experimental_rate = np.array([1.0, 1.0, 1.0]) + extent = np.array([0.1, 0.5, 0.9]) + + rss = opt.rss_increase_of_small_extents_impact( + [1.0], experimental_rate, simple_rate_law, + (t_data, temp_data), 1, + None, None, 0, None, None, None, None, None, + extent, 0.3, 2.0 + ) + + assert np.isclose(rss, 0.02, rtol=1e-10) + +def test_rss_small_rates_zones(): + """ + Test RSS with zones: + - Experimental rate: [0.1, 1.0, 2.0] + - Model rate: [0.2, 1.1, 2.1] + - Max rate = 2.0 + - Fraction to amplify = 4 (threshold = 0.5) + - Amplification factor = 2 + + First point amplified (rate < max/4): ((0.2-0.1)*2)² = 0.04 + Other points normal: (0.1)² + (0.1)² = 0.02 + Mean RSS = (0.04 + 0.02) / 3 = 0.02 + """ + def simple_rate_law(t, temp, k, *args): + return k * np.array([0.2, 1.1, 2.1]) + + t_data = np.array([1.0, 2.0, 3.0]) + temp_data = 298 * np.ones_like(t_data) + experimental_rate = np.array([0.1, 1.0, 2.0]) + + rss = opt.rss_increase_of_small_rates_impact_with_zones( + [1.0], experimental_rate, simple_rate_law, + (t_data, temp_data), 1, + None, None, 0, None, None, None, None, None, + 4.0, 2.0 + ) + + assert np.isclose(rss, 0.02, rtol=1e-10) \ No newline at end of file diff --git a/requirements.txt b/requirements.txt index 8868c32..7d771ba 100644 Binary files a/requirements.txt and b/requirements.txt differ