Calculator.cpp 27 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920
  1. #include "Calculator.h"
  2. #include "Container.h"
  3. #include <iostream>
  4. #include <fstream>
  5. #include <string>
  6. #include <sstream>
  7. #include <list>
  8. #include <tuple>
  9. #include <utility>
  10. #include <exception>
  11. #include <algorithm>
  12. #include <numeric>
  13. #include <vector>
  14. #include <cassert>
  15. #include <cstring>
  16. #include <cstdlib>
  17. #include <exception>
  18. #include "ActiveArty.h"
  19. using namespace std;
  20. Calculator::Calculator() : numBat(0), numMis(0), numTarg(0)
  21. {
  22. Container A;
  23. ifstream fileArty("arty.txt");
  24. try
  25. {
  26. if (fileArty.is_open())
  27. {
  28. string input;
  29. while (getline(fileArty, input))
  30. {
  31. stringstream sinput(input);
  32. char delim = ',';
  33. char delim2 = ';';
  34. string xcoord, ycoord, typearty, numinbattery, numshells, typeshells, typefuze, temperature, maxTime;
  35. getline(sinput, xcoord, delim);
  36. getline(sinput, ycoord, delim);
  37. getline(sinput, typearty, delim);
  38. getline(sinput, numinbattery, delim);
  39. getline(sinput, numshells, delim);
  40. getline(sinput, typeshells, delim);
  41. getline(sinput, typefuze, delim);
  42. getline(sinput, temperature, delim);
  43. getline(sinput, maxTime, delim2);
  44. numBat++;
  45. artyInfo.push_back({ xcoord, ycoord, typearty,numinbattery, numshells, typeshells, typefuze, temperature, maxTime });
  46. A.checkCorrectArty(typearty, static_cast<size_t>(stoi(numinbattery)), static_cast<size_t>(stoi(numshells)),
  47. typeshells, typefuze, stod(temperature));
  48. }
  49. }
  50. else
  51. {
  52. throw runtime_error("Problem while opening file to read");
  53. }
  54. fileArty.close();
  55. }
  56. catch (exception e)
  57. {
  58. cerr << e.what() << endl;
  59. exit(1);
  60. }
  61. ifstream fileTarg("targets.txt");
  62. try
  63. {
  64. if (fileTarg.is_open())
  65. {
  66. string input;
  67. while (getline(fileTarg, input))
  68. {
  69. stringstream sinput(input);
  70. char delim = ',';
  71. char delim2 = ';';
  72. string xcoord, ycoord, typetarg, mission,entrenchedArmoured, geometry, preparing;
  73. getline(sinput, xcoord, delim);
  74. getline(sinput, ycoord, delim);
  75. getline(sinput, typetarg, delim);
  76. getline(sinput, entrenchedArmoured, delim);
  77. getline(sinput, mission, delim);
  78. getline(sinput, geometry, delim);
  79. getline(sinput, preparing, delim2);
  80. numTarg++;
  81. targInfo.push_back({ xcoord, ycoord, typetarg, mission, entrenchedArmoured, geometry, preparing });
  82. A.checkCorrectTarget(typetarg, mission, entrenchedArmoured, geometry, preparing);
  83. }
  84. }
  85. else
  86. {
  87. throw runtime_error("Problem while opening file to read");
  88. }
  89. fileTarg.close();
  90. }
  91. catch (exception e)
  92. {
  93. cerr << e.what() << endl;
  94. exit(1);
  95. }
  96. }
  97. Calculator::~Calculator() {
  98. AMtx->clear();
  99. delete AMtx;
  100. BVec->clear();
  101. delete BVec;
  102. CVec->clear();
  103. delete CVec;
  104. }
  105. size_t Calculator::countRangesAndAmount(size_t whatTarget)
  106. {
  107. vector<ActiveArty> units;
  108. double xTarg = stod(targInfo[whatTarget][0]);
  109. double yTarg = stod(targInfo[whatTarget][1]);
  110. size_t numVariants = 0;
  111. Container A;
  112. size_t actualCharge = 0;
  113. for_each(artyInfo.begin(), artyInfo.end(), [&](auto& _n)
  114. {
  115. double xArty = stod(_n[0]);
  116. double yArty = stod(_n[1]);
  117. string name = _n[2];
  118. size_t numInBattery = static_cast<size_t>(stoi(_n[3]));
  119. size_t numShells = static_cast<size_t>(stoi(_n[4]));
  120. double range = pow(pow(xArty - xTarg, 2) + pow(yArty - yTarg, 2), 0.5);
  121. actualCharge = 1;
  122. //A.getCharge(name, range);
  123. if (range > 0 && range < A.getRange(name) && numShells > 0) {
  124. units.push_back({ range, stod(_n[7]), A.getType(name),A.getCaliber(name),
  125. numInBattery, numShells, _n[5], _n[6], name}); //ñîçäà¸ì ñïèñîê àêòèâíîé íà ýòó öåëü àðòèëëåðèè
  126. }
  127. });
  128. //range, temperatire, <howitzer, caliber>, <numinbatt, numshellscommon>, <typeshells, typefuze>
  129. string targName = targInfo[whatTarget][2];
  130. size_t move = calculatedNum.size();
  131. size_t j = 0;
  132. vector<pair<size_t, size_t>> filling = vector<pair<size_t, size_t>>(numBat, make_pair(UINT_MAX, UINT_MAX));
  133. calculatedNum.insert(calculatedNum.end(), filling.begin(), filling.end()); //ïî äåôîëòó àêòèâíàÿ àðòèëëåðèÿ - âñÿ, íî äëÿ íåàêòèâíîé ìàêñèìàëüíûå ÷èñëà â ïàðàõ
  134. switch (A.checkStationary(targName))
  135. {
  136. case (true):
  137. for_each(units.begin(), units.end(), [&](auto& _n)
  138. {
  139. bool isCluster = (_n.typeshells == "cluster") ? true : false;
  140. bool isShortened = (targInfo[whatTarget][5] == "shortened") ? true : false;
  141. bool isEntrenchedOrArmored = (targInfo[whatTarget][3] == "entrenched" || targInfo[whatTarget][3] == "armoured") ? true : false;
  142. bool isDestruction = (targInfo[whatTarget][4] == "destruction") ? true : false;
  143. double temperature = _n.temperature;
  144. size_t actualNumShells = A.getNormStationary(_n.typesCalibers, isCluster, false,
  145. isShortened, isEntrenchedOrArmored, isDestruction, false, _n.range, _n.numInBattery, targName);
  146. size_t actualMinutes = A.measureFireRate(_n.typesCalibers, actualCharge, temperature, actualNumShells, A.getSpecificType(_n.name));
  147. calculatedNum[move + j] = make_pair(actualNumShells, actualMinutes);
  148. numVariants++;
  149. j++;
  150. });
  151. break;
  152. case (false):
  153. for_each(units.begin(), units.end(), [&](auto& _n)
  154. {
  155. bool isCluster = (_n.typeshells == "cluster") ? true : false;
  156. bool isEntrenchedOrArmored = (targInfo[whatTarget][3] == "entrenched" || targInfo[whatTarget][3] == "armoured") ? true : false;
  157. double temperature = _n.temperature;
  158. size_t actualNumShells = A.getNormMoving(_n.typesCalibers, isCluster, isEntrenchedOrArmored, _n.numInBattery, targName);
  159. size_t actualMinutes = A.measureFireRate(_n.typesCalibers, actualCharge, temperature, actualNumShells, A.getSpecificType(_n.name));
  160. calculatedNum[move + j] = make_pair(actualNumShells, actualMinutes);
  161. numVariants++;
  162. j++;
  163. });
  164. }
  165. units.clear();
  166. return numVariants;
  167. }
  168. void Calculator::simplexInit(double a, double b)
  169. {
  170. vector<size_t> numVar;
  171. for (size_t i = 0; i < numTarg; i++)
  172. {
  173. numVar.push_back(countRangesAndAmount(i)); //ðàññ÷èòûâàåì âàðèàíòû íà êàæäóþ öåëü è äîáàâëÿåì èõ ÷èñëî ê âåêòîðó èíäåêñîâ
  174. }
  175. r = 2 * numBat + numTarg; //÷èñëî áàçèñíûõ ïåðåìåííûõ
  176. k = numTarg * numBat;
  177. //r - numTarg; //÷èñëî ñâîáîäíûõ ïåðåìåííûõ
  178. size_t i = 0;
  179. AMtx = new vector<vector<double>>(r, vector<double>(k, 0));
  180. auto it1 = AMtx->begin();
  181. size_t move = 0;
  182. advance(it1, numBat);
  183. for_each(AMtx->begin(), it1, [&](auto& _n)
  184. {
  185. i = 0;
  186. for_each(_n.begin(), _n.end(), [&](auto& _m)
  187. {
  188. if ((move <= i) && (i < move + numTarg))
  189. {
  190. _m = calculatedNum[i].first;
  191. }
  192. else
  193. {
  194. _m = 0;
  195. }
  196. i++;
  197. });
  198. move += numTarg;
  199. });
  200. auto it2 = it1;
  201. advance(it2, numBat);
  202. move = 0;
  203. for_each(it1, it2, [&](auto& _n)
  204. {
  205. i = 0;
  206. for_each(_n.begin(), _n.end(), [&](auto& _m)
  207. {
  208. if ((move <= i) && (i < move + numTarg))
  209. {
  210. _m = calculatedNum[i].second;
  211. }
  212. else
  213. {
  214. _m = 0;
  215. }
  216. i++;
  217. });
  218. move += numTarg;
  219. });
  220. i = 0;
  221. for_each(it2, AMtx->end(), [&](auto& _n)
  222. {
  223. size_t j = 0;
  224. auto it = _n.begin();
  225. advance(it, numTarg * numBat);
  226. for_each(_n.begin(), it, [&](auto& _m)
  227. {
  228. if ((j) % numTarg - i == 0)
  229. {
  230. _m = 1;
  231. }
  232. else
  233. {
  234. _m = 0;
  235. }
  236. j++;
  237. });
  238. i++;
  239. });
  240. // for (i = 0; i < 2 * numBat; i++)
  241. // {
  242. // (*AMtx)[i][i + numBat * numTarg] = 1.0; //äîáàâëÿåì áàçèñíûå ñëàãàåìûå
  243. // }
  244. i = 0;
  245. BVec = new vector<double>(r);
  246. auto it3 = BVec->begin();
  247. advance(it3, numBat);
  248. for_each(BVec->begin(), it3, [&](auto& _n)
  249. {
  250. _n = static_cast<double>(stof(artyInfo[i][4])); //÷èñëî ñíàðÿäîâ íà áàòàðåå
  251. i++;
  252. });
  253. i = 0;
  254. auto it4 = it3;
  255. advance(it4, numBat);
  256. for_each(it3, it4, [&](auto& _n)
  257. {
  258. _n = static_cast<long long int>(stof(artyInfo[i][8])); //âðåìÿ äî ïåðåäèñëîêàöèè áàòàðåè
  259. i++;
  260. });
  261. for_each(it4, BVec->end(), [&](auto& _n)
  262. {
  263. _n = 1; //ñêîëüêî áàòàðåé äîëæíî áûòü ïîäêëþ÷åíî ê ïîðàæåíèþ
  264. i++;
  265. });
  266. i = 0;
  267. CVec = new vector<double>(k);
  268. auto it5 = CVec->begin();
  269. advance(it5, numBat*numTarg);
  270. for_each(CVec->begin(), it5, [&](auto& _n)
  271. {
  272. _n = a * calculatedNum[i].first + b * calculatedNum[i].second;
  273. i++;
  274. });
  275. for_each(it5, CVec->end(), [&](auto& _n)
  276. {
  277. _n = 0;
  278. });
  279. val = 0;
  280. simplex(AMtx, BVec, CVec, *CVec, 0, k, 1, true, a, b);
  281. }
  282. void Calculator::output(vector<vector<double> >* a, vector<double>* b, vector<double>* c, vector<double> x, double opt, int numVars)
  283. {
  284. // open out.txt
  285. ofstream outfile;
  286. outfile.open("out.txt");
  287. // make sure it's open
  288. try
  289. {
  290. if (outfile.is_open())
  291. {
  292. // output the final simplex tableau seen
  293. outfile << "Final Tableau:" << endl;
  294. size_t i = 0;
  295. for_each(a->begin(), a->end(), [&](auto const& _n)
  296. {
  297. for_each(_n.begin(), _n.end(), [&](auto const& _m)
  298. {
  299. // output the A matrix
  300. outfile << _m << " ";
  301. });
  302. // output the b matrix
  303. outfile << "| " << (*b)[i] << endl;
  304. i++;
  305. });
  306. // output the c matrix
  307. for_each(c->begin(), c->end(), [&](auto const& _n)
  308. {
  309. outfile << _n << " ";
  310. });
  311. // output the optimal value found
  312. outfile << "| " << opt << endl;
  313. outfile << endl;
  314. outfile << "z* = " << opt << endl;
  315. // output the optimal x
  316. outfile << "x* = (";
  317. for (int i = 0; i < x.size(); i++)
  318. {
  319. if (i == x.size() - 1)
  320. outfile << x[i] << ")" << endl;
  321. else
  322. outfile << x[i] << ", ";
  323. }
  324. // close out.txt
  325. outfile.close();
  326. }
  327. // error
  328. else
  329. throw runtime_error("Could not open output file.");
  330. }
  331. catch (exception e)
  332. {
  333. cerr << e.what() << endl;
  334. exit(1);
  335. }
  336. }
  337. //! Outputs an unbounded solution to "out.txt"
  338. /*!
  339. \param a the A matrix
  340. \param b the b matrix
  341. \param c the c matrix
  342. \param minimize true if min, false if max
  343. \sa output(), outputInfeasible()
  344. */
  345. void Calculator::outputUnbounded(vector<vector<double> >* a, vector<double>* b, vector<double>* c, bool minimize)
  346. {
  347. // open out.txt
  348. ofstream outfile;
  349. outfile.open("out.txt");
  350. // make sure it's open
  351. try
  352. {
  353. if (outfile.is_open())
  354. {
  355. // output the final simplex tableau seen
  356. outfile << "Final Tableau:" << endl;
  357. size_t i = 0;
  358. for_each(a->begin(), a->end(), [&](auto const& _n)
  359. {
  360. for_each(_n.begin(), _n.end(), [&](auto const& _m)
  361. {
  362. // output the A matrix
  363. outfile << _m << " ";
  364. });
  365. // output the b matrix
  366. outfile << "| " << (*b)[i] << endl;
  367. i++;
  368. });
  369. // output the c matrix
  370. for_each(c->begin(), c->end(), [&](auto const& _n)
  371. {
  372. outfile << _n << " ";
  373. });
  374. outfile << endl;
  375. // optimal is unbounded
  376. if (minimize)
  377. outfile << "z* = -infinity" << endl;
  378. else
  379. outfile << "z* = infinity" << endl;
  380. outfile << "The program is unbounded." << endl;
  381. // close out.txt
  382. outfile.close();
  383. }
  384. else throw runtime_error("Could not open output file.");
  385. }
  386. catch (exception e)
  387. {
  388. cerr << e.what() << endl;
  389. exit(1);
  390. }
  391. }
  392. //! Outputs an infeasible solution to "out.txt"
  393. /*!
  394. \param a the A matrix
  395. \param b the b matrix
  396. \param c the c matrix
  397. \sa output(), outputUnbounded()
  398. */
  399. void Calculator::outputInfeasible(vector<vector<double> >* a, vector<double>* b, vector<double>* c)
  400. {
  401. // open out.txt
  402. ofstream outfile;
  403. outfile.open("out.txt");
  404. try
  405. {
  406. // make sure it's open
  407. if (outfile.is_open())
  408. {
  409. // output the final simplex tableau seen
  410. outfile << "Final Tableau:" << endl;
  411. size_t i = 0;
  412. for_each(a->begin(), a->end(), [&](auto const& _n)
  413. {
  414. for_each(_n.begin(), _n.end(), [&](auto const& _m)
  415. {
  416. // output the A matrix
  417. outfile << _m << " ";
  418. });
  419. // output the b matrix
  420. outfile << "| " << (*b)[i] << endl;
  421. i++;
  422. });
  423. // output the c matrix
  424. for_each(c->begin(), c->end(), [&](auto const& _n)
  425. {
  426. outfile << _n << " ";
  427. });
  428. outfile << endl;
  429. // problem is infeasible
  430. outfile << "The program is infeasible." << endl;
  431. // close out.txt
  432. outfile.close();
  433. }
  434. // error
  435. else throw runtime_error("Could not open output file.");
  436. }
  437. catch (exception e)
  438. {
  439. cerr << e.what() << endl;
  440. exit(1);
  441. }
  442. }
  443. void Calculator::reduce(vector<vector<double> >* a, vector<double>* b, vector<double>* c, double* opt, int row, int col)
  444. {
  445. // get the value needed to make the pivot element 1
  446. double pivdiv = 1 / (*a)[row][col];
  447. // reduce the pivot row with this value
  448. for (int i = 0; i < (*a)[0].size(); i++)
  449. {
  450. (*a)[row][i] = (*a)[row][i] * pivdiv;
  451. }
  452. (*b)[row] = (*b)[row] * pivdiv;
  453. // for every other row
  454. for (int i = 0; i < (*a).size(); i++)
  455. {
  456. if (i == row)
  457. continue;
  458. // get the value needed to make the pivot column element 0
  459. double coeff = (*a)[i][col] * (*a)[row][col];
  460. // reduce the row with this value
  461. for (int j = 0; j < (*a)[0].size(); j++)
  462. {
  463. (*a)[i][j] = (*a)[i][j] - ((*a)[row][j] * coeff);
  464. }
  465. (*b)[i] = (*b)[i] - ((*b)[row] * coeff);
  466. }
  467. // get the value needed to make the c element 0 in the pivot column
  468. double coeff = (*c)[col] * (*a)[row][col];
  469. // reduce the c row with this value
  470. for (int i = 0; i < (*c).size(); i++)
  471. {
  472. (*c)[i] = (*c)[i] - ((*a)[row][i] * coeff);
  473. }
  474. (*opt) = (*opt) - ((*b)[row] * coeff);
  475. }
  476. //! Performs row operations to make c elements in basic columns 0
  477. /*!
  478. \param a the A matrix
  479. \param b the b matrix
  480. \param c the c matrix
  481. \param opt the optimal value in the current tableau
  482. \sa reduce()
  483. */
  484. void Calculator::reduceC(vector<vector<double> >* a, vector<double>* b, vector<double>* c, double* opt)
  485. {
  486. // for each row
  487. for (int i = 0; i < (*a).size(); i++)
  488. {
  489. // default the row operation coefficient to 1
  490. double coeff = 1;
  491. for (int j = 0; j < (*a)[0].size(); j++)
  492. {
  493. // found a 1, check for a basic column
  494. if ((*a)[i][j] == 1)
  495. {
  496. bool basic = true;
  497. for (int k = 0; k < (*a).size(); k++)
  498. {
  499. // not 0 or 1, not basic
  500. if ((*a)[k][j] != 0 && k != i)
  501. basic = false;
  502. }
  503. // basic, make sure the c element becomes 0
  504. if (basic)
  505. coeff = (*c)[j] * (*a)[i][j];
  506. }
  507. // reduce the c element in this column
  508. (*c)[j] = (*c)[j] - ((*a)[i][j] * coeff);
  509. }
  510. // apply the reduction to the optimal value
  511. (*opt) = (*opt) - ((*b)[i] * coeff);
  512. }
  513. return;
  514. }
  515. //! Performs the Simplex algorithm on a given tableau
  516. /*!
  517. \param a the A matrix
  518. \param b the b matrix
  519. \param c the c matrix
  520. \param initialC the original c matrix (for calculating optimal)
  521. \param initialOpt the optimal value to begin the algorithm with
  522. \param numVars the number of variables in the program
  523. \param phase 1 or 2, corresponding to what phase of the Simplex Method we are on
  524. \param minimize true if min, false if max
  525. */
  526. void Calculator::simplex(vector<vector<double>>* a, vector<double>* b, vector<double>* c,
  527. vector<double> initialC, double initialOpt, int numVars, int phase, bool minimize, double _aa, double _bb)
  528. {
  529. // Phase 1
  530. if (phase == 1)
  531. {
  532. // make the phase 1 c row
  533. vector<double> p1c;
  534. // make original variables 0
  535. for (int i = 0; i < (*a)[0].size(); i++)
  536. {
  537. p1c.push_back(0.0);
  538. }
  539. // add as many 1s as rows
  540. for (int i = 0; i < (*a).size(); i++)
  541. {
  542. if (i < (*a).size() - numTarg) p1c.push_back(1.0);
  543. // add an identity matrix to A
  544. for (int j = 0; j < (*a).size() - numTarg; j++)
  545. {
  546. if (j == i)
  547. (*a)[i].push_back(1.0);
  548. else
  549. (*a)[i].push_back(0.0);
  550. }
  551. }
  552. double opt = 0.0;
  553. // reduce the new c row so basic variables have 0 cost
  554. reduceC(a, b, &p1c, &opt);
  555. // run the Simplex algorithm
  556. bool stop = false;
  557. while (!stop)
  558. {
  559. // get the pivot column
  560. double min_r = 0;
  561. int piv_col = 0;
  562. for (int i = 0; i < p1c.size(); i++)
  563. {
  564. // check for the smallest c < 0
  565. if (p1c[i] < min_r)
  566. {
  567. min_r = p1c[i];
  568. piv_col = i;
  569. }
  570. }
  571. // no c < 0, we're done
  572. if (min_r >= 0)
  573. {
  574. stop = true;
  575. continue;
  576. }
  577. // get the pivot row
  578. vector<double> ratios;
  579. // for each row
  580. for (int i = 0; i < (*a).size(); i++)
  581. {
  582. if ((*a)[i][piv_col] == 0)
  583. ratios.push_back(-1);
  584. else
  585. {
  586. // denegerate loop ahead, don't pivot here
  587. if ((*a)[i][piv_col] <= 0 && (*b)[i] == 0)
  588. ratios.push_back(-1);
  589. else
  590. ratios.push_back((*b)[i] / (*a)[i][piv_col]);
  591. }
  592. }
  593. double min_ratio = ratios[0];
  594. int piv_row = 0;
  595. for (int i = 1; i < ratios.size(); i++)
  596. {
  597. if (ratios[i] >= 0 && (ratios[i] < min_ratio || min_ratio < 0))
  598. {
  599. min_ratio = ratios[i];
  600. piv_row = i;
  601. }
  602. }
  603. // no positive ratios, stop here
  604. if (min_ratio < 0)
  605. {
  606. stop = true;
  607. continue;
  608. }
  609. // reduce the tableau on the pivot row and column
  610. reduce(a, b, &p1c, &opt, piv_row, piv_col);
  611. // return to step 1 of the Simplex algorithm
  612. long long int indFrac = -1;
  613. double maxFracPart = 0.0;
  614. bool isInteger = true;
  615. for (auto x_i : (*b))
  616. {
  617. if (x_i != floor(x_i))
  618. {
  619. isInteger = false;
  620. break;
  621. }
  622. }
  623. if (!isInteger)
  624. {
  625. for (size_t i = 0; i < (*b).size(); i++)
  626. {
  627. double fracPart = 0.;
  628. if ((*b)[i] > 0) fracPart = (*b)[i] - floor((*b)[i]);
  629. else fracPart = (*b)[i] - ceil((*b)[i]);
  630. if (fracPart > maxFracPart) {
  631. maxFracPart = fracPart;
  632. indFrac = i;
  633. }
  634. }
  635. r++;
  636. k++;
  637. numVars++;
  638. AMtx->push_back(vector<double>(k + 2 * numBat, 0.));
  639. for (size_t i = 0; i < r - 1; i++)
  640. {
  641. (*AMtx)[i].push_back(0.);
  642. }
  643. //åù¸ îãðàíè÷åíèå è ïåðåìåííàÿ
  644. for (size_t j = 0; j < k + 2 * numBat - 1; j++)
  645. {
  646. if ((*AMtx)[indFrac][j] > 0) (*AMtx)[r - 1][j] = (*AMtx)[indFrac][j] - floor((*AMtx)[indFrac][j]);
  647. else (*AMtx)[r - 1][j] = (*AMtx)[indFrac][j] - ceil((*AMtx)[indFrac][j]);
  648. }
  649. double bVal = 0;
  650. if ((*b)[indFrac] > 0) bVal = (*b)[indFrac] - floor((*b)[indFrac]);
  651. else bVal = (*b)[indFrac] - ceil((*b)[indFrac]);
  652. BVec->push_back(abs(bVal));
  653. if (bVal > 0)
  654. {
  655. (*AMtx)[r - 1][k + 2 * numBat - 1] = 1.;
  656. }
  657. else
  658. {
  659. (*AMtx)[r - 1][k + 2 * numBat - 1] = -1.;
  660. }
  661. CVec->push_back(0.);
  662. //gomoriInit(_aa, _bb);
  663. for (size_t j = 0; j < k; j++)
  664. {
  665. if (((*CVec)[j] > 0.) && ((*CVec)[j] < 1e-5)) (*CVec)[j] = 0.;
  666. }
  667. for (size_t i = 0; i < r; i++)
  668. {
  669. if (((*BVec)[i] > 0) && ((*BVec)[i] < 1e-5)) (*BVec)[i] = 0.;
  670. for (size_t j = 0; j < k + 2 * numBat; j++)
  671. {
  672. if (((*AMtx)[i][j] > 0) && ((*AMtx)[i][j] < 1e-5)) (*AMtx)[i][j] = 0.;
  673. }
  674. }
  675. }
  676. }
  677. // phase 1 done, cut off the added variables
  678. for (int i = 0; i < (*a).size(); i++)
  679. {
  680. (*a)[i].resize(numVars);
  681. }
  682. for (size_t j = initialC.size(); j < CVec->size(); j++)
  683. {
  684. initialC.push_back((*CVec)[j]);
  685. }
  686. // reset optimal value
  687. opt = 0;
  688. // reduce the new c row so basic variables have 0 cost
  689. reduceC(a, b, c, &opt);
  690. // go to phase 2
  691. simplex(a, b, c, initialC, opt, numVars, 2, minimize, _aa, _bb);
  692. }
  693. // Phase 2
  694. else
  695. {
  696. double opt = initialOpt;
  697. // run the Simplex algorithm
  698. bool stop = false;
  699. while (!stop)
  700. {
  701. // get the pivot column
  702. double min_r = 0;
  703. int piv_col = 0;
  704. for (int i = 0; i < (*c).size(); i++)
  705. {
  706. // check for the smallest c < 0
  707. if ((*c)[i] < min_r)
  708. {
  709. min_r = (*c)[i];
  710. piv_col = i;
  711. }
  712. }
  713. // no c < 0, we're done
  714. if (min_r >= 0)
  715. {
  716. stop = true;
  717. // check for infeasibility
  718. bool infeasible = true;
  719. for (int i = 0; i < (*c).size(); i++)
  720. {
  721. // any c > 0, feasible
  722. if ((*c)[i] > 0)
  723. {
  724. infeasible = false;
  725. break;
  726. }
  727. }
  728. if (infeasible)
  729. {
  730. outputInfeasible(a, b, c);
  731. continue;
  732. }
  733. // get x* from the matrix
  734. vector<double> x;
  735. // for each column
  736. for (int i = 0; i < numVars; i++)
  737. {
  738. // check if x_i is basic
  739. bool basic = true;
  740. int basic_i = 0;
  741. for (int j = 0; j < (*a).size(); j++)
  742. {
  743. // not 0 or 1, not basic
  744. if ((*a)[j][i] != 0 && (*a)[j][i] != 1)
  745. {
  746. basic = false;
  747. }
  748. // get the basic index
  749. if ((*a)[j][i] == 1)
  750. basic_i = j;
  751. }
  752. // basic x, get the b value
  753. if (basic)
  754. x.push_back((*b)[basic_i]);
  755. // nonbasic, 0
  756. else
  757. x.push_back(0.0);
  758. }
  759. // get optimal
  760. opt = 0;
  761. // multiply x values with cost function, sum up
  762. for (int i = 0; i < x.size(); i++) {
  763. opt += x[i] * initialC[i];
  764. }
  765. // min function, just output
  766. if (minimize)
  767. {
  768. output(a, b, c, x, opt, numVars);
  769. }
  770. // max function, output with negative optimal
  771. else if (!minimize)
  772. {
  773. output(a, b, c, x, -1 * opt, numVars);
  774. }
  775. continue;
  776. }
  777. // get the pivot row
  778. vector<double> ratios;
  779. // for each row
  780. for (int i = 0; i < (*a).size(); i++)
  781. {
  782. if ((*a)[i][piv_col] == 0)
  783. ratios.push_back(-1);
  784. else
  785. {
  786. // denegerate loop ahead, don't pivot here
  787. if ((*a)[i][piv_col] <= 0 && (*b)[i] == 0)
  788. ratios.push_back(-1);
  789. else
  790. ratios.push_back((*b)[i] / (*a)[i][piv_col]);
  791. }
  792. }
  793. double min_ratio = ratios[0];
  794. int piv_row = 0;
  795. for (int i = 1; i < ratios.size(); i++)
  796. {
  797. if (ratios[i] >= 0 && (ratios[i] < min_ratio || min_ratio < 0))
  798. {
  799. min_ratio = ratios[i];
  800. piv_row = i;
  801. }
  802. }
  803. // no positive ratios, problem is unbounded
  804. if (min_ratio < 0)
  805. {
  806. stop = true;
  807. outputUnbounded(a, b, c, minimize);
  808. continue;
  809. }
  810. // reduce the tableau on the pivot row and column
  811. reduce(a, b, c, &opt, piv_row, piv_col);
  812. // return to step 1 of the Simplex algorithm
  813. }
  814. }
  815. }
  816. void Calculator::gomoriInit(double a, double b)
  817. {
  818. }