1 module lbfgsd.math; 2 3 private static import std.math; 4 5 private static import std.numeric; 6 7 private import lbfgsd.autodiff; 8 9 /// optimized (x * x) 10 auto square(T)(const scope T x) @safe @nogc pure nothrow 11 { 12 static if (is(T : Variable!(U, N), U, size_t N)) 13 { 14 T y = void; 15 y.d[] = (x.a + x.a) * x.d[]; 16 y.a = x.a * x.a; 17 return y; 18 } 19 else 20 return x * x; 21 } 22 23 unittest 24 { 25 alias Var = Variable!(double, 2); 26 auto x = Var(1.0, 0); 27 auto y = Var(2.0, 1); 28 29 auto sx = square(x); 30 assert(sx.a == 1); 31 assert(sx.d[0] == 2); 32 assert(sx.d[1] == 0); 33 34 auto sy = square(y); 35 assert(sy.a == 4); 36 assert(sy.d[0] == 0); 37 assert(sy.d[1] == 4); 38 39 auto z = 3.0; 40 auto sz = square(z); 41 assert(sz == 9); 42 } 43 44 auto dotProduct(T, U)(const scope T[] a, const scope U[] b) 45 { 46 static if (is(T : Variable!(S, N), S, size_t N) || is(U : Variable!(W, M), W, size_t M)) 47 { 48 assert(a.length == b.length); 49 50 import std.traits; 51 52 alias Q = Unqual!(typeof(T.init * U.init)); 53 auto sum0 = Q(0), sum1 = Q(0); 54 55 const all_endp = a.length; 56 const smallblock_endp = all_endp & ~3; 57 const bigblock_endp = all_endp & ~15; 58 59 size_t i = 0; 60 auto ap = a[]; 61 auto bp = b[]; 62 for (; i != bigblock_endp; i += 16, ap = ap[16 .. $], bp = bp[16 .. $]) 63 { 64 sum0 += ap[0] * bp[0]; 65 sum1 += ap[1] * bp[1]; 66 sum0 += ap[2] * bp[2]; 67 sum1 += ap[3] * bp[3]; 68 sum0 += ap[4] * bp[4]; 69 sum1 += ap[5] * bp[5]; 70 sum0 += ap[6] * bp[6]; 71 sum1 += ap[7] * bp[7]; 72 sum0 += ap[8] * bp[8]; 73 sum1 += ap[9] * bp[9]; 74 sum0 += ap[10] * bp[10]; 75 sum1 += ap[11] * bp[11]; 76 sum0 += ap[12] * bp[12]; 77 sum1 += ap[13] * bp[13]; 78 sum0 += ap[14] * bp[14]; 79 sum1 += ap[15] * bp[15]; 80 } 81 82 for (; i != smallblock_endp; i += 4, ap = ap[4 .. $], bp = bp[4 .. $]) 83 { 84 sum0 += ap[0] * bp[0]; 85 sum1 += ap[1] * bp[1]; 86 sum0 += ap[2] * bp[2]; 87 sum1 += ap[3] * bp[3]; 88 } 89 90 for (; i != all_endp; ++i, ap = ap[1 .. $], bp = bp[1 .. $]) 91 { 92 sum0 += ap[0] * bp[0]; 93 } 94 95 return sum0 + sum1; 96 } 97 else 98 return std.numeric.dotProduct(a, b); 99 } 100 101 unittest 102 { 103 alias Var = Variable!(double, 3); 104 auto xs = new Var[3]; 105 auto ys = new Var[3]; 106 foreach (i; 0 .. xs.length) 107 { 108 ys[i] = xs[i] = Var(i, i); 109 } 110 auto z = dotProduct(xs, ys); 111 assert(z.a == 5); 112 assert(z.d[0] == 0); 113 assert(z.d[1] == 2); 114 assert(z.d[2] == 4); 115 } 116 117 unittest 118 { 119 auto xs = new double[3]; 120 auto ys = new double[3]; 121 foreach (i; 0 .. xs.length) 122 { 123 ys[i] = xs[i] = i; 124 } 125 auto z = dotProduct(xs, ys); 126 assert(z == 5); 127 } 128 129 unittest 130 { 131 alias Var = Variable!(double, 3); 132 auto xs = new Var[3]; 133 auto ys = new double[3]; 134 foreach (i; 0 .. xs.length) 135 { 136 xs[i] = Var(i, i); 137 ys[i] = i; 138 } 139 auto z = dotProduct(xs, ys); 140 assert(z.a == 5); 141 assert(z.d[0] == 0); 142 assert(z.d[1] == 1); 143 assert(z.d[2] == 2); 144 145 auto w = dotProduct(ys, xs); 146 assert(z.a == w.a); 147 assert(z.d[0] == w.d[0]); 148 assert(z.d[1] == w.d[1]); 149 assert(z.d[2] == w.d[2]); 150 } 151 152 unittest 153 { 154 alias Var = Variable!(double, 1000); 155 auto xs = new Var[1000]; 156 foreach (i; 0 .. xs.length) 157 xs[i] = Var(i, i); 158 159 auto y = dotProduct(xs, xs); 160 assert(y.a == 332833500); 161 foreach (i; 0 .. xs.length) 162 assert(y.d[i] == 2 * i); 163 } 164 165 T sum(T)(const scope T[] xs) @safe @nogc pure nothrow 166 { 167 assert(xs.length > 0); 168 T y = xs[0]; 169 foreach (i; 1 .. xs.length) 170 y += xs[i]; 171 return y; 172 } 173 174 @safe pure nothrow unittest 175 { 176 alias Var = Variable!(double, 3); 177 auto xs = new Var[3]; 178 xs[0] = Var(0.2, 0); 179 xs[1] = Var(0.4, 1); 180 xs[2] = Var(0.6, 2); 181 182 auto sx = sum(xs); 183 assert(std.math.approxEqual(sx.a, 1.2)); 184 assert(sx.d[0] == 1); 185 assert(sx.d[1] == 1); 186 assert(sx.d[2] == 1); 187 188 auto ys = new double[3]; 189 ys[0] = 0.2; 190 ys[1] = 0.4; 191 ys[2] = 0.6; 192 auto sy = sum(ys); 193 assert(std.math.approxEqual(sy, 1.2)); 194 } 195 196 T sumsq(T)(const scope T[] xs) @safe @nogc pure nothrow 197 { 198 assert(xs.length > 0); 199 T y = square(xs[0]); 200 foreach (i; 1 .. xs.length) 201 y += square(xs[i]); 202 return y; 203 } 204 205 @safe pure nothrow unittest 206 { 207 alias Var = Variable!(double, 3); 208 auto xs = new Var[3]; 209 xs[0] = Var(0.2, 0); 210 xs[1] = Var(0.4, 1); 211 xs[2] = Var(0.6, 2); 212 213 auto sx = sumsq(xs); 214 assert(sx.a == 0.56); 215 assert(sx.d[0] == 0.4); 216 assert(sx.d[1] == 0.8); 217 assert(sx.d[2] == 1.2); 218 219 auto ys = new double[3]; 220 ys[0] = 0.2; 221 ys[1] = 0.4; 222 ys[2] = 0.6; 223 224 auto sy = sumsq(ys); 225 assert(sy == 0.56); 226 } 227 228 T sumxmy2(T, U)(const scope T[] xs, const scope U[] ys) @safe @nogc pure nothrow 229 { 230 assert(xs.length > 0); 231 assert(xs.length == ys.length); 232 233 T sum = square(xs[0] - ys[0]); 234 foreach (i; 1 .. xs.length) 235 sum += square(xs[i] - ys[i]); 236 return sum; 237 } 238 239 unittest 240 { 241 alias Var = Variable!(double, 3); 242 auto xs = new Var[3]; 243 xs[0] = Var(0.2, 0); 244 xs[1] = Var(0.4, 1); 245 xs[2] = Var(0.6, 2); 246 auto ys = new Var[3]; 247 ys[0] = Var(0.1, 0); 248 ys[1] = Var(0.2, 1); 249 ys[2] = Var(0.3, 2); 250 251 auto s = sumxmy2(xs, ys); 252 assert(s == square(xs[0] - ys[0]) + square(xs[1] - ys[1]) + square(xs[2] - ys[2])); 253 } 254 255 T sqrt(T)(const scope T x) @safe @nogc pure nothrow 256 { 257 static if (is(T : Variable!(U, N), U, size_t N)) 258 { 259 const t = std.math.sqrt(x.a); 260 T y = void; 261 y.d[] = x.d[] * (0.5 / t); 262 y.a = t; 263 return y; 264 } 265 else 266 return std.math.sqrt(x); 267 } 268 269 @safe pure nothrow unittest 270 { 271 alias Var = Variable!(double, 2); 272 auto x = Var(2.0, 0); 273 double y = 2.0; 274 275 auto z = sqrt(x); 276 auto w = sqrt(y); 277 278 assert(z.a == w); 279 assert(std.math.approxEqual(z.d[0], 1 / (2 * w))); 280 assert(std.math.approxEqual(z.d[1], 0)); 281 } 282 283 T exp(T)(const scope T x) @safe @nogc pure nothrow 284 { 285 static if (is(T : Variable!(U, N), U, size_t N)) 286 { 287 const t = std.math.exp(x.a); 288 T y = void; 289 y.d[] = t * x.d[]; 290 y.a = t; 291 return y; 292 } 293 else 294 return std.math.exp(x); 295 } 296 297 @safe pure nothrow unittest 298 { 299 alias Var = Variable!(double, 2); 300 auto x = Var(1.0, 0); 301 double y = 1.0; 302 303 auto z = exp(x); 304 auto w = exp(y); 305 306 assert(z.a == w); 307 assert(std.math.approxEqual(z.d[0], w)); 308 assert(std.math.approxEqual(z.d[1], 0)); 309 } 310 311 T log(T)(const scope T x) @safe @nogc pure nothrow 312 { 313 static if (is(T : Variable!(U, N), U, size_t N)) 314 { 315 T y = void; 316 y.d[] = x.d[] / x.a; 317 y.a = std.math.log(x.a); 318 return y; 319 } 320 else 321 return std.math.log(x); 322 } 323 324 @safe pure nothrow unittest 325 { 326 alias Var = Variable!(double, 2); 327 auto x = Var(2.0, 0); 328 double y = 2.0; 329 330 auto z = log(x); 331 auto w = log(y); 332 333 assert(z.a == w); 334 assert(std.math.approxEqual(z.d[0], 0.5)); 335 assert(std.math.approxEqual(z.d[1], 0)); 336 } 337 338 T sin(T)(const scope T x) @safe @nogc pure nothrow 339 { 340 static if (is(T : Variable!(U, N), U, size_t N)) 341 { 342 T y = void; 343 y.d[] = std.math.cos(x.a) * x.d[]; 344 y.a = std.math.sin(x.a); 345 return y; 346 } 347 else 348 return std.math.sin(x); 349 } 350 351 @safe pure nothrow unittest 352 { 353 alias Var = Variable!(double, 2); 354 auto x = Var(2.0, 0); 355 double y = 2.0; 356 357 auto z = sin(x); 358 auto w = sin(y); 359 360 assert(std.math.approxEqual(z.a, 0.909297427)); 361 assert(std.math.approxEqual(z.d[0], -0.416146837)); 362 assert(std.math.approxEqual(z.d[1], 0)); 363 assert(std.math.approxEqual(w, 0.909297427)); 364 } 365 366 T cos(T)(const scope T x) @safe @nogc pure nothrow 367 { 368 static if (is(T : Variable!(U, N), U, size_t N)) 369 { 370 T y = void; 371 y.d[] = -std.math.sin(x.a) * x.d[]; 372 y.a = std.math.cos(x.a); 373 return y; 374 } 375 else 376 return std.math.cos(x); 377 } 378 379 @safe pure nothrow unittest 380 { 381 alias Var = Variable!(double, 2); 382 auto x = Var(2.0, 0); 383 double y = 2.0; 384 385 auto z = cos(x); 386 auto w = cos(y); 387 388 assert(std.math.approxEqual(z.a, -0.416146837)); 389 assert(std.math.approxEqual(z.d[0], -0.909297427)); 390 assert(std.math.approxEqual(z.d[1], 0)); 391 assert(std.math.approxEqual(w, -0.416146837)); 392 } 393 394 T tan(T)(const scope T x) @safe @nogc pure nothrow 395 { 396 static if (is(T : Variable!(U, N), U, size_t N)) 397 { 398 T y = void; 399 auto t = std.math.tan(x.a); 400 y.d[] = (1 + t * t) * x.d[]; 401 y.a = t; 402 return y; 403 } 404 else 405 return std.math.tan(x); 406 } 407 408 @safe pure nothrow unittest 409 { 410 alias Var = Variable!(double, 2); 411 auto x = Var(2.0, 0); 412 double y = 2.0; 413 414 auto z = tan(x); 415 auto w = tan(y); 416 417 assert(std.math.approxEqual(z.a, w)); 418 assert(std.math.approxEqual(z.d[0], 1 + w * w)); 419 assert(std.math.approxEqual(z.d[1], 0)); 420 } 421 422 T sinh(T)(const scope T x) @safe @nogc pure nothrow 423 { 424 static if (is(T : Variable!(U, N), U, size_t N)) 425 { 426 T y = void; 427 y.d[] = std.math.cosh(x.a) * x.d[]; 428 y.a = std.math.sinh(x.a); 429 return y; 430 } 431 else 432 return std.math.sinh(x); 433 } 434 435 @safe pure nothrow unittest 436 { 437 alias Var = Variable!(double, 2); 438 auto x = Var(1.0, 0); 439 double y = 1.0; 440 441 auto z = sinh(x); 442 auto w = sinh(y); 443 444 assert(z.a == w); 445 assert(std.math.approxEqual(z.d[0], std.math.cosh(1.0))); 446 assert(std.math.approxEqual(z.d[1], 0)); 447 } 448 449 T cosh(T)(const scope T x) @safe @nogc pure nothrow 450 { 451 static if (is(T : Variable!(U, N), U, size_t N)) 452 { 453 T y = void; 454 y.d[] = std.math.sinh(x.a) * x.d[]; 455 y.a = std.math.cosh(x.a); 456 return y; 457 } 458 else 459 return std.math.cosh(x); 460 } 461 462 @safe pure nothrow unittest 463 { 464 alias Var = Variable!(double, 2); 465 auto x = Var(1.0, 0); 466 double y = 1.0; 467 468 auto z = cosh(x); 469 auto w = cosh(y); 470 471 assert(z.a == w); 472 assert(std.math.approxEqual(z.d[0], std.math.sinh(1.0))); 473 assert(std.math.approxEqual(z.d[1], 0)); 474 } 475 476 T tanh(T)(const scope T x) @safe @nogc pure nothrow 477 { 478 static if (is(T : Variable!(U, N), U, size_t N)) 479 { 480 const t = std.math.tanh(x.a); 481 T y = void; 482 y.d[] = (1 - t * t) * x.d[]; 483 y.a = t; 484 return y; 485 } 486 else 487 return std.math.tanh(x); 488 } 489 490 @safe pure nothrow unittest 491 { 492 alias Var = Variable!(double, 2); 493 auto x = Var(1.0, 0); 494 double y = 1.0; 495 496 auto z = tanh(x); 497 auto w = tanh(y); 498 499 assert(z.a == w); 500 assert(std.math.approxEqual(z.d[0], 1 - w * w)); 501 assert(std.math.approxEqual(z.d[1], 0)); 502 } 503 504 T asinh(T)(const scope T x) @safe @nogc pure nothrow 505 { 506 static if (is(T : Variable!(U, N), U, size_t N)) 507 { 508 T y = void; 509 y.d[] = x.d[] / std.math.sqrt(x.a * x.a + 1); 510 y.a = std.math.asinh(x.a); 511 return y; 512 } 513 else 514 return std.math.asinh(x); 515 } 516 517 @safe pure nothrow unittest 518 { 519 alias Var = Variable!(double, 2); 520 auto x = Var(0.5, 0); 521 double y = 0.5; 522 523 auto z = asinh(x); 524 auto w = asinh(y); 525 526 assert(z.a == w); 527 assert(std.math.approxEqual(z.d[0], 0.894427)); 528 assert(std.math.approxEqual(z.d[1], 0)); 529 } 530 531 T acosh(T)(const scope T x) @safe @nogc pure nothrow 532 { 533 static if (is(T : Variable!(U, N), U, size_t N)) 534 { 535 T y = void; 536 y.d[] = x.d[] / std.math.sqrt(x.a * x.a - 1); 537 y.a = std.math.acosh(x.a); 538 return y; 539 } 540 else 541 return std.math.acosh(x); 542 } 543 544 @safe pure nothrow unittest 545 { 546 alias Var = Variable!(double, 2); 547 auto x = Var(1.5, 0); 548 double y = 1.5; 549 550 auto z = acosh(x); 551 auto w = acosh(y); 552 553 assert(z.a == w); 554 assert(std.math.approxEqual(z.d[0], 0.894427)); 555 assert(std.math.approxEqual(z.d[1], 0)); 556 } 557 558 T atanh(T)(const scope T x) @safe @nogc pure nothrow 559 { 560 static if (is(T : Variable!(U, N), U, size_t N)) 561 { 562 T y = void; 563 y.d[] = x.d[] / (1 - x.a * x.a); 564 y.a = std.math.atanh(x.a); 565 return y; 566 } 567 else 568 return std.math.atanh(x); 569 } 570 571 @safe pure nothrow unittest 572 { 573 alias Var = Variable!(double, 2); 574 auto x = Var(0.5, 0); 575 double y = 0.5; 576 577 auto z = atanh(x); 578 auto w = atanh(y); 579 580 assert(z.a == w); 581 assert(std.math.approxEqual(z.d[0], 1.33333)); 582 assert(std.math.approxEqual(z.d[1], 0)); 583 } 584 585 /// correlation coefficient 586 auto correl(T, U)(const scope T[] xs, const scope U[] ys) 587 { 588 assert(xs.length == ys.length); 589 import std.traits; 590 591 alias Q = Unqual!(typeof(T.init * U.init)); 592 593 immutable N = xs.length; 594 auto xa = sum(xs) / N; 595 auto ya = sum(ys) / N; 596 597 return correl(xs, ys, xa, ya); 598 } 599 600 @safe unittest 601 { 602 alias Var = Variable!(double, 2); 603 604 auto xa = Var(1, 0); 605 auto ya = Var(1, 1); 606 607 auto xs = new Var[3]; 608 auto ys = new Var[3]; 609 foreach (i, ref x; xs) 610 x = xa * i; 611 foreach (i, ref y; ys) 612 y = ya * i; 613 614 auto c = correl(xs, ys); 615 import std.algorithm, std.math; 616 617 assert(approxEqual(c.a, 1.0)); 618 assert(equal!approxEqual(c.d[], [0.0, 0.0])); 619 } 620 621 auto correl(T, U)(const scope T[] xs, const scope U[] ys, const scope T xa, const scope T ya) 622 { 623 assert(xs.length == ys.length); 624 import std.traits; 625 626 alias Q = Unqual!(typeof(T.init * U.init)); 627 628 auto xv = T(0); 629 auto yv = U(0); 630 auto cv = Q(0); 631 632 foreach (i; 0 .. xs.length) 633 { 634 auto tx = xs[i] - xa; 635 auto ty = ys[i] - ya; 636 cv += tx * ty; 637 xv += square(tx); 638 yv += square(ty); 639 } 640 return cv / (sqrt(xv) * sqrt(yv)); 641 } 642 643 T sigmoid(T)(const scope T x) 644 { 645 return T(1) / (1 + exp(-x)); 646 } 647 648 T swish(T)(const scope T x) 649 { 650 return x * sigmoid(x); 651 } 652 653 T softplus(T)(const scope T x) 654 { 655 return log(1 + exp(x)); 656 } 657 658 T mish(T)(const scope T x) 659 { 660 return x * tanh(softplus(x)); 661 }