1 #Region "Microsoft.VisualBasic::eb5faf8917db4d5f0508d8f368988149, Microsoft.VisualBasic.Core\Extensions\Math\Correlations\Correlations.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 Correlations
35     
36     '         Properties: PearsonDefault
37     
38     '         Function: __kldPart, (+2 OverloadsGetPearson, JaccardIndex, kendallTauBeta, KLD
39     '                   rankKendallTauBeta, SW
40     '         Structure Pearson
41     
42     '             Properties: P
43     
44     '             FunctionMeasure, RankPearson, ToString
45     
46     '         Delegate Function
47     
48     '             Function: __getOrder, CorrelationMatrix, Spearman
49     
50     '             Sub: throwNotAgree
51     '         Structure spcc
52     
53     
54     '             Structure __spccInner
55     
56     
57     
58     
59     
60     
61     
62     '     Module Beta
63     
64     '         Function: betacf, betai, erfcc, gammln
65     
66     
67     
68     ' /********************************************************************************/
69
70 #End Region
71
72 Imports System.Runtime.CompilerServices
73 Imports Microsoft.VisualBasic.CommandLine.Reflection
74 Imports Microsoft.VisualBasic.ComponentModel.DataStructures
75 Imports Microsoft.VisualBasic.Language
76 Imports Microsoft.VisualBasic.Language.Default
77 Imports Microsoft.VisualBasic.Linq.Extensions
78 Imports Microsoft.VisualBasic.Scripting.MetaData
79 Imports DataSet = Microsoft.VisualBasic.ComponentModel.DataSourceModel.NamedValue(Of System.Collections.Generic.Dictionary(Of StringDouble))
80 Imports sys = System.Math
81 Imports Vector = Microsoft.VisualBasic.ComponentModel.DataSourceModel.NamedValue(Of Double())
82
83 Namespace Math.Correlations
84
85     ''' <summary>
86     ''' 计算两个数据向量之间的相关度的大小
87     ''' </summary>
88     <Package("Correlations", Category:=APICategories.UtilityTools, Publisher:="amethyst.asuka@gcmodeller.org")>
89     Public Module Correlations
90
91         ''' <summary>
92         ''' The Jaccard index, also known as Intersection over Union and the Jaccard similarity coefficient 
93         ''' (originally coined coefficient de communauté by Paul Jaccard), is a statistic used for comparing 
94         ''' the similarity and diversity of sample sets. The Jaccard coefficient measures similarity between 
95         ''' finite sample sets, and is defined as the size of the intersection divided by the size of the 
96         ''' union of the sample sets.
97         ''' 
98         ''' https://en.wikipedia.org/wiki/Jaccard_index
99         ''' </summary>
100         ''' <typeparam name="T"></typeparam>
101         ''' <param name="a"></param>
102         ''' <param name="b"></param>
103         ''' <param name="equal"></param>
104         ''' <returns></returns>
105         Public Function JaccardIndex(Of T)(a As IEnumerable(Of T), b As IEnumerable(Of T), Optional equal As Func(Of ObjectObjectBoolean) = NothingAs Double
106             If equal Is Nothing Then
107                 equal = Function(x, y) x.Equals(y)
108             End If
109
110             Dim setA As New [Set](a, equal)
111             Dim setB As New [Set](b, equal)
112
113             ' (If A and B are both empty, we define J(A,B) = 1.)
114             If setA.Length = setB.Length AndAlso setA.Length = 0 Then
115                 Return 1
116             End If
117
118             Dim intersects = setA And setB
119             Dim union = setA + setB
120             Dim similarity# = intersects.Length / union.Length
121
122             Return similarity
123         End Function
124
125         ''' <summary>
126         ''' Sandelin-Wasserman similarity function.(假若所有的元素都是0-1之间的话,结果除以2可以得到相似度)
127         ''' </summary>
128         ''' <param name="x"></param>
129         ''' <param name="y"></param>
130         ''' <returns></returns>
131         <ExportAPI("SW"Info:="Sandelin-Wasserman similarity function")>
132         Public Function SW(x As Double(), y As Double()) As Double
133             Dim p = From i As Integer
134                     In x.Sequence
135                     Select x(i) - y(i)
136             Dim s# = Aggregate n As Double In p Into Sum(n ^ 2)
137             s = 2 - s
138             Return s
139         End Function
140
141         ''' <summary>
142         ''' Kullback-Leibler divergence
143         ''' </summary>
144         ''' <param name="x"></param>
145         ''' <param name="y"></param>
146         ''' <returns></returns>
147         <ExportAPI("KLD"Info:="Kullback-Leibler divergence")>
148         Public Function KLD(x As Double(), y As Double()) As Double
149             Dim index As Integer() = x.Sequence
150             Dim a As Double = (From i As Integer In index Select __kldPart(x(i), y(i))).Sum
151             Dim b As Double = (From i As Integer In index Select __kldPart(y(i), x(i))).Sum
152             Dim value As Double = (a + b) / 2
153             Return value
154         End Function
155
156         Private Function __kldPart(Xa#, Ya#) As Double
157             If Xa = 0R Then
158                 Return 0R
159             End If
160             Dim value As Double = Xa * sys.Log(Xa / Ya)  ' 0 * n = 0
161             Return value
162         End Function
163
164 #Region "https://en.wikipedia.org/wiki/Kendall_tau_distance"
165
166         ''' <summary>
167         ''' Provides rank correlation coefficient metrics Kendall tau
168         ''' </summary>
169         ''' <param name="x"></param>
170         ''' <param name="y"></param>
171         ''' <returns></returns>
172         ''' <remarks>
173         ''' https://github.com/felipebravom/RankCorrelation
174         ''' </remarks>
175         Public Function rankKendallTauBeta(x As Double(), y As Double()) As Double
176             Debug.Assert(x.Length = y.Length)
177             Dim x_n As Integer = x.Length
178             Dim y_n As Integer = y.Length
179             Dim x_rank As Double() = New Double(x_n - 1) {}
180             Dim y_rank As Double() = New Double(y_n - 1) {}
181             Dim sorted As New SortedDictionary(Of Double?, HashSet(Of Integer?))()
182
183             For i As Integer = 0 To x_n - 1
184                 Dim v As Double = x(i)
185                 If sorted.ContainsKey(v) = False Then
186                     sorted(v) = New HashSet(Of Integer?)()
187                 End If
188                 sorted(v).Add(i)
189             Next
190
191             Dim c As Integer = 1
192             For Each v As Double In sorted.Keys.OrderByDescending(Function(k) k)
193                 Dim r As Double = 0
194                 For Each i As Integer In sorted(v)
195                     r += c
196                     c += 1
197                 Next
198
199                 r /= sorted(v).Count
200
201                 For Each i As Integer In sorted(v)
202                     x_rank(i) = r
203                 Next
204             Next
205
206             sorted.Clear()
207             For i As Integer = 0 To y_n - 1
208                 Dim v As Double = y(i)
209                 If sorted.ContainsKey(v) = False Then
210                     sorted(v) = New HashSet(Of Integer?)()
211                 End If
212                 sorted(v).Add(i)
213             Next
214
215             c = 1
216             For Each v As Double In sorted.Keys.OrderByDescending(Function(k) k)
217                 Dim r As Double = 0
218                 For Each i As Integer In sorted(v)
219                     r += c
220                     c += 1
221                 Next
222
223                 r /= (sorted(v).Count)
224
225                 For Each i As Integer In sorted(v)
226                     y_rank(i) = r
227                 Next
228             Next
229
230             Return kendallTauBeta(x_rank, y_rank)
231         End Function
232
233         ''' <summary>
234         ''' Provides rank correlation coefficient metrics Kendall tau
235         ''' </summary>
236         ''' <param name="x"></param>
237         ''' <param name="y"></param>
238         ''' <returns></returns>
239         ''' <remarks>
240         ''' https://github.com/felipebravom/RankCorrelation
241         ''' </remarks>
242         Public Function kendallTauBeta(x As Double(), y As Double()) As Double
243             Debug.Assert(x.Length = y.Length)
244
245             Dim c As Integer = 0
246             Dim d As Integer = 0
247             Dim xTies As New Dictionary(Of Double?, HashSet(Of Integer?))()
248             Dim yTies As New Dictionary(Of Double?, HashSet(Of Integer?))()
249
250             For i As Integer = 0 To x.Length - 2
251                 For j As Integer = i + 1 To x.Length - 1
252                     If x(i) > x(j) AndAlso y(i) > y(j) Then
253                         c += 1
254                     ElseIf x(i) < x(j) AndAlso y(i) < y(j) Then
255                         c += 1
256                     ElseIf x(i) > x(j) AndAlso y(i) < y(j) Then
257                         d += 1
258                     ElseIf x(i) < x(j) AndAlso y(i) > y(j) Then
259                         d += 1
260                     Else
261                         If x(i) = x(j) Then
262                             If xTies.ContainsKey(x(i)) = False Then
263                                 xTies(x(i)) = New HashSet(Of Integer?)()
264                             End If
265                             xTies(x(i)).Add(i)
266                             xTies(x(i)).Add(j)
267                         End If
268
269                         If y(i) = y(j) Then
270                             If yTies.ContainsKey(y(i)) = False Then
271                                 yTies(y(i)) = New HashSet(Of Integer?)()
272                             End If
273                             yTies(y(i)).Add(i)
274                             yTies(y(i)).Add(j)
275                         End If
276                     End If
277                 Next
278             Next
279
280             Dim diff As Integer = c - d
281             Dim denom As Double = 0
282
283             Dim n0 As Double = (x.Length * (x.Length - 1)) / 2.0
284             Dim n1 As Double = 0
285             Dim n2 As Double = 0
286
287             For Each t As Double In xTies.Keys
288                 Dim s As Double = xTies(t).Count
289                 n1 += (s * (s - 1)) / 2
290             Next
291
292             For Each t As Double In yTies.Keys
293                 Dim s As Double = yTies(t).Count
294                 n2 += (s * (s - 1)) / 2
295             Next
296
297             denom = sys.Sqrt((n0 - n1) * (n0 - n2))
298
299             If denom = 0 Then
300                 denom += 0.000000001
301             End If
302
303             Dim td As Double = diff / (denom) ' 0.000..1 added on 11/02/2013 fixing NaN error
304
305             Debug.Assert(td >= -1 AndAlso td <= 1, td)
306
307             Return td
308         End Function
309 #End Region
310
311         ''' <summary>
312         ''' will regularize the unusual case of complete correlation
313         ''' </summary>
314         Const TINY As Double = 1.0E-20
315
316         ''' <summary>
317         '''
318         ''' </summary>
319         ''' <param name="x"></param>
320         ''' <param name="y"></param>
321         ''' <param name="prob"></param>
322         ''' <param name="prob2"></param>
323         ''' <param name="z">fisher's z trasnformation</param>
324         ''' <returns></returns>
325         ''' <remarks>
326         ''' checked by Excel
327         ''' </remarks>
328         <ExportAPI("Pearson")>
329         Public Function GetPearson(x#(), y#(), Optional ByRef prob# = 0, Optional ByRef prob2# = 0, Optional ByRef z# = 0) As Double
330             Dim t#, df#
331             Dim pcc As Double = GetPearson(x, y)
332             Dim n As Integer = x.Length
333
334             ' fisher's z trasnformation
335             z = 0.5 * sys.Log((1.0 + pcc + TINY) / (1.0 - pcc + TINY))
336
337             ' student's t probability
338             df = n - 2
339             t = pcc * sys.Sqrt(df / ((1.0 - pcc + TINY) * (1.0 + pcc + TINY)))
340
341             prob = Beta.betai(0.5 * df, 0.5, df / (df + t * t))
342             ' for a large n
343             prob2 = Beta.erfcc(Abs(z * sys.Sqrt(n - 1.0)) / 1.4142136)
344
345             Return pcc
346         End Function
347
348         Public Structure Pearson
349             Dim pearson#
350             Dim pvalue#
351             Dim pvalue2#
352             Dim Z#
353
354             Public ReadOnly Property P As Double
355                 <MethodImpl(MethodImplOptions.AggressiveInlining)>
356                 Get
357                     Return -Math.Log10(pvalue)
358                 End Get
359             End Property
360
361             Public Overrides Function ToString() As String
362                 Return $"{pearson} @ {pvalue}"
363             End Function
364
365             Public Shared Function Measure(x As IEnumerable(Of Double), y As IEnumerable(Of Double)) As Pearson
366                 Dim pvalue1, pvalue2, z As Double
367                 Dim pearson# = GetPearson(x.ToArray, y.ToArray, pvalue1, pvalue2, z)
368
369                 Return New Pearson With {
370                     .pearson = pearson,
371                     .pvalue = pvalue1,
372                     .pvalue2 = pvalue2,
373                     .Z = z
374                 }
375             End Function
376
377             Public Shared Function RankPearson(x As IEnumerable(Of Double), y As IEnumerable(Of Double)) As Pearson
378                 Dim r1 = x.Ranking(Strategies.FractionalRanking)
379                 Dim r2 = y.Ranking(Strategies.FractionalRanking)
380
381                 Return Measure(r1, r2)
382             End Function
383         End Structure
384
385         Public ReadOnly Property PearsonDefault As DefaultValue(Of ICorrelation) = New ICorrelation(AddressOf GetPearson).AsDefault
386
387         ''' <summary>
388         ''' Pearson correlations
389         ''' </summary>
390         ''' <param name="x#"></param>
391         ''' <param name="y#"></param>
392         ''' <returns></returns>
393         <ExportAPI("Pearson")> Public Function GetPearson(x#(), y#()) As Double
394             Dim j As Integer, n As Integer = x.Length
395             Dim yt As Double, xt As Double
396             Dim syy As Double = 0.0, sxy As Double = 0.0, sxx As Double = 0.0
397             Dim ay As Double = 0.0, ax As Double = 0.0
398
399             For j = 0 To n - 1
400                 ' finds the mean
401                 ax += x(j)
402                 ay += y(j)
403             Next
404
405             ax /= n
406             ay /= n
407
408             For j = 0 To n - 1
409                 ' compute correlation coefficient
410                 xt = x(j) - ax
411                 yt = y(j) - ay
412                 sxx += xt * xt
413                 syy += yt * yt
414                 sxy += xt * yt
415             Next
416
417             Return sxy / (Sqrt(sxx * syy) + TINY)
418         End Function
419
420         ''' <summary>
421         ''' 相关性的计算分析函数
422         ''' </summary>
423         ''' <param name="X"></param>
424         ''' <param name="Y"></param>
425         ''' <returns></returns>
426         Public Delegate Function ICorrelation(X#(), Y#()) As Double
427
428         Const VectorSizeMustAgree$ = "[X:={0}, Y:={1}] The vector length betwen the two samples is not agreed!!!"
429
430         <MethodImpl(MethodImplOptions.AggressiveInlining)>
431         Private Sub throwNotAgree(x#(), y#())
432             Dim message$ = String.Format(VectorSizeMustAgree, x.Length, y.Length)
433             Throw New DataException(message)
434         End Sub
435
436         ''' <summary>
437         ''' This method should not be used in cases where the data set is truncated; that is,
438         ''' when the Spearman correlation coefficient is desired for the top X records
439         ''' (whether by pre-change rank or post-change rank, or both), the user should use the
440         ''' Pearson correlation coefficient formula given above.
441         ''' (斯皮尔曼相关性)
442         ''' </summary>
443         ''' <param name="X"></param>
444         ''' <param name="Y"></param>
445         ''' <returns></returns>
446         ''' <remarks>
447         ''' https://en.wikipedia.org/wiki/Spearman%27s_rank_correlation_coefficient
448         ''' checked!
449         ''' </remarks>
450         '''
451         <ExportAPI("Spearman",
452                Info:="This method should not be used in cases where the data set is truncated; 
453                that is, when the Spearman correlation coefficient is desired for the top X records 
454                (whether by pre-change rank or post-change rank, or both), the user should use the 
455                Pearson correlation coefficient formula given above.")>
456         Public Function Spearman#(X#(), Y#())
457             If X.Length <> Y.Length Then
458                 Call throwNotAgree(X, Y)
459             ElseIf X.Length = 1 Then
460                 Throw New DataException(UnableMeasures)
461             End If
462
463             ' size n
464             Dim n As Integer = X.Length
465             Dim Xx As spcc() = __getOrder(X)
466             Dim Yy As spcc() = __getOrder(Y)
467
468             Dim deltaSum# = Aggregate i As Integer
469                             In n.Sequence
470                             Into Sum((Xx(i).rank - Yy(i).rank) ^ 2)
471             Dim spcc = 1 - 6 * deltaSum / (n ^ 3 - n)
472
473             Return spcc
474         End Function
475
476         Const UnableMeasures$ = "Samples number just equals 1, the function unable to measure the correlation!!!"
477
478         Private Function __getOrder(samples#()) As spcc()
479             Dim dat = (From i As Integer
480                        In samples.Sequence
481                        Select spcc = New spcc.__spccInner With {  ' 原有的顺序
482                            .i = i,
483                            .val = samples(i)
484                        }
485                        Order By spcc.val Ascending).ToArray  ' 从小到大排序
486             Dim buf = (From p As Integer  ' rank
487                        In dat.Sequence
488                        Select spcc = New spcc With {
489                            .rank = p,
490                            .data = dat(p)
491                        }
492                        Group spcc By spcc.data.val Into Group) _
493                         .ToDictionary(Function(x) x.val,
494                                       Function(x) x.Group.ToArray)
495
496             Dim rankList As New List(Of spcc)
497
498             For Each item As spcc() In buf.Values
499                 If item.Length = 1 Then
500                     Call rankList.Add(item(Scan0))
501                 Else
502                     Dim rank As Double = item.Select(Function(x) x.rank).Average
503                     Dim array As spcc() = item.Select(
504                         Function(x) New spcc With {
505                             .rank = rank,
506                             .data = x.data
507                         }).ToArray
508                     Call rankList.AddRange(array)
509                 End If
510             Next
511
512             ' 重新按照原有的顺序返回
513             Return (From x As spcc
514                     In rankList
515                     Select x
516                     Order By x.data.i Ascending).ToArray
517         End Function
518
519         ''' <summary>
520         ''' 计算所需要的临时变量类型
521         ''' </summary>
522         Private Structure spcc
523             ''' <summary>
524             ''' 排序之后得到的位置
525             ''' </summary>
526             Public rank As Double
527             ''' <summary>
528             ''' 原始数据
529             ''' </summary>
530             Public data As __spccInner
531
532             Public Structure __spccInner
533                 ''' <summary>
534                 ''' 在序列之中原有的位置
535                 ''' </summary>
536                 Public i As Integer
537                 Public val As Double
538             End Structure
539         End Structure
540
541         ''' <summary>
542         ''' 输入的数据为一个对象属性的集合,默认的<paramref name="compute"/>计算方法为<see cref="GetPearson"/>
543         ''' </summary>
544         ''' <param name="data">``[ID, properties]``</param>
545         ''' <param name="compute">
546         ''' Using pearson method as default if this parameter is nothing.
547         ''' (默认的计算形式为<see cref="GetPearson"/>)
548         ''' </param>
549         ''' <returns></returns>
550         <Extension>
551         Public Function CorrelationMatrix(data As IEnumerable(Of Vector), Optional compute As ICorrelation = NothingAs DataSet()
552             Dim array As Vector() = data.ToArray
553             Dim outMatrix As New List(Of DataSet)
554
555             compute = compute Or PearsonDefault
556
557             For Each a As Vector In array
558                 Dim ca As New Dictionary(Of StringDouble)
559                 Dim v#() = a.Value
560
561                 For Each b In array
562                     ca(b.Name) = compute(v, b.Value)
563                 Next
564
565                 outMatrix += New DataSet With {
566                     .Name = a.Name,
567                     .Value = ca
568                 }
569             Next
570
571             Return outMatrix
572         End Function
573     End Module
574
575     Public Module Beta
576
577         Const SWITCH As Integer = 3000, MAXIT As Integer = 1000
578         Const EPS As Double = 0.0000003, FPMIN As Double = 1.0E-30
579
580         Public Function betai(a As Double, b As Double, x As DoubleAs Double
581             Dim bt As Double
582
583             If x < 0.0 OrElse x > 1.0 Then
584                 Throw New ArgumentException($"Bad x:={x} in routine betai")
585             End If
586             If x = 0.0 OrElse x = 1.0 Then
587                 bt = 0.0
588             Else
589                 bt = sys.Exp(gammln(a + b) - gammln(a) - gammln(b) + a * sys.Log(x) + b * sys.Log(1.0 - x))
590             End If
591             If x < (a + 1.0) / (a + b + 2.0) Then
592                 Return bt * betacf(a, b, x) / a
593             Else
594                 Return 1.0 - bt * betacf(b, a, 1.0 - x) / b
595             End If
596         End Function
597
598         Private Function gammln(xx As DoubleAs Double
599             Dim x As Double = xx, y As Double, tmp As Double, ser As Double
600             Dim j As Integer
601
602             y = x
603             tmp = x + 5.5
604             tmp -= (x + 0.5) * sys.Log(tmp)
605             ser = 1.00000000019001
606             For j = 0 To 5
607                 y += 1 ' ++y
608                 ser += cof(j) / y
609             Next
610
611             Return -tmp + sys.Log(2.506628274631 * ser / x)
612         End Function
613
614         ReadOnly cof As Double() = {
615             76.1800917294715,
616             -86.5053203294168,
617             24.0140982408309,
618             -1.23173957245015,
619             0.00120865097386618,
620             -0.000005395239384953
621         }
622
623         Private Function betacf(a As Double, b As Double, x As DoubleAs Double
624             Dim m As Integer, m2 As Integer
625             Dim aa As Double,
626             c As Double,
627             d As Double,
628             del As Double,
629             h As Double,
630             qab As Double,
631             qam As Double,
632             qap As Double
633
634             qab = a + b
635             qap = a + 1.0
636             qam = a - 1.0
637             c = 1.0
638             d = 1.0 - qab * x / qap
639             If sys.Abs(d) < FPMIN Then
640                 d = FPMIN
641             End If
642             d = 1.0 / d
643             h = d
644
645             For m = 1 To MAXIT
646                 m2 = 2 * m
647                 aa = m * (b - m) * x / ((qam + m2) * (a + m2))
648                 d = 1.0 + aa * d
649                 If sys.Abs(d) < FPMIN Then
650                     d = FPMIN
651                 End If
652                 c = 1.0 + aa / c
653                 If sys.Abs(c) < FPMIN Then
654                     c = FPMIN
655                 End If
656                 d = 1.0 / d
657                 h *= d * c
658                 aa = -(a + m) * (qab + m) * x / ((a + m2) * (qap + m2))
659                 d = 1.0 + aa * d
660                 If sys.Abs(d) < FPMIN Then
661                     d = FPMIN
662                 End If
663                 c = 1.0 + aa / c
664                 If sys.Abs(c) < FPMIN Then
665                     c = FPMIN
666                 End If
667                 d = 1.0 / d
668                 del = d * c
669                 h *= del
670                 If sys.Abs(del - 1.0) < EPS Then
671                     Exit For
672                 End If
673             Next
674             If m > MAXIT Then
675                 Dim msg As String =
676                 $"a:={a} or b:={b} too big, or MAXIT too small in betacf"
677                 Throw New ArgumentException(msg)
678             End If
679             Return h
680         End Function
681
682         Public Function erfcc(x As DoubleAs Double
683             Dim t As Double, z As Double, ans As Double
684
685             z = sys.Abs(x)
686             t = 1.0 / (1.0 + 0.5 * z)
687             ans = t * sys.Exp(-z * z - 1.26551223 +
688                            t * (1.00002368 +
689                            t * (0.37409196 +
690                            t * (0.09678418 +
691                            t * (-0.18628806 +
692                            t * (0.27886807 +
693                            t * (-1.13520398 +
694                            t * (1.48851587 +
695                            t * (-0.82215223 +
696                            t * 0.17087277)))))))))
697             Return If(x >= 0.0, ans, 2.0 - ans)
698         End Function
699     End Module
700 End Namespace