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