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 }