Fix convergence check in Newton method (#1663)

Fixes a bug where the Newton solver and potentially simulation would run longer than necessary (introduced in #446 / edf7d5b). Also saves some memory (2 * nx_solver^2).

> It seems that the value of `x_newton_` in `SteadystateProblem::applyNewtonsMethod`([applyNewtonsMethod](https://github.com/AMICI-dev/AMICI/blob/82d2cbb9d00ac9c3746b4e85c0fb1d6ebbff2863/src/steadystateproblem.cpp#L490)) doesn't change, but it's used to compute weighted root mean square of `xdot` and, therefore, is used to check convergence. As a result, convergence of Newton's method is achieved only due to `xdot_` becoming sufficiently small. 
For two models that require postequilibration I've observed simulations that finished with `posteq_status = [-4 1 0]`and `posteq_numstep = [0 0 0]`, i.e. Newton's method didn't converge, but the check here https://github.com/AMICI-dev/AMICI/blob/82d2cbb9d00ac9c3746b4e85c0fb1d6ebbff2863/src/steadystateproblem.cpp#L614 returned True and the simulation was considered successful with zero steps. 
Changing `x_newton_` to `x_` fixes the issue. 


* use x_ instead of x_newton_ for residuals computation
* fix typos
* Remove unused `SteadystateProblem::x_newton_` and `SteadystateProblem::rel_x_newton_`

include/amici/solver.h

@@ -632,7 +632,7 @@ class Solver {
                        AmiVectorArray &sx, AmiVector &xQ) const;
 
     /**
-     * @brief write solution from forward simulation
+     * @brief write solution from backward simulation
      * @param t time
      * @param xB adjoint state
      * @param dxB adjoint derivative state

include/amici/steadystateproblem.h

@@ -353,10 +353,6 @@ class SteadystateProblem {
     AmiVector ewt_;
     /** error weights for backward quadratures, dimension nplist() */
     AmiVector ewtQB_;
-    /** container for relative error calculation? */
-    AmiVector rel_x_newton_;
-    /** container for absolute error calculation? */
-    AmiVector x_newton_;
     /** state vector */
     AmiVector x_;
     /** old state vector */

src/backwardproblem.cpp

@@ -31,7 +31,7 @@ BackwardProblem::BackwardProblem(const ForwardProblem &fwd,
     root_idx_(fwd.getRootIndexes()),
     dJydx_(fwd.getDJydx()),
     dJzdx_(fwd.getDJzdx()) {
-        /* complement dJydx from postequilibration. This should overwrite
+        /* complement dJydx from postequilibration. This shouldn't overwrite
          * anything but only fill in previously 0 values, as only non-inf
          * timepoints are filled from fwd.
          */

src/steadystateproblem.cpp

@@ -20,7 +20,6 @@ namespace amici {
 
 SteadystateProblem::SteadystateProblem(const Solver &solver, const Model &model)
     : delta_(model.nx_solver), ewt_(model.nx_solver), ewtQB_(model.nplist()),
-      rel_x_newton_(model.nx_solver), x_newton_(model.nx_solver),
       x_(model.nx_solver), x_old_(model.nx_solver), dx_(model.nx_solver),
       xdot_(model.nx_solver), xdot_old_(model.nx_solver),
       sx_(model.nx_solver, model.nplist()), sdx_(model.nx_solver, model.nplist()),
@@ -505,11 +504,10 @@ void SteadystateProblem::applyNewtonsMethod(Model *model,
 
     /* Check for relative error, but make sure not to divide by 0!
         Ensure positivity of the state */
-    x_newton_ = x_;
     x_old_ = x_;
     xdot_old_ = xdot_;
 
-    wrms_ = getWrmsNorm(x_newton_, xdot_, newtonSolver->atol_,
+    wrms_ = getWrmsNorm(x_, xdot_, newtonSolver->atol_,
                         newtonSolver->rtol_, ewt_);
     bool converged = newton_retry ? false : wrms_ < RCONST(1.0);
     while (!converged && i_newtonstep < newtonSolver->max_steps) {
@@ -530,7 +528,7 @@ void SteadystateProblem::applyNewtonsMethod(Model *model,
 
         /* Compute new xdot and residuals */
         model->fxdot(t_, x_, dx_, xdot_);
-        realtype wrms_tmp = getWrmsNorm(x_newton_, xdot_, newtonSolver->atol_,
+        realtype wrms_tmp = getWrmsNorm(x_, xdot_, newtonSolver->atol_,
                                         newtonSolver->rtol_, ewt_);
 
         if (wrms_tmp < wrms_) {
@@ -555,7 +553,7 @@ void SteadystateProblem::applyNewtonsMethod(Model *model,
                 }
                 if (recheck_convergence) {
                   model->fxdot(t_, x_, dx_, xdot_);
-                  wrms_ = getWrmsNorm(x_newton_, xdot_, newtonSolver->atol_,
+                  wrms_ = getWrmsNorm(x_, xdot_, newtonSolver->atol_,
                                      newtonSolver->rtol_, ewt_);
                   converged = wrms_ < RCONST(1.0);
                 }
@@ -611,6 +609,7 @@ void SteadystateProblem::runSteadystateSimulation(const Solver *solver,
     if (backward)
         sensitivityFlag = SensitivityMethod::adjoint;
 
+    /* If run after Newton's method checks again if it converged */
     bool converged = checkConvergence(solver, model, sensitivityFlag);
     int sim_steps = 0;