PureBasic Forums - English

Code: Select all

/*-------------------------------------------------------------------------*/
/* This function uses a simple O(N^2) implementation.  It probably has a
 * smaller constant and therefore is useful in the small N case, and is also
 * useful for testing the relatively complex O(N log N) implementation.
 */

float kendallSmallN( float *arr1, float *arr2, int len )
{
    int m1 = 0, m2 = 0, s = 0, nPair , i,j ;
    float cor ;

/* printf("enter kendallSmallN\n") ; */

    for(i = 0; i < len; i++) {
        for(j = i + 1; j < len; j++) {
            if(arr2[i] > arr2[j]) {
                if (arr1[i] > arr1[j]) {
                    s++;
                } else if(arr1[i] < arr1[j]) {
                    s--;
                } else {
                    m1++;
                }
            } else if(arr2[i] < arr2[j]) {
                if (arr1[i] > arr1[j]) {
                    s--;
                } else if(arr1[i] < arr1[j]) {
                    s++;
                } else {
                    m1++;
                }
            } else {
                m2++;

                if(arr1[i] == arr1[j]) {
                    m1++;
                }
            }
        }
    }

    nPair = len * (len - 1) / 2;

    if( m1 < nPair && m2 < nPair )
      cor = s / ( sqrtf((float)(nPair-m1)) * sqrtf((float)(nPair-m2)) ) ;
    else
      cor = 0.0f ;

    return cor ;
}

Code: Select all

;PB Code
EnableExplicit
;==============KendalTauSlow=======================================
Procedure.d kendallSmallN(  *arr1.double,  *arr2.double,  len.l )
    Protected m1.l = 0, m2.l = 0, s.l = 0, nPair.l , i.l,j.l ,qSOD.q
    Protected cor.d 

    qSOD=SizeOf(double)

    For i = 0 To len-1
        For j = i + 1 To len-1 
            If (*arr2+i*qSOD) > (*arr2+j*qSOD)
                If (*arr1+i*qSOD) > (*arr1+j*qSOD)
                    s=s +1
                ElseIf (*arr1+i*qSOD) < (*arr1+j*qSOD)
                    s=s-1
                Else 
                    m1=m1+1
                EndIf
           ElseIf (*arr2+i*qSOD) < (*arr2+j*qSOD) 
                If (*arr1+i*qSOD) > (*arr1+j*qSOD) 
                    s=s-1
                ElseIf (*arr1+i*qSOD) < (*arr1+j*qSOD) 
                    s=s+1;
                Else 
                    m1=m1+1
                EndIf
            Else
                m2=m2+1

                If (*arr1+i*qSOD) = (*arr1+j*qSOD) 
                    m1=m1+1;
                EndIf
            EndIf
        Next j
    Next i

    nPair = len * (len - 1) / 2

    If( m1 < nPair And m2 < nPair )
      cor = s / ( Sqr(nPair-m1) * Sqr(nPair-m2) ) 
    Else
      cor = 0
    EndIf
    
    ProcedureReturn cor 
      
EndProcedure
;========================Run======================================
OpenConsole()                                 
ConsoleTitle ("PureBasic - Console Example:")                                                  
EnableGraphicalConsole(1)

  Define i.l, n.l, tau.d, taufast.d
  Dim a.d(100)
  Dim b.d(100)
 
  a(0)=1 :a(1)=1:  a(2)=1:  a(3)=2:  a(4)=3:  a(5)=3:  a(6)=4:  a(7)=4 :
  b(0)=1 :b(1)=2:  b(2)=1:  b(3)=3:  b(4)=3:  b(5)=5:  b(6)=5:  b(7)=5
  
  tau=KendallSmallN(a(),b(),8)
  PrintN("KendallSmallN="+StrD(tau)+ " Should be =0.869565")
  Input()
 CloseConsole()
End

Code: Select all

;PB Code
EnableExplicit
;==============KendalTauSlow=======================================
Procedure.d kendallSmallN(  Array arr1.f(1),  Array arr2.f(1),  len.l )
    Protected m1.l = 0, m2.l = 0, s.l = 0, nPair.l , i.l,j.l ,qSOD.q
    Protected cor.f

    qSOD=SizeOf(float)

    For i = 0 To len-1
        For j = i + 1 To len-1 
            If (arr2(i)) > (arr2(j))
                If (arr1(i)) > (arr1(j))
                    s=s +1
                ElseIf (arr1(i)) < (arr1(j))
                    s=s-1
                Else 
                    m1=m1+1
                EndIf
           ElseIf (arr2(i)) < (arr2(j)) 
                If (arr1(i)) > (arr1(j)) 
                    s=s-1
                ElseIf (arr1(i)) < (arr1(j)) 
                    s=s+1;
                Else 
                    m1=m1+1
                EndIf
            Else
                m2=m2+1

                If (arr1(i)) = (arr1(j)) 
                    m1=m1+1;
                EndIf
            EndIf
        Next j
    Next i

    nPair = len * (len - 1) / 2

    If( m1 < nPair And m2 < nPair )
      cor = s / ( Sqr(nPair-m1) * Sqr(nPair-m2) ) 
    Else
      cor = 0
    EndIf
    
    ProcedureReturn cor 
      
EndProcedure
;========================Run======================================
OpenConsole()                                 
ConsoleTitle ("PureBasic - Console Example:")                                                  
EnableGraphicalConsole(1)

  Define i.l, n.l, tau.d, taufast.f
  Dim a.f(100)
  Dim b.f(100)
 
  a(0)=1 :a(1)=1:  a(2)=1:  a(3)=2:  a(4)=3:  a(5)=3:  a(6)=4:  a(7)=4 :
  b(0)=1 :b(1)=2:  b(2)=1:  b(3)=3:  b(4)=3:  b(5)=5:  b(6)=5:  b(7)=5
  
  tau=KendallSmallN(a(),b(),8)
  PrintN("KendallSmallN="+StrD(tau)+ " Should be =0.869565")
  Input()
 CloseConsole()
End

Code: Select all

/*-------------------------------------------------------------------------*/
/* This file contains code to calculate Kendall's Tau in O(N log N) time in
 * a manner similar to the following reference:
 *
 * A Computer Method for Calculating Kendall's Tau with Ungrouped Data
 * William R. Knight Journal of the American Statistical Association,
 * Vol. 61, No. 314, Part 1 (Jun., 1966), pp. 436-439
 *
 * Copyright 2010 David Simcha
 *
 * License:
 * Boost Software License - Version 1.0 - August 17th, 2003
 *
 * Permission is hereby granted, free of charge, to any person or organization
 * obtaining a copy of the software and accompanying documentation covered by
 * this license (the "Software") to use, reproduce, display, distribute,
 * execute, and transmit the Software, and to prepare derivative works of the
 * Software, and to permit third-parties to whom the Software is furnished to
 * do so, all subject to the following:
 *
 * The copyright notices in the Software and this entire statement, including
 * the above license grant, this restriction and the following disclaimer,
 * must be included in all copies of the Software, in whole or in part, and
 * all derivative works of the Software, unless such copies or derivative
 * works are solely in the form of machine-executable object code generated by
 * a source language processor.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE, TITLE AND NON-INFRINGEMENT. IN NO EVENT
 * SHALL THE COPYRIGHT HOLDERS OR ANYONE DISTRIBUTING THE SOFTWARE BE LIABLE
 * FOR ANY DAMAGES OR OTHER LIABILITY, WHETHER IN CONTRACT, TORT OR OTHERWISE,
 * ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
 * DEALINGS IN THE SOFTWARE.
 *
 */
/*-------------------------------------------------------------------------*/
/**---------------- Adapted for AFNI by RWCox - April 2010 ---------------**/
/*-------------------------------------------------------------------------*/

#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <stdio.h>
#include <assert.h>
#include <time.h>

#ifdef SOLARIS_OLD      /* 'if' might be overkill, but just to be minimal... */
                        /*                               10 May 2010 [rickr] */
#include "machdep.h"
#else
#include <stdint.h>
#endif


/*-------------------------------------------------------------------------*/
/* Sorts in place, returns the bubble sort distance between the input array
 * and the sorted array.
 */

static int insertionSort(float *arr, int len)
{
    int maxJ, i,j , swapCount = 0;

/* printf("enter insertionSort len=%d\n",len) ; */

    if(len < 2) { return 0; }

    maxJ = len - 1;
    for(i = len - 2; i >= 0; --i) {
        float  val = arr[i];
        for(j=i; j < maxJ && arr[j + 1] < val; ++j) {
            arr[j] = arr[j + 1];
        }

        arr[j] = val;
        swapCount += (j - i);
    }

    return swapCount;
}

/*-------------------------------------------------------------------------*/

static int merge(float *from, float *to, int middle, int len)
{
    int bufIndex, leftLen, rightLen , swaps ;
    float *left , *right;

/* printf("enter merge\n") ; */

    bufIndex = 0;
    swaps = 0;

    left = from;
    right = from + middle;
    rightLen = len - middle;
    leftLen = middle;

    while(leftLen && rightLen) {
        if(right[0] < left[0]) {
            to[bufIndex] = right[0];
            swaps += leftLen;
            rightLen--;
            right++;
        } else {
            to[bufIndex] = left[0];
            leftLen--;
            left++;
        }
        bufIndex++;
    }

    if(leftLen) {
#pragma omp critical (MEMCPY)
        memcpy(to + bufIndex, left, leftLen * sizeof(float));
    } else if(rightLen) {
#pragma omp critical (MEMCPY)
        memcpy(to + bufIndex, right, rightLen * sizeof(float));
    }

    return swaps;
}

/*-------------------------------------------------------------------------*/
/* Sorts in place, returns the bubble sort distance between the input array
 * and the sorted array.
 */

static int mergeSort(float *x, float *buf, int len)
{
    int swaps , half ;

/* printf("enter mergeSort\n") ; */

    if(len < 10) {
        return insertionSort(x, len);
    }

    swaps = 0;

    if(len < 2) { return 0; }

    half = len / 2;
    swaps += mergeSort(x, buf, half);
    swaps += mergeSort(x + half, buf + half, len - half);
    swaps += merge(x, buf, half, len);

#pragma omp critical (MEMCPY)
    memcpy(x, buf, len * sizeof(float));
    return swaps;
}

/*-------------------------------------------------------------------------*/

static int getMs(float *data, int len)  /* Assumes data is sorted */
{
    int Ms = 0, tieCount = 0 , i ;

/* printf("enter getMs\n") ; */

    for(i = 1; i < len; i++) {
        if(data[i] == data[i-1]) {
            tieCount++;
        } else if(tieCount) {
            Ms += (tieCount * (tieCount + 1)) / 2;
            tieCount = 0;
        }
    }
    if(tieCount) {
        Ms += (tieCount * (tieCount + 1)) / 2;
    }
    return Ms;
}

/*-------------------------------------------------------------------------*/
/* This function calculates the Kendall correlation tau_b.
 * The arrays arr1 should be sorted before this call, and arr2 should be
 * re-ordered in lockstep.  This can be done by calling
 *   qsort_floatfloat(len,arr1,arr2)
 * for example.
 * Note also that arr1 and arr2 will be modified, so if they need to
 * be preserved, do so before calling this function.
 */

float kendallNlogN( float *arr1, float *arr2, int len )
{
    int m1 = 0, m2 = 0, tieCount, swapCount, nPair, s,i ;
    float cor ;

/* printf("enter kendallNlogN\n") ; */

    if( len < 2 ) return (float)0 ;

    nPair = len * (len - 1) / 2;
    s = nPair;

    tieCount = 0;
    for(i = 1; i < len; i++) {
        if(arr1[i - 1] == arr1[i]) {
            tieCount++;
        } else if(tieCount > 0) {
            insertionSort(arr2 + i - tieCount - 1, tieCount + 1);
            m1 += tieCount * (tieCount + 1) / 2;
            s += getMs(arr2 + i - tieCount - 1, tieCount + 1);
            tieCount = 0;
        }
    }
    if(tieCount > 0) {
        insertionSort(arr2 + i - tieCount - 1, tieCount + 1);
        m1 += tieCount * (tieCount + 1) / 2;
        s += getMs(arr2 + i - tieCount - 1, tieCount + 1);
    }

    swapCount = mergeSort(arr2, arr1, len);

    m2 = getMs(arr2, len);
    s -= (m1 + m2) + 2 * swapCount;

    if( m1 < nPair && m2 < nPair )
      cor = s / ( sqrtf((float)(nPair-m1)) * sqrtf((float)(nPair-m2)) ) ;
    else
      cor = 0.0f ;

    return cor ;
}

/*-------------------------------------------------------------------------*/
/* This function uses a simple O(N^2) implementation.  It probably has a
 * smaller constant and therefore is useful in the small N case, and is also
 * useful for testing the relatively complex O(N log N) implementation.
 */

float kendallSmallN( float *arr1, float *arr2, int len )
{
    int m1 = 0, m2 = 0, s = 0, nPair , i,j ;
    float cor ;

/* printf("enter kendallSmallN\n") ; */

    for(i = 0; i < len; i++) {
        for(j = i + 1; j < len; j++) {
            if(arr2[i] > arr2[j]) {
                if (arr1[i] > arr1[j]) {
                    s++;
                } else if(arr1[i] < arr1[j]) {
                    s--;
                } else {
                    m1++;
                }
            } else if(arr2[i] < arr2[j]) {
                if (arr1[i] > arr1[j]) {
                    s--;
                } else if(arr1[i] < arr1[j]) {
                    s++;
                } else {
                    m1++;
                }
            } else {
                m2++;

                if(arr1[i] == arr1[j]) {
                    m1++;
                }
            }
        }
    }

    nPair = len * (len - 1) / 2;

    if( m1 < nPair && m2 < nPair )
      cor = s / ( sqrtf((float)(nPair-m1)) * sqrtf((float)(nPair-m2)) ) ;
    else
      cor = 0.0f ;

    return cor ;
}

#if 0

int main( int argc , char *argv[] )
{
    float  a[100], b[100];
    float  smallNCor, smallNCov, largeNCor, largeNCov;
    int i;
    int N ;
#define NBOT  10
#define NTOP  2000
#define NSTEP 10

    /* Test the small N version against a few values obtained from the old
     * version in R.  Only exercising the small N version because the large
     * N version requires the first array to be sorted and the second to be
     * reordered in lockstep before it's called.*/
    {
        float  a1[] = {1,2,3,5,4};
        float  a2[] = {1,2,3,3,5};
        float  b1[] = {8,6,7,5,3,0,9};
        float  b2[] = {3,1,4,1,5,9,2};
        float  c1[] = {1,1,1,2,3,3,4,4};
        float  c2[] = {1,2,1,3,3,5,5,5};

        printf("kSmall 5 = %g\n", kendallSmallN(a1, a2, 5) ) ;
        printf("kSmall 7 = %g\n", kendallSmallN(b1, b2, 7) ) ;
        printf("kSmall 8 = %g\n", kendallSmallN(c1, c2, 8) ) ;
    }

    /* Now that we're confident that the simple, small N version works,
     * extensively test it against the much more complex and bug-prone
     * O(N log N) version.
     */
    for(i = 0; i < 100; i++) {
        int j, len;
        for(j = 0; j < 100; j++) {
            a[j] = rand() % 30;
            b[j] = rand() % 30;
        }

        len = rand() % 50 + 50;

        /* The large N version assumes that the first array is sorted.  This
         * will usually be made true in R before passing the arrays to the
         * C functions.
         */
        insertionSort(a, len);

        smallNCor = kendallSmallN(a, b, len);
        largeNCor = kendallNlogN (a, b, len);
        printf("Cor: kSmall = %g  kNlogN = %g\n",smallNCor,largeNCor) ;
    }

    printf("-- short tests done --\n") ;

    /* Speed test.  Compare the O(N^2) version, which is very similar to
     * R's current impl, to my O(N log N) version.
     */
    {
        float  *foo, *bar, *buf;
        int i;
        float  startTime, stopTime;

        foo = (float *) malloc(NTOP * sizeof(float));
        bar = (float *) malloc(NTOP * sizeof(float));
        buf = (float *) malloc(NTOP * sizeof(float));

     for( N=NBOT ; N <= NTOP ; N += NSTEP ){
        for(i = 0; i < N; i++) {
            foo[i] = rand();
            bar[i] = rand();
        }

        startTime = clock();
        smallNCor = kendallSmallN(foo, bar, N);
        stopTime = clock();
        printf("N=%d: slow = %f ms  val = %g\n", N,stopTime - startTime,smallNCor);

        startTime = clock();

        /* Only sorting first array.  Normally the second one would be
         * reordered in lockstep.
         */
#if 1
        qsort_floatfloat( N , foo , bar ) ;
#else
        mergeSort(foo, buf, N);
#endif
        largeNCor = kendallNlogN(foo, bar, N);
        stopTime = clock();
        printf("N=%d: fast = %f ms  val = %g\n", N,stopTime - startTime,largeNCor);
     }
    }

    return 0;
}

#endif

Code: Select all

;==============KendalTauSlow=======================================
Procedure.d kendallSmallN(  *arr1.double,  *arr2.double,  len.l )
    Protected m1.l = 0, m2.l = 0, s.l = 0, nPair.l , i.l,j.l ,qSOD.q
    Protected cor.d

    qSOD=SizeOf(double)

    For i = 0 To len-1
        For j = i + 1 To len-1
            If PeekD(*arr2+i*qSOD) > PeekD(*arr2+j*qSOD)
                If PeekD(*arr1+i*qSOD) > PeekD(*arr1+j*qSOD)
                    s=s +1
                ElseIf PeekD(*arr1+i*qSOD) < PeekD(*arr1+j*qSOD)
                    s=s-1
                Else
                    m1=m1+1
                EndIf
           ElseIf PeekD(*arr2+i*qSOD) < PeekD(*arr2+j*qSOD)
                If PeekD(*arr1+i*qSOD) > PeekD(*arr1+j*qSOD)
                    s=s-1
                ElseIf PeekD(*arr1+i*qSOD) < PeekD(*arr1+j*qSOD)
                    s=s+1;
                Else
                    m1=m1+1
                EndIf
            Else
                m2=m2+1

                If PeekD(*arr1+i*qSOD) = PeekD(*arr1+j*qSOD)
                    m1=m1+1;
                EndIf
            EndIf
        Next j
    Next i

    nPair = len * (len - 1) / 2

    If( m1 < nPair And m2 < nPair )
      cor = s / ( Sqr(nPair-m1) * Sqr(nPair-m2) )
    Else
      cor = 0
    EndIf
   
    ProcedureReturn cor
     
EndProcedure

Code: Select all

Structure DOUBLEARRAY
  d.d[0]
EndStructure
Procedure.d kendallSmallN(  *arr1.DOUBLEARRAY,  *arr2.DOUBLEARRAY,  len.l )
    Protected m1.l = 0, m2.l = 0, s.l = 0, nPair.l , i.l,j.l
    Protected cor.d

    For i = 0 To len-1
        For j = i + 1 To len-1
            If *arr2\d[i] > *arr2\d[j]
                If *arr1\d[i] > *arr1\d[j]
                    s=s +1
                ElseIf *arr1\d[i] < *arr1\d[j]
                    s=s-1
                Else
                    m1=m1+1
                EndIf
           ElseIf *arr2\d[i] < *arr2\d[j]
                If *arr1\d[i] > *arr1\d[j]
                    s=s-1
                ElseIf *arr1\d[i] < *arr1\d[j]
                    s=s+1;
                Else
                    m1=m1+1
                EndIf
            Else
                m2=m2+1

                If *arr1\d[i] = *arr1\d[j]
                    m1=m1+1;
                EndIf
            EndIf
        Next j
    Next i

    nPair = len * (len - 1) / 2

    If( m1 < nPair And m2 < nPair )
      cor = s / ( Sqr(nPair-m1) * Sqr(nPair-m2) )
    Else
      cor = 0
    EndIf
   
    ProcedureReturn cor
     
EndProcedure

Code: Select all

EnableExplicit

Structure DD
  d.d[0]
EndStructure
;================KendallTauSlow from WFME=================
Procedure.d KendallTaux(Array dArrIn.d(1), lNum.l)

; Kendall tau algorithm from  numerical receipes in Fortran  p493 , C++ p 647
  Protected n.l, ii.l, iis.l, n2.l, n1.l
  Protected svar.d, aa.d, a2.d, a1.d, z.d
  Dim ddata1.d(1)
  Dim ddata2.d(1)
  Protected j.l, k.l, tau.d

  iis=0
  n2=0
  n1=0

  n=lnum
  ReDim ddata1( n)
  ReDim ddata2 ( n)

  For ii= 1 To n
    ddata1(ii)=ii
    ddata2(ii)=dArrIn(ii)
  Next ii

  For j=1 To n-1
    For k=j+1 To n
      a1=ddata1(j)-ddata1(k)
      a2=ddata2(j)-ddata2(k)
      aa=a1*a2
      If aa <> 0.0 
        n1=n1+1
        n2=n2+1
        If aa>0.0
          iis=iis+1
        EndIf  
        If aa<=0.0
          iis=iis-1
        EndIf  
      Else
        If a1<>0.0
          n1=n1+1
        EndIf  
        If a2<>0.0
          n2=n2+1
        EndIf  
      EndIf ;aa<>0
    Next k
  Next j

  tau=iis/Sqr(n1*n2)
  svar=(4.0*n+10.0)/(9.0*n*(n-1.0))
  z=tau/Sqr(svar)
  ProcedureReturn tau
EndProcedure
;===============InsertionSort with pointers====================  
Procedure.l insertionSortx(*arr.double, len.l)
    Protected maxJ.l, i.l,j.l , swapCount.l = 0, val.d
    Protected  qSOD.q
    
    qSOD=SizeOf(double)
    
    If(len < 2) : ProcedureReturn 0 : EndIf
    maxJ = len - 1
    For i = len - 2 To  0 Step -1
        val = PeekD(*arr +i*qSOD)
        j=i
        While j<maxj And PeekD(*arr+(j+1)*qSOD )< val 
          PokeD(*arr+j*qSOD,PeekD((*arr+(j+1)*qSOD)))
          j=j+1
        Wend 
        PokeD(*arr+j*qSOD,val)
        swapCount= swapcount + (j - i);
    Next i

    ProcedureReturn swapCount
EndProcedure  
;===========InsertionSort with pointer & stucture=================
Procedure.l insertionSort(*arr.DD, len.l)
    Protected maxJ.l, i.l,j.l , swapCount.l = 0, val.d
    Protected  qSOD.q
    
    qSOD=SizeOf(double)
    
    If(len < 2) : ProcedureReturn 0 : EndIf
    maxJ = len - 1
    For i = len - 2 To  0 Step -1
        val = *arr\d[i] ;PeekD(*arr +i*qSOD)
        j=i
        ;While j<maxj And PeekD(*arr+(j+1)*qSOD )< val
        While j<maxj And *arr\d[j+1] < val  
           *arr\d[j] = *arr\d[j + 1]
          ;PokeD(*arr+j*qSOD,PeekD((*arr+(j+1)*qSOD)))
          j=j+1
        Wend ;}
        *arr\d[j] = val
        ;PokeD(*arr+j*qSOD,val)
        swapCount= swapcount + (j - i);
    Next i

    ProcedureReturn swapCount
EndProcedure  
 
;=================Merge====================================
Procedure.l merge( *from.DD, *too.DD,  middle.l,  len.l)
    Protected bufIndex.l, leftLen.l, rightLen.l , swaps.l, qSOD.q
    Protected *xleft.DD, *xright.DD, *temp

    qSOD=SizeOf(double)
    bufIndex = 0
    swaps = 0

    *xLeft = *from
    *xright = *from + middle*qSOD
    rightLen = len - middle
    leftLen = middle;

    While(leftLen And rightLen) 
      If *xright\d[0] < *xLeft\d[0]
            *too\d[bufindex]=*xright\d[0]
            swaps =swaps+ leftLen;
            rightLen=rightlen-1;
            *xright=*xright+qSOD
          Else 
            *too\d[bufIndex]=*xleft\d[0]
            leftlen=leftLen-1;
            *xleft=*xleft+qSOD
        EndIf
        bufindex=bufIndex+1
    Wend

    If leftLen 
        CopyMemory(*xleft, *too + bufIndex*qSOD, leftLen*qSOD) ;note souurce and dest opposite
    ElseIf rightLen 
        CopyMemory(*xright, *too + bufIndex*qSOD, rightLen*qSOD)
    EndIf

    ProcedureReturn swaps
  EndProcedure

;========================MergeSort===================================  
Procedure.l mergeSort(*x.DD, *buf.DD, len.l)

    Protected swaps.l , half.l, qSOD.q

    qSOD=SizeOf(double)

    If(len < 10) 
        ProcedureReturn insertionSort(*x, len);
    EndIf

    swaps = 0;

    If(len < 2) : ProcedureReturn  0 : EndIf

    half = len / 2;
    swaps = swaps + mergeSort(*x, *buf, half);
    swaps = swaps + mergeSort(*x + half*qSOD, *buf + half*qSOD, len - half);
    swaps = swaps + merge(*x, *buf, half, len);

    CopyMemory(*buf, *x, len * qSOD)
    ;note source and dest opposite from C memcpy(dest,source,sizeof(dest)
    ;CopyMemory(source,Dest,Size)
    ProcedureReturn swaps              
    
  EndProcedure 
;=======================GetMs====================================  
Procedure.l getMs(*data.DD, len.l)  ;/* Assumes Data is sorted */
    Protected Ms.l = 0, tieCount.l = 0 , i.l, qSOD.q 

    qSOD=SizeOf(double)

    For i = 1 To len-1  
      If *data\d[i] = *data\d[(i-1)]
            tieCount + 1
        ElseIf tieCount
            Ms =ms+ (tieCount * (tieCount + 1)) / 2;
            tieCount = 0;
        EndIf
    Next i
    If(tieCount) 
      Ms = ms + (tieCount * (tieCount + 1)) / 2;
      ;==================================================
      tiecount +1 ; note in orginial but not RWCox version
      ;===================================================
    EndIf
    ProcedureReturn Ms;
  EndProcedure  
;=======================KendallTauFast================================
; * Note also that arr1 And arr2 will be modified, so If they need To
; * be preserved, do so before calling this function.
Procedure.d kendallNlogN( *arr1.DD, *arr2.DD, len.l)

    Protected  m1.l = 0, m2.l = 0, tieCount.l, swapCount.l, nPair.l, s.l,i.l 
    Protected cor.d,qSOD.q, ktau.d, den1.d, den2.d

    qSOD=SizeOf(double)
    If len < 2 : ProcedureReturn 0 : EndIf

    nPair = len * (len - 1) / 2
    s = nPair

    tieCount = 0;
    For i = 1 To  len-1
        If *arr1\d[i - 1] = *arr1\d[i]
            tieCount +1
        ElseIf tieCount > 0
            insertionSort(*arr2 + (i - tieCount - 1)*qSOD, tieCount + 1);
            m1 = m1 + tieCount * (tieCount + 1) / 2;
            s = s + getMs(*arr2 + (i - tieCount - 1)*qSOD, tieCount + 1);
            tieCount +1
            tieCount = 0                                                         ;
        EndIf    
   Next i

    If(tieCount > 0) 
        insertionSort(*arr2 + (i - tieCount - 1)*qSOD, tieCount + 1);
        m1 = m1 + tieCount * (tieCount + 1) / 2;
        s = s + getMs(*arr2 + (i - tieCount - 1)*qSOD, tieCount + 1);
    EndIf

    swapCount = mergeSort(*arr2, *arr1, len)

    m2 = getMs(*arr2, len)
    s = s - ((m1 + m2) + 2 * swapCount)
    
    den1=nPair-m1
    den2=nPair-m2

    If( den1>0 And Den2>0 )
      ktau= s /Sqr(den1)/Sqr(den2)
    Else
      ktau = 0 
    EndIf
    
    ProcedureReturn ktau 
EndProcedure
;==============KendalTauSlow=======================================
Procedure.d kendallSmallN(  *arr1.DD,  *arr2.DD,  len.l )
    Protected m1.l = 0, m2.l = 0, s.l = 0, nPair.l , i.l,j.l
    Protected ktau.d, den1.d, den2.d

    For i = 0 To len-1
        For j = i + 1 To len-1
            If *arr2\d[i] > *arr2\d[j]
                If *arr1\d[i] > *arr1\d[j]
                    s=s +1
                ElseIf *arr1\d[i] < *arr1\d[j]
                    s=s-1
                Else
                    m1=m1+1
                EndIf
           ElseIf *arr2\d[i] < *arr2\d[j]
                If *arr1\d[i] > *arr1\d[j]
                    s=s-1
                ElseIf *arr1\d[i] < *arr1\d[j]
                    s=s+1;
                Else
                    m1=m1+1
                EndIf
            Else
                m2=m2+1

                If *arr1\d[i] = *arr1\d[j]
                    m1=m1+1;
                EndIf
            EndIf
        Next j
    Next i

    nPair = len * (len - 1) / 2
    ktau=0
      den1= nPair-m1
      den2= nPair-m2
      If( den1>0 And den2>0 )
        ktau = s /Sqr(den1)/Sqr(den2)
      EndIf
  
    ProcedureReturn ktau
     
EndProcedure
;========================Run======================================
OpenConsole()                                 
ConsoleTitle ("PureBasic - Console Example:")                                                  
EnableGraphicalConsole(1)

  Define i.l, n.l, tau.d,taux.d, taufast.d, s.s
  Dim c1.d(1000)
  Dim c2.d(1000)
 
    i=0
  If ReadFile(0, "c:\Text.txt")   
    While Eof(0) = 0           
      s=ReadString(0)
      c2(i)=ValD(s)
      c1(i)=i
      ;PrintN(Str(c1(i))+" "+StrD(c2(i)))
      i=i+1
    Wend
    CloseFile(0)              
  Else
    MessageRequester("Information","Couldn't open the file!")
  EndIf

  taux=KendallTaux(c2(),i-1)
  tau=KendallSmallN(c1(),c2(),i-1)
  
  ;insertionSort(c1(),i-1) ;if c1() not sorted then uncomment this line
  taufast=kendallNlogN(c1(),c2(),i-1)
  PrintN("KendallTaux="+StrD(taux))
  PrintN("KendallSmallN="+StrD(tau))
  PrintN("KendallTauFast="+StrD(taufast))
 
  Input()
 CloseConsole()
End