1 #Region "Microsoft.VisualBasic::8dc43eb21684515ed1a94093f961dbda, Microsoft.VisualBasic.Core\Language\Language\Java\MathUtils.vb"
2
3     ' Author:
4     ' 
5     '       asuka (amethyst.asuka@gcmodeller.org)
6     '       xie (genetics@smrucc.org)
7     '       xieguigang (xie.guigang@live.com)
8     ' 
9     ' Copyright (c) 2018 GPL3 Licensed
10     ' 
11     ' 
12     ' GNU GENERAL PUBLIC LICENSE (GPL3)
13     ' 
14     ' 
15     ' This program is free software: you can redistribute it and/or modify
16     ' it under the terms of the GNU General Public License as published by
17     ' the Free Software Foundation, either version 3 of the License, or
18     ' (at your option) any later version.
19     ' 
20     ' This program is distributed in the hope that it will be useful,
21     ' but WITHOUT ANY WARRANTY; without even the implied warranty of
22     ' MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
23     ' GNU General Public License for more details.
24     ' 
25     ' You should have received a copy of the GNU General Public License
26     ' along with this program. If not, see <http://www.gnu.org/licenses/>.
27
28
29
30     ' /********************************************************************************/
31
32     ' Summaries:
33
34     '     Module MathUtils
35     ' 
36     '         Properties: Seed
37     ' 
38     '         Function: getNormalized, (+2 Overloads) getTotal, hypot, nextBoolean, nextByte
39     '                   nextChar, nextDouble, nextExponential, nextFloat, nextGamma
40     '                   nextGaussian, (+2 Overloads) nextInt, nextInverseGaussian, nextLong, nextShort
41     '                   permuted, randomChoice, randomChoicePDF, randomLogDouble, sampleIndicesWithReplacement
42     '                   shuffled, uniform
43     ' 
44     '         Sub: nextBytes, permute, (+2 Overloads) shuffle
45     ' 
46     ' 
47     ' /********************************************************************************/
48
49 #End Region
50
51 Imports sys = System.Math
52
53 ' ' * MathUtils.java
54 ' *
55 ' * Copyright (C) 2002-2006 Alexei Drummond and Andrew Rambaut
56 ' *
57 ' * This file is part of BEAST.
58 ' * See the NOTICE file distributed with this work for additional
59 ' * information regarding copyright ownership and licensing.
60 ' *
61 ' * BEAST is free software; you can redistribute it and/or modify
62 ' * it under the terms of the GNU Lesser General Public License as
63 ' * published by the Free Software Foundation; either version 2
64 ' * of the License, or (at your option) any later version.
65 ' *
66 ' *  BEAST is distributed in the hope that it will be useful,
67 ' *  but WITHOUT ANY WARRANTY; without even the implied warranty of
68 ' *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
69 ' *  GNU Lesser General Public License for more details.
70 ' *
71 ' * You should have received a copy of the GNU Lesser General Public
72 ' * License along with BEAST; if not, write to the
73 ' * Free Software Foundation, Inc., 51 Franklin St, Fifth Floor,
74 ' * Boston, MA  02110-1301  USA
75 ' 
76
77 Namespace Language.Java
78
79     ''' <summary>
80     ''' Handy utility functions which have some Mathematical relavance.
81     ''' 
82     ''' @author Matthew Goode
83     ''' @author Alexei Drummond
84     ''' @author Gerton Lunter
85     ''' @version $Id: MathUtils.java,v 1.13 2006/08/31 14:57:24 rambaut Exp $
86     ''' </summary>
87     Public Module MathUtils
88
89         ''' <summary>
90         ''' A random number generator that is initialized with the clock when this
91         ''' class is loaded into the JVM. Use this for all random numbers. Note: This
92         ''' method or getting random numbers in not thread-safe. Since
93         ''' MersenneTwisterFast is currently (as of 9/01) not synchronized using this
94         ''' function may cause concurrency issues. Use the static get methods of the
95         ''' MersenneTwisterFast class for access to a single instance of the class,
96         ''' that has synchronization.
97         ''' </summary>
98         ReadOnly random As MersenneTwisterFast = MersenneTwisterFast.DEFAULT_INSTANCE
99
100         ''' <summary>
101         ''' Chooses one category if a cumulative probability distribution is given 
102         ''' </summary>
103         ''' <param name="cf"></param>
104         ''' <returns></returns>
105         Public Function randomChoice(cf As Double()) As Integer
106
107             Dim U As Double = random.nextDouble()
108
109             Dim s As Integer
110             If U <= cf(0) Then
111                 s = 0
112             Else
113                 For s = 1 To cf.Length - 1
114                     If U <= cf(s) AndAlso U > cf(s - 1) Then Exit For
115                 Next s
116             End If
117
118             Return s
119         End Function
120
121         ''' <param name="pdf">
122         '''            array of unnormalized probabilities </param>
123         ''' <returns> a sample according to an unnormalized probability distribution </returns>
124         Public Function randomChoicePDF(pdf As Double()) As Integer
125
126             Dim U As Double = random.nextDouble() * getTotal(pdf)
127             For i As Integer = 0 To pdf.Length - 1
128
129                 U -= pdf(i)
130                 If U < 0.0 Then Return i
131
132             Next i
133             For i As Integer = 0 To pdf.Length - 1
134                 Console.WriteLine(i & vbTab & pdf(i))
135             Next i
136             Throw New Exception("randomChoiceUnnormalized falls through -- negative components in input distribution?")
137         End Function
138
139         ''' <param name="array">
140         '''            to normalize </param>
141         ''' <returns> a new double array where all the values sum to 1. Relative ratios
142         '''         are preserved. </returns>
143         Public Function getNormalized(array As Double()) As Double()
144             Dim newArray As Double() = New Double(array.Length - 1) {}
145             Dim ___total As Double = getTotal(array)
146             For i As Integer = 0 To array.Length - 1
147                 newArray(i) = array(i) / ___total
148             Next i
149             Return newArray
150         End Function
151
152         ''' <param name="array">
153         '''            entries to be summed </param>
154         ''' <param name="start">
155         '''            start position </param>
156         ''' <param name="end">
157         '''            the index of the element after the last one to be included </param>
158         ''' <returns> the total of a the values in a range of an array </returns>
159         Public Function getTotal(array As Double(), start As Integer, [end] As Integer) As Double
160             Dim ___total As Double = 0.0
161             For i As Integer = start To [end] - 1
162                 ___total += array(i)
163             Next i
164             Return ___total
165         End Function
166
167         ''' <param name="array">
168         '''            to sum over </param>
169         ''' <returns> the total of the values in an array </returns>
170         Public Function getTotal(array As Double()) As Double
171             Return getTotal(array, 0, array.Length)
172
173         End Function
174
175         ' ===================== (Synchronized) Static access methods to the private
176         ' random instance ===========
177
178         ''' <summary>
179         ''' Access a default instance of this class, access is synchronized
180         ''' </summary>
181         Public Property Seed As Long
182             Get
183                 SyncLock random
184                     Return random.Seed
185                 End SyncLock
186             End Get
187             Set(seed As Long)
188                 SyncLock random
189                     random.Seed = seed
190                 End SyncLock
191             End Set
192         End Property
193
194         ''' <summary>
195         ''' Access a default instance of this class, access is synchronized
196         ''' </summary>
197         Public Function nextByte() As SByte
198             SyncLock random
199                 Return random.nextByte()
200             End SyncLock
201         End Function
202
203         ''' <summary>
204         ''' Access a default instance of this class, access is synchronized
205         ''' </summary>
206         Public Function nextBoolean() As Boolean
207             SyncLock random
208                 Return random.nextBoolean()
209             End SyncLock
210         End Function
211
212         ''' <summary>
213         ''' Access a default instance of this class, access is synchronized
214         ''' </summary>
215         Public Sub nextBytes(bs As SByte())
216             SyncLock random
217                 random.nextBytes(bs)
218             End SyncLock
219         End Sub
220
221         ''' <summary>
222         ''' Access a default instance of this class, access is synchronized
223         ''' </summary>
224         Public Function nextChar() As Char
225             SyncLock random
226                 Return random.nextChar()
227             End SyncLock
228         End Function
229
230         ''' <summary>
231         ''' Access a default instance of this class, access is synchronized
232         ''' </summary>
233         Public Function nextGaussian() As Double
234             SyncLock random
235                 Return random.nextGaussian()
236             End SyncLock
237         End Function
238
239         ' Mean = alpha / lambda
240         ' Variance = alpha / (lambda*lambda)
241
242         Public Function nextGamma(alpha As Double, lambda As Double) As Double
243             SyncLock random
244                 Return random.nextGamma(alpha, lambda)
245             End SyncLock
246         End Function
247
248         ''' <summary>
249         ''' Access a default instance of this class, access is synchronized
250         ''' </summary>
251         ''' <returns> a pseudo random double precision floating point number in [01) </returns>
252         Public Function nextDouble() As Double
253             SyncLock random
254                 Return random.nextDouble()
255             End SyncLock
256         End Function
257
258         ''' <returns> log of random variable in [0,1] </returns>
259         Public Function randomLogDouble() As Double
260             Return sys.Log(nextDouble())
261         End Function
262
263         ''' <summary>
264         ''' Access a default instance of this class, access is synchronized
265         ''' </summary>
266         Public Function nextExponential(lambda As Double) As Double
267             SyncLock random
268                 Return -1.0 * sys.Log(1 - random.nextDouble()) / lambda
269             End SyncLock
270         End Function
271
272         ''' <summary>
273         ''' Access a default instance of this class, access is synchronized
274         ''' </summary>
275         Public Function nextInverseGaussian(mu As Double, lambda As Double) As Double
276             SyncLock random
277                 '
278                 '  * CODE TAKEN FROM WIKIPEDIA. TESTING DONE WITH RESULTS GENERATED IN
279                 '  * R AND LOOK COMPARABLE
280                 '  
281                 Dim v As Double = random.nextGaussian() ' sample from a normal
282                 ' distribution with a mean of 0
283                 ' and 1 standard deviation
284                 Dim y As Double = v * v
285                 Dim x As Double = mu + (mu * mu * y) / (2 * lambda) - (mu / (2 * lambda)) * sys.Sqrt(4 * mu * lambda * y + mu * mu * y * y)
286                 Dim test As Double = MathUtils.nextDouble() ' sample from a uniform
287                 ' distribution between 0
288                 ' and 1
289                 If test <= (mu) / (mu + x) Then
290                     Return x
291                 Else
292                     Return (mu * mu) / x
293                 End If
294             End SyncLock
295         End Function
296
297         ''' <summary>
298         ''' Access a default instance of this class, access is synchronized
299         ''' </summary>
300         Public Function nextFloat() As Single
301             SyncLock random
302                 Return random.nextFloat()
303             End SyncLock
304         End Function
305
306         ''' <summary>
307         ''' Access a default instance of this class, access is synchronized
308         ''' </summary>
309         Public Function nextLong() As Long
310             SyncLock random
311                 Return random.nextLong()
312             End SyncLock
313         End Function
314
315         ''' <summary>
316         ''' Access a default instance of this class, access is synchronized
317         ''' </summary>
318         Public Function nextShort() As Short
319             SyncLock random
320                 Return random.nextShort()
321             End SyncLock
322         End Function
323
324         ''' <summary>
325         ''' Access a default instance of this class, access is synchronized
326         ''' </summary>
327         Public Function nextInt() As Integer
328             SyncLock random
329                 Return random.nextInt()
330             End SyncLock
331         End Function
332
333         ''' <summary>
334         ''' Access a default instance of this class, access is synchronized
335         ''' </summary>
336         Public Function nextInt(n As Integer) As Integer
337             SyncLock random
338                 Return random.nextInt(n)
339             End SyncLock
340         End Function
341
342         ''' 
343         ''' <param name="low"> </param>
344         ''' <param name="high"> </param>
345         ''' <returns> uniform between low and high </returns>
346         Public Function uniform(low As Double, high As Double) As Double
347             Return low + nextDouble() * (high - low)
348         End Function
349
350         ''' <summary>
351         ''' Shuffles an array.
352         ''' </summary>
353         Public Sub shuffle(array As Integer())
354             SyncLock random
355                 random.shuffle(array)
356             End SyncLock
357         End Sub
358
359         ''' <summary>
360         ''' Shuffles an array. Shuffles numberOfShuffles times
361         ''' </summary>
362         Public Sub shuffle(array As Integer(), numberOfShuffles As Integer)
363             SyncLock random
364                 random.shuffle(array, numberOfShuffles)
365             End SyncLock
366         End Sub
367
368         ''' <summary>
369         ''' Returns an array of shuffled indices of length l.
370         ''' </summary>
371         ''' <param name="l">
372         '''            length of the array required. </param>
373         Public Function shuffled(l As Integer) As Integer()
374             SyncLock random
375                 Return random.shuffled(l)
376             End SyncLock
377         End Function
378
379         Public Function sampleIndicesWithReplacement(length As Integer) As Integer()
380             SyncLock random
381                 Dim result As Integer() = New Integer(length - 1) {}
382                 For i As Integer = 0 To length - 1
383                     result(i) = random.nextInt(length)
384                 Next i
385                 Return result
386             End SyncLock
387         End Function
388
389         ''' <summary>
390         ''' Permutes an array.
391         ''' </summary>
392         Public Sub permute(array As Integer())
393             SyncLock random
394                 random.permute(array)
395             End SyncLock
396         End Sub
397
398         ''' <summary>
399         ''' Returns a uniform random permutation of 0,...,l-1
400         ''' </summary>
401         ''' <param name="l">
402         '''            length of the array required. </param>
403         Public Function permuted(l As Integer) As Integer()
404             SyncLock random
405                 Return random.permuted(l)
406             End SyncLock
407         End Function
408
409         'Public Function logHyperSphereVolume(dimension As Integer, radius As Double) As Double
410         '    Return dimension * (0.57236494292470008 + sys.Log(radius)) + -jebl.math.GammaFunction.lnGamma(dimension / 2.0 + 1.0)
411         'End Function
412
413         ''' <summary>
414         ''' Returns sqrt(a^2 + b^2) without under/overflow.
415         ''' </summary>
416         Public Function hypot(a As Double, b As Double) As Double
417             Dim r As Double
418             If sys.Abs(a) > sys.Abs(b) Then
419                 r = b / a
420                 r = sys.Abs(a) * sys.Sqrt(1 + r * r)
421             ElseIf b <> 0 Then
422                 r = a / b
423                 r = sys.Abs(b) * sys.Sqrt(1 + r * r)
424             Else
425                 r = 0.0
426             End If
427             Return r
428         End Function
429     End Module
430 End Namespace