1 module nxt.numeric_ex;
2 
3 import std.range: isInputRange, ElementType;
4 import std.traits: CommonType, isFloatingPoint;
5 import std.algorithm: reduce;
6 
7 /* import std.numeric: sum; */
8 
9 /** TODO Issue 4725: Remove when sum is standard in Phobos.
10     See_Also: http://d.puremagic.com/issues/show_bug.cgi?id=4725
11     See_Also: https://github.com/D-Programming-Language/phobos/pull/1205
12     See_Also: http://forum.dlang.org/thread/bug-4725-3@http.d.puremagic.com%2Fissues%2F
13  */
14 auto sum(Range, SumType = ElementType!Range)(Range range)
15     @safe pure nothrow if (isInputRange!Range)
16 {
17     return reduce!"a+b"(0, range);
18 }
19 
20 ///
21 unittest { assert([1, 2, 3, 4].sum == 10); }
22 
23 /** TODO Remove when product is standard in Phobos. */
24 auto product(Range)(Range range)
25     @safe pure nothrow if (isInputRange!Range)
26 {
27     return reduce!"a*b"(1, range);
28 }
29 
30 ///
31 unittest { assert([1, 2, 3, 4].product == 24); }
32 
33 // ==============================================================================================
34 
35 version(none)
36 {
37 
38     /** Computes $(LUCKY Discrete Signal Entropy) of input range $(D range).
39      */
40     auto signalEntropy(Range, RequestedBinType = double)(in Range range) @safe pure
41         if (isInputRange!Range &&
42             !is(CommonType!(ElementType!Range, F, G) == void))
43     {
44         enum normalized = true; // we need normalized histogram
45         import nxt.ngram: histogram, Kind, Storage, Symmetry;
46         auto hist = range.histogram!(Kind.saturated,
47                                      Storage.denseDynamic,
48                                      Symmetry.ordered,
49                                      RequestedBinType);
50         import std.numeric: entropy;
51         hist.normalize;
52         return hist[].entropy;
53     }
54 
55     ///
56     unittest
57     {
58         const ubyte[] p1 = [ 0, 1, 0, 1 ];
59         assert(p1.signalEntropy == 1);
60 
61         const ubyte[] p2 = [ 0, 1, 2, 3 ];
62         assert(p2.signalEntropy == 2);
63 
64         const ubyte[] p3 = [ 0, 255, 0, 255];
65         assert(p3.signalEntropy == 1);
66     }
67 
68     /** Returns: Element in $(D r ) that minimizes $(D fun).
69         LaTeX: \underset{x}{\arg\min} */
70     auto argmin_(alias fun, Range)(in Range r) @safe pure
71         if (isInputRange!Range &&
72             is(typeof(fun(r.front) <
73                       fun(r.front)) == bool))
74     {
75         auto bestX = r.front;
76         auto bestY = fun(bestX);
77         r.popFront();
78         foreach (e; r) {
79             auto candY = fun(e);     // candidate
80             if (candY > bestY) continue;
81             bestX = e;
82             bestY = candY;
83         }
84         return bestX;
85     }
86 
87     ///
88     unittest
89     {
90         assert(argmin!(x => x*x)([1, 2, 3]) == 1);
91         assert(argmin!(x => x*x)([3, 2, 1]) == 1);
92     }
93 
94   /** Returns: Element in $(D r ) that maximizes $(D fun).
95       LaTeX: \underset{x}{\arg\max} */
96     auto argmax_(alias fun, Range)(in Range r) @safe pure
97         if (isInputRange!Range &&
98             is(typeof(fun(r.front) >
99                       fun(r.front)) == bool))
100     {
101         auto bestX = r.front;
102         auto bestY = fun(bestX);
103         r.popFront();
104         foreach (e; r) {
105             auto candY = fun(e);     // candidate
106             if (candY < bestY) continue;
107             bestX = e;
108             bestY = candY;
109         }
110         return bestX;
111     }
112 
113     ///
114     unittest
115     {
116         assert(argmax!(x => x*x)([1, 2, 3]) == 3);
117         assert(argmax!(x => x*x)([3, 2, 1]) == 3);
118     }
119 
120 }
121 
122 /** Returns: Element in $(D r ) that minimizes $(D fun).
123     LaTeX: \underset{x}{\arg\min} */
124 auto argmin(alias fun, Range)(in Range r)
125     @safe pure if (isInputRange!Range &&
126                    is(typeof(fun(r.front) < fun(r.front)) == bool))
127 {
128     import std.front: front;
129     return typeof(r.front).max.reduce!((a,b) => fun(a) < fun(b) ? a : b)(r);
130 }
131 
132 ///
133 unittest
134 {
135     /* assert(argmin!(x => x*x)([1, 2, 3]) == 1); */
136     /* assert(argmin!(x => x*x)([3, 2, 1]) == 1); */
137 }
138 
139 /** Returns: Element in $(D r ) that maximizes $(D fun).
140     LaTeX: \underset{x}{\arg\max} */
141 auto argmax(alias fun, Range)(in Range r) @safe pure
142     if (isInputRange!Range &&
143         is(typeof(fun(r.front) >
144                   fun(r.front)) == bool))
145 {
146     import std.front: front;
147     return typeof(r.front).min.reduce!((a,b) => fun(a) > fun(b) ? a : b)(r);
148 }
149 
150 ///
151 unittest
152 {
153     /* assert(argmax!(x => x*x)([1, 2, 3]) == 3); */
154     /* assert(argmax!(x => x*x)([3, 2, 1]) == 3); */
155 }
156 
157 // ==============================================================================================
158 
159 /// Returns: 0.0 if x < edge, otherwise it returns 1.0.
160 CommonType!(T1, T2) step(T1, T2)(T1 edge, T2 x)
161 {
162     return x < edge ? 0 : 1;
163 }
164 
165 ///
166 unittest
167 {
168     assert(step(0, 1) == 1.0f);
169     assert(step(0, 10) == 1.0f);
170     assert(step(1, 0) == 0.0f);
171     assert(step(10, 0) == 0.0f);
172     assert(step(1, 1) == 1.0f);
173 }
174 
175 // ==============================================================================================
176 
177 /** Smoothstep from $(D edge0) to $(D edge1) at $(D x).
178     Returns: 0.0 if x <= edge0 and 1.0 if x >= edge1 and performs smooth
179     hermite interpolation between 0 and 1 when edge0 < x < edge1.
180     This is useful in cases where you would want a threshold function with a smooth transition.
181 */
182 CommonType!(T1,T2,T3) smoothstep(T1, T2, T3) (T1 edge0, T2 edge1, T3 x) @safe pure nothrow
183     if (isFloatingPoint!(CommonType!(T1,T2,T3)))
184 {
185     import std.algorithm: clamp;
186     x = clamp((x - edge0) / (edge1 - edge0), 0, 1);
187     return x * x * (3 - 2 * x);
188 }
189 
190 ///
191 unittest
192 {
193     //  assert(smoothstep(1, 0, 2) == 0);
194     assert(smoothstep(1.0, 0.0, 2.0) == 0);
195     assert(smoothstep(1.0, 0.0, 0.5) == 0.5);
196     // assert(almost_equal(smoothstep(0.0, 2.0, 0.5), 0.15625, 0.00001));
197 }
198 
199 /** Smootherstep from $(D edge0) to $(D edge1) at $(D x). */
200 @safe pure nothrow E smootherstep(E)(E edge0, E edge1, E x)
201     if (isFloatingPoint!(E))
202 {
203     // Scale, and clamp x to 0..1 range
204     import std.algorithm: clamp;
205     x = clamp((x - edge0) / (edge1 - edge0), 0.0, 1.0);
206     // Evaluate polynomial
207     return x*x*x*(x*(x*6 - 15) + 10); // evaluate polynomial
208 }
209 
210 ///
211 unittest
212 {
213     assert(smootherstep(1.0, 0.0, 2.0) == 0);
214 }