diff --git a/FEATURES.md b/FEATURES.md
index eb17f2ed97..f6ff74c98c 100644
--- a/FEATURES.md
+++ b/FEATURES.md
@@ -17,3 +17,5 @@ Fixed bug in presolve when pointers stored in HighsMatrixSlice get invalidated w
 Primal and dual residual tolerances - applied following IPM or PDLP solution - now documented as options
 
 Highs::getCols (Highs::getRows) now runs in linear time if the internal constraint matrix is stored column-wise (row-wise). Added ensureColwise/Rowwise to the Highs class, the C API and highspy so that users can set the internal constraint matrix storage orientation
+
+When columns and rows are deleted from the incumbent LP after a basic solution has been found, HiGHS no longer invalidates the basis. Now it maintains the basic and nonbasic status of the remaining variables and constraints. When the model is re-solved, this information is used to construct a starting basis.
\ No newline at end of file
diff --git a/app/HighsRuntimeOptions.h b/app/HighsRuntimeOptions.h
index 04d3133e4b..2892a2fa86 100644
--- a/app/HighsRuntimeOptions.h
+++ b/app/HighsRuntimeOptions.h
@@ -145,7 +145,7 @@ bool loadOptions(const HighsLogOptions& report_log_options, int argc,
         case HighsLoadOptionsStatus::kError:
           return false;
         case HighsLoadOptionsStatus::kEmpty:
-          writeOptionsToFile(stdout, options.records);
+          writeOptionsToFile(stdout, options.log_options, options.records);
           return false;
         default:
           break;
diff --git a/app/RunHighs.cpp b/app/RunHighs.cpp
index 8813755b4c..b5830bc370 100644
--- a/app/RunHighs.cpp
+++ b/app/RunHighs.cpp
@@ -47,6 +47,7 @@ int main(int argc, char** argv) {
   // settings defined in any options file.
   highs.passOptions(loaded_options);
   //  highs.writeOptions("Options.md");
+  highs.writeOptions("", true);
 
   // Load the model from model_file
   HighsStatus read_status = highs.readModel(model_file);
diff --git a/check/TestLpModification.cpp b/check/TestLpModification.cpp
index 30bdac0f06..9a443ea43f 100644
--- a/check/TestLpModification.cpp
+++ b/check/TestLpModification.cpp
@@ -2274,3 +2274,204 @@ TEST_CASE("row-wise-get-row-avgas", "[highs_data]") {
   testAvgasGetRow(h);
   testAvgasGetCol(h);
 }
+
+TEST_CASE("hot-start-after-delete", "[highs_data]") {
+  Highs h;
+  h.setOptionValue("output_flag", dev_run);
+  const HighsLp& lp = h.getLp();
+  const HighsInfo& info = h.getInfo();
+  const HighsBasis& basis = h.getBasis();
+  const HighsSolution& solution = h.getSolution();
+  std::string model = "avgas";
+  std::string model_file =
+      std::string(HIGHS_DIR) + "/check/instances/" + model + ".mps";
+  h.readModel(model_file);
+  h.run();
+  HighsInt ieration_count0 = info.simplex_iteration_count;
+  if (dev_run)
+    printf("Initial solve takes %d iterations and yields objective = %g\n",
+           int(info.simplex_iteration_count), info.objective_function_value);
+
+  HighsInt max_dim = std::max(lp.num_col_, lp.num_row_);
+  std::vector<double> cost(1);
+  std::vector<double> lower(1);
+  std::vector<double> upper(1);
+  HighsInt nnz;
+  std::vector<HighsInt> start(1);
+  std::vector<HighsInt> index(max_dim);
+  std::vector<double> value(max_dim);
+  HighsInt get_num;
+  HighsInt use_col, use_row;
+  for (HighsInt k = 0; k < 2; k++) {
+    if (dev_run) {
+      for (HighsInt iCol = 0; iCol < lp.num_col_; iCol++)
+        printf("Col %2d is %s\n", int(iCol),
+               basis.col_status[iCol] == HighsBasisStatus::kBasic ? "basic"
+                                                                  : "nonbasic");
+      printf("\n");
+    }
+    if (k == 0) {
+      use_col = 1;  // Nonbasic
+    } else {
+      use_col = 4;  // Basic
+    }
+    if (dev_run)
+      printf(
+          "\nDeleting and adding column %1d with status \"%s\" and value %g\n",
+          int(use_col),
+          h.basisStatusToString(basis.col_status[use_col]).c_str(),
+          solution.col_value[use_col]);
+
+    h.getCols(use_col, use_col, get_num, cost.data(), lower.data(),
+              upper.data(), nnz, start.data(), index.data(), value.data());
+
+    h.deleteCols(use_col, use_col);
+    if (dev_run) basis.printScalars();
+
+    h.addCol(cost[0], lower[0], upper[0], nnz, index.data(), value.data());
+
+    h.run();
+    if (dev_run)
+      printf(
+          "After deleting and adding column %1d, solve takes %d iterations and "
+          "yields objective = %g\n",
+          int(use_col), int(info.simplex_iteration_count),
+          info.objective_function_value);
+    REQUIRE(info.simplex_iteration_count < ieration_count0);
+  }
+
+  for (HighsInt k = 0; k < 2; k++) {
+    if (dev_run) {
+      for (HighsInt iRow = 0; iRow < lp.num_row_; iRow++)
+        printf("Row %2d is %s\n", int(iRow),
+               basis.row_status[iRow] == HighsBasisStatus::kBasic ? "basic"
+                                                                  : "nonbasic");
+    }
+    if (k == 0) {
+      use_row = 1;  // Nonbasic
+    } else {
+      use_row = 8;  // Basic
+    }
+    if (dev_run)
+      printf("\nDeleting and adding row %1d with status \"%s\" and value %g\n",
+             int(use_row),
+             h.basisStatusToString(basis.row_status[use_row]).c_str(),
+             solution.row_value[use_row]);
+
+    h.getRows(use_row, use_row, get_num, lower.data(), upper.data(), nnz,
+              start.data(), index.data(), value.data());
+
+    h.deleteRows(use_row, use_row);
+    if (dev_run) basis.printScalars();
+
+    h.addRow(lower[0], upper[0], nnz, index.data(), value.data());
+
+    h.run();
+    if (dev_run)
+      printf(
+          "After deleting and adding row %1d, solve takes %d iterations and "
+          "yields objective = %g\n",
+          int(use_row), int(info.simplex_iteration_count),
+          info.objective_function_value);
+    REQUIRE(info.simplex_iteration_count < ieration_count0);
+  }
+  std::vector<HighsInt> set = {1, 3, 4};
+  HighsInt num_set_en = set.size();
+  cost.resize(num_set_en);
+  lower.resize(num_set_en);
+  upper.resize(num_set_en);
+  start.resize(num_set_en);
+  index.resize(num_set_en * max_dim);
+  value.resize(num_set_en * max_dim);
+
+  h.getCols(num_set_en, set.data(), get_num, cost.data(), lower.data(),
+            upper.data(), nnz, start.data(), index.data(), value.data());
+
+  h.deleteCols(num_set_en, set.data());
+  if (dev_run) basis.printScalars();
+
+  h.addCols(get_num, cost.data(), lower.data(), upper.data(), nnz, start.data(),
+            index.data(), value.data());
+
+  h.run();
+  if (dev_run)
+    printf(
+        "After deleting and adding %d columns in set, solve takes %d "
+        "iterations and yields objective = %g\n",
+        int(get_num), int(info.simplex_iteration_count),
+        info.objective_function_value);
+  //  REQUIRE(info.simplex_iteration_count < ieration_count0);
+
+  h.getRows(num_set_en, set.data(), get_num, lower.data(), upper.data(), nnz,
+            start.data(), index.data(), value.data());
+
+  h.deleteRows(num_set_en, set.data());
+  if (dev_run) basis.printScalars();
+
+  h.addRows(get_num, lower.data(), upper.data(), nnz, start.data(),
+            index.data(), value.data());
+
+  h.run();
+  if (dev_run)
+    printf(
+        "After deleting and adding %d rows in set, solve takes %d iterations "
+        "and yields objective = %g\n",
+        int(get_num), int(info.simplex_iteration_count),
+        info.objective_function_value);
+  //  REQUIRE(info.simplex_iteration_count < ieration_count0);
+  std::vector<HighsInt> mask;
+  mask.assign(max_dim, 0);
+  mask[1] = 1;
+  mask[4] = 1;
+  mask[5] = 1;
+
+  h.getCols(mask.data(), get_num, cost.data(), lower.data(), upper.data(), nnz,
+            start.data(), index.data(), value.data());
+
+  h.deleteCols(mask.data());
+  if (dev_run) basis.printScalars();
+
+  h.addCols(get_num, cost.data(), lower.data(), upper.data(), nnz, start.data(),
+            index.data(), value.data());
+
+  h.run();
+  if (dev_run)
+    printf(
+        "After deleting and adding %d columns in mask, solve takes %d "
+        "iterations and yields objective = %g\n",
+        int(get_num), int(info.simplex_iteration_count),
+        info.objective_function_value);
+  //  REQUIRE(info.simplex_iteration_count < ieration_count0);
+
+  mask.assign(max_dim, 0);
+  mask[1] = 1;
+  mask[4] = 1;
+  mask[5] = 1;
+  mask[8] = 1;
+  mask[9] = 1;
+  HighsInt num_mask_en = mask.size();
+  cost.resize(num_mask_en);
+  lower.resize(num_mask_en);
+  upper.resize(num_mask_en);
+  start.resize(num_mask_en);
+  index.resize(num_mask_en * max_dim);
+  value.resize(num_mask_en * max_dim);
+
+  h.getRows(mask.data(), get_num, lower.data(), upper.data(), nnz, start.data(),
+            index.data(), value.data());
+
+  h.deleteRows(mask.data());
+  if (dev_run) basis.printScalars();
+
+  h.addRows(get_num, lower.data(), upper.data(), nnz, start.data(),
+            index.data(), value.data());
+
+  h.run();
+  if (dev_run)
+    printf(
+        "After deleting and adding %d rows in mask, solve takes %d iterations "
+        "and yields objective = %g\n",
+        int(get_num), int(info.simplex_iteration_count),
+        info.objective_function_value);
+  //  REQUIRE(info.simplex_iteration_count < ieration_count0);
+}
diff --git a/check/TestOptions.cpp b/check/TestOptions.cpp
index f511f8caef..b2a4142710 100644
--- a/check/TestOptions.cpp
+++ b/check/TestOptions.cpp
@@ -95,7 +95,7 @@ TEST_CASE("internal-options", "[highs_options]") {
   REQUIRE(options.small_matrix_value == 0.001);
   REQUIRE(options.mps_parser_type_free);
 
-  if (dev_run) reportOptions(stdout, options.records, true);
+  if (dev_run) reportOptions(stdout, report_log_options, options.records, true);
 
   return_status = checkOptions(report_log_options, options.records);
   REQUIRE(return_status == OptionStatus::kOk);
@@ -157,7 +157,7 @@ TEST_CASE("internal-options", "[highs_options]") {
 
   if (dev_run) {
     printf("\nAfter setting allowed_matrix_scale_factor to 1\n");
-    reportOptions(stdout, options.records);
+    reportOptions(stdout, report_log_options, options.records);
   }
 
   double allowed_matrix_scale_factor_double = 1e-7;
@@ -174,7 +174,7 @@ TEST_CASE("internal-options", "[highs_options]") {
 
   if (dev_run) {
     printf("\nAfter testing HighsInt options\n");
-    reportOptions(stdout, options.records);
+    reportOptions(stdout, report_log_options, options.records);
   }
 
   // Check setting double options
@@ -231,7 +231,7 @@ TEST_CASE("internal-options", "[highs_options]") {
                           options.log_options, options.records, model_file);
   REQUIRE(return_status == OptionStatus::kUnknownOption);
 
-  if (dev_run) reportOptions(stdout, options.records);
+  if (dev_run) reportOptions(stdout, report_log_options, options.records);
 
   bool get_mps_parser_type_free;
   return_status =
diff --git a/src/Highs.h b/src/Highs.h
index 8480b3db35..2986d292ff 100644
--- a/src/Highs.h
+++ b/src/Highs.h
@@ -312,11 +312,11 @@ class Highs {
 
   /**
    * @brief Write (deviations from default values of) the options to a
-   * file, with the extension ".html" producing HTML, otherwise using
-   * the standard format used to read options from a file.
+   * file, using the standard format used to read options from a file.
+   * Possible to write only deviations from default values.
    */
   HighsStatus writeOptions(const std::string& filename,  //!< The filename
-                           const bool report_only_deviations = false) const;
+                           const bool report_only_deviations = false);
 
   /**
    * @brief Returns the number of user-settable options
diff --git a/src/lp_data/HStruct.h b/src/lp_data/HStruct.h
index 77b9e6e11b..d69c36dede 100644
--- a/src/lp_data/HStruct.h
+++ b/src/lp_data/HStruct.h
@@ -58,14 +58,27 @@ struct HotStart {
 };
 
 struct HighsBasis {
+  // Logical flags for a HiGHS basis:
+  //
+  // valid: has been factored by HiGHS
+  //
+  // alien: a basis that's been set externally, so cannot be assumed
+  // to even have the right number of basic and nonbasic variables
+  //
+  // useful: a basis that may be useful
+  //
+  // Need useful since, by default, a basis is alien but not useful
   bool valid = false;
   bool alien = true;
+  bool useful = false;
   bool was_alien = true;
   HighsInt debug_id = -1;
   HighsInt debug_update_count = -1;
   std::string debug_origin_name = "None";
   std::vector<HighsBasisStatus> col_status;
   std::vector<HighsBasisStatus> row_status;
+  void print(std::string message = "") const;
+  void printScalars(std::string message = "") const;
   void invalidate();
   void clear();
 };
diff --git a/src/lp_data/Highs.cpp b/src/lp_data/Highs.cpp
index df38dad281..94ed199839 100644
--- a/src/lp_data/Highs.cpp
+++ b/src/lp_data/Highs.cpp
@@ -139,7 +139,8 @@ HighsStatus Highs::resetOptions() {
 }
 
 HighsStatus Highs::writeOptions(const std::string& filename,
-                                const bool report_only_deviations) const {
+                                const bool report_only_deviations) {
+  this->logHeader();
   HighsStatus return_status = HighsStatus::kOk;
   FILE* file;
   HighsFileType file_type;
@@ -148,15 +149,16 @@ HighsStatus Highs::writeOptions(const std::string& filename,
       openWriteFile(filename, "writeOptions", file, file_type), return_status,
       "openWriteFile");
   if (return_status == HighsStatus::kError) return return_status;
+  if (filename == "") file_type = HighsFileType::kMinimal;
   // Report to user that options are being written to a file
   if (filename != "")
     highsLogUser(options_.log_options, HighsLogType::kInfo,
                  "Writing the option values to %s\n", filename.c_str());
-  return_status =
-      interpretCallStatus(options_.log_options,
-                          writeOptionsToFile(file, options_.records,
-                                             report_only_deviations, file_type),
-                          return_status, "writeOptionsToFile");
+  return_status = interpretCallStatus(
+      options_.log_options,
+      writeOptionsToFile(file, options_.log_options, options_.records,
+                         report_only_deviations, file_type),
+      return_status, "writeOptionsToFile");
   if (file != stdout) fclose(file);
   return return_status;
 }
@@ -716,6 +718,7 @@ HighsStatus Highs::readBasis(const std::string& filename) {
   // Update the HiGHS basis and invalidate any simplex basis for the model
   basis_ = read_basis;
   basis_.valid = true;
+  basis_.useful = true;
   // Follow implications of a new HiGHS basis
   newHighsBasis();
   // Can't use returnFromHighs since...
@@ -1265,8 +1268,16 @@ HighsStatus Highs::solve() {
   // is available, simplex should surely be chosen.
   const bool solver_will_use_basis = options_.solver == kSimplexString ||
                                      options_.solver == kHighsChooseString;
-  if ((basis_.valid || options_.presolve == kHighsOffString ||
-       unconstrained_lp) &&
+  const bool has_basis = basis_.useful;
+  if (has_basis) {
+    assert(basis_.col_status.size() ==
+           static_cast<size_t>(incumbent_lp.num_col_));
+    assert(basis_.row_status.size() ==
+           static_cast<size_t>(incumbent_lp.num_row_));
+  }
+  if (basis_.valid) assert(basis_.useful);
+
+  if ((has_basis || options_.presolve == kHighsOffString || unconstrained_lp) &&
       solver_will_use_basis) {
     // There is a valid basis for the problem, presolve is off, or LP
     // has no constraint matrix, and the solver will use the basis
@@ -1274,7 +1285,7 @@ HighsStatus Highs::solve() {
         "LP without presolve, or with basis, or unconstrained";
     // If there is a valid HiGHS basis, refine any status values that
     // are simply HighsBasisStatus::kNonbasic
-    if (basis_.valid) refineBasis(incumbent_lp, solution_, basis_);
+    if (basis_.useful) refineBasis(incumbent_lp, solution_, basis_);
     solveLp(incumbent_lp,
             "Solving LP without presolve, or with basis, or unconstrained",
             this_solve_original_lp_time);
@@ -1416,6 +1427,7 @@ HighsStatus Highs::solve() {
         basis_.debug_origin_name = "Presolve to empty";
         basis_.valid = true;
         basis_.alien = false;
+        basis_.useful = true;
         basis_.was_alien = false;
         solution_.value_valid = true;
         solution_.dual_valid = true;
@@ -1551,6 +1563,7 @@ HighsStatus Highs::solve() {
           solution_.dual_valid = true;
           // Set basis and its status
           basis_.valid = true;
+          basis_.useful = true;
           basis_.col_status = presolve_.data_.recovered_basis_.col_status;
           basis_.row_status = presolve_.data_.recovered_basis_.row_status;
           basis_.debug_origin_name += ": after postsolve";
@@ -2339,6 +2352,7 @@ HighsStatus Highs::setBasis(const HighsBasis& basis,
     basis_ = basis;
   }
   basis_.valid = true;
+  basis_.useful = true;
   if (origin != "") basis_.debug_origin_name = origin;
   assert(basis_.debug_origin_name != "");
   assert(!basis_.alien);
@@ -4157,6 +4171,7 @@ HighsStatus Highs::callRunPostsolve(const HighsSolution& solution,
       // Set basis and its status
       //
       // basis_.valid = true;
+      // basis_.useful = true;
       // basis_.col_status = presolve_.data_.recovered_basis_.col_status;
       // basis_.row_status = presolve_.data_.recovered_basis_.row_status;
       basis_.debug_origin_name += ": after postsolve";
@@ -4270,10 +4285,12 @@ void Highs::forceHighsSolutionBasisSize() {
   if (basis_.col_status.size() != static_cast<size_t>(model_.lp_.num_col_)) {
     basis_.col_status.resize(model_.lp_.num_col_);
     basis_.valid = false;
+    basis_.useful = false;
   }
   if (basis_.row_status.size() != static_cast<size_t>(model_.lp_.num_row_)) {
     basis_.row_status.resize(model_.lp_.num_row_);
     basis_.valid = false;
+    basis_.useful = false;
   }
 }
 
diff --git a/src/lp_data/HighsInterface.cpp b/src/lp_data/HighsInterface.cpp
index e522e46db9..627338c7e1 100644
--- a/src/lp_data/HighsInterface.cpp
+++ b/src/lp_data/HighsInterface.cpp
@@ -305,7 +305,7 @@ HighsStatus Highs::addColsInterface(
   HighsLp& lp = model_.lp_;
   HighsBasis& basis = basis_;
   HighsScale& scale = lp.scale_;
-  bool& valid_basis = basis.valid;
+  bool& useful_basis = basis.useful;
   bool& lp_has_scaling = lp.scale_.has_scaling;
 
   // Check that if nonzeros are to be added then the model has a positive number
@@ -415,7 +415,7 @@ HighsStatus Highs::addColsInterface(
                                       &scale.col[lp.num_col_]);
   }
   // Update the basis corresponding to new nonbasic columns
-  if (valid_basis) appendNonbasicColsToBasisInterface(ext_num_new_col);
+  if (useful_basis) appendNonbasicColsToBasisInterface(ext_num_new_col);
 
   // Possibly add column names
   lp.addColNames("", ext_num_new_col);
@@ -471,7 +471,8 @@ HighsStatus Highs::addRowsInterface(HighsInt ext_num_new_row,
   HighsLp& lp = model_.lp_;
   HighsBasis& basis = basis_;
   HighsScale& scale = lp.scale_;
-  bool& valid_basis = basis.valid;
+  bool& useful_basis = basis.useful;
+
   bool& lp_has_scaling = lp.scale_.has_scaling;
 
   // Check that if nonzeros are to be added then the model has a positive number
@@ -559,7 +560,7 @@ HighsStatus Highs::addRowsInterface(HighsInt ext_num_new_row,
                                        &scale.row[lp.num_row_]);
   }
   // Update the basis corresponding to new basic rows
-  if (valid_basis) appendBasicRowsToBasisInterface(ext_num_new_row);
+  if (useful_basis) appendBasicRowsToBasisInterface(ext_num_new_row);
 
   // Possibly add row names
   lp.addRowNames("", ext_num_new_row);
@@ -576,6 +577,69 @@ HighsStatus Highs::addRowsInterface(HighsInt ext_num_new_row,
   return return_status;
 }
 
+void deleteBasisEntries(std::vector<HighsBasisStatus>& status,
+                        bool& deleted_basic, bool& deleted_nonbasic,
+                        const HighsIndexCollection& index_collection,
+                        const HighsInt entry_dim) {
+  assert(ok(index_collection));
+  assert(static_cast<size_t>(entry_dim) == status.size());
+  HighsInt from_k;
+  HighsInt to_k;
+  limits(index_collection, from_k, to_k);
+  if (from_k > to_k) return;
+
+  HighsInt delete_from_entry;
+  HighsInt delete_to_entry;
+  HighsInt keep_from_entry;
+  HighsInt keep_to_entry = -1;
+  HighsInt current_set_entry = 0;
+  HighsInt new_num_entry = 0;
+  deleted_basic = false;
+  deleted_nonbasic = false;
+  for (HighsInt k = from_k; k <= to_k; k++) {
+    updateOutInIndex(index_collection, delete_from_entry, delete_to_entry,
+                     keep_from_entry, keep_to_entry, current_set_entry);
+    // Account for the initial entries being kept
+    if (k == from_k) new_num_entry = delete_from_entry;
+    // Identify whether a basic or a nonbasic entry has been deleted
+    for (HighsInt entry = delete_from_entry; entry <= delete_to_entry;
+         entry++) {
+      if (status[entry] == HighsBasisStatus::kBasic) {
+        deleted_basic = true;
+      } else {
+        deleted_nonbasic = true;
+      }
+    }
+    if (delete_to_entry >= entry_dim - 1) break;
+    for (HighsInt entry = keep_from_entry; entry <= keep_to_entry; entry++) {
+      status[new_num_entry] = status[entry];
+      new_num_entry++;
+    }
+    if (keep_to_entry >= entry_dim - 1) break;
+  }
+  status.resize(new_num_entry);
+}
+
+void deleteBasisCols(HighsBasis& basis,
+                     const HighsIndexCollection& index_collection,
+                     const HighsInt original_num_col) {
+  bool deleted_basic;
+  bool deleted_nonbasic;
+  deleteBasisEntries(basis.col_status, deleted_basic, deleted_nonbasic,
+                     index_collection, original_num_col);
+  if (deleted_basic) basis.valid = false;
+}
+
+void deleteBasisRows(HighsBasis& basis,
+                     const HighsIndexCollection& index_collection,
+                     const HighsInt original_num_row) {
+  bool deleted_basic;
+  bool deleted_nonbasic;
+  deleteBasisEntries(basis.row_status, deleted_basic, deleted_nonbasic,
+                     index_collection, original_num_row);
+  if (deleted_nonbasic) basis.valid = false;
+}
+
 void Highs::deleteColsInterface(HighsIndexCollection& index_collection) {
   HighsLp& lp = model_.lp_;
   HighsBasis& basis = basis_;
@@ -587,13 +651,24 @@ void Highs::deleteColsInterface(HighsIndexCollection& index_collection) {
 
   lp.deleteCols(index_collection);
   model_.hessian_.deleteCols(index_collection);
-  assert(lp.num_col_ <= original_num_col);
-  if (lp.num_col_ < original_num_col) {
-    // Nontrivial deletion so reset the model_status and invalidate
-    // the Highs basis
-    model_status_ = HighsModelStatus::kNotset;
-    basis.valid = false;
+  // Bail out if no columns were actually deleted
+  if (lp.num_col_ == original_num_col) return;
+
+  assert(lp.num_col_ < original_num_col);
+
+  // Nontrivial deletion so reset the model_status and update any
+  // Highs basis
+  model_status_ = HighsModelStatus::kNotset;
+  if (basis_.useful) {
+    assert(basis_.col_status.size() == static_cast<size_t>(original_num_col));
+    // Have a full set of column basis status values, so maintain
+    // them, and only invalidate the basis if a basic column has been
+    // deleted
+    deleteBasisCols(basis_, index_collection, original_num_col);
+  } else {
+    assert(!basis.valid);
   }
+
   if (lp.scale_.has_scaling) {
     deleteScale(lp.scale_.col, index_collection);
     lp.scale_.col.resize(lp.num_col_);
@@ -632,13 +707,24 @@ void Highs::deleteRowsInterface(HighsIndexCollection& index_collection) {
   HighsInt original_num_row = lp.num_row_;
 
   lp.deleteRows(index_collection);
-  assert(lp.num_row_ <= original_num_row);
-  if (lp.num_row_ < original_num_row) {
-    // Nontrivial deletion so reset the model_status and invalidate
-    // the Highs basis
-    model_status_ = HighsModelStatus::kNotset;
-    basis.valid = false;
+  // Bail out if no rows were actually deleted
+  if (lp.num_row_ == original_num_row) return;
+
+  assert(lp.num_row_ < original_num_row);
+
+  // Nontrivial deletion so reset the model_status and update any
+  // Highs basis
+  model_status_ = HighsModelStatus::kNotset;
+  if (basis_.useful) {
+    assert(basis_.row_status.size() == static_cast<size_t>(original_num_row));
+    // Have a full set of row basis status values, so maintain them,
+    // and only invalidate the basis if a nonbasic row has been
+    // deleted
+    deleteBasisRows(basis_, index_collection, original_num_row);
+  } else {
+    assert(!basis.valid);
   }
+
   if (lp.scale_.has_scaling) {
     deleteScale(lp.scale_.row, index_collection);
     lp.scale_.row.resize(lp.num_row_);
@@ -1169,14 +1255,17 @@ void Highs::setNonbasicStatusInterface(
 }
 
 void Highs::appendNonbasicColsToBasisInterface(const HighsInt ext_num_new_col) {
+  if (ext_num_new_col == 0) return;
   HighsBasis& highs_basis = basis_;
-  if (!highs_basis.valid) return;
+  if (!highs_basis.useful) return;
   const bool has_simplex_basis = ekk_instance_.status_.has_basis;
   SimplexBasis& simplex_basis = ekk_instance_.basis_;
   HighsLp& lp = model_.lp_;
 
+  assert(highs_basis.col_status.size() == static_cast<size_t>(lp.num_col_));
+  assert(highs_basis.row_status.size() == static_cast<size_t>(lp.num_row_));
+
   // Add nonbasic structurals
-  if (ext_num_new_col == 0) return;
   HighsInt newNumCol = lp.num_col_ + ext_num_new_col;
   HighsInt newNumTot = newNumCol + lp.num_row_;
   highs_basis.col_status.resize(newNumCol);
@@ -1243,13 +1332,17 @@ void Highs::appendNonbasicColsToBasisInterface(const HighsInt ext_num_new_col) {
 }
 
 void Highs::appendBasicRowsToBasisInterface(const HighsInt ext_num_new_row) {
+  if (ext_num_new_row == 0) return;
   HighsBasis& highs_basis = basis_;
-  if (!highs_basis.valid) return;
+  if (!highs_basis.useful) return;
   const bool has_simplex_basis = ekk_instance_.status_.has_basis;
   SimplexBasis& simplex_basis = ekk_instance_.basis_;
   HighsLp& lp = model_.lp_;
+
+  assert(highs_basis.col_status.size() == static_cast<size_t>(lp.num_col_));
+  assert(highs_basis.row_status.size() == static_cast<size_t>(lp.num_row_));
+
   // Add basic logicals
-  if (ext_num_new_row == 0) return;
   // Add the new rows to the Highs basis
   HighsInt newNumRow = lp.num_row_ + ext_num_new_row;
   highs_basis.row_status.resize(newNumRow);
@@ -1538,6 +1631,7 @@ HighsStatus Highs::setHotStartInterface(const HotStart& hot_start) {
     nonbasicMove[num_col + iRow] = move;
   }
   basis_.valid = true;
+  basis_.useful = true;
   ekk_instance_.status_.has_basis = true;
   ekk_instance_.setNlaRefactorInfo();
   ekk_instance_.updateStatus(LpAction::kHotStart);
@@ -1944,6 +2038,12 @@ HighsStatus Highs::elasticityFilterReturn(
 
   run_status = this->deleteCols(original_num_col, lp.num_col_ - 1);
   assert(run_status == HighsStatus::kOk);
+  //
+  // Now that deleteRows and deleteCols may yield a valid basis, the
+  // lack of dual values triggers an assert in
+  // getKktFailures. Ultimately (#2081) the dual values will be
+  // available but, for now, make the basis invalid.
+  basis_.valid = false;
 
   run_status =
       this->changeColsCost(0, original_num_col - 1, original_col_cost.data());
diff --git a/src/lp_data/HighsOptions.cpp b/src/lp_data/HighsOptions.cpp
index 29ddd37629..15c3fb6204 100644
--- a/src/lp_data/HighsOptions.cpp
+++ b/src/lp_data/HighsOptions.cpp
@@ -811,15 +811,17 @@ void resetLocalOptions(std::vector<OptionRecord*>& option_records) {
   }
 }
 
-HighsStatus writeOptionsToFile(FILE* file,
+HighsStatus writeOptionsToFile(FILE* file, const HighsLogOptions& log_options,
                                const std::vector<OptionRecord*>& option_records,
                                const bool report_only_deviations,
                                const HighsFileType file_type) {
-  reportOptions(file, option_records, report_only_deviations, file_type);
+  reportOptions(file, log_options, option_records, report_only_deviations,
+                file_type);
   return HighsStatus::kOk;
 }
 
-void reportOptions(FILE* file, const std::vector<OptionRecord*>& option_records,
+void reportOptions(FILE* file, const HighsLogOptions& log_options,
+                   const std::vector<OptionRecord*>& option_records,
                    const bool report_only_deviations,
                    const HighsFileType file_type) {
   HighsInt num_options = option_records.size();
@@ -831,22 +833,27 @@ void reportOptions(FILE* file, const std::vector<OptionRecord*>& option_records,
       if (!kAdvancedInDocumentation) continue;
     }
     if (type == HighsOptionType::kBool) {
-      reportOption(file, ((OptionRecordBool*)option_records[index])[0],
+      reportOption(file, log_options,
+                   ((OptionRecordBool*)option_records[index])[0],
                    report_only_deviations, file_type);
     } else if (type == HighsOptionType::kInt) {
-      reportOption(file, ((OptionRecordInt*)option_records[index])[0],
+      reportOption(file, log_options,
+                   ((OptionRecordInt*)option_records[index])[0],
                    report_only_deviations, file_type);
     } else if (type == HighsOptionType::kDouble) {
-      reportOption(file, ((OptionRecordDouble*)option_records[index])[0],
+      reportOption(file, log_options,
+                   ((OptionRecordDouble*)option_records[index])[0],
                    report_only_deviations, file_type);
     } else {
-      reportOption(file, ((OptionRecordString*)option_records[index])[0],
+      reportOption(file, log_options,
+                   ((OptionRecordString*)option_records[index])[0],
                    report_only_deviations, file_type);
     }
   }
 }
 
-void reportOption(FILE* file, const OptionRecordBool& option,
+void reportOption(FILE* file, const HighsLogOptions& log_options,
+                  const OptionRecordBool& option,
                   const bool report_only_deviations,
                   const HighsFileType file_type) {
   if (!report_only_deviations || option.default_value != *option.value) {
@@ -865,40 +872,53 @@ void reportOption(FILE* file, const OptionRecordBool& option,
       fprintf(file, "%s = %s\n", option.name.c_str(),
               highsBoolToString(*option.value).c_str());
     } else {
-      fprintf(file, "%s = %s\n", option.name.c_str(),
-              highsBoolToString(*option.value).c_str());
+      std::string line =
+          highsFormatToString("Set option %s to %s\n", option.name.c_str(),
+                              highsBoolToString(*option.value).c_str());
+      if (file == stdout) {
+        highsLogUser(log_options, HighsLogType::kInfo, "%s", line.c_str());
+      } else {
+        fprintf(file, "%s", line.c_str());
+      }
     }
   }
 }
 
-void reportOption(FILE* file, const OptionRecordInt& option,
+void reportOption(FILE* file, const HighsLogOptions& log_options,
+                  const OptionRecordInt& option,
                   const bool report_only_deviations,
                   const HighsFileType file_type) {
   if (!report_only_deviations || option.default_value != *option.value) {
     if (file_type == HighsFileType::kMd) {
-      fprintf(file,
-              "## %s\n- %s\n- Type: integer\n- Range: {%" HIGHSINT_FORMAT
-              ", %" HIGHSINT_FORMAT "}\n- Default: %" HIGHSINT_FORMAT "\n\n",
-              highsInsertMdEscapes(option.name).c_str(),
-              highsInsertMdEscapes(option.description).c_str(),
-              option.lower_bound, option.upper_bound, option.default_value);
+      fprintf(
+          file,
+          "## %s\n- %s\n- Type: integer\n- Range: {%d, %d}\n- Default: %d\n\n",
+          highsInsertMdEscapes(option.name).c_str(),
+          highsInsertMdEscapes(option.description).c_str(),
+          int(option.lower_bound), int(option.upper_bound),
+          int(option.default_value));
     } else if (file_type == HighsFileType::kFull) {
       fprintf(file, "\n# %s\n", option.description.c_str());
       fprintf(file,
               "# [type: integer, advanced: %s, range: {%" HIGHSINT_FORMAT
-              ", %" HIGHSINT_FORMAT "}, default: %" HIGHSINT_FORMAT "]\n",
+              ", %d}, default: %d]\n",
               highsBoolToString(option.advanced).c_str(), option.lower_bound,
-              option.upper_bound, option.default_value);
-      fprintf(file, "%s = %" HIGHSINT_FORMAT "\n", option.name.c_str(),
-              *option.value);
+              int(option.upper_bound), int(option.default_value));
+      fprintf(file, "%s = %d\n", option.name.c_str(), int(*option.value));
     } else {
-      fprintf(file, "%s = %" HIGHSINT_FORMAT "\n", option.name.c_str(),
-              *option.value);
+      std::string line = highsFormatToString(
+          "Set option %s to %d\n", option.name.c_str(), int(*option.value));
+      if (file == stdout) {
+        highsLogUser(log_options, HighsLogType::kInfo, "%s", line.c_str());
+      } else {
+        fprintf(file, "%s", line.c_str());
+      }
     }
   }
 }
 
-void reportOption(FILE* file, const OptionRecordDouble& option,
+void reportOption(FILE* file, const HighsLogOptions& log_options,
+                  const OptionRecordDouble& option,
                   const bool report_only_deviations,
                   const HighsFileType file_type) {
   if (!report_only_deviations || option.default_value != *option.value) {
@@ -917,12 +937,19 @@ void reportOption(FILE* file, const OptionRecordDouble& option,
               option.upper_bound, option.default_value);
       fprintf(file, "%s = %g\n", option.name.c_str(), *option.value);
     } else {
-      fprintf(file, "%s = %g\n", option.name.c_str(), *option.value);
+      std::string line = highsFormatToString(
+          "Set option %s to %g\n", option.name.c_str(), *option.value);
+      if (file == stdout) {
+        highsLogUser(log_options, HighsLogType::kInfo, "%s", line.c_str());
+      } else {
+        fprintf(file, "%s", line.c_str());
+      }
     }
   }
 }
 
-void reportOption(FILE* file, const OptionRecordString& option,
+void reportOption(FILE* file, const HighsLogOptions& log_options,
+                  const OptionRecordString& option,
                   const bool report_only_deviations,
                   const HighsFileType file_type) {
   // Don't report for the options file if writing to an options file
@@ -943,7 +970,14 @@ void reportOption(FILE* file, const OptionRecordString& option,
               option.default_value.c_str());
       fprintf(file, "%s = %s\n", option.name.c_str(), (*option.value).c_str());
     } else {
-      fprintf(file, "%s = %s\n", option.name.c_str(), (*option.value).c_str());
+      std::string line =
+          highsFormatToString("Set option %s to \"%s\"\n", option.name.c_str(),
+                              (*option.value).c_str());
+      if (file == stdout) {
+        highsLogUser(log_options, HighsLogType::kInfo, "%s", line.c_str());
+      } else {
+        fprintf(file, "%s", line.c_str());
+      }
     }
   }
 }
diff --git a/src/lp_data/HighsOptions.h b/src/lp_data/HighsOptions.h
index 842f432b00..ceeabbb1cc 100644
--- a/src/lp_data/HighsOptions.h
+++ b/src/lp_data/HighsOptions.h
@@ -233,22 +233,28 @@ OptionStatus getLocalOptionType(
 void resetLocalOptions(std::vector<OptionRecord*>& option_records);
 
 HighsStatus writeOptionsToFile(
-    FILE* file, const std::vector<OptionRecord*>& option_records,
+    FILE* file, const HighsLogOptions& report_log_options,
+    const std::vector<OptionRecord*>& option_records,
     const bool report_only_deviations = false,
     const HighsFileType file_type = HighsFileType::kFull);
-void reportOptions(FILE* file, const std::vector<OptionRecord*>& option_records,
+void reportOptions(FILE* file, const HighsLogOptions& report_log_options,
+                   const std::vector<OptionRecord*>& option_records,
                    const bool report_only_deviations = false,
                    const HighsFileType file_type = HighsFileType::kFull);
-void reportOption(FILE* file, const OptionRecordBool& option,
+void reportOption(FILE* file, const HighsLogOptions& report_log_options,
+                  const OptionRecordBool& option,
                   const bool report_only_deviations,
                   const HighsFileType file_type);
-void reportOption(FILE* file, const OptionRecordInt& option,
+void reportOption(FILE* file, const HighsLogOptions& report_log_options,
+                  const OptionRecordInt& option,
                   const bool report_only_deviations,
                   const HighsFileType file_type);
-void reportOption(FILE* file, const OptionRecordDouble& option,
+void reportOption(FILE* file, const HighsLogOptions& report_log_options,
+                  const OptionRecordDouble& option,
                   const bool report_only_deviations,
                   const HighsFileType file_type);
-void reportOption(FILE* file, const OptionRecordString& option,
+void reportOption(FILE* file, const HighsLogOptions& report_log_options,
+                  const OptionRecordString& option,
                   const bool report_only_deviations,
                   const HighsFileType file_type);
 
diff --git a/src/lp_data/HighsSolution.cpp b/src/lp_data/HighsSolution.cpp
index 5375a3b6b9..a28060861a 100644
--- a/src/lp_data/HighsSolution.cpp
+++ b/src/lp_data/HighsSolution.cpp
@@ -677,7 +677,7 @@ double computeObjectiveValue(const HighsLp& lp, const HighsSolution& solution) {
 // and any solution values
 void refineBasis(const HighsLp& lp, const HighsSolution& solution,
                  HighsBasis& basis) {
-  assert(basis.valid);
+  assert(basis.useful);
   assert(isBasisRightSize(lp, basis));
   const bool have_highs_solution = solution.value_valid;
 
@@ -1319,6 +1319,7 @@ HighsStatus ipxBasicSolutionToHighsBasicSolution(
   highs_solution.value_valid = true;
   highs_solution.dual_valid = true;
   highs_basis.valid = true;
+  highs_basis.useful = true;
   return HighsStatus::kOk;
 }
 
@@ -1356,9 +1357,14 @@ HighsStatus formSimplexLpBasisAndFactor(HighsLpSolverObject& solver_object,
   // If new scaling is performed, the hot start information is
   // no longer valid
   if (new_scaling) ekk_instance.clearHotStart();
-  if (basis.alien) {
-    // An alien basis needs to be checked for rank deficiency, and
-    // possibly completed if it is rectangular
+  const bool check_basis = basis.alien || (!basis.valid && basis.useful);
+  if (check_basis) {
+    // The basis needs to be checked for rank deficiency, and possibly
+    // completed if it is rectangular
+    //
+    // If it's not valid but useful, but not alien,
+    // accommodateAlienBasis will assert, so make the basis alien
+    basis.alien = true;
     assert(!only_from_known_basis);
     accommodateAlienBasis(solver_object);
     basis.alien = false;
@@ -1605,9 +1611,32 @@ void HighsSolution::clear() {
 
 void HighsObjectiveSolution::clear() { this->col_value.clear(); }
 
+void HighsBasis::print(std::string message) const {
+  if (!this->useful) return;
+  this->printScalars(message);
+  for (HighsInt iCol = 0; iCol < HighsInt(this->col_status.size()); iCol++)
+    printf("Basis: col_status[%2d] = %d\n", int(iCol),
+           int(this->col_status[iCol]));
+  for (HighsInt iRow = 0; iRow < HighsInt(this->row_status.size()); iRow++)
+    printf("Basis: row_status[%2d] = %d\n", int(iRow),
+           int(this->row_status[iRow]));
+}
+
+void HighsBasis::printScalars(std::string message) const {
+  printf("\nBasis: %s\n", message.c_str());
+  printf(" valid = %d\n", this->valid);
+  printf(" alien = %d\n", this->alien);
+  printf(" useful = %d\n", this->useful);
+  printf(" was_alien = %d\n", this->was_alien);
+  printf(" debug_id = %d\n", int(this->debug_id));
+  printf(" debug_update_count = %d\n", int(this->debug_update_count));
+  printf(" debug_origin_name = %s\n", this->debug_origin_name.c_str());
+}
+
 void HighsBasis::invalidate() {
   this->valid = false;
   this->alien = true;
+  this->useful = false;
   this->was_alien = true;
   this->debug_id = -1;
   this->debug_update_count = -1;
diff --git a/src/lp_data/HighsSolve.cpp b/src/lp_data/HighsSolve.cpp
index 1062b35e5f..acc81abf43 100644
--- a/src/lp_data/HighsSolve.cpp
+++ b/src/lp_data/HighsSolve.cpp
@@ -390,6 +390,7 @@ HighsStatus solveUnconstrainedLp(const HighsOptions& options, const HighsLp& lp,
   solution.value_valid = true;
   solution.dual_valid = true;
   basis.valid = true;
+  basis.useful = true;
   highs_info.basis_validity = kBasisValidityValid;
   setSolutionStatus(highs_info);
   if (highs_info.num_primal_infeasibilities) {
diff --git a/src/mip/HighsMipSolverData.cpp b/src/mip/HighsMipSolverData.cpp
index 928897c980..4c411f8ffc 100644
--- a/src/mip/HighsMipSolverData.cpp
+++ b/src/mip/HighsMipSolverData.cpp
@@ -1265,6 +1265,7 @@ void HighsMipSolverData::performRestart() {
     root_basis.row_status.resize(postSolveStack.getOrigNumRow(),
                                  HighsBasisStatus::kBasic);
     root_basis.valid = true;
+    root_basis.useful = true;
 
     for (HighsInt i = 0; i < mipsolver.model_->num_col_; ++i)
       root_basis.col_status[postSolveStack.getOrigColIndex(i)] =
@@ -1379,6 +1380,7 @@ void HighsMipSolverData::basisTransfer() {
     firstrootbasis.row_status.assign(numRow, HighsBasisStatus::kNonbasic);
     firstrootbasis.valid = true;
     firstrootbasis.alien = true;
+    firstrootbasis.useful = true;
 
     for (HighsInt i = 0; i < numRow; ++i) {
       HighsBasisStatus status =
@@ -1929,6 +1931,7 @@ void HighsMipSolverData::evaluateRootNode() {
     firstrootbasis.row_status.assign(mipsolver.numRow(),
                                      HighsBasisStatus::kBasic);
     firstrootbasis.valid = true;
+    firstrootbasis.useful = true;
   }
 
   if (cutpool.getNumCuts() != 0) {
diff --git a/src/presolve/HPresolve.cpp b/src/presolve/HPresolve.cpp
index eb075e53fb..6c72a06046 100644
--- a/src/presolve/HPresolve.cpp
+++ b/src/presolve/HPresolve.cpp
@@ -6307,8 +6307,7 @@ void HPresolve::debug(const HighsLp& lp, const HighsOptions& options) {
 
     sol = reducedsol;
     basis = reducedbasis;
-    postsolve_stack.undoUntil(options, flagRow, flagCol, sol, basis,
-                              tmp.numReductions());
+    postsolve_stack.undoUntil(options, sol, basis, tmp.numReductions());
 
     HighsBasis temp_basis;
     HighsSolution temp_sol;
@@ -6330,6 +6329,7 @@ void HPresolve::debug(const HighsLp& lp, const HighsOptions& options) {
     temp_sol.row_value.resize(model.num_row_);
     calculateRowValuesQuad(model, sol);
     temp_basis.valid = true;
+    temp_basis.useful = true;
     refineBasis(model, temp_sol, temp_basis);
     Highs highs;
     highs.passOptions(options);
@@ -6353,8 +6353,7 @@ void HPresolve::debug(const HighsLp& lp, const HighsOptions& options) {
                          ARstart, ARindex, ARvalue);
     sol = reducedsol;
     basis = reducedbasis;
-    postsolve_stack.undoUntil(options, flagRow, flagCol, sol, basis,
-                              reductionLim);
+    postsolve_stack.undoUntil(options, sol, basis, reductionLim);
 
     calculateRowValuesQuad(model, sol);
     kktinfo = dev_kkt_check::initInfo();
diff --git a/src/presolve/HighsPostsolveStack.h b/src/presolve/HighsPostsolveStack.h
index bb1b1f785b..e324af710b 100644
--- a/src/presolve/HighsPostsolveStack.h
+++ b/src/presolve/HighsPostsolveStack.h
@@ -784,9 +784,7 @@ class HighsPostsolveStack {
   */
 
   // Only used for debugging
-  void undoUntil(const HighsOptions& options,
-                 const std::vector<HighsInt>& flagRow,
-                 const std::vector<HighsInt>& flagCol, HighsSolution& solution,
+  void undoUntil(const HighsOptions& options, HighsSolution& solution,
                  HighsBasis& basis, size_t numReductions) {
     reductionValues.resetPosition();
 
diff --git a/src/qpsolver/a_quass.cpp b/src/qpsolver/a_quass.cpp
index 56ff800596..d60bc999e5 100644
--- a/src/qpsolver/a_quass.cpp
+++ b/src/qpsolver/a_quass.cpp
@@ -116,6 +116,7 @@ static QpAsmStatus quass2highs(Instance& instance, Settings& settings,
   }
   highs_basis.valid = true;
   highs_basis.alien = false;
+  highs_basis.useful = true;
   return qp_asm_return_status;
 }
 
diff --git a/src/simplex/HApp.h b/src/simplex/HApp.h
index 9fc35589b4..178593cbc6 100644
--- a/src/simplex/HApp.h
+++ b/src/simplex/HApp.h
@@ -136,6 +136,22 @@ inline HighsStatus solveLpSimplex(HighsLpSolverObject& solver_object) {
   // If new scaling is performed, the hot start information is
   // no longer valid
   if (new_scaling) ekk_instance.clearHotStart();
+  //
+  if (!status.has_basis && !basis.valid && basis.useful) {
+    // There is no simplex basis, but there is a useful HiGHS basis
+    // that is not validated
+    assert(basis.col_status.size() ==
+           static_cast<size_t>(incumbent_lp.num_col_));
+    assert(basis.row_status.size() ==
+           static_cast<size_t>(incumbent_lp.num_row_));
+    HighsStatus return_status = formSimplexLpBasisAndFactor(solver_object);
+    if (return_status != HighsStatus::kOk)
+      return returnFromSolveLpSimplex(solver_object, HighsStatus::kError);
+    // formSimplexLpBasisAndFactor may introduce variables with
+    // HighsBasisStatus::kNonbasic, so refine it
+    refineBasis(incumbent_lp, solution, basis);
+    basis.valid = true;
+  }
   // Move the LP to EKK, updating other EKK pointers and any simplex
   // NLA pointers, since they may have moved if the LP has been
   // modified
diff --git a/src/simplex/HEkk.cpp b/src/simplex/HEkk.cpp
index 2dd2ff85ae..c6da49da49 100644
--- a/src/simplex/HEkk.cpp
+++ b/src/simplex/HEkk.cpp
@@ -1465,6 +1465,7 @@ HighsBasis HEkk::getHighsBasis(HighsLp& use_lp) const {
   }
   highs_basis.valid = true;
   highs_basis.alien = false;
+  highs_basis.useful = true;
   highs_basis.was_alien = false;
   highs_basis.debug_id =
       (HighsInt)(build_synthetic_tick_ + total_synthetic_tick_);