1 /**
2 Authors: [Ilya Yaroshenko](http://9il.github.io)
3 
4 Copyright: © 2014-2015 [Ilya Yaroshenko](http://9il.github.io)
5 
6 License: MIT
7 */
8 module atmosphere.utilities;
9 
10 import core.stdc.tgmath;
11 
12 import core.stdc.string;
13 import std.traits;
14 import std.range;
15 import std.compiler;
16 
17 import std.math : isFinite, isIdentical, approxEqual, isNaN, NaN;
18 
19 package:
20 
21 version(LDC)
22 {
23 	import ldc.intrinsics;
24 	pragma(LDC_inline_ir)
25 		R inlineIR(string s, R, P...)(P);
26 }
27 
28 package:
29 
30 import cblas;
31 import simple_matrix;
32 
33 version(LDC)
34 {
35 	T sum(T)(in T[] a)
36 	{
37 		T ret = 0;
38 		foreach(j; 0..a.length)
39 			static if(is(Unqual!T == double))
40 			ret = inlineIR!(`
41 				%r = fadd fast double %0, %1
42 				ret double %r`, double)(ret, a[j]);
43 			else
44 			static if(is(Unqual!T == float))
45 			ret = inlineIR!(`
46 				%r = fadd fast float %0, %1
47 				ret float %r`, float)(ret, a[j]);
48 			else
49 			ret += a[j];
50 		return ret;
51 	}
52 }
53 else
54 {
55 	T sum(T)(in T[] a)
56 	{
57 		T ret = 0;
58 		foreach(j; 0..a.length)
59 			ret += a[j];
60 		return ret;
61 	}
62 }
63 
64 unittest {
65 	import std.range : iota, array;
66 	static import std.algorithm;
67 	foreach(i; 0.0..30.0)
68 		assert(std.algorithm.sum(iota(i)) == iota(i).array.sum);
69 }
70 
71 
72 ///
73 auto avg(Range)(Range range)
74 {
75 	return range.sum / range.length;
76 }
77 
78 ///
79 void normalize(F)(F[] range)
80 {
81 	immutable s = range.sum;
82 	assert(s.isFinite);
83 	assert(s > 0);
84 	foreach(ref elem; range)
85 		elem /= s;
86 }
87 
88 ///
89 void gemv(M, F)(M m, in F[] a, F[] b)
90 in {
91 	assert (m.width == a.length);
92 	assert (m.height == b.length);
93 }
94 body {
95 	static if(is(M : Matrix!(T), T) && (is(Unqual!T == double) || is(Unqual!T == float)))
96 	{
97 		assert(m.ptr);
98 		assert(m.shift >= m.width);
99 		cblas.gemv(
100 			Order.RowMajor,
101 			Transpose.NoTrans,
102 			cast(blasint)b.length,
103 		 	cast(blasint)a.length,
104 			1,//F
105 			m.ptr,
106 			cast(blasint)m.shift,
107 			a.ptr,
108 			1,
109 			0,//F
110 			b.ptr,
111 			1);
112 	}
113 	else
114 	static if(is(M : TransposedMatrix!T, T) && is(Unqual!T == double) || is(Unqual!T == float))
115 	{
116 		assert(m.matrix.ptr);
117 		assert(m.matrix.shift >= m.matrix.width);
118 		cblas.gemv(
119 			Order.RowMajor,
120 			Transpose.Trans,
121 			cast(blasint)a.length,
122 		 	cast(blasint)b.length,
123 			1.0,
124 			m.matrix.ptr,
125 			cast(blasint)m.matrix.shift,
126 			a.ptr,
127 			1,
128 			0.0,
129 			b.ptr,
130 			1);
131 	}
132 	else
133 	{
134 		foreach(ref e; b)
135 		{
136 			assert(!m.empty);
137 			e = dotProduct(a, m.front);
138 			m.popFront;
139 		}
140 	}
141 }
142 
143 unittest
144 {
145 	const ar = [
146 	 1.000,  6.000,   2.000,
147 	 8.000,  3.000,   7.000,
148 	 3.000,  5.000,   2.000,
149 	53.000, 23.000, 123.000,
150 	];
151 	auto m = Matrix!(const double)(ar.ptr, 4, 3);
152 	const a = [
153 	42.000,
154 	35.000,
155 	12.000,
156 	];
157 	auto b = new double[4];
158 	gemv(m, a, b);
159 	assert(b == [ 
160 	 276.000,
161 	 525.000,
162 	 325.000,
163 	4507.000,
164 	]);
165 
166 }
167 
168 unittest
169 {
170 	const ar = [
171   	1.000,   8.000,  3.000,  53.000,
172   	6.000,   3.000,  5.000,  23.000,
173   	2.000,   7.000,  2.000, 123.000,
174 	];
175 	auto m = Matrix!(const double)(ar.ptr, 3, 4);
176 	const a = [
177 	42.000,
178 	35.000,
179 	12.000,
180 	];
181 	auto b = new double[4];
182 	gemv(m.transposed, a, b);
183 	assert(b == [ 
184 	 276.000,
185 	 525.000,
186 	 325.000,
187 	4507.000,
188 	]);
189 
190 }
191 
192 
193 version(LDC)
194 {
195 	auto dotProduct(Range1, Range2)(Range1 a, Range2 b)
196 		if(isInputRange!Range1
197 		&& isInputRange!Range2 
198 		&& (!isArray!Range1 || !isArray!Range2)
199 		&& is(Unqual!(ElementType!Range1) == Unqual!(ElementType!Range2))
200 		)
201 	{
202 		alias T = Unqual!(ElementType!Range1);
203 		T ret = 0;
204 		while(!a.empty)
205 		{
206 			static if(is(Unqual!T == double))
207 				ret = inlineIR!(`
208 					%d = fmul fast double %1, %2
209 					%r = fadd fast double %0, %d
210 					ret double %r`, double)(ret, a.front, b.front);
211 			else
212 			static if(is(Unqual!T == float))
213 				ret = inlineIR!(`
214 					%d = fmul fast float %1, %2
215 					%r = fadd fast float %0, %d
216 					ret float %r`, float)(ret, a.front, b.front);
217 			else
218 				ret += a.front * b.front;
219 			a.popFront;
220 			b.popFront;
221 		}
222 		assert(b.empty);
223 		return ret;
224 	}
225 }
226 else
227 {
228 	auto dotProduct(Range1, Range2)(Range1 a, Range2 b)
229 		if(isInputRange!Range1
230 		&& isInputRange!Range2 
231 		&& (!isArray!Range1 || !isArray!Range2)
232 		&& is(Unqual!(ElementType!Range1) : Unqual!(ElementType!Range2))
233 		)
234 	{
235 		alias T = Unqual!(CommonType!(ElementType!Range1, ElementType!Range2));
236 		T ret = 0;
237 		while(!a.empty)
238 		{
239 			ret += a.front * b.front;
240 			a.popFront;
241 			b.popFront;
242 		}
243 		assert(b.empty);
244 		return ret;
245 	}
246 }
247 
248 version(LDC)
249 {
250 	T dotProduct(T)(in T[] a, in T[] b)
251 	{
252 		T ret = 0;
253 		foreach(j; 0..a.length)
254 			static if(is(Unqual!T == double))
255 			ret = inlineIR!(`
256 				%d = fmul fast double %1, %2
257 				%r = fadd fast double %0, %d
258 				ret double %r`, double)(ret, a[j], b[j]);
259 			else
260 			static if(is(Unqual!T == float))
261 			ret = inlineIR!(`
262 				%d = fmul fast float %1, %2
263 				%r = fadd fast float %0, %d
264 				ret float %r`, float)(ret, a[j], b[j]);
265 			else
266 			ret += a[j] * b[j];
267 		return ret;
268 	}
269 }
270 else
271 {
272 	T dotProduct(T)(in T[] a, in T[] b)
273 	{
274 		T ret = 0;
275 		foreach(j; 0..a.length)
276 			ret += a[j] * b[j];
277 		return ret;
278 	}
279 }
280 
281 version(LDC)
282 {
283 	T dotProductInverse(T)(in T[] a, in T[] b)
284 	{
285 		T ret = 0;
286 		foreach(j; 0..a.length)
287 			static if(is(Unqual!T == double))
288 			ret = inlineIR!(`
289 				%d = fdiv fast double %1, %2
290 				%r = fadd fast double %0, %d
291 				ret double %r`, double)(ret, a[j], b[j]);
292 			else
293 			static if(is(Unqual!T == float))
294 			ret = inlineIR!(`
295 				%d = fdiv fast float %1, %2
296 				%r = fadd fast float %0, %d
297 				ret float %r`, float)(ret, a[j], b[j]);
298 			else
299 			ret += a[j] / b[j];
300 		return ret;
301 	}
302 }
303 else
304 {
305 	T dotProductInverse(T)(in T[] a, in T[] b)
306 	{
307 		T ret = 0;
308 		foreach(j; 0..a.length)
309 			ret += a[j] / b[j];
310 		return ret;
311 	}
312 }
313 
314 
315 version(LDC)
316 {
317 	T dotProductInverse2(T)(in T[] a, in T[] b, T[] c)
318 	{
319 		T ret = 0;
320 		foreach(j; 0..a.length)
321 		{
322 			static if(is(Unqual!T == double))
323 			{
324 				ret = inlineIR!(`
325 					%d = fdiv fast double %1, %2
326 					%r = fadd fast double %0, %d
327 					ret double %r`, double)(ret, b[j], a[j]);
328 				c[j] = a[j] - b[j];
329 			}
330 			else
331 			static if(is(Unqual!T == float))
332 			{
333 				ret = inlineIR!(`
334 					%d = fdiv fast float %1, %2
335 					%r = fadd fast float %0, %d
336 					ret float %r`, float)(ret, b[j], a[j]);
337 				c[j] = a[j] - b[j];
338 			}
339 			else
340 			{
341 				ret += b[j] / a[j];
342 				c[j] = a[j] - b[j];
343 			}
344 		}
345 		return ret;
346 	}
347 }
348 else
349 {
350 	T dotProductInverse2(T)(in T[] a, in T[] b, T[] c)
351 	{
352 		T ret = 0;
353 		foreach(j; 0..a.length)
354 		{
355 			ret += b[j] / a[j];
356 			c[j] = a[j] - b[j];
357 		}
358 		return ret;
359 	}
360 }
361 
362 /**
363 Struct that represent flat matrix.
364 Useful for sliding windows.
365 */
366 struct MatrixColumnsSlider(F)
367 {
368 	Matrix!F _matrix;
369 	Matrix!F matrix;
370 
371 	this(size_t maxHeight, size_t maxWidth, size_t height)
372 	{
373 		_matrix = Matrix!F(maxHeight, maxWidth);
374 		_matrix.width = _matrix.shift;
375 		matrix.ptr = _matrix.ptr;
376 		matrix.shift = _matrix.shift;
377 		matrix.height = height;
378 	}
379 
380 	void popFrontN(size_t n)
381 	in 
382 	{
383 		assert(n <= matrix.width, "n > matrix.width");
384 	}
385 	body 
386 	{
387 		if(n < matrix.width)
388 		{
389 			matrix.width -= n;
390 			matrix.ptr += n;
391 		}
392 		else
393 		{ 
394 			reset;
395 		}
396 	}
397 
398 	void popFront()
399 	{
400 		popFrontN(1);
401 	}
402 
403 	void reset()
404 	{
405 		matrix.ptr = _matrix.ptr;
406 		matrix.width = 0;
407 	}
408 
409 	void putBackN(size_t n)
410 	in
411 	{
412 		assert(matrix.shift >= matrix.width+n);
413 	}
414 	body 
415 	{
416 		if(n > _matrix.ptrEnd-matrix.ptrEnd)
417 		{
418 			bringToFront();
419 		}
420 		matrix.width += n;
421 	}
422 
423 	void putBack()
424 	{
425 		putBackN(1);
426 	}
427 
428 	void bringToFront()
429 	{
430 		if(matrix.width)
431 		{
432 			memmove(_matrix.ptr, matrix.ptr, (matrix.shift*matrix.height)*F.sizeof);					
433 		}
434 		matrix.ptr = _matrix.ptr;
435 	}
436 }
437 
438 ///
439 template convertTo(alias InterfaceTemp)
440 {
441 	InterfaceTemp!F convertTo(Fun, F = ReturnType!Fun)(Fun fun)
442 	{
443 		static assert(isFloatingPoint!F);
444 		return new class (fun) InterfaceTemp!F {
445 			Fun fun;
446 			this(Fun fun) { this.fun = fun; }
447 			F opCall(F x) { return fun(x); }
448 		};
449 	}
450 }
451 
452 ///
453 unittest
454 {
455 	import std.math;
456 	import atmosphere.pdf;
457 
458 	real fun(real x)
459 	{
460 		// 1/sqrt(2 PI)
461 		enum c = 0.398942280401432677939946L;
462 		return c * exp(-0.5f * x * x);
463 	}
464 
465 	PDF!real pdf = convertTo!PDF(&fun);
466 }
467 
468 /++
469 Generates random permutation
470 +/
471 size_t[] randomPermutation(size_t length)
472 {
473 	import core.memory;
474 	import std.random : uniform;
475 	import std.algorithm : makeIndex;
476 	auto indexesR = new size_t[length];
477 	scope(exit) 
478 		GC.free(indexesR.ptr);
479 	auto indexesS = new size_t[length];
480 	foreach(j, ref index; indexesR)
481 	{
482 		index = uniform!"[]"(0, size_t.max);
483 		indexesS[j] = j;
484 	}
485 	makeIndex(indexesR, indexesS);
486 	return indexesS;
487 }