1	#Region "Microsoft.VisualBasic::b1b479ece4f79a113f702c4ab2b1db29, 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 Overloads) GetPearson, JaccardIndex, kendallTauBeta, KLD
39	'                   rankKendallTauBeta, SW
40	'         Structure Pearson
41	'
42	'             Properties: P
43	'
44	'             Function: Measure, 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 String, Double))
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 Object, Object, Boolean) = Nothing) As 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	''' <summary>
386	''' 默认使用Pearson相似度
387	''' </summary>
388	''' <returns></returns>
389	Public ReadOnly Property PearsonDefault As DefaultValue(Of ICorrelation) = New ICorrelation(AddressOf GetPearson)
390
391	''' <summary>
392	''' Pearson correlations
393	''' </summary>
394	''' <param name="x#"></param>
395	''' <param name="y#"></param>
396	''' <returns></returns>
397	<ExportAPI("Pearson")> Public Function GetPearson(x#(), y#()) As Double
398	Dim j As Integer, n As Integer = x.Length
399	Dim yt As Double, xt As Double
400	Dim syy As Double = 0.0, sxy As Double = 0.0, sxx As Double = 0.0
401	Dim ay As Double = 0.0, ax As Double = 0.0
402
403	For j = 0 To n - 1
404	' finds the mean
405	ax += x(j)
406	ay += y(j)
407	Next
408
409	ax /= n
410	ay /= n
411
412	For j = 0 To n - 1
413	' compute correlation coefficient
414	xt = x(j) - ax
415	yt = y(j) - ay
416	sxx += xt * xt
417	syy += yt * yt
418	sxy += xt * yt
419	Next
420
421	Return sxy / (Sqrt(sxx * syy) + TINY)
422	End Function
423
424	''' <summary>
425	''' 相关性的计算分析函数
426	''' </summary>
427	''' <param name="X"></param>
428	''' <param name="Y"></param>
429	''' <returns></returns>
430	Public Delegate Function ICorrelation(X#(), Y#()) As Double
431
432	Const VectorSizeMustAgree$ = "[X:={0}, Y:={1}] The vector length betwen the two samples is not agreed!!!"
433
434	<MethodImpl(MethodImplOptions.AggressiveInlining)>
435	Private Sub throwNotAgree(x#(), y#())
436	Dim message$ = String.Format(VectorSizeMustAgree, x.Length, y.Length)
437	Throw New DataException(message)
438	End Sub
439
440	''' <summary>
441	''' This method should not be used in cases where the data set is truncated; that is,
442	''' when the Spearman correlation coefficient is desired for the top X records
443	''' (whether by pre-change rank or post-change rank, or both), the user should use the
444	''' Pearson correlation coefficient formula given above.
445	''' (斯皮尔曼相关性)
446	''' </summary>
447	''' <param name="X"></param>
448	''' <param name="Y"></param>
449	''' <returns></returns>
450	''' <remarks>
451	''' https://en.wikipedia.org/wiki/Spearman%27s_rank_correlation_coefficient
452	''' checked!
453	''' </remarks>
454	'''
455	<ExportAPI("Spearman",
456	Info:="This method should not be used in cases where the data set is truncated;
457	that is, when the Spearman correlation coefficient is desired for the top X records
458	(whether by pre-change rank or post-change rank, or both), the user should use the
459	Pearson correlation coefficient formula given above.")>
460	Public Function Spearman#(X#(), Y#())
461	If X.Length <> Y.Length Then
462	Call throwNotAgree(X, Y)
463	ElseIf X.Length = 1 Then
464	Throw New DataException(UnableMeasures)
465	End If
466
467	' size n
468	Dim n As Integer = X.Length
469	Dim Xx As spcc() = __getOrder(X)
470	Dim Yy As spcc() = __getOrder(Y)
471
472	Dim deltaSum# = Aggregate i As Integer
473	In n.Sequence
474	Into Sum((Xx(i).rank - Yy(i).rank) ^ 2)
475	Dim spcc = 1 - 6 * deltaSum / (n ^ 3 - n)
476
477	Return spcc
478	End Function
479
480	Const UnableMeasures$ = "Samples number just equals 1, the function unable to measure the correlation!!!"
481
482	Private Function __getOrder(samples#()) As spcc()
483	Dim dat = (From i As Integer
484	In samples.Sequence
485	Select spcc = New spcc.__spccInner With {  ' 原有的顺序
486	.i = i,
487	.val = samples(i)
488	}
489	Order By spcc.val Ascending).ToArray  ' 从小到大排序
490	Dim buf = (From p As Integer  ' rank
491	In dat.Sequence
492	Select spcc = New spcc With {
493	.rank = p,
494	.data = dat(p)
495	}
496	Group spcc By spcc.data.val Into Group) _
497	.ToDictionary(Function(x) x.val,
498	Function(x) x.Group.ToArray)
499
500	Dim rankList As New List(Of spcc)
501
502	For Each item As spcc() In buf.Values
503	If item.Length = 1 Then
504	Call rankList.Add(item(Scan0))
505	Else
506	Dim rank As Double = item.Select(Function(x) x.rank).Average
507	Dim array As spcc() = item.Select(
508	Function(x) New spcc With {
509	.rank = rank,
510	.data = x.data
511	}).ToArray
512	Call rankList.AddRange(array)
513	End If
514	Next
515
516	' 重新按照原有的顺序返回
517	Return (From x As spcc
518	In rankList
519	Select x
520	Order By x.data.i Ascending).ToArray
521	End Function
522
523	''' <summary>
524	''' 计算所需要的临时变量类型
525	''' </summary>
526	Private Structure spcc
527	''' <summary>
528	''' 排序之后得到的位置
529	''' </summary>
530	Public rank As Double
531	''' <summary>
532	''' 原始数据
533	''' </summary>
534	Public data As __spccInner
535
536	Public Structure __spccInner
537	''' <summary>
538	''' 在序列之中原有的位置
539	''' </summary>
540	Public i As Integer
541	Public val As Double
542	End Structure
543	End Structure
544
545	''' <summary>
546	''' 输入的数据为一个对象属性的集合，默认的<paramref name="compute"/>计算方法为<see cref="GetPearson"/>
547	''' </summary>
548	''' <param name="data">``[ID, properties]``</param>
549	''' <param name="compute">
550	''' Using pearson method as default if this parameter is nothing.
551	''' (默认的计算形式为<see cref="GetPearson"/>)
552	''' </param>
553	''' <returns></returns>
554	<Extension>
555	Public Function CorrelationMatrix(data As IEnumerable(Of Vector), Optional compute As ICorrelation = Nothing) As DataSet()
556	Dim array As Vector() = data.ToArray
557	Dim outMatrix As New List(Of DataSet)
558
559	compute = compute Or PearsonDefault
560
561	For Each a As Vector In array
562	Dim ca As New Dictionary(Of String, Double)
563	Dim v#() = a.Value
564
565	For Each b In array
566	ca(b.Name) = compute(v, b.Value)
567	Next
568
569	outMatrix += New DataSet With {
570	.Name = a.Name,
571	.Value = ca
572	}
573	Next
574
575	Return outMatrix
576	End Function
577	End Module
578
579	Public Module Beta
580
581	Const SWITCH As Integer = 3000, MAXIT As Integer = 1000
582	Const EPS As Double = 0.0000003, FPMIN As Double = 1.0E-30
583
584	Public Function betai(a As Double, b As Double, x As Double) As Double
585	Dim bt As Double
586
587	If x < 0.0 OrElse x > 1.0 Then
588	Throw New ArgumentException($"Bad x:={x} in routine betai")
589	End If
590	If x = 0.0 OrElse x = 1.0 Then
591	bt = 0.0
592	Else
593	bt = sys.Exp(gammln(a + b) - gammln(a) - gammln(b) + a * sys.Log(x) + b * sys.Log(1.0 - x))
594	End If
595	If x < (a + 1.0) / (a + b + 2.0) Then
596	Return bt * betacf(a, b, x) / a
597	Else
598	Return 1.0 - bt * betacf(b, a, 1.0 - x) / b
599	End If
600	End Function
601
602	Private Function gammln(xx As Double) As Double
603	Dim x As Double = xx, y As Double, tmp As Double, ser As Double
604	Dim j As Integer
605
606	y = x
607	tmp = x + 5.5
608	tmp -= (x + 0.5) * sys.Log(tmp)
609	ser = 1.00000000019001
610	For j = 0 To 5
611	y += 1 ' ++y
612	ser += cof(j) / y
613	Next
614
615	Return -tmp + sys.Log(2.506628274631 * ser / x)
616	End Function
617
618	ReadOnly cof As Double() = {
619	76.1800917294715,
620	-86.5053203294168,
621	24.0140982408309,
622	-1.23173957245015,
623	0.00120865097386618,
624	-0.000005395239384953
625	}
626
627	Private Function betacf(a As Double, b As Double, x As Double) As Double
628	Dim m As Integer, m2 As Integer
629	Dim aa As Double,
630	c As Double,
631	d As Double,
632	del As Double,
633	h As Double,
634	qab As Double,
635	qam As Double,
636	qap As Double
637
638	qab = a + b
639	qap = a + 1.0
640	qam = a - 1.0
641	c = 1.0
642	d = 1.0 - qab * x / qap
643	If sys.Abs(d) < FPMIN Then
644	d = FPMIN
645	End If
646	d = 1.0 / d
647	h = d
648
649	For m = 1 To MAXIT
650	m2 = 2 * m
651	aa = m * (b - m) * x / ((qam + m2) * (a + m2))
652	d = 1.0 + aa * d
653	If sys.Abs(d) < FPMIN Then
654	d = FPMIN
655	End If
656	c = 1.0 + aa / c
657	If sys.Abs(c) < FPMIN Then
658	c = FPMIN
659	End If
660	d = 1.0 / d
661	h *= d * c
662	aa = -(a + m) * (qab + m) * x / ((a + m2) * (qap + m2))
663	d = 1.0 + aa * d
664	If sys.Abs(d) < FPMIN Then
665	d = FPMIN
666	End If
667	c = 1.0 + aa / c
668	If sys.Abs(c) < FPMIN Then
669	c = FPMIN
670	End If
671	d = 1.0 / d
672	del = d * c
673	h *= del
674	If sys.Abs(del - 1.0) < EPS Then
675	Exit For
676	End If
677	Next
678	If m > MAXIT Then
679	Dim msg As String =
680	$"a:={a} or b:={b} too big, or MAXIT too small in betacf"
681	Throw New ArgumentException(msg)
682	End If
683	Return h
684	End Function
685
686	Public Function erfcc(x As Double) As Double
687	Dim t As Double, z As Double, ans As Double
688
689	z = sys.Abs(x)
690	t = 1.0 / (1.0 + 0.5 * z)
691	ans = t * sys.Exp(-z * z - 1.26551223 +
692	t * (1.00002368 +
693	t * (0.37409196 +
694	t * (0.09678418 +
695	t * (-0.18628806 +
696	t * (0.27886807 +
697	t * (-1.13520398 +
698	t * (1.48851587 +
699	t * (-0.82215223 +
700	t * 0.17087277)))))))))
701	Return If(x >= 0.0, ans, 2.0 - ans)
702	End Function
703	End Module
704	End Namespace