1
0

metaviscon.cpp 8.4 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279
  1. #include <iostream>
  2. #include <vector>
  3. #include <string>
  4. #include <algorithm>
  5. #include <random>
  6. struct r2 {
  7. double x = 0, y = 0;
  8. };
  9. std::string convertToStdString(const std::vector<bool>& vec);
  10. double euclidenDistance(const std::vector<bool>& vecA, const std::vector<bool>& vecB, int amountIfBits);
  11. double neighbourProbability(const std::vector<std::vector<bool>>& solutions, const int j, const int i, int amountOfBits, double omega);
  12. double euclidenDistance(const r2& y1, const r2& y2);
  13. double calculate_tDistributed_q(const std::vector<r2>& solutions, int j, int i);
  14. double calculate_thing(const r2& y1, r2 y2);
  15. double cost(const std::vector<double>& p, std::vector<double>& q, size_t amountOfSolutions);
  16. r2 computeGradient(int i, const std::vector<double>& q, const std::vector<double>& p, const std::vector<r2>& ySolutions, int amountOfSolutions);
  17. void updateQ(const size_t& amountOfSolutions, std::vector<double>& matrixQ, std::vector<r2>& ySoltuion);
  18. int main23(int argc, char * argv[] )
  19. {
  20. if (argc > 1) {
  21. std::cout << "Hello World! " << argv[1] << std::endl;
  22. }
  23. //Generate Solutions;
  24. const int amountOfBits = 100;
  25. size_t amountOfSolutions = 5;
  26. const double omega = 5;
  27. /*size_t amountOfSolutions = (size_t)std::pow(2, amountOfBits);
  28. std::vector<std::vector<bool> > solutions(amountOfSolutions);
  29. for (int i = 0; i < amountOfSolutions; i++) {
  30. std::vector<bool> newVec(amountOfBits);
  31. for (int k = 0; k < amountOfBits; k++) {
  32. int maxHalf = (int)std::pow(2, k);
  33. int max = maxHalf * 2;
  34. newVec[k] = i % max < maxHalf;
  35. }
  36. solutions[i] = newVec;
  37. }
  38. auto rng = std::default_random_engine{};
  39. std::shuffle(solutions.begin(), solutions.end(), rng);
  40. solutions.erase(solutions.begin() + solutions.size() / 2);
  41. amountOfSolutions /= 2;*/
  42. std::vector<std::vector<bool> > solutions(amountOfSolutions);
  43. std::random_device rd; //Will be used to obtain a seed for the random number engine
  44. std::mt19937 gen(rd()); //Standard mersenne_twister_engine seeded with rd()
  45. std::uniform_real_distribution<double> doubleDistr(0.0, 1.0);
  46. for (int i = 0; i < amountOfSolutions; i++) {
  47. std::vector<bool> newVec(amountOfBits);
  48. double rnd = doubleDistr(gen);
  49. for (int k = 0; k < amountOfBits; k++) {
  50. newVec[k] = doubleDistr(gen) < rnd;
  51. }
  52. solutions[i] = newVec;
  53. }
  54. for (int i = 0; i < amountOfSolutions; i++) {
  55. std::cout << "[" << convertToStdString(solutions[i]) << "]" << std::endl;
  56. }
  57. //Normal
  58. std::vector<double> matrixProbability(amountOfSolutions * amountOfSolutions);
  59. for (size_t i = 0; i < amountOfSolutions; i++) {
  60. for (size_t j = 0; j < amountOfSolutions; j++) {
  61. if (j == i) {
  62. matrixProbability[i * amountOfSolutions + j] = 0.0;
  63. }
  64. else {
  65. matrixProbability[i * amountOfSolutions + j] = neighbourProbability(solutions, j, i, amountOfBits, omega);
  66. }
  67. }
  68. }
  69. std::vector<double> matrixP(amountOfSolutions * amountOfSolutions);
  70. double sum = 0.0;
  71. //Symmettric
  72. for (size_t x = 0; x < amountOfSolutions; x++) {
  73. for (size_t y = 0; y < amountOfSolutions; y++) {
  74. double value = (matrixProbability[x * amountOfSolutions + y] + matrixProbability[y * amountOfSolutions + x]) / 2.0;
  75. matrixP[x * amountOfSolutions + y] = value;
  76. sum += value;
  77. }
  78. }
  79. //Normalized
  80. for (size_t x = 0; x < amountOfSolutions; x++) {
  81. for (size_t y = 0; y < amountOfSolutions; y++) {
  82. matrixP[x * amountOfSolutions + y] = matrixP[x * amountOfSolutions + y] / sum;
  83. }
  84. }
  85. //sample random solution
  86. std::vector<r2> ySolutions(amountOfSolutions);
  87. std::uniform_real_distribution<double> zeroDistr(-1.0, 1.0);
  88. std::for_each(ySolutions.begin(), ySolutions.end(), [&gen, &zeroDistr](r2& r) {r.x = zeroDistr(gen);r.y = zeroDistr(gen);});
  89. for (const r2& r : ySolutions) {
  90. std::cout << r.x << ", " << r.y << std::endl;
  91. }
  92. //generate q vec
  93. std::vector<double> matrixQ(amountOfSolutions * amountOfSolutions);
  94. updateQ(amountOfSolutions, matrixQ, ySolutions);
  95. /*std::cout << "MatrixP:" << std::endl;
  96. for (size_t x = 0; x < amountOfSolutions; x++) {
  97. for (size_t y = 0; y < amountOfSolutions; y++) {
  98. std::cout << matrixP[x * amountOfSolutions + y] << " ";
  99. }
  100. std::cout << std::endl;
  101. }
  102. std::cout << "MatrixQ:" << std::endl;
  103. for (size_t x = 0; x < amountOfSolutions; x++) {
  104. for (size_t y = 0; y < amountOfSolutions; y++) {
  105. std::cout << matrixQ[x * amountOfSolutions + y] << " ";
  106. }
  107. std::cout << std::endl;
  108. }
  109. std::cout << "Cost: " << cost(matrixP, matrixQ, amountOfSolutions) << std::endl;*/
  110. //update
  111. int iterations = 100;
  112. double learningRate = 10;
  113. std::vector<r2> update(ySolutions);
  114. for (int iter = 0; iter < iterations; iter++) {
  115. //updateQ
  116. updateQ(amountOfSolutions, matrixQ, ySolutions);
  117. for (size_t i = 0; i < amountOfSolutions; i++) {
  118. //computeGradient
  119. const r2 gradientFromI= computeGradient(i, matrixQ, matrixP, ySolutions, amountOfSolutions);
  120. r2& yi = ySolutions[i];
  121. //update Y
  122. r2& ui = update[i];
  123. ui.x = learningRate * gradientFromI.x + 0.5 * ui.x;
  124. ui.y = learningRate * gradientFromI.y + 0.5 * ui.y;
  125. yi.x = yi.x + ui.x;
  126. yi.y = yi.y + ui.y;
  127. }
  128. //
  129. //std::cout << "Cost: " << cost(matrixP, matrixQ, amountOfSolutions) << std::endl;
  130. std::cout << "Iter: " << iter << std::endl;
  131. for (const r2& r : ySolutions) {
  132. std::cout << r.x << ", " << r.y << std::endl;
  133. }
  134. getchar();
  135. }
  136. return 0;
  137. }
  138. void updateQ(const size_t& amountOfSolutions, std::vector<double>& matrixQ, std::vector<r2>& ySoltuion)
  139. {
  140. double sum = 0.0;
  141. for (size_t i = 0; i < amountOfSolutions; i++) {
  142. for (size_t j = 0; j < amountOfSolutions; j++) {
  143. if (j == i) {
  144. matrixQ[i * amountOfSolutions + j] = 0.0;
  145. }
  146. else {
  147. double value = calculate_tDistributed_q(ySoltuion, j, i);
  148. matrixQ[i * amountOfSolutions + j] = value;
  149. sum += value;
  150. }
  151. }
  152. }
  153. for (size_t i = 0; i < matrixQ.size(); i++) matrixQ[i] /= sum;
  154. }
  155. std::string convertToStdString(const std::vector<bool>& vec) {
  156. std::string result("");
  157. int count = 0;
  158. for (const bool& boolean : vec) {
  159. result += boolean ? "1" : "0";
  160. if (boolean) count++;
  161. }
  162. return result + ":" + std::to_string(count);
  163. }
  164. double euclidenDistance(const std::vector<bool>& vecA, const std::vector<bool>& vecB, int amountOfBits) {
  165. int count = 0;
  166. for (int i = 0; i < amountOfBits; i++) {
  167. if (vecA[i] ^ vecB[i]) count++;
  168. }
  169. return std::sqrt(count);
  170. }
  171. double neighbourProbability(const std::vector<std::vector<bool>>& solutions,const int j,const int i, int amountOfBits, double omega) {
  172. double numerator = std::exp(-euclidenDistance(solutions[i], solutions[j], amountOfBits)/(2 * omega * omega));
  173. double denominator = 0.0;
  174. for (int k = 0; k < solutions.size(); k++) {
  175. if (k == i) continue;
  176. denominator += std::exp(-euclidenDistance(solutions[i], solutions[k], amountOfBits) / (2 * omega * omega));
  177. }
  178. return numerator/ denominator;
  179. }
  180. double euclidenDistance(const r2& y1, const r2& y2) {
  181. return std::sqrt(std::pow(y1.x - y2.x, 2) + std::pow(y1.y - y2.y, 2));
  182. }
  183. double calculate_tDistributed_q(const std::vector<r2>& solutions, int j, int i) {
  184. double numerator = calculate_thing(solutions[i], solutions[j]);
  185. double denominator = 0.0;
  186. for (int k = 0; k < solutions.size(); k++) {
  187. if (k == i) continue;
  188. denominator += calculate_thing(solutions[i], solutions[k]);
  189. }
  190. return numerator / denominator;
  191. }
  192. double calculate_thing(const r2& y1, r2 y2) {
  193. return 1 / (1 + euclidenDistance(y1, y2));
  194. }
  195. double cost(const std::vector<double>& p, std::vector<double>& q, size_t amountOfSolutions) {
  196. double sum = 0.0;
  197. for (size_t x = 0; x < amountOfSolutions; x++) {
  198. for (size_t y = 0; y < amountOfSolutions; y++) {
  199. sum += q[x * amountOfSolutions + y];
  200. }
  201. }
  202. for (size_t x = 0; x < amountOfSolutions; x++) {
  203. for (size_t y = 0; y < amountOfSolutions; y++) {
  204. q[x * amountOfSolutions + y] /= sum;
  205. }
  206. }
  207. double cost = 0.0;
  208. for (size_t i = 0; i < amountOfSolutions; i++) {
  209. for (size_t j = 0; j < amountOfSolutions; j++) {
  210. if (i == j) continue;
  211. double pji = p[j * amountOfSolutions + i];
  212. double value = pji * std::log(pji / q[j * amountOfSolutions + i]);
  213. cost += value;
  214. }
  215. }
  216. return cost;
  217. }
  218. r2 computeGradient(int i, const std::vector<double>& q, const std::vector<double>& p, const std::vector<r2>& ySolutions, int amountOfSolutions) {
  219. r2 value;
  220. const r2& yi = ySolutions[i];
  221. for (size_t j = 0; j < amountOfSolutions; j++) {
  222. if (j == i) continue;
  223. double scalar = p[j * amountOfSolutions + i] - q[j * amountOfSolutions + i];
  224. const r2& yj = ySolutions[j];
  225. value.x += scalar * (yi.x - yj.y);
  226. value.y += scalar * (yi.y - yj.y);
  227. }
  228. value.x *= 4;
  229. value.y *= 4;
  230. return value;
  231. }