//
//    Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
//    File:    midual1.cc
//
//    Purpose:  Demonstrates how to compute dual values with
//              respect to a fixed integer solution for a MIO problem.
//
//              minimize    500 s1 + 300 s2 + 10 x1 + 14 x2
//              subject to  x1 + x2 >= 100
//                          0 <= x1 <= 70 s1
//                          0 <= x2 <= 80 s2
//                          s1, s2 - binary
//
#include <iostream>
#include <iomanip>
#include "fusion.h"

using namespace mosek::fusion;
using namespace monty;

int main(int argc, char ** argv)
{
  /* M is the initial mixed-integer model */
  Model::t M = new Model("midual1"); auto _M = finally([&]() { M->dispose(); });

  auto x = M->variable("x", 2, Domain::greaterThan(0));
  auto s = M->variable("s", 2, Domain::binary());

  M->objective(ObjectiveSense::Minimize,
               Expr::add(Expr::dot(new_array_ptr<double, 1>({10,14}), x),
                         Expr::dot(new_array_ptr<double, 1>({500,300}), s)));

  M->constraint("demand", Expr::sum(x), Domain::greaterThan(100));
  M->constraint("production",
                Expr::sub(x, Expr::mulElm(new_array_ptr<double, 1>({70,80}), s)),
                Domain::lessThan(0));

  M->solve();

  if (M->getProblemStatus() != ProblemStatus::PrimalFeasible)
    return -1;  // Unsuitable problem status, exiting

  std::cout << std::setprecision(2)
            << "x = " << (*(x->level()))[0] << ", " << (*(x->level()))[1] << std::endl ;

  /* F is the continuous fixed model */
  Model::t F = M->getFixedModel(); auto _F = finally([&]() { F->dispose(); });

  F->solve();

  if (F->getProblemStatus() != ProblemStatus::PrimalAndDualFeasible)
    return -1;  // Unsuitable problem status, exiting

  auto demand = F->getConstraint("demand");
  auto production = F->getConstraint("production");
  auto xfix = F->getVariable("x");

  std::cout << std::setprecision(2)
            << "xfix = " << (*(xfix->level()))[0] << ", " << (*(xfix->level()))[1] << std::endl
            << "demand dual = " << (*(demand->dual()))[0] << std::endl
            << "production dual = " << (*(production->dual()))[0] << ", " << (*(production->dual()))[1]
            << std::endl;
  return 0;
}