#include <iostream>
#include <Eigen/Dense>
#include <vector>
#include <algorithm>
#include <set>
#include <numeric>
#include <list>
using namespace std;
using Eigen::MatrixXd;
template< typename T >
std::vector<size_t> sort_indexes(const std::vector<T>& v) {
std::vector<size_t> indx(v.size());
iota(indx.begin(), indx.end(), 0);
stable_sort(indx.begin(), indx.end(), [&v](size_t ia, size_t ib) {
return v[ia] > v[ib];
});
return indx;
}
bool is_exist(size_t col, const std::set<size_t>& cols, size_t N) {
if (cols.find(col) != cols.end()) return true;
return false;
}
bool find_new(
const Eigen::MatrixXd& input,
size_t row,
std::set<size_t> usedCols,
std::list<std::pair<size_t, size_t>>& result)
{
// fmt::print("row: {}, new: {}\n", row, usedCols);
if (row == input.rows() && usedCols.size() == input.cols()) return true;
for (size_t col = 0; col < input.cols(); ++col) {
if (usedCols.count(col)) continue;
if (input(row, col) == 0) {
auto status = usedCols.insert(col);
if (find_new(input, row + 1, usedCols, result)) {
// fmt::print("get it: {}, {}; {}\n", row, col, usedCols);
result.push_front({ row, col });
return true;
}
else {
usedCols.erase(status.first);
// fmt::print("try: {}, {}\n", row, col);
}
}
}
return false;
}
std::list<std::pair<size_t, size_t>> get_valid_indexes(const Eigen::MatrixXd& input) {
std::set<size_t> usedCols;
std::list<std::pair<size_t, size_t>> result;
find_new(input, 0, usedCols, result);
return result;
}
int solve(const Eigen::MatrixXd& input, std::list<std::pair<size_t, size_t>>& matched) {
// step 0: weight must positive
Eigen::MatrixXd m(input);
const auto minRow = m.rowwise().minCoeff();
Eigen::MatrixXd m1 = m.colwise() - minRow;
const auto minCol = m1.colwise().minCoeff();
Eigen::MatrixXd result = m1.rowwise() - minCol;
do {
// get zero count of each row and col, then sort it by count
size_t total = std::count(result.data(), result.data() + result.size(), 0);
// fmt::print("step 2:\ntotal {}\n", total);
size_t times = 0;
Eigen::MatrixXd mTempInf = result;
Eigen::MatrixXd mTempLine = Eigen::MatrixXd::Zero(result.rows(), result.cols());
Eigen::MatrixXd mTempLineCross = Eigen::MatrixXd::Zero(result.rows(), result.cols());
do {
size_t ridx = 0, cidx = 0;
std::vector<size_t> vRowCounts, vColCounts;
for (size_t row = 0; row < mTempInf.rows(); ++row) {
const auto& r = mTempInf.row(row);
size_t c = r.size() - r.count();
vRowCounts.push_back(c);
}
for (size_t col = 0; col < mTempInf.cols(); ++col) {
const auto& c = mTempInf.col(col);
size_t value = c.size() - c.count();
vColCounts.push_back(value);
}
// sort
std::vector<size_t> vRows = sort_indexes(vRowCounts);
std::vector<size_t> vCols = sort_indexes(vColCounts);
size_t rValue = vRowCounts[vRows[ridx]];
size_t cValue = vColCounts[vCols[cidx]];
if (rValue > cValue) {
auto rows = mTempInf.row(vRows[ridx]);
rows = (rows.array() == 0).select(INFINITY, rows);
mTempLine.row(vRows[ridx]) = result.row(vRows[ridx]);
mTempLineCross.row(vRows[ridx]) += Eigen::VectorXd::Ones(result.cols());
ridx += 1;
}
else {
auto cols = mTempInf.col(vCols[cidx]);
cols = (cols.array() == 0).select(INFINITY, cols);
mTempLine.col(vCols[cidx]) = result.col(vCols[cidx]);
mTempLineCross.col(vCols[cidx]) += Eigen::VectorXd::Ones(result.cols());
cidx += 1;
}
times += 1;
} while (mTempInf.count() < mTempInf.size());
if (times >= input.rows()) {
break;
}
double minimal = mTempInf.minCoeff();
double maximal = mTempLineCross.maxCoeff();
mTempInf -= Eigen::MatrixXd::Ones(mTempInf.rows(), mTempInf.cols()) * minimal;
mTempInf = (mTempInf.array() == INFINITY).select(0, mTempInf);
auto pluss = (mTempLineCross.array() == maximal).select(minimal, Eigen::MatrixXd::Zero(result.rows(), result.cols()));
result = mTempInf + mTempLine + pluss;
} while (true);
// get position of matched
auto indexes = get_valid_indexes(result);
auto itr = indexes.begin();
for (size_t offset = 0; offset < input.rows(); ++itr, ++offset);
matched.assign(indexes.begin(), itr);
return 1;
}
int main()
{
MatrixXd m(6, 6);
m << 62, 75, 80, 93, 95, 97,
75, 80, 82, 85, 71, 97,
80, 75, 81, 98, 90, 97,
78, 82, 84, 80, 50, 98,
90, 85, 85, 80, 85, 99,
65, 75, 80, 75, 68, 96;
cout << m << endl;
std::list<std::pair<size_t, size_t>> indexes;
int ret = solve(m, indexes);
int weight;
if (ret == 1) {
weight = 0;
for (auto indx = indexes.begin(); indx != indexes.end(); ++indx) {
weight += m(indx->first, indx->second);
cout << indx->first << "-->" << indx->second << endl;
}
}
cout << weight << endl;
return ret;
}