## 2 Materials and Methods  
### Non-Linear Least Square Fitting Algorithm  
Consider a set of data \((v_{j},p_{j}),j=1,\ldots,m\) and a function  
\[p_{x}(\
u):\mathbb{R}\longrightarrow\mathbb{R},x\in\mathbb{R}^{n}.\]  
The residuals \(r_{j}(x):\mathbb{R}^{n}\longrightarrow\mathbb{R}\) are vector-valued functions corresponding to the set of data \((v_{j},p_{j})\) possibly defined as follows:  
* \(r_{j}(x)=|p_{x}(v_{j})-p_{j}|\), \(j=1,\ldots,m\), (absolute error) or
* \(r_{j}(x)=\frac{|p_{x}(v_{j})-p_{j}|}{|p_{j}|}\), \(j=1,\ldots,m\) (relative error)Let \(R(x):\mathbb{R}^{n}\longrightarrow\mathbb{R}^{m}\) be  
\[R(x)=(r_{1}(x),r_{2}(x),\ldots,r_{m}(x)). \tag{1}\]  
Thus, \(R\) is a differentiable function, where \(n\) represents the number of unknowns (the dimension of \(x\)) and \(m\) is the number of data points, with \(m>n\). The general box-constrained nonlinear least squares (NLLS) problem is defined as  
\[\underset{l\leq x\leq u}{\operatorname{argmin}}f(x), \tag{2}\]  
where the vectors \(l,u\in\mathbb{R}^{n}\) represent the lower and upper bounds, respectively, with \(l_{i}<u_{i},i=1,\ldots,n\), and the objective function is defined by the Euclidean norm of \(R(x)\):  
\[f(x)=\frac{1}{2}(\|R(x)\|_{2})^{2}, \tag{3}\]  
whose gradient \(\
abla f\) is given by  
\[\
abla f(x)=J(x)^{T}\cdot R(x), \tag{4}\]  
with \(J\) being the Jacobian matrix of \(R(x)\). Given the constraints of the box-bounded domain, a diagonal scaling matrix \(D\) with non-zero elements \(d_{i}\) is defined as follows:  
\[d_{i}(x)= \begin{cases}x_{i}-u_{i}&\text{ if }\
abla_{i}f(x)<0\\ x_{i}-l_{i}&\text{ if }\
abla_{i}f(x)\geq 0,\end{cases} \]  
in which \(\
abla_{i},i=1,\ldots,n\) denotes the \(i-\)th derivative of the function \(f\). The first-order necessary conditions for the optimality of problem (2), in the case of minima within the box or on its boundaries, are expressed as [11]  
\[D(x)\cdot\
abla f(x)=0. \tag{5}\]  
To solve problem (2), HEMSim implements a Gauss-Newton method enhanced with a trust-region strategy, combining the rapid convergence of [PERSON]'s method with the robustness of the steepest descent approach. The algorithm generates a sequence of feasible iterates \(\{x_{k}\}\) defined by  
\[x_{k+1}=x_{k}+s_{k},\quad\forall k\in\mathbb{N}. \tag{6}\]  
At each iteration \(k\), the trust-region subproblem centered at \(x_{k}\) with radius \(\Delta_{k}\) involves defining a quadratic approximate model \(m_{k}\) using a second-order [PERSON] expansion of \(f\) around \(x_{k}\) with step \(s\):  
\[m_{k}(s)=\frac{1}{2}\|I_{k}\cdot s+R_{k}\|_{2}^{2},\qquad\forall s\in\mathbb{ R}^{n}, \tag{7}\]  
where \(J_{k}=J(x_{k})\) and \(R_{k}=R(x_{k})\). The trust-region strategy requires computing \(s\) by solving  
\[\min_{s}\{m_{k}(s):\quad\|s\|_{2}\leq\Delta_{k}\}. \tag{8}\]  
To determine a feasible step \(s_{k}\), a trial solution \(s_{\text{trial}}\) to the trust-region subproblem is found. By applying the stationary point condition to (7), the overdetermined system  
\[J_{k}\cdot s_{k}=-R_{k}, \tag{9}\]  
is derived, and the minimum norm solution \(s_{k}^{N}\) is given by  
\[s_{k}^{N}=-J_{k}^{+}\cdot R_{k}, \tag{10}\]  
where \(J_{k}^{+}\) is the Moore-Penrose pseudoinverse of \(J_{k}\), computed via singular value decomposition (SVD) [14]. This ensures the possibility of solving the problem even if the rank of \(J_{k}\) is not maximum. If \(\|s_{k}^{N}\|<\Delta_{k}\), then  
\[s_{\text{trial}}=s_{k}^{N} \tag{11}\]  
is used. Otherwise, if the absolute minimum of \(m_{k}(s)\) lies outside the trust region, the solution to subproblem (8) is chosen as the projection of the solution of the unbounded problem onto the trust region boundary:  
\[\min_{s}\{m_{k}(s):\quad\|s\|_{2}=\Delta_{k}\}. \tag{12}\]  
This problem can be solved exactly using Lagrange multipliers [17] or approximately by considering solutions in the subspace generated by the gradient of \(f\), \(\text{span}\{-D_{k}\
abla f\}\), to find the Cauchy point, which is the minimum of \(m_{k}(s)\) in the gradient direction within the trust region.