API tools faq
paste
Advertisement
Guest User

Untitled

a guest
Aug 23rd, 2017
78
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 17.89 KB | None | 0 0
raw download clone embed print report
  1. // functions from numeric.js that my solution needs
  2. function base0to1(A) {
  3. if(typeof A !== "object") { return A; }
  4. var ret = [], i,n=A.length;
  5. for(i=0;i<n;i++) ret[i+1] = base0to1(A[i]);
  6. return ret;
  7. }
  8. function base1to0(A) {
  9. if(typeof A !== "object") { return A; }
  10. var ret = [], i,n=A.length;
  11. for(i=1;i<n;i++) ret[i-1] = base1to0(A[i]);
  12. return ret;
  13. }
  14.  
  15. function dpofa(a, lda, n, info) {
  16. var i, j, jm1, k, t, s;
  17.  
  18. for (j = 1; j <= n; j = j + 1) {
  19. info[1] = j;
  20. s = 0;
  21. jm1 = j - 1;
  22. if (jm1 < 1) {
  23. s = a[j][j] - s;
  24. if (s <= 0) {
  25. break;
  26. }
  27. a[j][j] = Math.sqrt(s);
  28. } else {
  29. for (k = 1; k <= jm1; k = k + 1) {
  30. //~ t = a[k][j] - ddot(k - 1, a[1][k], 1, a[1][j], 1);
  31. t = a[k][j];
  32. for (i = 1; i < k; i = i + 1) {
  33. t = t - (a[i][j] * a[i][k]);
  34. }
  35. t = t / a[k][k];
  36. a[k][j] = t;
  37. s = s + t * t;
  38. }
  39. s = a[j][j] - s;
  40. if (s <= 0) {
  41. break;
  42. }
  43. a[j][j] = Math.sqrt(s);
  44. }
  45. info[1] = 0;
  46. }
  47. }
  48.  
  49. function dposl(a, lda, n, b) {
  50. var i, k, kb, t;
  51.  
  52. for (k = 1; k <= n; k = k + 1) {
  53. //~ t = ddot(k - 1, a[1][k], 1, b[1], 1);
  54. t = 0;
  55. for (i = 1; i < k; i = i + 1) {
  56. t = t + (a[i][k] * b[i]);
  57. }
  58.  
  59. b[k] = (b[k] - t) / a[k][k];
  60. }
  61.  
  62. for (kb = 1; kb <= n; kb = kb + 1) {
  63. k = n + 1 - kb;
  64. b[k] = b[k] / a[k][k];
  65. t = -b[k];
  66. //~ daxpy(k - 1, t, a[1][k], 1, b[1], 1);
  67. for (i = 1; i < k; i = i + 1) {
  68. b[i] = b[i] + (t * a[i][k]);
  69. }
  70. }
  71. }
  72.  
  73. function dpori(a, lda, n) {
  74. var i, j, k, kp1, t;
  75.  
  76. for (k = 1; k <= n; k = k + 1) {
  77. a[k][k] = 1 / a[k][k];
  78. t = -a[k][k];
  79. //~ dscal(k - 1, t, a[1][k], 1);
  80. for (i = 1; i < k; i = i + 1) {
  81. a[i][k] = t * a[i][k];
  82. }
  83.  
  84. kp1 = k + 1;
  85. if (n < kp1) {
  86. break;
  87. }
  88. for (j = kp1; j <= n; j = j + 1) {
  89. t = a[k][j];
  90. a[k][j] = 0;
  91. //~ daxpy(k, t, a[1][k], 1, a[1][j], 1);
  92. for (i = 1; i <= k; i = i + 1) {
  93. a[i][j] = a[i][j] + (t * a[i][k]);
  94. }
  95. }
  96. }
  97.  
  98. }
  99.  
  100. function qpgen2(dmat, dvec, fddmat, n, sol, crval, amat,
  101. bvec, fdamat, q, meq, iact, nact, iter, work, ierr) {
  102.  
  103. var i, j, l, l1, info, it1, iwzv, iwrv, iwrm, iwsv, iwuv, nvl, r, iwnbv,
  104. temp, sum, t1, tt, gc, gs, nu,
  105. t1inf, t2min,
  106. vsmall, tmpa, tmpb,
  107. go;
  108.  
  109. r = Math.min(n, q);
  110. l = 2 * n + (r * (r + 5)) / 2 + 2 * q + 1;
  111.  
  112. vsmall = 1.0e-60;
  113. do {
  114. vsmall = vsmall + vsmall;
  115. tmpa = 1 + 0.1 * vsmall;
  116. tmpb = 1 + 0.2 * vsmall;
  117. } while (tmpa <= 1 || tmpb <= 1);
  118.  
  119. for (i = 1; i <= n; i = i + 1) {
  120. work[i] = dvec[i];
  121. }
  122. for (i = n + 1; i <= l; i = i + 1) {
  123. work[i] = 0;
  124. }
  125. for (i = 1; i <= q; i = i + 1) {
  126. iact[i] = 0;
  127. }
  128.  
  129. info = [];
  130.  
  131. if (ierr[1] === 0) {
  132. dpofa(dmat, fddmat, n, info);
  133. if (info[1] !== 0) {
  134. ierr[1] = 2;
  135. return;
  136. }
  137. dposl(dmat, fddmat, n, dvec);
  138. dpori(dmat, fddmat, n);
  139. } else {
  140. for (j = 1; j <= n; j = j + 1) {
  141. sol[j] = 0;
  142. for (i = 1; i <= j; i = i + 1) {
  143. sol[j] = sol[j] + dmat[i][j] * dvec[i];
  144. }
  145. }
  146. for (j = 1; j <= n; j = j + 1) {
  147. dvec[j] = 0;
  148. for (i = j; i <= n; i = i + 1) {
  149. dvec[j] = dvec[j] + dmat[j][i] * sol[i];
  150. }
  151. }
  152. }
  153.  
  154. crval[1] = 0;
  155. for (j = 1; j <= n; j = j + 1) {
  156. sol[j] = dvec[j];
  157. crval[1] = crval[1] + work[j] * sol[j];
  158. work[j] = 0;
  159. for (i = j + 1; i <= n; i = i + 1) {
  160. dmat[i][j] = 0;
  161. }
  162. }
  163. crval[1] = -crval[1] / 2;
  164. ierr[1] = 0;
  165.  
  166. iwzv = n;
  167. iwrv = iwzv + n;
  168. iwuv = iwrv + r;
  169. iwrm = iwuv + r + 1;
  170. iwsv = iwrm + (r * (r + 1)) / 2;
  171. iwnbv = iwsv + q;
  172.  
  173. for (i = 1; i <= q; i = i + 1) {
  174. sum = 0;
  175. for (j = 1; j <= n; j = j + 1) {
  176. sum = sum + amat[j][i] * amat[j][i];
  177. }
  178. work[iwnbv + i] = Math.sqrt(sum);
  179. }
  180. nact = 0;
  181. iter[1] = 0;
  182. iter[2] = 0;
  183.  
  184. function fn_goto_50() {
  185. iter[1] = iter[1] + 1;
  186.  
  187. l = iwsv;
  188. for (i = 1; i <= q; i = i + 1) {
  189. l = l + 1;
  190. sum = -bvec[i];
  191. for (j = 1; j <= n; j = j + 1) {
  192. sum = sum + amat[j][i] * sol[j];
  193. }
  194. if (Math.abs(sum) < vsmall) {
  195. sum = 0;
  196. }
  197. if (i > meq) {
  198. work[l] = sum;
  199. } else {
  200. work[l] = -Math.abs(sum);
  201. if (sum > 0) {
  202. for (j = 1; j <= n; j = j + 1) {
  203. amat[j][i] = -amat[j][i];
  204. }
  205. bvec[i] = -bvec[i];
  206. }
  207. }
  208. }
  209.  
  210. for (i = 1; i <= nact; i = i + 1) {
  211. work[iwsv + iact[i]] = 0;
  212. }
  213.  
  214. nvl = 0;
  215. temp = 0;
  216. for (i = 1; i <= q; i = i + 1) {
  217. if (work[iwsv + i] < temp * work[iwnbv + i]) {
  218. nvl = i;
  219. temp = work[iwsv + i] / work[iwnbv + i];
  220. }
  221. }
  222. if (nvl === 0) {
  223. return 999;
  224. }
  225.  
  226. return 0;
  227. }
  228.  
  229. function fn_goto_55() {
  230. for (i = 1; i <= n; i = i + 1) {
  231. sum = 0;
  232. for (j = 1; j <= n; j = j + 1) {
  233. sum = sum + dmat[j][i] * amat[j][nvl];
  234. }
  235. work[i] = sum;
  236. }
  237.  
  238. l1 = iwzv;
  239. for (i = 1; i <= n; i = i + 1) {
  240. work[l1 + i] = 0;
  241. }
  242. for (j = nact + 1; j <= n; j = j + 1) {
  243. for (i = 1; i <= n; i = i + 1) {
  244. work[l1 + i] = work[l1 + i] + dmat[i][j] * work[j];
  245. }
  246. }
  247.  
  248. t1inf = true;
  249. for (i = nact; i >= 1; i = i - 1) {
  250. sum = work[i];
  251. l = iwrm + (i * (i + 3)) / 2;
  252. l1 = l - i;
  253. for (j = i + 1; j <= nact; j = j + 1) {
  254. sum = sum - work[l] * work[iwrv + j];
  255. l = l + j;
  256. }
  257. sum = sum / work[l1];
  258. work[iwrv + i] = sum;
  259. if (iact[i] < meq) {
  260. // continue;
  261. break;
  262. }
  263. if (sum < 0) {
  264. // continue;
  265. break;
  266. }
  267. t1inf = false;
  268. it1 = i;
  269. }
  270.  
  271. if (!t1inf) {
  272. t1 = work[iwuv + it1] / work[iwrv + it1];
  273. for (i = 1; i <= nact; i = i + 1) {
  274. if (iact[i] < meq) {
  275. // continue;
  276. break;
  277. }
  278. if (work[iwrv + i] < 0) {
  279. // continue;
  280. break;
  281. }
  282. temp = work[iwuv + i] / work[iwrv + i];
  283. if (temp < t1) {
  284. t1 = temp;
  285. it1 = i;
  286. }
  287. }
  288. }
  289.  
  290. sum = 0;
  291. for (i = iwzv + 1; i <= iwzv + n; i = i + 1) {
  292. sum = sum + work[i] * work[i];
  293. }
  294. if (Math.abs(sum) <= vsmall) {
  295. if (t1inf) {
  296. ierr[1] = 1;
  297. // GOTO 999
  298. return 999;
  299. } else {
  300. for (i = 1; i <= nact; i = i + 1) {
  301. work[iwuv + i] = work[iwuv + i] - t1 * work[iwrv + i];
  302. }
  303. work[iwuv + nact + 1] = work[iwuv + nact + 1] + t1;
  304. // GOTO 700
  305. return 700;
  306. }
  307. } else {
  308. sum = 0;
  309. for (i = 1; i <= n; i = i + 1) {
  310. sum = sum + work[iwzv + i] * amat[i][nvl];
  311. }
  312. tt = -work[iwsv + nvl] / sum;
  313. t2min = true;
  314. if (!t1inf) {
  315. if (t1 < tt) {
  316. tt = t1;
  317. t2min = false;
  318. }
  319. }
  320.  
  321. for (i = 1; i <= n; i = i + 1) {
  322. sol[i] = sol[i] + tt * work[iwzv + i];
  323. if (Math.abs(sol[i]) < vsmall) {
  324. sol[i] = 0;
  325. }
  326. }
  327.  
  328. crval[1] = crval[1] + tt * sum * (tt / 2 + work[iwuv + nact + 1]);
  329. for (i = 1; i <= nact; i = i + 1) {
  330. work[iwuv + i] = work[iwuv + i] - tt * work[iwrv + i];
  331. }
  332. work[iwuv + nact + 1] = work[iwuv + nact + 1] + tt;
  333.  
  334. if (t2min) {
  335. nact = nact + 1;
  336. iact[nact] = nvl;
  337.  
  338. l = iwrm + ((nact - 1) * nact) / 2 + 1;
  339. for (i = 1; i <= nact - 1; i = i + 1) {
  340. work[l] = work[i];
  341. l = l + 1;
  342. }
  343.  
  344. if (nact === n) {
  345. work[l] = work[n];
  346. } else {
  347. for (i = n; i >= nact + 1; i = i - 1) {
  348. if (work[i] === 0) {
  349. // continue;
  350. break;
  351. }
  352. gc = Math.max(Math.abs(work[i - 1]), Math.abs(work[i]));
  353. gs = Math.min(Math.abs(work[i - 1]), Math.abs(work[i]));
  354. if (work[i - 1] >= 0) {
  355. temp = Math.abs(gc * Math.sqrt(1 + gs * gs / (gc * gc)));
  356. } else {
  357. temp = -Math.abs(gc * Math.sqrt(1 + gs * gs / (gc * gc)));
  358. }
  359. gc = work[i - 1] / temp;
  360. gs = work[i] / temp;
  361.  
  362. if (gc === 1) {
  363. // continue;
  364. break;
  365. }
  366. if (gc === 0) {
  367. work[i - 1] = gs * temp;
  368. for (j = 1; j <= n; j = j + 1) {
  369. temp = dmat[j][i - 1];
  370. dmat[j][i - 1] = dmat[j][i];
  371. dmat[j][i] = temp;
  372. }
  373. } else {
  374. work[i - 1] = temp;
  375. nu = gs / (1 + gc);
  376. for (j = 1; j <= n; j = j + 1) {
  377. temp = gc * dmat[j][i - 1] + gs * dmat[j][i];
  378. dmat[j][i] = nu * (dmat[j][i - 1] + temp) - dmat[j][i];
  379. dmat[j][i - 1] = temp;
  380.  
  381. }
  382. }
  383. }
  384. work[l] = work[nact];
  385. }
  386. } else {
  387. sum = -bvec[nvl];
  388. for (j = 1; j <= n; j = j + 1) {
  389. sum = sum + sol[j] * amat[j][nvl];
  390. }
  391. if (nvl > meq) {
  392. work[iwsv + nvl] = sum;
  393. } else {
  394. work[iwsv + nvl] = -Math.abs(sum);
  395. if (sum > 0) {
  396. for (j = 1; j <= n; j = j + 1) {
  397. amat[j][nvl] = -amat[j][nvl];
  398. }
  399. bvec[nvl] = -bvec[nvl];
  400. }
  401. }
  402. // GOTO 700
  403. return 700;
  404. }
  405. }
  406.  
  407. return 0;
  408. }
  409.  
  410. function fn_goto_797() {
  411. l = iwrm + (it1 * (it1 + 1)) / 2 + 1;
  412. l1 = l + it1;
  413. if (work[l1] === 0) {
  414. // GOTO 798
  415. return 798;
  416. }
  417. gc = Math.max(Math.abs(work[l1 - 1]), Math.abs(work[l1]));
  418. gs = Math.min(Math.abs(work[l1 - 1]), Math.abs(work[l1]));
  419. if (work[l1 - 1] >= 0) {
  420. temp = Math.abs(gc * Math.sqrt(1 + gs * gs / (gc * gc)));
  421. } else {
  422. temp = -Math.abs(gc * Math.sqrt(1 + gs * gs / (gc * gc)));
  423. }
  424. gc = work[l1 - 1] / temp;
  425. gs = work[l1] / temp;
  426.  
  427. if (gc === 1) {
  428. // GOTO 798
  429. return 798;
  430. }
  431. if (gc === 0) {
  432. for (i = it1 + 1; i <= nact; i = i + 1) {
  433. temp = work[l1 - 1];
  434. work[l1 - 1] = work[l1];
  435. work[l1] = temp;
  436. l1 = l1 + i;
  437. }
  438. for (i = 1; i <= n; i = i + 1) {
  439. temp = dmat[i][it1];
  440. dmat[i][it1] = dmat[i][it1 + 1];
  441. dmat[i][it1 + 1] = temp;
  442. }
  443. } else {
  444. nu = gs / (1 + gc);
  445. for (i = it1 + 1; i <= nact; i = i + 1) {
  446. temp = gc * work[l1 - 1] + gs * work[l1];
  447. work[l1] = nu * (work[l1 - 1] + temp) - work[l1];
  448. work[l1 - 1] = temp;
  449. l1 = l1 + i;
  450. }
  451. for (i = 1; i <= n; i = i + 1) {
  452. temp = gc * dmat[i][it1] + gs * dmat[i][it1 + 1];
  453. dmat[i][it1 + 1] = nu * (dmat[i][it1] + temp) - dmat[i][it1 + 1];
  454. dmat[i][it1] = temp;
  455. }
  456. }
  457.  
  458. return 0;
  459. }
  460.  
  461. function fn_goto_798() {
  462. l1 = l - it1;
  463. for (i = 1; i <= it1; i = i + 1) {
  464. work[l1] = work[l];
  465. l = l + 1;
  466. l1 = l1 + 1;
  467. }
  468.  
  469. work[iwuv + it1] = work[iwuv + it1 + 1];
  470. iact[it1] = iact[it1 + 1];
  471. it1 = it1 + 1;
  472. if (it1 < nact) {
  473. // GOTO 797
  474. return 797;
  475. }
  476.  
  477. return 0;
  478. }
  479.  
  480. function fn_goto_799() {
  481. work[iwuv + nact] = work[iwuv + nact + 1];
  482. work[iwuv + nact + 1] = 0;
  483. iact[nact] = 0;
  484. nact = nact - 1;
  485. iter[2] = iter[2] + 1;
  486.  
  487. return 0;
  488. }
  489.  
  490. go = 0;
  491. while (true) {
  492. go = fn_goto_50();
  493. if (go === 999) {
  494. return;
  495. }
  496. while (true) {
  497. go = fn_goto_55();
  498. if (go === 0) {
  499. break;
  500. }
  501. if (go === 999) {
  502. return;
  503. }
  504. if (go === 700) {
  505. if (it1 === nact) {
  506. fn_goto_799();
  507. } else {
  508. while (true) {
  509. fn_goto_797();
  510. go = fn_goto_798();
  511. if (go !== 797) {
  512. break;
  513. }
  514. }
  515. fn_goto_799();
  516. }
  517. }
  518. }
  519. }
  520.  
  521. }
  522.  
  523. function solveQP(Dmat, dvec, Amat, bvec, meq, factorized) {
  524. Dmat = base0to1(Dmat);
  525. dvec = base0to1(dvec);
  526. Amat = base0to1(Amat);
  527. var i, n, q,
  528. nact, r,
  529. crval = [], iact = [], sol = [], work = [], iter = [],
  530. message;
  531.  
  532. meq = meq || 0;
  533. factorized = factorized ? base0to1(factorized) : [undefined, 0];
  534. bvec = bvec ? base0to1(bvec) : [];
  535.  
  536. // In Fortran the array index starts from 1
  537. n = Dmat.length - 1;
  538. q = Amat[1].length - 1;
  539.  
  540. if (!bvec) {
  541. for (i = 1; i <= q; i = i + 1) {
  542. bvec[i] = 0;
  543. }
  544. }
  545. for (i = 1; i <= q; i = i + 1) {
  546. iact[i] = 0;
  547. }
  548. nact = 0;
  549. r = Math.min(n, q);
  550. for (i = 1; i <= n; i = i + 1) {
  551. sol[i] = 0;
  552. }
  553. crval[1] = 0;
  554. for (i = 1; i <= (2 * n + (r * (r + 5)) / 2 + 2 * q + 1); i = i + 1) {
  555. work[i] = 0;
  556. }
  557. for (i = 1; i <= 2; i = i + 1) {
  558. iter[i] = 0;
  559. }
  560.  
  561. qpgen2(Dmat, dvec, n, n, sol, crval, Amat,
  562. bvec, n, q, meq, iact, nact, iter, work, factorized);
  563.  
  564. message = "";
  565. if (factorized[1] === 1) {
  566. message = "constraints are inconsistent, no solution!";
  567. }
  568. if (factorized[1] === 2) {
  569. message = "matrix D in quadratic function is not positive definite!";
  570. }
  571.  
  572. return {
  573. solution: base1to0(sol),
  574. value: base1to0(crval),
  575. unconstrained_solution: base1to0(dvec),
  576. iterations: base1to0(iter),
  577. iact: base1to0(iact),
  578. message: message
  579. };
  580. }
  581.  
  582. // main
  583. // solve QP with 1/2 (xT) D (x) - (dT)(x) s.t. Ax >= b
  584.  
  585. var h = [2, 5, 5]; // vector of heights
  586. var a = [8, 15, 16]; // input vector
  587. var N = h.length;
  588. var eta = 2;
  589. var k = 5;
  590.  
  591. // print h
  592. document.write("h...n");
  593. for (i = 0; i < N; i++){
  594. document.write(h[i] + ", ");
  595. }
  596.  
  597.  
  598. // print a
  599. document.write("na...n");
  600. for (i = 0; i < N; i++){
  601. document.write(a[i] + ", ");
  602. }
  603.  
  604. // print N
  605. document.write("nN..." + N + "n");
  606.  
  607. // print eta
  608. document.write("neta..." + eta + "n");
  609.  
  610. // print k
  611. document.write("nk..." + k + "n");
  612.  
  613. // construct DMat = 2I, and dvec = 2a
  614. var DMat = [];
  615. for (i = 0; i < N; i++){
  616. DMat[i] = [];
  617. for (j = 0; j < N; j++){
  618. if (i != j){
  619. DMat[i][j] = 0;
  620. }
  621. else{
  622. DMat[i][j] = 2;
  623. }
  624. }
  625. }
  626.  
  627. // print DMat
  628. document.write("nDMat...n");
  629. for (i = 0; i < N; i++){
  630. for (j = 0; j < N; j++){
  631. document.write(DMat[i][j] + ", ");
  632. }
  633. document.write("n");
  634. }
  635.  
  636. // now dvec
  637. var dvec = [];
  638. for (i = 0; i < N; i++){
  639. dvec[i] = a[i] * 2;
  640. }
  641.  
  642. // print dvec
  643. document.write("ndvec...n");
  644. for (i = 0; i < N; i++){
  645. document.write(dvec[i] + ", ");
  646. }
  647.  
  648. // construct AMat and bvec for constraints
  649. var AMat = [];
  650. for (i = 0; i < N; i++){
  651. AMat[i] = [];
  652. for (j = 0; j < N; j++){
  653. if (i == j){
  654. AMat[i][j] = 1;
  655. }
  656. else if ((j+1) == i){
  657. AMat[i][j] = -1;
  658. }
  659. else{
  660. AMat[i][j] = 0;
  661. }
  662. }
  663. }
  664.  
  665. // print AMat
  666. document.write("nAMat...n");
  667. for (i = 0; i < N; i++){
  668. for (j = 0; j < N; j++){
  669. document.write(AMat[i][j] + ", ");
  670. }
  671. document.write("n");
  672. }
  673.  
  674. // now bvec
  675. var bvec = [];
  676. for (i = 0; i < N; i++){
  677. if (i == 0){
  678. bvec[i] = k;
  679. }
  680. else{
  681. bvec[i] = h[i - 1] + eta;
  682. }
  683. }
  684.  
  685. // print bvec
  686. document.write("nbvec...n");
  687. for (i = 0; i < N; i++){
  688. document.write(bvec[i] + ", ");
  689. }
  690.  
  691. // solve QP with 1/2 (xT) D (x) - (dT)(x) s.t. Ax >= b
  692. var res = solveQP(DMat, dvec, AMat, bvec, 0, false);
  693. document.write("nSolution...n");
  694.  
  695. for (i = 0; i < res.solution.length; i++){
  696. document.write(res.solution[i] + ", ");
  697. }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!