Support iterative solve of the same solver_model (when using highs) about linopy HOT 4 OPEN

class Solver:
    """A persistent instance of a solver model which can be updated and resolved
    
    Concrete versions of this class replace f"run_{solver_name}"
    """
    def __init__(self, model, sanitize_zeros=True, problem_fn=None, log_fn=None, **solver_options):
        self.sanitize_zeros = sanitize_zeros
        self.problem_fn = problem_fn
        self.log_fn = log_fn
        self.solver_model = self.create_solver_model(model)
        self.solver_options = solver_options
        
    @staticmethod
    def get_solver(solver_name, **kwargs):
        if solver_name == "highs":
            return HighsSolver(**kwargs)
        if solver_name == "gurobi":
            return GurobiSolver(**kwargs)
        ...
    
    def run(self) -> SolverStatus:
        """Solver the solver_model"""
        raise NotImplementedError()
    
    def compute_IIS(self):
        raise NotImplementedError()
    
    def update_cost_coeffs(self, c) -> None:
        """Set the cost vector on the solver_model"""
        raise NotImplementedError()
    
    def update_constraint(self, constraint_name: str, values: xr.DataArray) -> None:
        """Set the cost vector on the solver_model"""
        pass
    
    def update_rhs(self, b) -> None:
        """Set the b vector"""
        pass
    
    def get_result(self) -> Result: 
        pass

    def close(self) -> None:
        """Free any resources (eg release the license if using gurobi)""""
        pass

@dataclass
class Result:
    ...
       
    def update_model(self, model):
        model.objective._value = self.solution.objective
        model.status = self.status.status.value
        model.termination_condition = self.status.termination_condition.value
        model.solver_model = self.solver_model

        if not self.status.is_ok:
            return
        self.assign_solution(self, model)
        self.assign_dual(self, model)

    def assign_solution(self, model):
        # map solution and dual to original shape which includes missing values
        sol = self.solution.primal.copy()
        sol.loc[-1] = nan

        for name, var in model.variables.items():
            idx = np.ravel(var.labels)
            try:
                vals = sol[idx].values.reshape(var.labels.shape)
            except KeyError:
                vals = sol.reindex(idx).values.reshape(var.labels.shape)
            var.solution = xr.DataArray(vals, var.coords)

    def assign_dual(self, model):
        # map solution and dual to original shape which includes missing values
        if not self.solution.dual.empty:
            dual = self.solution.dual.copy()
            dual.loc[-1] = nan

            for name, con in model.constraints.items():
                idx = np.ravel(con.labels)
                try:
                    vals = dual[idx].values.reshape(con.labels.shape)
                except KeyError:
                    vals = dual.reindex(idx).values.reshape(con.labels.shape)
                con.dual = xr.DataArray(vals, con.labels.coords)

def solve(
        self,
        ...
    ):
        if remote:
            ...
            return self.status, self.termination_condition

        solver = Solver.get_solver(
            solver_name=solver_name, 
            io_api=io_api, 
            env=env, 
            sanitize_zeros=sanitize_zeros, 
            problem_fn=problem_fn,
            solution_fn=solution_fn,
            log_fn=log_fn,
            basis_fn=basis_fn,
            warmstart_fn=warmstart_fn,
            **solver_options
        )
        with solver:
          result = solver.run()
          result.info()
          
          result.update_model(self)
          
        return result