public static int medoid(boolean parallel, Collection<Integer> indecies, double tol, List<? extends Vec> X, DistanceMetric dm, List<Double> accel)
    {
        final int N = indecies.size();
        
        if(tol <= 0 || N < SystemInfo.LogicalCores)//Really just not enough points, lets simplify
            return PAM.medoid(parallel, indecies, X, dm, accel);
        
        
        final double log2d = Math.log(1)-Math.log(tol);
        
        /**
         * Online estimate of the standard deviation that will be used
         */
        final OnLineStatistics distanceStats;
        /**
         * This array contains the current sum of all distance computations done
         * for each index. Corresponds to mu in the paper.
         */
        AtomicDoubleArray totalDistSum = new AtomicDoubleArray(N);
        /**
         * This array contains the current number of distance computations that
         * have been done for each feature index. Corresponds to T_i in the
         * paper.
         */
        AtomicIntegerArray totalDistCount = new AtomicIntegerArray(N);
        final int[] indx_map = indecies.stream().mapToInt(i->i).toArray();
        final boolean symetric = dm.isSymmetric();
        final double[] lower_bound_est = new double[N];
        final double[] upper_bound_est = new double[N];

        ThreadLocal<Random> localRand = ThreadLocal.withInitial(RandomUtil::getRandom);

        //First pass, lets pull every "arm" (compute a dsitance) for each datumn at least once, so that we have estiamtes to work with. 
        distanceStats = ParallelUtils.run(parallel, N, (start, end)->
        {
            Random rand = localRand.get();
            OnLineStatistics localStats = new OnLineStatistics();
            for(int i = start; i < end; i++)
            {
                int j = rand.nextInt(N);
                while(j == i)
                    j = rand.nextInt(N);
                

                double d_ij = dm.dist(indx_map[i], indx_map[j], X, accel);
                localStats.add(d_ij);
                totalDistSum.addAndGet(i, d_ij);
                totalDistCount.incrementAndGet(i);
                if(symetric)
                {
                    totalDistSum.addAndGet(j, d_ij);
                    totalDistCount.incrementAndGet(j);
                }
            }
            
            return localStats;
        }, (a,b)-> OnLineStatistics.add(a, b));
        
        //Now lets prepare the lower and upper bound estimates
        ConcurrentSkipListSet<Integer> lowerQ = new ConcurrentSkipListSet<>((Integer o1, Integer o2) -> 
        {
            int cmp = Double.compare(lower_bound_est[o1], lower_bound_est[o2]);
            if(cmp == 0)//same bounds, but sort by identity to avoid issues
                cmp = o1.compareTo(o2);
            return cmp;
        });
        
        ConcurrentSkipListSet<Integer> upperQ = new ConcurrentSkipListSet<>((Integer o1, Integer o2) -> 
        {
            int cmp = Double.compare(upper_bound_est[o1], upper_bound_est[o2]);
            if(cmp == 0)//same bounds, but sort by identity to avoid issues
                cmp = o1.compareTo(o2);
            return cmp;
        });

        ParallelUtils.run(parallel, N, (start, end)->
        {
            double v = distanceStats.getVarance();
            for(int i = start; i < end; i++)
            {
                int T_i = totalDistCount.get(i);
                double c_i = Math.sqrt(2*v*log2d/T_i);
                lower_bound_est[i] = totalDistSum.get(i)/T_i - c_i;
                upper_bound_est[i] = totalDistSum.get(i)/T_i + c_i;
                lowerQ.add(i);
                upperQ.add(i);
            }
        });
        
        
        //Now lets start sampling! 
        
        //how many points should we pick and sample? Not really discussed in paper- but a good idea for efficency (dont want to pay that Q cost as much as possible)
        /**
         * to-pull is how many arms we will select per iteration
         */
        int num_to_pull;
        /**
         * to sample is how many random pairs we will pick for each pulled arm
         */
        int samples;
        
        if(parallel)
        {
            num_to_pull = Math.max(SystemInfo.LogicalCores, 32);
            samples = Math.min(32, N-1);
        }
        else
        {
            num_to_pull = Math.min(32, N);
            samples = Math.min(32, N-1);
        }
        
        /**
         * The levers we will pull this iteration, and then add back in
         */
        IntList to_pull = new IntList();
        /**
         * the levers we must add back in but not update b/c they hit max evaluations and the confidence bound is tight
         */
        IntList toAddBack = new IntList();
        boolean[] isExact = new boolean[N];
        Arrays.fill(isExact, false);
        int numExact = 0;
        
        
        while(numExact < N)//loop should break out before this ever happens
        {
            to_pull.clear();
            toAddBack.clear();

            //CONVERGENCE CEHCK
            if(upper_bound_est[upperQ.first()] < lower_bound_est[lowerQ.first()])
            {
                //WE are done!
                return indx_map[upperQ.first()];
            }
            

            while(to_pull.size() < num_to_pull)
            {
                
                if(lowerQ.isEmpty())
                    break;//we've basically evaluated everyone
                int i = lowerQ.pollFirst();
                
                
                if(totalDistCount.get(i) >= N-1 && !isExact[i])//Lets just replace with exact value
                {
                    double avg_d_i = ParallelUtils.run(parallel, N, (start, end)->
                    {
                        double d = 0;
                        for (int j = start; j < end; j++)
                            if (i != j)
                                d += dm.dist(indx_map[i], indx_map[j], X, accel);
                        return d;
                    }, (a, b)->a+b);
                    avg_d_i /= N-1;
                    
                    upperQ.remove(i);
                    lower_bound_est[i] = upper_bound_est[i] = avg_d_i;
                    totalDistSum.set(i, avg_d_i);
                    totalDistCount.set(i, N);
                    isExact[i] = true;
                    numExact++;
//                    System.out.println("Num Exact: " + numExact);
                    //OK, exavt value for datumn I is set. 
                    toAddBack.add(i);
                }
                

                if(!isExact[i])
                    to_pull.add(i);
            }
            
            //OK, lets now pull a bunch of levers / measure distances
            
            OnLineStatistics changeInStats = ParallelUtils.run(parallel, to_pull.size(), (start, end)->
            {
                Random rand = localRand.get();
                OnLineStatistics localStats = new OnLineStatistics();
                for(int i_count = start; i_count < end; i_count++)
                {
                    int i = to_pull.get(i_count);
                    for(int j_count = 0; j_count < samples; j_count++)
                    {
                        int j = rand.nextInt(N);
                        while(j == i)
                            j = rand.nextInt(N);
                        
                        double d_ij = dm.dist(indx_map[i], indx_map[j], X, accel);
                        localStats.add(d_ij);
                        totalDistSum.addAndGet(i, d_ij);
                        totalDistCount.incrementAndGet(i);
                        if(symetric && !isExact[j])
                        {
                            totalDistSum.addAndGet(j, d_ij);
                            totalDistCount.incrementAndGet(j);
                        }
                    }
                }
                
                return localStats;
            }, (a,b) -> OnLineStatistics.add(a, b));
            
            if(!to_pull.isEmpty())//might be empty if everyone went over the threshold
                distanceStats.add(changeInStats);
            
            //update bounds and re-insert
            double v = distanceStats.getVarance();
            //we are only updating the bounds on the levers we pulled
            //that may mean some old bounds are stale
            //these values are exact
            lowerQ.addAll(toAddBack);
            upperQ.addAll(toAddBack);
            upperQ.removeAll(to_pull);
            for(int i : to_pull)
            {
                int T_i = totalDistCount.get(i);
                double c_i = Math.sqrt(2*v*log2d/T_i);
                lower_bound_est[i] = totalDistSum.get(i)/T_i - c_i;
                upper_bound_est[i] = totalDistSum.get(i)/T_i + c_i;
                lowerQ.add(i);
                upperQ.add(i);
            }
        }
        
        //We can reach this point on small N or low D datasets. Iterate and return the correct value
        int bestIndex = 0;
        for(int i = 1; i < N; i++)
            if(lower_bound_est[i] < lower_bound_est[bestIndex])
                bestIndex = i;
        
        return bestIndex;
    }
public static int medoid ( boolean parallel , Collection < Integer > indecies , double tol , List < ? extends Vec > X , DistanceMetric dm , List < Double > accel ) { final int N = indecies . size ( ) ; if ( tol <= 0 || N < SystemInfo . LogicalCores ) return PAM . medoid ( parallel , indecies , X , dm , accel ) ; final double log2d = Math . log ( 1 ) - Math . log ( tol ) ; final OnLineStatistics distanceStats ; AtomicDoubleArray totalDistSum = new AtomicDoubleArray ( N ) ; AtomicIntegerArray totalDistCount = new AtomicIntegerArray ( N ) ; final int [ ] indx_map = indecies . stream ( ) . mapToInt ( i -> i ) . toArray ( ) ; final boolean symetric = dm . isSymmetric ( ) ; final double [ ] lower_bound_est = new double [ N ] ; final double [ ] upper_bound_est = new double [ N ] ; ThreadLocal < Random > localRand = ThreadLocal . withInitial ( RandomUtil :: getRandom ) ; distanceStats = ParallelUtils . run ( parallel , N , ( start , end ) -> { Random rand = localRand . get ( ) ; OnLineStatistics localStats = new OnLineStatistics ( ) ; for ( int i = start ; i < end ; i ++ ) { int j = rand . nextInt ( N ) ; while ( j == i ) j = rand . nextInt ( N ) ; double d_ij = dm . dist ( indx_map [ i ] , indx_map [ j ] , X , accel ) ; localStats . add ( d_ij ) ; totalDistSum . addAndGet ( i , d_ij ) ; totalDistCount . incrementAndGet ( i ) ; if ( symetric ) { totalDistSum . addAndGet ( j , d_ij ) ; totalDistCount . incrementAndGet ( j ) ; } } return localStats ; } , ( a , b ) -> OnLineStatistics . add ( a , b ) ) ; ConcurrentSkipListSet < Integer > lowerQ = new ConcurrentSkipListSet < > ( ( Integer o1 , Integer o2 ) -> { int cmp = Double . compare ( lower_bound_est [ o1 ] , lower_bound_est [ o2 ] ) ; if ( cmp == 0 ) cmp = o1 . compareTo ( o2 ) ; return cmp ; } ) ; ConcurrentSkipListSet < Integer > upperQ = new ConcurrentSkipListSet < > ( ( Integer o1 , Integer o2 ) -> { int cmp = Double . compare ( upper_bound_est [ o1 ] , upper_bound_est [ o2 ] ) ; if ( cmp == 0 ) cmp = o1 . compareTo ( o2 ) ; return cmp ; } ) ; ParallelUtils . run ( parallel , N , ( start , end ) -> { double v = distanceStats . getVarance ( ) ; for ( int i = start ; i < end ; i ++ ) { int T_i = totalDistCount . get ( i ) ; double c_i = Math . sqrt ( 2 * v * log2d / T_i ) ; lower_bound_est [ i ] = totalDistSum . get ( i ) / T_i - c_i ; upper_bound_est [ i ] = totalDistSum . get ( i ) / T_i + c_i ; lowerQ . add ( i ) ; upperQ . add ( i ) ; } } ) ; int num_to_pull ; int samples ; if ( parallel ) { num_to_pull = Math . max ( SystemInfo . LogicalCores , 32 ) ; samples = Math . min ( 32 , N - 1 ) ; } else { num_to_pull = Math . min ( 32 , N ) ; samples = Math . min ( 32 , N - 1 ) ; } IntList to_pull = new IntList ( ) ; IntList toAddBack = new IntList ( ) ; boolean [ ] isExact = new boolean [ N ] ; Arrays . fill ( isExact , false ) ; int numExact = 0 ; while ( numExact < N ) { to_pull . clear ( ) ; toAddBack . clear ( ) ; if ( upper_bound_est [ upperQ . first ( ) ] < lower_bound_est [ lowerQ . first ( ) ] ) { return indx_map [ upperQ . first ( ) ] ; } while ( to_pull . size ( ) < num_to_pull ) { if ( lowerQ . isEmpty ( ) ) break ; int i = lowerQ . pollFirst ( ) ; if ( totalDistCount . get ( i ) >= N - 1 && ! isExact [ i ] ) { double avg_d_i = ParallelUtils . run ( parallel , N , ( start , end ) -> { double d = 0 ; for ( int j = start ; j < end ; j ++ ) if ( i != j ) d += dm . dist ( indx_map [ i ] , indx_map [ j ] , X , accel ) ; return d ; } , ( a , b ) -> a + b ) ; avg_d_i /= N - 1 ; upperQ . remove ( i ) ; lower_bound_est [ i ] = upper_bound_est [ i ] = avg_d_i ; totalDistSum . set ( i , avg_d_i ) ; totalDistCount . set ( i , N ) ; isExact [ i ] = true ; numExact ++ ; toAddBack . add ( i ) ; } if ( ! isExact [ i ] ) to_pull . add ( i ) ; } OnLineStatistics changeInStats = ParallelUtils . run ( parallel , to_pull . size ( ) , ( start , end ) -> { Random rand = localRand . get ( ) ; OnLineStatistics localStats = new OnLineStatistics ( ) ; for ( int i_count = start ; i_count < end ; i_count ++ ) { int i = to_pull . get ( i_count ) ; for ( int j_count = 0 ; j_count < samples ; j_count ++ ) { int j = rand . nextInt ( N ) ; while ( j == i ) j = rand . nextInt ( N ) ; double d_ij = dm . dist ( indx_map [ i ] , indx_map [ j ] , X , accel ) ; localStats . add ( d_ij ) ; totalDistSum . addAndGet ( i , d_ij ) ; totalDistCount . incrementAndGet ( i ) ; if ( symetric && ! isExact [ j ] ) { totalDistSum . addAndGet ( j , d_ij ) ; totalDistCount . incrementAndGet ( j ) ; } } } return localStats ; } , ( a , b ) -> OnLineStatistics . add ( a , b ) ) ; if ( ! to_pull . isEmpty ( ) ) distanceStats . add ( changeInStats ) ; double v = distanceStats . getVarance ( ) ; lowerQ . addAll ( toAddBack ) ; upperQ . addAll ( toAddBack ) ; upperQ . removeAll ( to_pull ) ; for ( int i : to_pull ) { int T_i = totalDistCount . get ( i ) ; double c_i = Math . sqrt ( 2 * v * log2d / T_i ) ; lower_bound_est [ i ] = totalDistSum . get ( i ) / T_i - c_i ; upper_bound_est [ i ] = totalDistSum . get ( i ) / T_i + c_i ; lowerQ . add ( i ) ; upperQ . add ( i ) ; } } int bestIndex = 0 ; for ( int i = 1 ; i < N ; i ++ ) if ( lower_bound_est [ i ] < lower_bound_est [ bestIndex ] ) bestIndex = i ; return bestIndex ; }
Computes the medoid of a sub-set of data
@param parallel whether or not the computation should be done using multiple cores
@param indecies the indexes of the points to get the medoid of
@param tol
@param X the list of all data
@param dm the distance metric to get the medoid with respect to
@param accel the acceleration cache for the distance metric
@return the index value contained within indecies that is the medoid