refactor(exp/inverse): Simplify activation parameter to a single scalar stretch

The parameter `q` for the inverse problem, representing muscle activation, is changed from a 6-DoF Voigt-notation vector to a single scalar per active cell. This scalar represents the principal stretch factor, gamma.

The full activation matrix is now derived from this scalar, assuming volume preservation, with diagonal components `[1/gamma, sqrt(gamma), sqrt(gamma)]`. This significantly reduces the dimensionality of the optimization problem from `6 * n_active_cells` to `n_active_cells`, making it better-posed and easier to solve.

Additionally, the maximum number of iterations for the conjugate gradient linear solver in the adjoint method is increased to improve the accuracy of the gradient calculation.

Author: liblaf

exp/2025/09/24/inverse-grin/src/33-inverse-no-reg.py

@@ -182,7 +182,9 @@ def muscle_orientation(self) -> Float[Array, "c 3 3"]:
     def make_params(self, q: Float[Array, "ca 6"]) -> Params:
         activation: Float[Array, "c 6"] = sim_jax.rest_activation(self.input.n_cells)
         # q = q.at[:, :3].set(jnp.exp(q[:, :3]))
-        activation = activation.at[self.active_mask].set(q)
+        activation = activation.at[self.active_mask, 0].set(jnp.reciprocal(q))
+        activation = activation.at[self.active_mask, 1].set(jnp.sqrt(q))
+        activation = activation.at[self.active_mask, 2].set(jnp.sqrt(q))
         # activation = sim_jax.transform_activation(activation, self.muscle_orientation)
         return Params(activation=activation)
 
@@ -240,7 +242,7 @@ def fun_and_jac(self, q: Array) -> tuple[Scalar, Array]:
                 -dLdu_free,
                 tol=1e-5,
                 atol=1e-15,
-                maxiter=ic(u_free.size),
+                maxiter=ic(10 * u_free.size),
                 M=lambda x: P_free * x,
             )
             logger.info("linear solve > info: {}", info)
@@ -395,12 +397,10 @@ def callback(intermediate_result: optim.Solution) -> None:
         # result: pv.UnstructuredGrid = mesh.warp_by_vector("solution")  # pyright: ignore[reportAssignmentType]
         writer.append(mesh)
 
+    q_init: Float[Array, " ca"] = jnp.ones((inverse.n_active_cells,))
     inverse.reg_mean_weight = 0.0
     inverse.reg_shear_weight = 0.0
     inverse.reg_volume_weight = 0.0
-
-    q_init: Float[Array, "ca 6"] = sim_jax.rest_activation(inverse.n_active_cells)
-    # q_init = q_init.at[:, :3].set(jnp.log(1.0))
     callback(optim.Solution({"x": q_init}))
     optimizer = optim.MinimizerScipy(
         jit=False, method="L-BFGS-B", tol=1e-15, options={}