diff --git a/.gitignore b/.gitignore
index 88fa19b4c..b4a858689 100644
--- a/.gitignore
+++ b/.gitignore
@@ -12,5 +12,9 @@ Testing
 *.pyc
 *.user
 *~
+
+.idea
+*.csv
+
 .vscode
 /src/expr-info.cc
diff --git a/solvers/mosek/main.cc b/solvers/mosek/main.cc
index 8aadb02b2..7fabdc46b 100644
--- a/solvers/mosek/main.cc
+++ b/solvers/mosek/main.cc
@@ -4,8 +4,7 @@
 std::unique_ptr<mp::BasicBackend> CreateMosekBackend();
 
 extern "C" int main1(int, char **argv) {
-  return
-      mp::RunBackendApp(argv, CreateMosekBackend);
+  return mp::RunBackendApp(argv, CreateMosekBackend);
 }
 
 extern "C" int main2(int, char** argv, CCallbacks cb) {
diff --git a/solvers/mosek/mosekbackend.cc b/solvers/mosek/mosekbackend.cc
index 400b1d950..e49548754 100644
--- a/solvers/mosek/mosekbackend.cc
+++ b/solvers/mosek/mosekbackend.cc
@@ -15,7 +15,7 @@ namespace {
 
 
 bool InterruptMosek(void* prob) {
-  // TODO 
+  // TODO
   return true;
 }
 
@@ -50,8 +50,7 @@ MosekBackend::~MosekBackend() {
   CloseSolver();
 }
 
-static void MSKAPI printstr(void* handle,
-  const char* str)
+static void MSKAPI printstr(void* handle, const char* str)
 {
   fmt::print("{}", str);
   fflush(stdout);
@@ -63,7 +62,7 @@ void MosekBackend::OpenSolver() {
 
   MSKenv_t env = NULL;
   MSKtask_t task;
-  if (initialize) { // If an initialization callback is provided, 
+  if (initialize) { // If an initialization callback is provided,
     // use it to create the environment
     MOSEK_CCALL(MSK_makeenv(&env, NULL));
     env = (MSKenv_t)initialize();
@@ -98,13 +97,15 @@ std::string MosekBackend::GetSolverVersion() {
   return fmt::format("{}.{}.{}", MSK_VERSION_MAJOR,
     MSK_VERSION_MINOR, MSK_VERSION_REVISION);
 }
+
 void MosekBackend::InitOptionParsing() {
   OpenSolver();
 }
 
 bool MosekBackend::IsMIP() const {
-  return getIntAttr(MSK_IINF_ANA_PRO_NUM_VAR_BIN)+
-    getIntAttr(MSK_IINF_ANA_PRO_NUM_VAR_INT);
+  int is_int_sol = 0;
+  MOSEK_CCALL(MSK_solutiondef(lp(), MSK_SOL_ITG, &is_int_sol));
+  return is_int_sol;
 }
 
 bool MosekBackend::IsQCP() const {
@@ -173,45 +174,32 @@ void MosekBackend::SetInterrupter(mp::Interrupter *inter) {
 }
 
 MSKsoltypee MosekBackend::GetSolutionTypeToFetch() {
+  int b = 0;
+
+  // Handle MIP
+  MOSEK_CCALL(MSK_solutiondef(lp(), MSK_SOL_ITG, &b));
+  if (b) return MSK_SOL_ITG;
+
+  // Handle LP
   MSKproblemtypee type;
-  MSK_getprobtype(lp(), &type);
-  int optimizer;
-  GetSolverOption(MSK_IPAR_OPTIMIZER, optimizer);
-  MSKoptimizertypee optype = (MSKoptimizertypee)optimizer;
-
-  // TODO Implement proper logic here
-  if (IsMIP())
-    return MSK_SOL_ITG;
-
-  if (optype == MSK_OPTIMIZER_FREE) { // decide by problem type
-    switch (type) {
-      case MSK_PROBTYPE_LO:
-        return MSK_SOL_BAS;
-      default:
-        return MSK_SOL_ITR;
-    }
-  }
-  else { // if algorithm is specifically chosen
-    switch(optype)
-    {
-    case MSK_OPTIMIZER_CONIC:
-    case MSK_OPTIMIZER_INTPNT:
-      return MSK_SOL_ITR;
-    case MSK_OPTIMIZER_MIXED_INT:
-      return MSK_SOL_ITG;
-    default: // all simplex
-      return MSK_SOL_BAS;
-    }
+  MOSEK_CCALL(MSK_getprobtype(lp(), &type));
+  if (type == MSK_PROBTYPE_LO)
+  {
+    MOSEK_CCALL(MSK_solutiondef(lp(), MSK_SOL_BAS, &b));
+    if (b) return MSK_SOL_BAS;
   }
 
+  // All other case
+  return MSK_SOL_ITR;
 }
+
 void MosekBackend::Solve() {
   if (!storedOptions_.exportFile_.empty()) {
     ExportModel(storedOptions_.exportFile_);
   }
   MOSEK_CCALL(MSK_optimizetrm(lp(), &termCode_));
   solToFetch_ = GetSolutionTypeToFetch();
-  MSK_getsolsta(lp(), solToFetch_, &solSta_); 
+  MOSEK_CCALL(MSK_getsolsta(lp(), solToFetch_, &solSta_));
   WindupMOSEKSolve();
 }
 
@@ -228,18 +216,21 @@ void MosekBackend::ReportMOSEKResults() {
   if (need_multiple_solutions())
     ReportMOSEKPool();
 }
+
 std::vector<double> MosekBackend::getPoolSolution(int i)
 {
   std::vector<double> vars(NumVars());
   // TODO get solutions in the pool 
   return vars;
 }
+
 double MosekBackend::getPoolObjective(int i)
 {
   double obj;
   // TODO get objective value of solution i
   return obj;
 }
+
 void MosekBackend::ReportMOSEKPool() {
   if (!IsMIP())
     return;
@@ -284,8 +275,9 @@ std::pair<int, std::string> MosekBackend::ConvertMOSEKStatus() {
     case MSK_SOL_STA_PRIM_AND_DUAL_FEAS:
       return { sol::UNCERTAIN, "feasible solution" };
     case MSK_SOL_STA_PRIM_INFEAS_CER:
+      return { sol::INFEASIBLE, "primal infeasible solution" };
     case MSK_SOL_STA_DUAL_INFEAS_CER:
-      return { sol::INFEASIBLE, "unfeasible solution" };
+      return { sol::INF_OR_UNB, "dual infeasible solution" };
     case MSK_SOL_STA_UNKNOWN:
     case MSK_SOL_STA_PRIM_ILLPOSED_CER:
     case MSK_SOL_STA_DUAL_ILLPOSED_CER:
@@ -320,7 +312,7 @@ static const mp::OptionValueInfo alg_values_method[] = {
   { "1", "Dual simplex", 1},
   { "2", "Automatic (default)", 2},
   { "3", "Free simplex", 3},
-  { "4", "Interior-poit method", 4},
+  { "4", "Interior-point method", 4},
   { "5", "Mixed-integer optimizer", 5},
   { "6", "Primal simplex", 6},
 };
@@ -350,163 +342,392 @@ void MosekBackend::InitCustomOptions() {
       "Default = \"\" (don't export the model).",
       storedOptions_.exportFile_);
 
-  // Example of direct solver option (set directly by the framework)
+  AddStoredOption("mip:constructsol mipconstructsol",
+      "Sets MSK_IPAR_MIO_CONSTRUCT_SOL. If set to MSK_ON and all integer variables "
+      "have been given a value for which a feasible mixed integer solution exists, "
+      "then MOSEK generates an initial solution to the mixed integer problem by "
+      "fixing all integer values and solving the remaining problem."
+      "Default = OFF",
+      storedOptions_.MIPConstructSol_);
+
   AddSolverOption("tech:threads threads",
     "Controls the number of threads employed by the optimizer. "
     "Default 0 ==> number of threads used will be equal to the number "
     "of cores detected on the machine.",
     MSK_IPAR_NUM_THREADS, 0, INT_MAX);
 
+  AddSolverOption("mip:relgapconst miorelgapconst",
+    "This value is used to compute the relative gap for the solution "
+    "to an integer optimization problem."
+    "Default = 1.0e-10",
+    MSK_DPAR_MIO_REL_GAP_CONST, 0.0, DBL_MAX);
+
   AddSolverOption("tech:outlev outlev",
     "0*/1: Whether to write mosek log lines to stdout.",
     MSK_IPAR_LOG, 0, 1);
 }
 
 double MosekBackend::MIPGap() {
-  return getDblAttr(MSK_DINF_MIO_OBJ_ABS_GAP);
+  return getDblAttr(MSK_DINF_MIO_OBJ_REL_GAP);
 }
+
 double MosekBackend::BestDualBound() {
-  // TODO
-  return 0;
-  //return getDblAttr(MOSEK_DBLATTR_BESTBND);
+  return getDblAttr(MSK_DINF_MIO_OBJ_BOUND);
 }
 
 double MosekBackend::MIPGapAbs() {
-  return std::fabs(
-    ObjectiveValue() - BestDualBound());
+  return getDblAttr(MSK_DINF_MIO_OBJ_ABS_GAP);
 }
 
+SensRanges MosekBackend::GetSensRanges()
+{
+  SensRanges sensr;
+
+  // Set sensitivity range parameters for constraints
+  int lencon = NumLinCons();
+  std::vector<MSKint32t> cindex(lencon);
+  std::vector<MSKmarke> cmarklb(lencon);
+  std::vector<MSKmarke> cmarkub(lencon);
+  for (int i=0; i<lencon; i++)
+  {
+    cindex[i] = i;
+    cmarklb[i] = MSK_MARK_LO;
+    cmarkub[i] = MSK_MARK_UP;
+  }
+  std::vector<MSKrealt> cpricelblo(lencon);
+  std::vector<MSKrealt> cpricelbhi(lencon);
+  std::vector<MSKrealt> cpriceublo(lencon);
+  std::vector<MSKrealt> cpriceubhi(lencon);
+  std::vector<MSKrealt> crangelblo(lencon);
+  std::vector<MSKrealt> crangelbhi(lencon);
+  std::vector<MSKrealt> crangeublo(lencon);
+  std::vector<MSKrealt> crangeubhi(lencon);
+
+  // Set sensitivity range parameters for variables
+  int lenvar = NumVars();
+  std::vector<MSKint32t> vindex(lenvar);
+  std::vector<MSKmarke> vmarklb(lenvar);
+  std::vector<MSKmarke> vmarkub(lenvar);
+  for (int i=0; i<lenvar; i++)
+  {
+    vindex[i] = i;
+    vmarklb[i] = MSK_MARK_LO;
+    vmarkub[i] = MSK_MARK_UP;
+  }
+  std::vector<MSKrealt> vpricelblo(lenvar);
+  std::vector<MSKrealt> vpricelbhi(lenvar);
+  std::vector<MSKrealt> vpriceublo(lenvar);
+  std::vector<MSKrealt> vpriceubhi(lenvar);
+  std::vector<MSKrealt> vrangelblo(lenvar);
+  std::vector<MSKrealt> vrangelbhi(lenvar);
+  std::vector<MSKrealt> vrangeublo(lenvar);
+  std::vector<MSKrealt> vrangeubhi(lenvar);
+
+  std::vector<MSKrealt> opricelo(lenvar);
+  std::vector<MSKrealt> opricehi(lenvar);
+  std::vector<MSKrealt> orangelo(lenvar);
+  std::vector<MSKrealt> orangehi(lenvar);
+
+  // Compute sensitivity ranges for lower bounds
+  MOSEK_CCALL(
+    MSK_primalsensitivity(
+        lp(),
+        lencon,
+        cindex.data(),
+        cmarklb.data(),
+        lenvar,
+        vindex.data(),
+        vmarklb.data(),
+        cpricelblo.data(),
+        cpricelbhi.data(),
+        crangelblo.data(),
+        crangelbhi.data(),
+        vpricelblo.data(),
+        vpricelbhi.data(),
+        vrangelblo.data(),
+        vrangelbhi.data()
+    )
+  );
+
+  // Compute sensitivity ranges for upper bounds
+  MOSEK_CCALL(
+    MSK_primalsensitivity(
+        lp(),
+        lencon,
+        cindex.data(),
+        cmarkub.data(),
+        lenvar,
+        vindex.data(),
+        vmarkub.data(),
+        cpriceublo.data(),
+        cpriceubhi.data(),
+        crangeublo.data(),
+        crangeubhi.data(),
+        vpriceublo.data(),
+        vpriceubhi.data(),
+        vrangeublo.data(),
+        vrangeubhi.data()
+    )
+  );
+
+  // Compute sensitivity ranges for objective coeffs
+  MOSEK_CCALL(
+    MSK_dualsensitivity(
+        lp(),
+        lenvar,
+        vindex.data(),
+        opricelo.data(),
+        opricehi.data(),
+        orangelo.data(),
+        orangehi.data()
+    )
+  );
+
+  // Get bounds and add them to sensitivity ranges
+  std::vector<MSKboundkeye> vbk(lenvar);
+  std::vector<MSKrealt> vbl(lenvar), vbu(lenvar);
+  MOSEK_CCALL(MSK_getvarboundslice(lp(), (MSKint32t)0, (MSKint32t)lenvar, vbk.data(), vbl.data(), vbu.data()));
+  std::transform(vrangelblo.begin(), vrangelblo.end(), vbl.begin(), vrangelblo.begin(), std::plus<MSKrealt>());
+  std::transform(vrangelbhi.begin(), vrangelbhi.end(), vbl.begin(), vrangelbhi.begin(), std::plus<MSKrealt>());
+  std::transform(vrangeublo.begin(), vrangeublo.end(), vbu.begin(), vrangeublo.begin(), std::plus<MSKrealt>());
+  std::transform(vrangeubhi.begin(), vrangeubhi.end(), vbu.begin(), vrangeubhi.begin(), std::plus<MSKrealt>());
+  for (int i=0; i<lenvar; i++)
+  {
+    if (vbk[i] == MSK_BK_UP || vbk[i] == MSK_BK_FR)
+    {
+      vrangelblo[i] = -MSK_INFINITY;
+      vrangelbhi[i] = MSK_INFINITY;
+    }
+    if (vbk[i] == MSK_BK_LO || vbk[i] == MSK_BK_FR)
+    {
+      vrangeublo[i] = -MSK_INFINITY;
+      vrangeubhi[i] = MSK_INFINITY;
+    }
+  }
+
+  std::vector<MSKboundkeye> cbk(lencon);
+  std::vector<MSKrealt> cbl(lencon), cbu(lencon);
+  MOSEK_CCALL(MSK_getconboundslice(lp(), (MSKint32t)0, (MSKint32t)lencon, cbk.data(), cbl.data(), cbu.data()));
+  std::transform(crangelblo.begin(), crangelblo.end(), cbl.begin(), crangelblo.begin(), std::plus<MSKrealt>());
+  std::transform(crangelbhi.begin(), crangelbhi.end(), cbl.begin(), crangelbhi.begin(), std::plus<MSKrealt>());
+  std::transform(crangeublo.begin(), crangeublo.end(), cbu.begin(), crangeublo.begin(), std::plus<MSKrealt>());
+  std::transform(crangeubhi.begin(), crangeubhi.end(), cbu.begin(), crangeubhi.begin(), std::plus<MSKrealt>());
+  for (int i=0; i<lencon; i++)
+  {
+    if (cbk[i] == MSK_BK_UP || cbk[i] == MSK_BK_FR)
+    {
+      crangelblo[i] = -MSK_INFINITY;
+      crangelbhi[i] = MSK_INFINITY;
+    }
+    if (cbk[i] == MSK_BK_LO || cbk[i] == MSK_BK_FR)
+    {
+      crangeublo[i] = -MSK_INFINITY;
+      crangeubhi[i] = MSK_INFINITY;
+    }
+  }
+
+  std::vector<MSKrealt> c(lenvar);
+  MOSEK_CCALL(MSK_getc(lp(), c.data()));
+  std::transform(orangelo.begin(), orangelo.end(), c.begin(), orangelo.begin(), std::plus<MSKrealt>());
+  std::transform(orangehi.begin(), orangehi.end(), c.begin(), orangehi.begin(), std::plus<MSKrealt>());
+
+  // Return sensitivity ranges
+  auto mvlbhi = GetValuePresolver().PostsolveGenericDbl({ {vrangelbhi} });
+  auto mvlblo = GetValuePresolver().PostsolveGenericDbl({ {vrangelblo} });
+  auto mvubhi = GetValuePresolver().PostsolveGenericDbl({ {vrangeubhi} });
+  auto mvublo = GetValuePresolver().PostsolveGenericDbl({ {vrangeublo} });
+  auto mvobjhi = GetValuePresolver().PostsolveGenericDbl({ {orangehi} });
+  auto mvobjlo = GetValuePresolver().PostsolveGenericDbl({ {orangelo} });
+  auto mvconlbhi = GetValuePresolver().PostsolveGenericDbl({ {}, {{{CG_Linear, crangelbhi}}} });
+  auto mvconlblo = GetValuePresolver().PostsolveGenericDbl({ {}, {{{CG_Linear, crangelblo}}} });
+  auto mvconubhi = GetValuePresolver().PostsolveGenericDbl({ {}, {{{CG_Linear, crangeubhi}}} });
+  auto mvconublo = GetValuePresolver().PostsolveGenericDbl({ {}, {{{CG_Linear, crangeublo}}} });
+  std::vector<MSKrealt> rhs(lencon, 0.0);
+  auto mvrhshi = GetValuePresolver().PostsolveGenericDbl({ {}, {{{CG_Linear, rhs}}} });
+  auto mvrhslo = GetValuePresolver().PostsolveGenericDbl({ {}, {{{CG_Linear, rhs}}} });
+
+  sensr.varlbhi = mvlbhi.GetVarValues()();
+  sensr.varlblo = mvlblo.GetVarValues()();
+  sensr.varubhi = mvubhi.GetVarValues()();
+  sensr.varublo = mvublo.GetVarValues()();
+  sensr.varobjhi = mvobjhi.GetVarValues()();
+  sensr.varobjlo = mvobjlo.GetVarValues()();
+  sensr.conlbhi = mvconlbhi.GetConValues()();
+  sensr.conlblo = mvconlblo.GetConValues()();
+  sensr.conubhi = mvconubhi.GetConValues()();
+  sensr.conublo = mvconublo.GetConValues()();
+  sensr.conrhshi = mvrhshi.GetConValues()();
+  sensr.conrhslo = mvrhslo.GetConValues()();
+
+  return sensr;
+}
 
 ArrayRef<int> MosekBackend::VarStatii() {
-  // TODO 
   std::vector<int> vars(NumVars());
-  /*
-  MOSEK_GetBasis(lp(), vars.data(), NULL);
-  for (auto& s : vars) {
-    switch (s) {
-    case MOSEK_BASIS_BASIC:
-      s = (int)BasicStatus::bas;
-      break;
-    case MOSEK_BASIS_LOWER:
-      s = (int)BasicStatus::low;
-      break;
-    case MOSEK_BASIS_UPPER:
-      s = (int)BasicStatus::upp;
-      break;
-    case MOSEK_BASIS_SUPERBASIC:
-      s = (int)BasicStatus::sup;
-      break;
-    case MOSEK_BASIS_FIXED:
-      s = (int)BasicStatus::equ;
-      break;
-    default:
-      MP_RAISE(fmt::format("Unknown Mosek VBasis value: {}", s));
+  MOSEK_CCALL(MSK_getskx(lp(), solToFetch_, (MSKstakeye *)vars.data()));
+
+  for (auto& s : vars)
+  {
+    switch (s)
+    {
+      case MSK_SK_BAS:
+        s = (int)BasicStatus::bas;
+        break;
+      case MSK_SK_LOW:
+        s = (int)BasicStatus::low;
+        break;
+      case MSK_SK_UPR:
+        s = (int)BasicStatus::upp;
+        break;
+      case MSK_SK_SUPBAS:
+        s = (int)BasicStatus::sup;
+        break;
+      case MSK_SK_FIX:
+        s = (int)BasicStatus::equ;
+        break;
+      default:
+        MP_RAISE(fmt::format("Unknown Mosek VBasis value: {}", s));
     }
   }
-  */
+
   return vars;
 }
 
 ArrayRef<int> MosekBackend::ConStatii() {
-  // TODO
   std::vector<int> cons(NumLinCons());
-  /*
-  MOSEK_GetBasis(lp(), NULL, cons.data());
-  for (auto& s : cons) {
-    switch (s) {
-    case MOSEK_BASIS_BASIC:
-      s = (int)BasicStatus::bas;
-      break;
-    case MOSEK_BASIS_LOWER:
-      s = (int)BasicStatus::low;
-      break;
-    case MOSEK_BASIS_UPPER:
-      s = (int)BasicStatus::upp;
-      break;
-    case MOSEK_BASIS_SUPERBASIC:
-      s = (int)BasicStatus::sup;
-      break;
-    case MOSEK_BASIS_FIXED:
-      s = (int)BasicStatus::equ;
-      break;
-    default:
-      MP_RAISE(fmt::format("Unknown Mosek VBasis value: {}", s));
+  MOSEK_CCALL(MSK_getskc(lp(), solToFetch_, (MSKstakeye *)cons.data()));
+
+  for (auto& s : cons)
+  {
+    switch (s)
+    {
+      case MSK_SK_BAS:
+        s = (int)BasicStatus::bas;
+        break;
+      case MSK_SK_LOW:
+        s = (int)BasicStatus::low;
+        break;
+      case MSK_SK_UPR:
+        s = (int)BasicStatus::upp;
+        break;
+      case MSK_SK_SUPBAS:
+        s = (int)BasicStatus::sup;
+        break;
+      case MSK_SK_FIX:
+        s = (int)BasicStatus::equ;
+        break;
+      default:
+        MP_RAISE(fmt::format("Unknown Mosek VBasis value: {}", s));
     }
-  }*/
+  }
   return cons;
 }
 
 void MosekBackend::VarStatii(ArrayRef<int> vst) {
-  // TODO
-  int index[1];
-  std::vector<int> stt(vst.data(), vst.data() + vst.size());
-  /*
-  for (auto j = stt.size(); j--; ) {
-    auto& s = stt[j];
-    switch ((BasicStatus)s) {
-    case BasicStatus::bas:
-      s = MOSEK_BASIS_BASIC;
-      break;
-    case BasicStatus::low:
-      s = MOSEK_BASIS_LOWER;
-      break;
-    case BasicStatus::equ:
-      s = MOSEK_BASIS_FIXED;
-      break;
-    case BasicStatus::upp:
-      s = MOSEK_BASIS_UPPER;
-      break;
-    case BasicStatus::sup:
-    case BasicStatus::btw:
-      s = MOSEK_BASIS_SUPERBASIC;
-      break;
-    case BasicStatus::none:
-      /// 'none' is assigned to new variables. Compute low/upp/sup:
-      /// Depending on where 0.0 is between bounds
-      double lb, ub;
-      index[0] = (int)j;
-      if(!MOSEK_GetColInfo(lp(), MOSEK_DBLINFO_LB, 1, index, &lb) && 
-        !MOSEK_GetColInfo(lp(), MOSEK_DBLINFO_UB, 1, index, &ub))
-      { 
-        if (lb >= -1e-6)
-          s = -1;
-        else if (ub <= 1e-6)
-          s = -2;
-        else
-          s = -3;  // or, leave at 0?
-      }
-      break;
-    default:
-      MP_RAISE(fmt::format("Unknown AMPL var status value: {}", s));
+  // BasicStatus enum values: from 0 to 6: none, bas, sup, low, upp, equ, btw
+  std::vector<int> skx(vst.data(), vst.data() + vst.size());
+
+  // TODO: remove this if we are sure that it is the same.
+  MSKint32t numvar;
+  MOSEK_CCALL(MSK_getnumvar(lp(), &numvar));
+  assert(numvar == vst.size());
+
+  for (auto j = 0; j < skx.size(); j++)
+  {
+    switch ((BasicStatus)skx[j])
+    {
+      case BasicStatus::bas:
+        skx[j] = MSK_SK_BAS;
+        break;
+      case BasicStatus::low:
+        skx[j] = MSK_SK_LOW;
+        break;
+      case BasicStatus::equ:
+        skx[j] = MSK_SK_FIX;
+        break;
+      case BasicStatus::upp:
+        skx[j] = MSK_SK_UPR;
+        break;
+      case BasicStatus::sup:
+      case BasicStatus::btw:
+        skx[j] = MSK_SK_SUPBAS;
+        break;
+      case BasicStatus::none:
+        /// 'none' is assigned to new variables (added after solving).
+        /// Compute low/upp/fix/sup, depending on where initial value 0.0 is between the bounds
+        MSKboundkeye bk;
+        MSKrealt bl, bu, tolx;
+        MOSEK_CCALL(MSK_getvarbound(lp(), (MSKint32t)j, &bk, &bl, &bu));
+        MOSEK_CCALL(MSK_getdouparam(lp(), MSK_DPAR_DATA_TOL_X, &tolx));
+
+        if (0.0 < bl + tolx) skx[j] = MSK_SK_LOW;
+        else if (bu - tolx < 0.0) skx[j] = MSK_SK_UPR;
+        else if (bu - bl < tolx) skx[j] = MSK_SK_FIX;
+        else skx[j] = MSK_SK_SUPBAS;
+
+        break;
+      default:
+        MP_RAISE(fmt::format("Unknown AMPL var status value: {}", skx[j]));
     }
   }
-  MOSEK_SetBasis(lp(), stt.data(), NULL);
-  */
+
+  // Input the variable status keys
+  MOSEK_CCALL(MSK_putskx(lp(), solToFetch_, (MSKstakeye *)skx.data()));
 }
 
 void MosekBackend::ConStatii(ArrayRef<int> cst) {
-  // TODO
-  /*
-  std::vector<int> stt(cst.data(), cst.data() + cst.size());
-  for (auto& s : stt) {
-    switch ((BasicStatus)s) {
-    case BasicStatus::bas:
-      s = MOSEK_BASIS_BASIC;
-      break;
-    case BasicStatus::none:   // for 'none', which is the status
-    case BasicStatus::upp:    // assigned to new rows, it seems good to guess
-    case BasicStatus::sup:    // a valid status.
-    case BasicStatus::low:    // 
-    case BasicStatus::equ:    // For active constraints, it is usually 'sup'.
-    case BasicStatus::btw:    // We could compute slack to decide though.
-      s = MOSEK_BASIS_SUPERBASIC;
-      break;
-    default:
-      MP_RAISE(fmt::format("Unknown AMPL con status value: {}", s));
+  std::vector<int> skc(cst.data(), cst.data() + cst.size());
+
+  // TODO: remove this if we are sure that it is the same.
+  MSKint32t numcon;
+  MOSEK_CCALL(MSK_getnumcon(lp(), &numcon));
+  assert(numcon == cst.size());
+
+  // The status key of a constraint is the status key of the logical (slack) variable assigned to it.
+  // The slack is defined as: l <= a'x <= u rewritten as a'x - s = 0, l <= s <= u.
+  for (auto j = 0; j < skc.size(); j++)
+  {
+    skc[j] = MSK_SK_UNK;
+    switch ((BasicStatus)skc[j])
+    {
+      case BasicStatus::bas:
+        skc[j] = MSK_SK_BAS;
+        break;
+      case BasicStatus::sup:
+      case BasicStatus::btw:
+        skc[j] = MSK_SK_SUPBAS;
+        break;
+      case BasicStatus::low:
+        skc[j] = MSK_SK_LOW;
+        break;
+      case BasicStatus::upp:
+        skc[j] = MSK_SK_UPR;
+        break;
+      case BasicStatus::equ:
+        skc[j] = MSK_SK_FIX;
+        break;
+      case BasicStatus::none:
+        /// 'none' is assigned to new constraints (added after solving).
+        /// Compute low/upp/fix/sup, depending on where initial value 0.0 is between the bounds
+        MSKboundkeye bk;
+        MSKrealt bl, bu, tolx;
+        MOSEK_CCALL(MSK_getconbound(lp(), (MSKint32t)j, &bk, &bl, &bu));
+        // MSK_DPAR_DATA_TOL_X is both for variables and constraints
+        MOSEK_CCALL(MSK_getdouparam(lp(), MSK_DPAR_DATA_TOL_X, &tolx));
+
+        if (0.0 < bl + tolx) skc[j] = MSK_SK_LOW;
+        else if (bu - tolx < 0.0) skc[j] = MSK_SK_UPR;
+        else if (bu - bl < tolx) skc[j] = MSK_SK_FIX;
+        else skc[j] = MSK_SK_SUPBAS;
+
+        break;
+      default:
+        MP_RAISE(fmt::format("Unknown AMPL var status value: {}", skc[j]));
     }
   }
-  MOSEK_SetBasis(lp(), NULL, stt.data());
-  */
+
+  // Input the variable status keys
+  MOSEK_CCALL(MSK_putskc(lp(), solToFetch_, (MSKstakeye *)skc.data()));
 }
 
 SolutionBasis MosekBackend::GetBasis() {
@@ -535,13 +756,66 @@ void MosekBackend::SetBasis(SolutionBasis basis) {
   ConStatii(constt);
 }
 
+void MosekBackend::AddPrimalDualStart(Solution sol)
+{
+  auto mv = GetValuePresolver().PresolveSolution(
+        { sol.primal, sol.dual } );
+  auto x0 = mv.GetVarValues()();
+  auto pi0 = mv.GetConValues()(CG_Linear);
+  MOSEK_CCALL(MSK_putxx(lp(), solToFetch_, (MSKrealt *)x0.data()));
+  MOSEK_CCALL(MSK_puty(lp(), solToFetch_, (MSKrealt *)pi0.data()));
+}
+
+void MosekBackend::AddMIPStart(ArrayRef<double> x0_unpres)
+{
+  auto mv = GetValuePresolver().PresolveSolution( { x0_unpres } );
+  auto x0 = mv.GetVarValues()();
+
+  if (Mosek_mip_construct_sol())
+  {
+    // Use integer part of given solution to construct initial solution by solving for continuous part.
+    MOSEK_CCALL(MSK_putintparam(lp(), MSK_IPAR_MIO_CONSTRUCT_SOL, MSK_ON));
+  }
+  MOSEK_CCALL(MSK_putxx(lp(), solToFetch_, (MSKrealt *)x0.data()));
+}
+
+ArrayRef<double> MosekBackend::Ray()
+{
+  // Problem is checked to be dual infeasible at this point, so the ray is the primal solution variable.
+  // (Primal can be unbounded or infeasible.)
+  std::vector<double> xx(NumVars());
+  MOSEK_CCALL(MSK_getxx(lp(), solToFetch_, (MSKrealt *)xx.data()));
+  // Argument is a ModelValues<ValueMap>, which can be constructed from a ValueMap for only the variables.
+  auto mv = GetValuePresolver().PostsolveSolution({xx});
+  // We get the ValueMap for variables, and by calling that, the value itself.
+  // (Similar to .MoveOut(0) for single key maps)
+  auto uray = mv.GetVarValues()();
+  return uray;
+}
+
+ArrayRef<double> MosekBackend::DRay()
+{
+  // Problem is checked to be (primal) infeasible at this point, so the ray is the dual solution variable.
+  // (Dual can be unbounded or infeasible.)
+  std::vector<double> y(NumLinCons());
+  MOSEK_CCALL(MSK_gety(lp(), solToFetch_, (MSKrealt *)y.data()));
+  // Argument is a ModelValues<ValueMap>, which is constructed now from two ValueMaps, an empty for variables, and
+  // a second one for constraints. The constraints ValueMap is given as a std::map for only linear constraints
+  auto mv = GetValuePresolver().PostsolveSolution({
+                                                    {},
+                                                    { { {CG_Linear, std::move(y)} } }
+                                                  });
+  // We get the ValueMap for constraints, and by calling MoveOut(), the value itself.
+  // (Similar to .MoveOut(0) for single key maps)
+  return mv.GetConValues().MoveOut();
+}
+
 } // namespace mp
 
 
 // AMPLs
 
-AMPLS_MP_Solver* AMPLSOpenMosek(
-  const char* slv_opt) {
+AMPLS_MP_Solver* AMPLSOpenMosek(const char* slv_opt) {
   return AMPLS__internal__Open(std::unique_ptr<mp::BasicBackend>{new mp::MosekBackend()},
     slv_opt);
 }
@@ -551,6 +825,5 @@ void AMPLSCloseMosek(AMPLS_MP_Solver* slv) {
 }
 
 MSKtask_t GetMosekmodel(AMPLS_MP_Solver* slv) {
-  return
-    dynamic_cast<mp::MosekBackend*>(AMPLSGetBackend(slv))->lp();
+  return dynamic_cast<mp::MosekBackend*>(AMPLSGetBackend(slv))->lp();
 }
diff --git a/solvers/mosek/mosekbackend.h b/solvers/mosek/mosekbackend.h
index 3e33840b1..8f8a22a14 100644
--- a/solvers/mosek/mosekbackend.h
+++ b/solvers/mosek/mosekbackend.h
@@ -21,12 +21,23 @@ class MosekBackend :
   MosekBackend();
   ~MosekBackend();
 
+  /*----------------------------------------------------*\
+  | Standard and optional methods to provide or retrieve |
+  | information to/from or manipulate the solver. Most   |
+  | of them override placeholders from base classes.     |
+  \*----------------------------------------------------*/
+
   /// Name displayed in messages
   static const char* GetSolverName() { return "MOSEK"; }
   std::string GetSolverVersion();
-  
+
+  /// AMPL solver name is used to parse solver options
+  /// for the [name]_options environment variable.
+  /// This is only done if the [executable_name]_options
+  /// variable is not provided.
   static const char* GetAMPLSolverName() { return "mosek"; }
   static const char* GetAMPLSolverLongName() { return "AMPL-MOSEK"; }
+
   static const char* GetBackendName();
   static const char* GetBackendLongName() { return nullptr; }
 
@@ -36,7 +47,6 @@ class MosekBackend :
   void FinishOptionParsing() override;
 
 
-
   ////////////////////////////////////////////////////////////
   /////////////// OPTIONAL STANDARD FEATURES /////////////////
   ////////////////////////////////////////////////////////////
@@ -44,28 +54,41 @@ class MosekBackend :
   // that may or may not need additional functions. 
   USING_STD_FEATURES;
 
+  // LP basis info, status keys
   ALLOW_STD_FEATURE(BASIS, true)
-  // TODO If getting/setting a basis is supported, implement the 
-  // accessor and the setter below
   SolutionBasis GetBasis() override;
   void SetBasis(SolutionBasis) override;
 
- /**
-  * Get MIP Gap
-  **/
-  // TODO Implement to return MIP gap
+  // LP hotstart
+  // Set primal/dual initial guesses for continuous case
+  ALLOW_STD_FEATURE(WARMSTART, true)
+  void AddPrimalDualStart(Solution) override;
+
+  // MIP hotstart
+  ALLOW_STD_FEATURE( MIPSTART, true )
+  void AddMIPStart(ArrayRef<double> x0) override;
+
+  // Obtain inf/unbounded rays
+  ALLOW_STD_FEATURE(RAYS, true)
+  ArrayRef<double> Ray() override;
+  ArrayRef<double> DRay() override;
+
+  // MIP gap
   // (adds option mip:return_gap)
   ALLOW_STD_FEATURE(RETURN_MIP_GAP, true)
   double MIPGap() override;
   double MIPGapAbs() override;
-  /**
-  * Get MIP dual bound
-  **/
-  // TODO Implement to return the best dual bound value
+
+  // MIP dual bound
   // (adds option mip:bestbound)
   ALLOW_STD_FEATURE(RETURN_BEST_DUAL_BOUND, true)
   double BestDualBound() override;
 
+  // Sensitivity analysis
+  // Report sensitivity analysis suffixes (postsolved)
+  ALLOW_STD_FEATURE(SENSITIVITY_ANALYSIS, true)
+  SensRanges GetSensRanges() override;
+
   /////////////////////////// Model attributes /////////////////////////
   bool IsMIP() const override;
   bool IsQCP() const override;
@@ -82,7 +105,9 @@ class MosekBackend :
 
   Solution GetSolution() override;
   ArrayRef<double> GetObjectiveValues() override
-  { return std::vector<double>{ObjectiveValue()}; } 
+  {
+    return std::vector<double>{ObjectiveValue()};
+  }
 
 
   //////////////////// [[ Implementation details ]] //////////////////////
@@ -126,8 +151,8 @@ class MosekBackend :
   void VarStatii(ArrayRef<int>);
   void ConStatii(ArrayRef<int>);
 
-  ArrayRef<int> VarsIIS();
-  pre::ValueMapInt ConsIIS();
+  //ArrayRef<int> VarsIIS();
+  //pre::ValueMapInt ConsIIS();
 
   // utility function used to decide which solution to fetch after optimization
   // sets member solToFetch_ as side effect
@@ -141,9 +166,15 @@ class MosekBackend :
   /// These options are stored in the class
   struct Options {
     std::string exportFile_;
+
+    // Whether to set MSK_IPAR_MIO_CONSTRUCT_SOL
+    int MIPConstructSol_=0;
   };
   Options storedOptions_;
 
+protected:
+  /**** Option accessors ****/
+  int Mosek_mip_construct_sol() const { return storedOptions_.MIPConstructSol_; }
 
 };
 
diff --git a/solvers/mosek/mosekmodelapi.cc b/solvers/mosek/mosekmodelapi.cc
index bbd6e6421..bd9e82216 100644
--- a/solvers/mosek/mosekmodelapi.cc
+++ b/solvers/mosek/mosekmodelapi.cc
@@ -9,8 +9,7 @@ namespace mp {
 /// Defining the function in ...modelapi.cc
 /// for recompilation speed
 std::unique_ptr<BasicModelManager>
-CreateMosekModelMgr(MosekCommon& cc, Env& e,
-                     pre::BasicValuePresolver*& pPre) {
+CreateMosekModelMgr(MosekCommon& cc, Env& e, pre::BasicValuePresolver*& pPre) {
   return CreateModelMgrWithFlatConverter<
       MosekModelAPI, MIPFlatConverter >(cc, e, pPre);
 }
@@ -82,16 +81,11 @@ void MosekModelAPI::SetQuadraticObjective(int iobj, const QuadraticObjective& qo
     throw std::runtime_error("Multiple quadratic objectives not supported");
   }
 }
+
 void MosekModelAPI::AddLinearConstraint(
     MSKtask_t lp, size_t size, MSKboundkey_enum key,
     double lb, double ub,
     const int* vindex, const double* values, const char* name) {
-
-  if (lb == -std::numeric_limits<double>::infinity())
-    lb = -MSK_DPAR_DATA_TOL_BOUND_INF;
-  if (ub == std::numeric_limits<double>::infinity())
-    ub = MSK_DPAR_DATA_TOL_BOUND_INF;
-
   /* Linear + quadratic constraints are preallocated in
      InitProblemModificationPhase() */
   MOSEK_CCALL(MSK_putarow(lp, n_alg_cons_, size, vindex, values));
@@ -101,21 +95,45 @@ void MosekModelAPI::AddLinearConstraint(
   ++n_alg_cons_;
 }
 
-void MosekModelAPI::AddConstraint(const LinConRange& lc) {
-  AddLinearConstraint(lp(), lc.size(), MSK_BK_RA, lc.lb(), lc.ub(),
-    lc.pvars(), lc.pcoefs(), lc.name());
+// Some linear constraints are processed as Range, so we set the bound keys here.
+// Sometimes the separate AddConstraint functions are also called (e.g. logical constraints), see below.
+void MosekModelAPI::AddConstraint(const LinConRange& lc)
+{
+  double lb = lc.lb();
+  double ub = lc.ub();
+  MSKboundkey_enum key;
+
+  if (lb == -std::numeric_limits<double>::infinity())
+  {
+    if (ub == std::numeric_limits<double>::infinity())
+      key = MSK_BK_FR;
+    else
+      key = MSK_BK_UP;
+  }
+  else
+  {
+    if (ub == std::numeric_limits<double>::infinity())
+      key = MSK_BK_LO;
+    else if (lb == ub)
+      key = MSK_BK_FX;
+    else
+      key = MSK_BK_RA;
+  }
+
+  AddLinearConstraint(lp(), lc.size(), key, lb, ub, lc.pvars(), lc.pcoefs(), lc.name());
 }
-void MosekModelAPI::AddConstraint(const LinConLE& lc) {
-  AddLinearConstraint(lp(), lc.size(), MSK_BK_UP, lc.lb(), lc.ub(),
-    lc.pvars(), lc.pcoefs(), lc.name());
+
+void MosekModelAPI::AddConstraint(const LinConLE& lc)
+{
+  AddLinearConstraint(lp(), lc.size(), MSK_BK_UP, lc.lb(), lc.ub(), lc.pvars(), lc.pcoefs(), lc.name());
 }
-void MosekModelAPI::AddConstraint(const LinConEQ& lc) {
-  AddLinearConstraint(lp(), lc.size(), MSK_BK_FX, lc.lb(), lc.ub(),
-    lc.pvars(), lc.pcoefs(), lc.name());
+void MosekModelAPI::AddConstraint(const LinConEQ& lc)
+{
+  AddLinearConstraint(lp(), lc.size(), MSK_BK_FX, lc.lb(), lc.ub(), lc.pvars(), lc.pcoefs(), lc.name());
 }
-void MosekModelAPI::AddConstraint(const LinConGE& lc) {
-  AddLinearConstraint(lp(), lc.size(), MSK_BK_LO, lc.lb(), lc.ub(),
-    lc.pvars(), lc.pcoefs(), lc.name());
+void MosekModelAPI::AddConstraint(const LinConGE& lc)
+{
+  AddLinearConstraint(lp(), lc.size(), MSK_BK_LO, lc.lb(), lc.ub(), lc.pvars(), lc.pcoefs(), lc.name());
 }
 
 void MosekModelAPI::AddConstraint(const IndicatorConstraintLinLE &ic)  {
diff --git a/test/end2end/cases/categorized/fast/suf_common/modellist.json b/test/end2end/cases/categorized/fast/suf_common/modellist.json
index 4f8191b26..74ffcccef 100644
--- a/test/end2end/cases/categorized/fast/suf_common/modellist.json
+++ b/test/end2end/cases/categorized/fast/suf_common/modellist.json
@@ -38,8 +38,6 @@
       "Buy['MCH'].senslblo": -7.274725274725274,
       "Buy['CHK'].sensubhi": 10.97922077922078,
       "Buy['MTL'].sensublo": 6.235955056179775,
-      "Diet_Max['A'].sensrhslo": 0.0,
-      "Diet_Max['A'].sensrhshi": 0.0,
       "Diet_Max['A'].senslbhi": 0.0,
       "if Diet_Max['CAL'].senslblo<=-1e20 then 1": 1,
       "Diet_Min['CAL'].senslblo": 15246.0,
@@ -57,12 +55,8 @@
     "tags" : ["linear", "continuous", "sens"],
     "objective" : 1128.642,
     "values": {
-      "Diet['A'].sensrhslo": 7850,
-      "Diet['A'].sensrhshi": 11600,
-      "Diet['C'].sensrhslo": 9466.66666666666666,
-      "Diet['C'].sensrhshi": 11500,
-      "if Diet['A'].senslblo<=-1e20 then 1": 1,
-      "Diet['A'].senslbhi": 10000,
+      "Diet['A'].senslblo": -1450,
+      "Diet['A'].senslbhi": 2300,
       "Diet['A'].sensublo": 7850,
       "Diet['A'].sensubhi": 11600,
       "Diet['C'].sensublo": 9466.66666666666666,
@@ -81,10 +75,6 @@
     "objective" : 1128.642,
     "comment": "Diet_Min{} constraints are converted to <= by AMPL",
     "values": {
-      "Diet_Max['A'].sensrhslo": 7850,
-      "Diet_Max['A'].sensrhshi": 11600,
-      "Diet_Max['C'].sensrhslo": 9466.66666666666666,
-      "Diet_Max['C'].sensrhshi": 11500,
       "if Diet_Max['A'].senslblo<=-1e20 then 1": 1,
       "if Diet_Max['A'].senslbhi>=1e20 then 1": 1,
       "Diet_Max['A'].sensublo": 7850,
diff --git a/test/end2end/scripts/python/Solver.py b/test/end2end/scripts/python/Solver.py
index 8c220f8ee..9602af9c6 100644
--- a/test/end2end/scripts/python/Solver.py
+++ b/test/end2end/scripts/python/Solver.py
@@ -639,7 +639,6 @@ def __init__(self, exeName, timeout=None, nthreads=None,
                  ModelTags.nonlinear,
                  ModelTags.log,
                  # ModelTags.trigonometric
-
                  ModelTags.relax,
                  ModelTags.multiobj
                  }
@@ -708,7 +707,8 @@ def _getAMPLOptionsName(self):
         return "mosek"
 
     def __init__(self, exeName, timeout=None, nthreads=None, otherOptions=None):
-        stags = {ModelTags.continuous, ModelTags.integer, ModelTags.binary,
-                 ModelTags.quadratic
-                 }
+        stags = {ModelTags.continuous, ModelTags.linear, ModelTags.integer,
+                 ModelTags.binary, ModelTags.warmstart, ModelTags.mipstart,
+                 ModelTags.return_mipgap, ModelTags.sens, ModelTags.sstatus}
+        # ModelTags.quadratic
         super().__init__(exeName, timeout, nthreads, otherOptions, stags)