Use SciPy orthogonal distance regression

appPy/src/api/fit_api.py

@@ -57,13 +57,12 @@ def fit_Vapor_api():
         'p_data': True,
         'T_data': True,
         'nparams0': False,
-        'skip_T_p_optimization': False
     }
     params = unpack_request_schema(request, params_schema)
 
     fit = Fit_Vapor_plot(*params.values())
     fit.optimize_p()
-    if not params['skip_T_p_optimization']: fit.optimize_T_p()
+    fit.optimize_T_p()
     payload = fit.serialize()
     fit.tabulate()
     payload['plot_p'] = fit.plot_p(mode='svg')

appPy/src/fit/Fit_Vapor.py

@@ -1,18 +1,19 @@
 import numpy as np
-from scipy.optimize import least_squares, root
+from scipy.optimize import least_squares
+from scipy.odr import Model, Data, ODR
 from src.config import cst
 from src.TD.vapor_models.supported_models import supported_models
-from src.utils.errors import AppException
 from src.utils.io.echo import echo, underline_echo
 from .Fit import Fit
 from .utils import RMS, AAD, const_param_wrappers
 
 
 class Fit_Vapor(Fit):
 
-    def __init__(self, compound, model_name, p_data, T_data, params0, skip_T_p_optimization=False):
+    def __init__(self, compound, model_name, p_data, T_data, params0):
         """
         Create non-linear regression problem for given vapor pressure datasets and a selected model.
+        Performs both normal non-linear regression, and orthogonal distance regression (ODR).
 
         compound (str): name of the compound.
         model_name (str): name of the model to be fitted.
@@ -22,84 +23,70 @@ def __init__(self, compound, model_name, p_data, T_data, params0, skip_T_p_optim
         skip_T_p_optimization (bool): if True, only p residuals will be optimized.
         """
         super().__init__(supported_models, model_name, params0)
-        self.keys_to_serialize.extend(['is_T_p_optimized', 'RMS_inter', 'AAD_inter', 'nparams_inter', 'T_min', 'T_max'])
+        self.keys_to_serialize.extend(
+            ['is_T_p_optimized', 'odr_messages', 'RMS_inter', 'AAD_inter', 'nparams_inter', 'T_min', 'T_max'])
         self.compound = compound
         self.p_data = np.array(p_data)
         self.T_data = np.array(T_data)
+        # will be offered as model T_min & T_max
         self.T_min = np.min(T_data)
         self.T_max = np.max(T_data)
 
-        # calculate weight factor for T residuals
-        p_span = np.max(self.p_data) - np.min(self.p_data)
-        T_span = self.T_max - self.T_min
-        self.T_wt_factor = p_span / T_span
-
         self.T_tab = self.p_tab_inter = self.p_tab_final = None
         self.nparams_inter = self.params_inter = None
         self.RMS_inter = self.AAD_inter = None
         self.is_T_p_optimized = False  # in case of Fit_Vapor, Fit.is_optimized means is_p_optimized
+        self.odr_messages = []  # scipy.ODR has an unusual output format, as array of messages, but no status
 
-        # in order to be consistent: if optimizing T,p, ALWAYS aim for T,p-residual, and if not, then always p-residual
-        self.resid_fun = self.get_p_residuals if skip_T_p_optimization else self.get_T_p_residuals
-        self.RMS_init = RMS(self.resid_fun(self.params0))
-        self.AAD_init = AAD(self.resid_fun(self.params0))
+        # even for ODR, only p residuals are calculated for simplicity, so they may actually increase after ODR
+        self.RMS_init = RMS(self.get_p_residuals(self.params0))
+        self.AAD_init = AAD(self.get_p_residuals(self.params0))
 
     def get_p_residuals(self, params):
-        """Calculate residuals for given params."""
-        # calculating p residuals is easy, as p is dependent variable
+        """With given params, calculate residuals only of the dependent variable (p)."""
         p_calc = self.model.fun(self.T_data, *params)
         return p_calc - self.p_data
 
-    def get_T_residuals(self, params):
-        # model is generally transcendental for T, so we need to solve to get T residuals
-        # as if we were calculating boiling point at a given pressure
-        T_calc = np.zeros(len(self.T_data))
-        for i in range(len(self.T_data)):
-            # pylint: disable=cell-var-from-loop
-            T_boil_resid = lambda T: self.model.fun(T, *params) - self.p_data[i]
-            sol = root(fun=T_boil_resid, x0=self.T_data[i], tol=cst.T_boil_tol)
-            if not sol.success:
-                raise AppException(f'Cannot get T,p residual – failed to find boiling point at {self.p_data[i]} kPa')
-            T_calc[i] = sol.x[0]
-        return T_calc - self.T_data
-
-    def get_T_p_residuals(self, params):
-        return np.concatenate((self.get_T_residuals(params) * self.T_wt_factor, self.get_p_residuals(params)))
-
     def optimize_p(self):
-        """Optimize the model params, considering only error of dependent variable (p). When T,p-optimization is skipped, we are done, and the inter results stay as final."""
+        """Optimize the model params, considering only error of dependent variable (p)."""
         var_params0, wrapped_fun, merge_params = const_param_wrappers(self.get_p_residuals, self.params0,
                                                                       self.const_param_names, self.model.param_names)
         result = least_squares(wrapped_fun, var_params0, method='lm')
-        if result.status <= 0: raise AppException(f'p-optimization failed! Status {result.status}: {result.message}')
+        if result.status <= 0:
+            self.warn(f'p-optimization failed! Status {result.status}: {result.message}')
+            return
         self.is_optimized = True
-        self.params_inter = self.params = merge_params(result.x)
-        self.nparams_inter = self.nparams = self.set_named_params(self.params_inter)
-        self.RMS_inter = self.RMS_final = RMS(self.resid_fun(self.params))
-        self.AAD_inter = self.AAD_final = AAD(self.resid_fun(self.params))
+        self.params_inter = merge_params(result.x)
+        self.nparams_inter = self.set_named_params(self.params_inter)
+        self.RMS_inter = RMS(self.get_p_residuals(self.params_inter))
+        self.AAD_inter = AAD(self.get_p_residuals(self.params_inter))
 
     def optimize_T_p(self):
-        """Optimize the model params, considering errors of both dependent (p) & independent (T) variables. Overwrite the final results."""
-        if not self.is_optimized: raise AppException('Optimization of p residuals must be performed first')
-        var_params_inter, wrapped_fun, merge_params = const_param_wrappers(self.get_T_p_residuals, self.params_inter,
-                                                                           self.const_param_names,
-                                                                           self.model.param_names)
+        """Optimize the model params, considering errors of both dependent (p) & independent (T) variables."""
+        full_params0 = self.params_inter if self.is_optimized else self.params  # use intermediate params if available
+        fun2wrap = lambda params, x: self.model.fun(x, *params)  # swap the order of params and x
+        var_params0, wrapped_fun, merge_params = const_param_wrappers(fun2wrap, full_params0, self.const_param_names,
+                                                                      self.model.param_names)
         try:
-            result = least_squares(wrapped_fun, var_params_inter, method='lm')
-            if result.status <= 0:
-                raise AppException(f'T,p-optimization failed! Status {result.status}: {result.message}')
-        except AppException as e:
+            odr_model = Model(wrapped_fun)
+            odr_data = Data(self.T_data, self.p_data)
+            odr = ODR(odr_data, odr_model, beta0=var_params0)
+            odr_result = odr.run()
+        except Exception as e:  # pylint: disable=broad-exception-caught
             self.warn(f'T,p-optimization failed with error: {e}')
             return
         self.is_T_p_optimized = True
-        self.params = merge_params(result.x)
+        self.params = merge_params(odr_result.beta)
         self.nparams = self.set_named_params(self.params)
-        self.RMS_final = RMS(self.resid_fun(self.params))
-        self.AAD_final = AAD(self.resid_fun(self.params))
+        self.RMS_final = RMS(self.get_p_residuals(self.params))
+        self.AAD_final = AAD(self.get_p_residuals(self.params))
+        self.odr_messages = odr_result.stopreason
+        self.odr_messages = self.odr_messages if isinstance(self.odr_messages, list) else [self.odr_messages]
 
     def tabulate(self):
         self.T_tab = np.linspace(self.T_min, self.T_max, cst.x_points_smooth_plot)
-        self.p_tab_inter = self.model.fun(self.T_tab, *self.params_inter)
+        if self.is_optimized:
+            self.p_tab_inter = self.model.fun(self.T_tab, *self.params_inter)
         if self.is_T_p_optimized:
             self.p_tab_final = self.model.fun(self.T_tab, *self.params)
 
@@ -111,8 +98,8 @@ def report(self):
         echo('')
 
         echo('Optimized parameters as intermediate & final value:')
-        for (name, p_inter, p) in zip(self.model.param_names, self.params_inter, self.params):
-            echo(f'  {name} = {p_inter:7.4g} {p:9.4g}')
+        for (name, p_inter, p_final) in zip(self.model.param_names, self.params_inter, self.params):
+            echo(f'  {name} = {p_inter:7.4g} {p_final:9.4g}')
 
         echo('')
         echo(f'Initial RMS = {self.RMS_init:.3g}')

appPy/src/fit/utils.py

@@ -41,7 +41,7 @@ def const_param_wrappers(fun, params0, const_param_names, model_param_names):
     merge_params = lambda new_params: overlay_vectors(const_params, const_param_idxs, new_params)
 
     # wrapped residual as function of variable model parameters, which are merged with const parameters
-    wrapped_fun = lambda var_params: fun(overlay_vectors(const_params, const_param_idxs, var_params))
+    wrapped_fun = lambda var_params, *args: fun(overlay_vectors(const_params, const_param_idxs, var_params), *args)
 
     # expose everything required to start and finish the problem
     return var_params0, wrapped_fun, merge_params

appPy/src/plot/Fit_Vapor_plot.py

@@ -3,15 +3,10 @@
 from src.utils.UoM import convert_T, convert_p
 from src.fit.Fit_Vapor import Fit_Vapor
 from src.plot.plot_io import init_plot, finish_plot
-from src.utils.errors import AppException
 
 
 class Fit_Vapor_plot(Fit_Vapor):
 
-    def __check_is_tabulated(self):
-        if self.T_tab is None or self.p_tab_inter is None:
-            raise AppException('No tabulated data to plot, must call tabulate() first!')
-
     def draw_data(self):
         T_data_disp = convert_T(self.T_data)
         p_data_disp = convert_p(self.p_data)
@@ -25,7 +20,7 @@ def decorate(self):
 
     def plot_p(self, mode):
         """Overlay p-fitted model over vapor pressure data."""
-        self.__check_is_tabulated()
+        if not self.is_optimized: return None
         init_plot(mode)
         self.draw_data()
 
@@ -38,8 +33,7 @@ def plot_p(self, mode):
 
     def plot_T_p(self, mode):
         """Overlay T,p-fitted model over vapor pressure data."""
-        if self.p_tab_final is None: return None
-        self.__check_is_tabulated()
+        if not self.is_T_p_optimized: return None
         init_plot(mode)
         self.draw_data()
 

appUI/src/actions/FitVapor/FitVaporResultsDialog.tsx

@@ -1,15 +1,17 @@
-import { Dispatch, FC, useMemo } from 'react';
+import { Dispatch, FC, useMemo, useState } from 'react';
 import Spreadsheet from 'react-spreadsheet';
-import { Button, DialogContent, Stack, styled } from '@mui/material';
+import { Box, Button, Collapse, DialogContent, Stack, styled } from '@mui/material';
 import { ResponsiveDialog } from '../../components/Mui/ResponsiveDialog.tsx';
 import { DialogTitleWithX } from '../../components/Mui/DialogTitle.tsx';
 import { DialogProps } from '../../adapters/types/DialogProps.ts';
 import { VaporFitRequest, VaporFitResponse } from '../../adapters/api/types/fitTypes.ts';
 import { AnalysisWarnings } from '../../components/AnalysisResults/AnalysisWarnings.tsx';
 import { PlotWithDownload } from '../../components/charts/PlotWithDownload.tsx';
+import { ExpandHelpButton } from '../../components/Mui/ExpandHelpButton.tsx';
 import { makeReadOnly, matrixToSpreadsheetData, spreadsheetToSigDgts } from '../../adapters/logic/spreadsheet.ts';
 import { fromNamedParams } from '../../adapters/logic/nparams.ts';
 import { sigDgtsMetrics, sigDgtsParams } from '../../adapters/logic/numbers.ts';
+import { FitVaporResultsHelp } from './FitVaporResultsHelp.tsx';
 
 const fitQualityMetrics = ['Root mean square', 'Average abs. deviation'];
 
@@ -30,6 +32,8 @@ export const FitVaporResultsDialog: FC<FitVaporResultsDialogProps> = ({
     data,
     setFitResults,
 }) => {
+    const [infoOpen, setInfoOpen] = useState(false);
+
     const optimizedP = data.is_optimized;
     const optimizedTP = data.is_T_p_optimized;
     const paramNames = useMemo(() => fromNamedParams(data.nparams0)[0], [data]);
@@ -80,27 +84,42 @@ export const FitVaporResultsDialog: FC<FitVaporResultsDialogProps> = ({
                 <h4 className="h-margin">Fitted model parameters</h4>
                 <Spreadsheet data={paramsSpreadsheetData} columnLabels={paramNames} rowLabels={rowLabels} />
                 <h4 className="h-margin">What to do with the results?</h4>
-                <Stack direction="row" gap={2} mb={4}>
+                <Stack direction="row" gap={2} mb={2}>
                     <NormalCaseButton variant="contained" color="error" onClick={handleClose}>
                         Reject
                     </NormalCaseButton>
-                    <NormalCaseButton variant="contained" onClick={acceptP}>
-                        Accept {optimizedTP && 'p-optimized'}
-                    </NormalCaseButton>
+                    {optimizedP && (
+                        <NormalCaseButton variant="contained" onClick={acceptP}>
+                            Accept {optimizedTP && 'p-optimized'}
+                        </NormalCaseButton>
+                    )}
                     {optimizedTP && (
                         <NormalCaseButton variant="contained" onClick={acceptTP}>
                             Accept T,p-optimized
                         </NormalCaseButton>
                     )}
+                    <ExpandHelpButton infoOpen={infoOpen} setInfoOpen={setInfoOpen} />
                 </Stack>
-                <h4 className="h-margin">p-optimized model plot</h4>
-                <PlotWithDownload svgContent={data.plot_p} fileName={`fit chart ${req.compound} ${req.model_name}`} />
+                <Box mb={3}>
+                    <Collapse in={infoOpen}>
+                        <FitVaporResultsHelp />
+                    </Collapse>
+                </Box>
+                {data.plot_p && (
+                    <>
+                        <h4 className="h-margin">p-optimized model plot</h4>
+                        <PlotWithDownload
+                            svgContent={data.plot_p}
+                            fileName={`p-fit chart ${req.compound} ${req.model_name}`}
+                        />
+                    </>
+                )}
                 {data.plot_T_p && (
                     <>
                         <h4 className="h-margin">T,p-optimized model plot</h4>
                         <PlotWithDownload
                             svgContent={data.plot_T_p}
-                            fileName={`fit chart ${req.compound} ${req.model_name}`}
+                            fileName={`T,p-fit chart ${req.compound} ${req.model_name}`}
                         />
                     </>
                 )}

appUI/src/actions/FitVapor/FitVaporResultsHelp.tsx

@@ -0,0 +1,29 @@
+import { FC } from 'react';
+import { Box } from '@mui/material';
+import { spacingN } from '../../contexts/MUITheme.tsx';
+
+export const FitVaporResultsHelp: FC = () => (
+    <Box ml={4}>
+        Two optimization methods are employed:
+        <ul style={{ marginBottom: 0, marginTop: spacingN(1) }}>
+            <li>
+                <u>p-optimization</u>: standard curve fitting, which only considers residuals of the dependent variable
+                (<i>p</i>).
+            </li>
+            <li>
+                <u>T,p-optimization</u>: orthogonal distance regression, which considers residuals of both variables (
+                <i>T</i>, <i>p</i>).
+                <ul>
+                    <li>
+                        This is <b>recommended</b>, as both variables have significant uncertainties.
+                    </li>
+                    <li>But it may be unstable, so both methods are offered.</li>
+                    <li>
+                        Note that RMS & AAD are always calculated only from <i>p</i> residuals, so they{' '}
+                        <i>may actually increase</i> with T,p-optimization.
+                    </li>
+                </ul>
+            </li>
+        </ul>
+    </Box>
+);

appUI/src/adapters/api/types/fitTypes.ts

@@ -74,18 +74,25 @@ export type VaporFitRequest = CompoundIdentifier & {
     p_data: number[];
     T_data: number[];
     nparams0?: NamedParams;
-    skip_T_p_optimization?: boolean;
 };
 
-export type VaporFitResponse = FitMetrics & {
+export type VaporFitResponse = AnalysisResult & {
+    is_optimized: boolean;
     is_T_p_optimized: boolean;
-    RMS_inter: number | null;
-    AAD_inter: number | null;
-    nparams0: NamedParams;
-    nparams_inter: NamedParams;
-    nparams: NamedParams;
+    odr_messages: string[];
     T_min: number;
     T_max: number;
-    plot_p: string;
+    plot_p?: string;
     plot_T_p?: string;
+
+    nparams0: NamedParams;
+    nparams_inter?: NamedParams;
+    nparams?: NamedParams;
+
+    RMS_init: number;
+    RMS_inter: number | null;
+    RMS_final: number | null;
+    AAD_init: number;
+    AAD_inter: number | null;
+    AAD_final: number | null;
 };

appUI/src/components/Mui/ExpandHelpButton.tsx

@@ -0,0 +1,15 @@
+import { Dispatch, FC, SetStateAction } from 'react';
+import { IconButton } from '@mui/material';
+import { HelpOutline, KeyboardArrowDown, KeyboardArrowUp } from '@mui/icons-material';
+
+export type ExpandHelpButtonProps = {
+    infoOpen: boolean;
+    setInfoOpen: Dispatch<SetStateAction<boolean>>;
+};
+
+export const ExpandHelpButton: FC<ExpandHelpButtonProps> = ({ infoOpen, setInfoOpen }) => (
+    <IconButton onClick={() => setInfoOpen((prev) => !prev)}>
+        <HelpOutline />
+        {infoOpen ? <KeyboardArrowUp /> : <KeyboardArrowDown />}
+    </IconButton>
+);