Concours: moyen le plus rapide de trier un grand nombre de données gaussiennes

Suite à l'intérêt que suscite cette question , j'ai pensé qu'il serait intéressant de proposer des réponses un peu plus objectives et quantitatives en proposant un concours.

L'idée est simple: j'ai généré un fichier binaire contenant 50 millions de doubles distribués en gauss (moyenne: 0, stdev 1). Le but est de créer un programme qui les triera en mémoire le plus rapidement possible. Une implémentation de référence très simple en python prend 1m4 à compléter. Jusqu'où pouvons-nous aller?

Les règles sont les suivantes: répondez avec un programme qui ouvre le fichier "gaussian.dat" et trie les nombres en mémoire (il n'est pas nécessaire de les sortir), ainsi que des instructions pour la construction et l'exécution du programme. Le programme doit pouvoir fonctionner sur ma machine Linux Arch (ce qui signifie que vous pouvez utiliser n’importe quel langage ou bibliothèque de programmation facilement installable sur ce système).

Le programme doit être raisonnablement lisible, afin que je puisse m'assurer qu'il peut être lancé en toute sécurité (pas de solution uniquement pour l'assembleur, s'il vous plaît!).

Je vais exécuter les réponses sur ma machine (quad core, 4 gigaoctets de RAM). La solution la plus rapide obtiendra la réponse acceptée et une prime de 100 points :)

Le programme utilisé pour générer les nombres:

#!/usr/bin/env python
import random
from array import array
from sys import argv
count=int(argv[1])
a=array('d',(random.gauss(0,1) for x in xrange(count)))
f=open("gaussian.dat","wb")
a.tofile(f)

L'implémentation de référence simple:

#!/usr/bin/env python
from array import array
from sys import argv
count=int(argv[1])
a=array('d')
a.fromfile(open("gaussian.dat"),count)
print "sorting..."
b=sorted(a)

EDIT: seulement 4 Go de RAM, désolé

EDIT # 2: Notez que le but du concours est de voir si nous pouvons utiliser des informations préalables sur les données . ce n'est pas supposé être un match nul entre différentes implémentations de langage de programmation!

Voici une solution en C ++ qui partitionne d'abord les nombres en compartiments avec le même nombre d'éléments attendu, puis trie chaque compartiment séparément. Il précalcule une table de la fonction de distribution cumulative basée sur certaines formules de Wikipedia, puis interpole les valeurs de cette table pour obtenir une approximation rapide.

Plusieurs étapes s'exécutent dans plusieurs threads pour utiliser les quatre cœurs.

#include <cstdlib>
#include <math.h>
#include <stdio.h>
#include <algorithm>

#include <tbb/parallel_for.h>

using namespace std;

typedef unsigned long long ull;

double signum(double x) {
    return (x<0) ? -1 : (x>0) ? 1 : 0;
}

const double fourOverPI = 4 / M_PI;

double erf(double x) {
    double a = 0.147;
    double x2 = x*x;
    double ax2 = a*x2;
    double f1 = -x2 * (fourOverPI + ax2) / (1 + ax2);
    double s1 = sqrt(1 - exp(f1));
    return signum(x) * s1;
}

const double sqrt2 = sqrt(2);

double cdf(double x) {
    return 0.5 + erf(x / sqrt2) / 2;
}

const int cdfTableSize = 200;
const double cdfTableLimit = 5;
double* computeCdfTable(int size) {
    double* res = new double[size];
    for (int i = 0; i < size; ++i) {
        res[i] = cdf(cdfTableLimit * i / (size - 1));
    }
    return res;
}
const double* const cdfTable = computeCdfTable(cdfTableSize);

double cdfApprox(double x) {
    bool negative = (x < 0);
    if (negative) x = -x;
    if (x > cdfTableLimit) return negative ? cdf(-x) : cdf(x);
    double p = (cdfTableSize - 1) * x / cdfTableLimit;
    int below = (int) p;
    if (p == below) return negative ? -cdfTable[below] : cdfTable[below];
    int above = below + 1;
    double ret = cdfTable[below] +
            (cdfTable[above] - cdfTable[below])*(p - below);
    return negative ? 1 - ret : ret;
}

void print(const double* arr, int len) {
    for (int i = 0; i < len; ++i) {
        printf("%e; ", arr[i]);
    }
    puts("");
}

void print(const int* arr, int len) {
    for (int i = 0; i < len; ++i) {
        printf("%d; ", arr[i]);
    }
    puts("");
}

void fillBuckets(int N, int bucketCount,
        double* data, int* partitions,
        double* buckets, int* offsets) {
    for (int i = 0; i < N; ++i) {
        ++offsets[partitions[i]];
    }

    int offset = 0;
    for (int i = 0; i < bucketCount; ++i) {
        int t = offsets[i];
        offsets[i] = offset;
        offset += t;
    }
    offsets[bucketCount] = N;

    int next[bucketCount];
    memset(next, 0, sizeof(next));
    for (int i = 0; i < N; ++i) {
        int p = partitions[i];
        int j = offsets[p] + next[p];
        ++next[p];
        buckets[j] = data[i];
    }
}

class Sorter {
public:
    Sorter(double* data, int* offsets) {
        this->data = data;
        this->offsets = offsets;
    }

    static void radixSort(double* arr, int len) {
        ull* encoded = (ull*)arr;
        for (int i = 0; i < len; ++i) {
            ull n = encoded[i];
            if (n & signBit) {
                n ^= allBits;
            } else {
                n ^= signBit;
            }
            encoded[i] = n;
        }

        const int step = 11;
        const ull mask = (1ull << step) - 1;
        int offsets[8][1ull << step];
        memset(offsets, 0, sizeof(offsets));

        for (int i = 0; i < len; ++i) {
            for (int b = 0, j = 0; b < 64; b += step, ++j) {
                int p = (encoded[i] >> b) & mask;
                ++offsets[j][p];
            }
        }

        int sum[8] = {0};
        for (int i = 0; i <= mask; i++) {
            for (int b = 0, j = 0; b < 64; b += step, ++j) {
                int t = sum[j] + offsets[j][i];
                offsets[j][i] = sum[j];
                sum[j] = t;
            }
        }

        ull* copy = new ull[len];
        ull* current = encoded;
        for (int b = 0, j = 0; b < 64; b += step, ++j) {
            for (int i = 0; i < len; ++i) {
                int p = (current[i] >> b) & mask;
                copy[offsets[j][p]] = current[i];
                ++offsets[j][p];
            }

            ull* t = copy;
            copy = current;
            current = t;
        }

        if (current != encoded) {
            for (int i = 0; i < len; ++i) {
                encoded[i] = current[i];
            }
        }

        for (int i = 0; i < len; ++i) {
            ull n = encoded[i];
            if (n & signBit) {
                n ^= signBit;
            } else {
                n ^= allBits;
            }
            encoded[i] = n;
        }
    }

    void operator() (tbb::blocked_range<int>& range) const {
        for (int i = range.begin(); i < range.end(); ++i) {
            double* begin = &data[offsets[i]];
            double* end = &data[offsets[i+1]];
            //std::sort(begin, end);
            radixSort(begin, end-begin);
        }
    }

private:
    double* data;
    int* offsets;
    static const ull signBit = 1ull << 63;
    static const ull allBits = ~0ull;
};

void sortBuckets(int bucketCount, double* data, int* offsets) {
    Sorter sorter(data, offsets);
    tbb::blocked_range<int> range(0, bucketCount);
    tbb::parallel_for(range, sorter);
    //sorter(range);
}

class Partitioner {
public:
    Partitioner(int bucketCount, double* data, int* partitions) {
        this->data = data;
        this->partitions = partitions;
        this->bucketCount = bucketCount;
    }

    void operator() (tbb::blocked_range<int>& range) const {
        for (int i = range.begin(); i < range.end(); ++i) {
            double d = data[i];
            int p = (int) (cdfApprox(d) * bucketCount);
            partitions[i] = p;
        }
    }

private:
    double* data;
    int* partitions;
    int bucketCount;
};

const int bucketCount = 512;
int offsets[bucketCount + 1];

int main(int argc, char** argv) {
    if (argc != 2) {
        printf("Usage: %s N\n N = the size of the input\n", argv[0]);
        return 1;
    }

    puts("initializing...");
    int N = atoi(argv[1]);
    double* data = new double[N];
    double* buckets = new double[N];
    memset(offsets, 0, sizeof(offsets));
    int* partitions = new int[N];

    puts("loading data...");
    FILE* fp = fopen("gaussian.dat", "rb");
    if (fp == 0 || fread(data, sizeof(*data), N, fp) != N) {
        puts("Error reading data");
        return 1;
    }
    //print(data, N);

    puts("assigning partitions...");
    tbb::parallel_for(tbb::blocked_range<int>(0, N),
            Partitioner(bucketCount, data, partitions));

    puts("filling buckets...");
    fillBuckets(N, bucketCount, data, partitions, buckets, offsets);
    data = buckets;

    puts("sorting buckets...");
    sortBuckets(bucketCount, data, offsets);

    puts("done.");

    /*
    for (int i = 0; i < N-1; ++i) {
        if (data[i] > data[i+1]) {
            printf("error at %d: %e > %e\n", i, data[i], data[i+1]);
        }
    }
    */

    //print(data, N);

    return 0;
}

Pour le compiler et l'exécuter, utilisez cette commande:

g++ -O3 -ltbb -o gsort gsort.cpp && time ./gsort 50000000

ÉDITER: Tous les compartiments sont maintenant placés dans le même tableau afin d’éviter de les copier dans le tableau. La taille de la table avec des valeurs précalculées a également été réduite, car les valeurs sont suffisamment précises. Néanmoins, si je change le nombre de compartiments supérieurs à 256, l'exécution du programme prend plus de temps qu'avec ce nombre de compartiments.

EDIT: Le même algorithme, langage de programmation différent. J'ai utilisé C ++ à la place de Java et le temps d'exécution a été réduit de ~ 3.2s à ~ 2.35s sur ma machine. Le nombre optimal de compartiments est toujours autour de 256 (encore une fois, sur mon ordinateur).

À propos, TBB est vraiment génial.

EDIT: Je me suis inspiré de l'excellente solution d'Alexandru et j'ai remplacé la dernière fois std :: sort par une version modifiée de sa sorte radix. J'ai utilisé une méthode différente pour traiter les nombres positifs / négatifs, même s'il faut plus de passages dans le tableau. J'ai également décidé de trier le tableau exactement et de supprimer le tri d'insertion. Je passerai plus tard un peu de temps à tester l’influence de ces modifications sur les performances et éventuellement leur annulation. Cependant, en utilisant le tri de base, le temps a été réduit d’environ 2,35 à environ 1,63.

Sans être malin, pour vous fournir un trieur naïf beaucoup plus rapide, en voici un en C qui devrait être à peu près équivalent à votre trieur Python:

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
int cmp(const void* av, const void* bv) {
    double a = *(const double*)av;
    double b = *(const double*)bv;
    return a < b ? -1 : a > b ? 1 : 0;
}
int main(int argc, char** argv) {
    if (argc <= 1)
        return puts("No argument!");
    unsigned count = atoi(argv[1]);

    double *a = malloc(count * sizeof *a);

    FILE *f = fopen("gaussian.dat", "rb");
    if (fread(a, sizeof *a, count, f) != count)
        return puts("fread failed!");
    fclose(f);

    puts("sorting...");
    double *b = malloc(count * sizeof *b);
    memcpy(b, a, count * sizeof *b);
    qsort(b, count, sizeof *b, cmp);
    return 0;
}

Compilé avec gcc -O3, sur ma machine, cela prend plus d’une minute de moins que le Python: environ 11 secondes contre 87 secondes.

J'ai partitionné en segments en fonction de l'écart type qui devrait le mieux être divisé en 4ème. Edit: réécrit sur la partition en fonction de la valeur x dans http://en.wikipedia.org/wiki/Error_function#Table_of_values

http://www.wolframalpha.com/input/?i=percentages+by++normal+distribution

J'ai essayé d'utiliser des seaux plus petits, mais cela semblait avoir peu d'effet une fois 2 * sur le nombre de cœurs disponibles. Sans collections parallèles, cela prendrait 37 secondes sur ma boîte et 24 avec les collections parallèles. En cas de partitionnement via une distribution, vous ne pouvez pas simplement utiliser un tableau, il y a donc un surplus de temps système. Je ne suis pas clair sur le moment où une valeur serait boxed / unboxed dans scala.

J'utilise scala 2.9 pour la collection parallèle. Vous pouvez simplement télécharger sa distribution tar.gz.

Pour compiler: scalac SortFile.scala (je viens de le copier directement dans le dossier scala / bin.

Pour exécuter: JAVA_OPTS = "- Xmx4096M" ./scala SortFile (je l'ai exécuté avec 2 Go de RAM et j'ai à peu près le même temps)

Edit: supprimé allocateDirect, plus lent que juste allouer. Suppression de l’amorçage de la taille initiale des tampons de tableau. En fait, il a lu toutes les valeurs 50000000. Réécrit pour éviter, espérons-le, les problèmes de substitution automatique (encore plus lent que naïf c)

import java.io.FileInputStream;
import java.nio.ByteBuffer
import java.nio.ByteOrder
import scala.collection.mutable.ArrayBuilder


object SortFile {

//used partition numbers from Damascus' solution
val partList = List(0, 0.15731, 0.31864, 0.48878, 0.67449, 0.88715, 1.1503, 1.5341)

val listSize = partList.size * 2;
val posZero = partList.size;
val neg = partList.map( _ * -1).reverse.zipWithIndex
val pos = partList.map( _ * 1).zipWithIndex.reverse

def partition(dbl:Double): Int = { 

//for each partition, i am running through the vals in order
//could make this a binary search to be more performant... but our list size is 4 (per side)

  if(dbl < 0) { return neg.find( dbl < _._1).get._2  }
  if(dbl > 0) { return posZero  + pos.find( dbl > _._1).get._2  }
      return posZero; 

}

  def main(args: Array[String])
    { 

    var l = 0
    val dbls = new Array[Double](50000000)
    val partList = new Array[Int](50000000)
    val pa = Array.fill(listSize){Array.newBuilder[Double]}
    val channel = new FileInputStream("../../gaussian.dat").getChannel()
    val bb = ByteBuffer.allocate(50000000 * 8)
    bb.order(ByteOrder.LITTLE_ENDIAN)
    channel.read(bb)
    bb.rewind
    println("Loaded" + System.currentTimeMillis())
    var dbl = 0.0
    while(bb.hasRemaining)
    { 
      dbl = bb.getDouble
      dbls.update(l,dbl) 

      l+=1
    }
    println("Beyond first load" + System.currentTimeMillis());

    for( i <- (0 to 49999999).par) { partList.update(i, partition(dbls(i)))}

    println("Partition computed" + System.currentTimeMillis() )
    for(i <- (0 to 49999999)) { pa(partList(i)) += dbls(i) }
    println("Partition completed " + System.currentTimeMillis())
    val toSort = for( i <- pa) yield i.result()
    println("Arrays Built" + System.currentTimeMillis());
    toSort.par.foreach{i:Array[Double] =>scala.util.Sorting.quickSort(i)};

    println("Read\t" + System.currentTimeMillis());

  }
}

Il suffit de mettre cela dans un fichier cs et de le compiler avec csc en théorie: (nécessite mono)

using System;
using System.IO;
using System.Threading;

namespace Sort
{
    class Program
    {
        const int count = 50000000;
        static double[][] doubles;
        static WaitHandle[] waiting = new WaitHandle[4];
        static AutoResetEvent[] events = new AutoResetEvent[4];

        static double[] Merge(double[] left, double[] right)
        {
            double[] result = new double[left.Length + right.Length];
            int l = 0, r = 0, spot = 0;
            while (l < left.Length && r < right.Length)
            {
                if (right[r] < left[l])
                    result[spot++] = right[r++];
                else
                    result[spot++] = left[l++];
            }
            while (l < left.Length) result[spot++] = left[l++];
            while (r < right.Length) result[spot++] = right[r++];
            return result;
        }

        static void ThreadStart(object data)
        {
            int index = (int)data;
            Array.Sort(doubles[index]);
            events[index].Set();
        }

        static void Main(string[] args)
        {
            System.Diagnostics.Stopwatch watch = new System.Diagnostics.Stopwatch();
            watch.Start();
            byte[] bytes = File.ReadAllBytes(@"..\..\..\SortGuassian\Data.dat");
            doubles = new double[][] { new double[count / 4], new double[count / 4], new double[count / 4], new double[count / 4] };
            for (int i = 0; i < 4; i++)
            {
                for (int j = 0; j < count / 4; j++)
                {
                    doubles[i][j] = BitConverter.ToDouble(bytes, i * count/4 + j * 8);
                }
            }
            Thread[] threads = new Thread[4];
            for (int i = 0; i < 4; i++)
            {
                threads[i] = new Thread(ThreadStart);
                waiting[i] = events[i] = new AutoResetEvent(false);
                threads[i].Start(i);
            }
            WaitHandle.WaitAll(waiting);
            double[] left = Merge(doubles[0], doubles[1]);
            double[] right = Merge(doubles[2], doubles[3]);
            double[] result = Merge(left, right);
            watch.Stop();
            Console.WriteLine(watch.Elapsed.ToString());
            Console.ReadKey();
        }
    }
}

Puisque vous savez quelle est la distribution, vous pouvez utiliser un tri O (N) à indexation directe. (Si vous vous demandez ce que c'est, supposons que vous avez un paquet de 52 cartes et que vous voulez le trier. Ayez juste 52 bacs et jetez chaque carte dans son propre bac.)

Vous avez 5e7 doubles. Allouer un tableau de résultats R de 5e7 doubles. Prenez chaque nombre xet obtenez i = phi(x) * 5e7. Fondamentalement faire R[i] = x. Avoir un moyen de gérer les collisions, par exemple en déplaçant le nombre avec lequel elle entre en collision (comme dans le codage de hachage simple). Alternativement, vous pouvez rendre R plusieurs fois plus grand, avec une valeur vide unique . À la fin, vous venez de balayer les éléments de R.

phiest juste la fonction de distribution cumulative gaussienne. Il convertit un nombre distribué gaussien entre +/- infini en un nombre distribué uniforme entre 0 et 1. Un moyen simple de le calculer consiste à rechercher dans une table et à effectuer une interpolation.

Voici une autre solution séquentielle:

#include <stdio.h>
#include <stdlib.h>
#include <algorithm>
#include <ctime>

typedef unsigned long long ull;

int size;
double *dbuf, *copy;
int cnt[8][1 << 16];

void sort()
{
  const int step = 10;
  const int start = 24;
  ull mask = (1ULL << step) - 1;

  ull *ibuf = (ull *) dbuf;
  for (int i = 0; i < size; i++) {
    for (int w = start, v = 0; w < 64; w += step, v++) {
      int p = (~ibuf[i] >> w) & mask;
      cnt[v][p]++;
    }
  }

  int sum[8] = { 0 };
  for (int i = 0; i <= mask; i++) {
    for (int w = start, v = 0; w < 64; w += step, v++) {
      int tmp = sum[v] + cnt[v][i];
      cnt[v][i] = sum[v];
      sum[v] = tmp;
    }
  }

  for (int w = start, v = 0; w < 64; w += step, v++) {
    ull *ibuf = (ull *) dbuf;
    for (int i = 0; i < size; i++) {
      int p = (~ibuf[i] >> w) & mask;
      copy[cnt[v][p]++] = dbuf[i];
    }

    double *tmp = copy;
    copy = dbuf;
    dbuf = tmp;
  }

  for (int p = 0; p < size; p++)
    if (dbuf[p] >= 0.) {
      std::reverse(dbuf + p, dbuf + size);
      break;
    }

  // Insertion sort
  for (int i = 1; i < size; i++) {
    double value = dbuf[i];
    if (value < dbuf[i - 1]) {
      dbuf[i] = dbuf[i - 1];
      int p = i - 1;
      for (; p > 0 && value < dbuf[p - 1]; p--)
        dbuf[p] = dbuf[p - 1];
      dbuf[p] = value;
    }
  }
}

int main(int argc, char **argv) {
  size = atoi(argv[1]);
  dbuf = new double[size];
  copy = new double[size];

  FILE *f = fopen("gaussian.dat", "r");
  fread(dbuf, size, sizeof(double), f);
  fclose(f);

  clock_t c0 = clock();
  sort();
  printf("Finished after %.3f\n", (double) ((clock() - c0)) / CLOCKS_PER_SEC);
  return 0;
}

Je doute que cela bat la solution multi-thread, mais les temps sur mon ordinateur portable i7 sont (stdsort est la solution C ++ fournie dans une autre réponse):

$ g++ -O3 mysort.cpp -o mysort && ./mysort 50000000
Finished after 2.10
$ g++ -O3 stdsort.cpp -o stdsort && ./stdsort
Finished after 7.12

Notez que cette solution a une complexité temporelle linéaire (car elle utilise la représentation spéciale des doubles).

EDIT : Correction de l'ordre des éléments à augmenter.

EDIT : amélioration de la vitesse de presque une demi-seconde.

EDIT : amélioration de la vitesse de 0,7 seconde. A rendu l'algorithme plus convivial en cache.

EDIT : amélioration de la vitesse d'une seconde. Comme il n’ya que 50.000.000 d’éléments, je peux partiellement trier la mantisse et utiliser le type insert (qui respecte le cache) pour réparer les éléments déplacés. Cette idée supprime environ deux itérations de la dernière boucle de tri de base.

EDIT : 0,16 secondes en moins. Premièrement, std :: reverse peut être éliminé si l'ordre de tri est inversé.

Prendre la solution de Christian Ammer et la mettre en parallèle avec les blocs de construction filetés d'Intel

#include <iostream>
#include <fstream>
#include <algorithm>
#include <vector>
#include <ctime>
#include <tbb/parallel_sort.h>

int main(void)
{
    std::ifstream ifs("gaussian.dat", std::ios::binary | std::ios::in);
    std::vector<double> values;
    values.reserve(50000000);
    double d;
    while (ifs.read(reinterpret_cast<char*>(&d), sizeof(double)))
    values.push_back(d);
    clock_t c0 = clock();
    tbb::parallel_sort(values.begin(), values.end());
    std::cout << "Finished after "
              << static_cast<double>((clock() - c0)) / CLOCKS_PER_SEC
              << std::endl;
}

Si vous avez accès à la bibliothèque Performance Primitives (IPP) d'Intel, vous pouvez utiliser son type de base. Il suffit de remplacer

#include <tbb/parallel_sort.h>

avec

#include "ipps.h"

et

tbb::parallel_sort(values.begin(), values.end());

avec

std::vector<double> copy(values.size());
ippsSortRadixAscend_64f_I(&values[0], &copy[0], values.size());

Sur mon portable dual core, les timings sont

C               16.4 s
C#              20 s
C++ std::sort   7.2 s
C++ tbb         5 s
C++ ipp         4.5 s
python          too long

Que diriez-vous d’une implémentation de type quicksort parallèle qui choisisse ses valeurs de pivot en fonction des statistiques de la distribution, assurant ainsi des partitions de taille égale? Le premier pivot serait à la moyenne (zéro dans ce cas), la paire suivante aux 25ème et 75ème centiles (+/- -0,67499 écarts-types), et ainsi de suite, chaque partition divisant par deux le jeu de données restant plus ou moins moins parfaitement.

Très moche (pourquoi utiliser des tableaux quand je peux utiliser des variables finissant par des nombres), mais un code rapide (mon premier essai avec std :: threads), temps entier (temps réel) sur mon système 1,8 s (en comparaison avec std :: sort () 4,8 s), compilez avec g ++ -std = c ++ 0x -O3 -march = natif -pthread Passez simplement des données via stdin (fonctionne uniquement pour 50M).

#include <cstdio>
#include <iostream>
#include <algorithm>
#include <thread>
using namespace std;
const size_t size=50000000;

void pivot(double* start,double * end, double middle,size_t& koniec){
    double * beg=start;
    end--;
    while (start!=end){
        if (*start>middle) swap (*start,*end--);
        else start++;
    }
    if (*end<middle) start+=1;
    koniec= start-beg;
}
void s(double * a, double* b){
    sort(a,b);
}
int main(){
    double *data=new double[size];
    FILE *f = fopen("gaussian.dat", "rb");
    fread(data,8,size,f);
    size_t end1,end2,end3,temp;
    pivot(data, data+size,0,end2);
    pivot(data, data+end2,-0.6745,end1);
    pivot(data+end2,data+size,0.6745,end3);
    end3+=end2;
    thread ts1(s,data,data+end1);
    thread ts2(s,data+end1,data+end2);
    thread ts3(s,data+end2,data+end3);
    thread ts4(s,data+end3,data+size);
    ts1.join(),ts2.join(),ts3.join(),ts4.join();
    //for (int i=0; i<size-1; i++){
    //  if (data[i]>data[i+1]) cerr<<"BLAD\n";
    //}
    fclose(f);
    //fwrite(data,8,size,stdout);
}

// Edit modifié pour lire le fichier gaussian.dat.

Une solution C ++ utilisant std::sort(éventuellement plus rapide que qsort, en ce qui concerne les performances de qsort vs std :: sort )

#include <iostream>
#include <fstream>
#include <algorithm>
#include <vector>
#include <ctime>

int main(void)
{
    std::ifstream ifs("C:\\Temp\\gaussian.dat", std::ios::binary | std::ios::in);
    std::vector<double> values;
    values.reserve(50000000);
    double d;
    while (ifs.read(reinterpret_cast<char*>(&d), sizeof(double)))
        values.push_back(d);
    clock_t c0 = clock();
    std::sort(values.begin(), values.end());
    std::cout << "Finished after "
              << static_cast<double>((clock() - c0)) / CLOCKS_PER_SEC
              << std::endl;
}

Je ne peux pas dire gaussian.datavec certitude combien de temps cela prend car je n'ai que 1 Go sur ma machine et avec le code Python donné, je ne pouvais créer qu'un fichier avec seulement 25 millions de doublons (sans erreur de mémoire). Mais je suis très intéressé par la durée d'exécution de l'algorithme std :: sort.

Voici un mélange de la sorte de radix d'Alexandru avec le pivot intelligent de Zjarek. Compilez-le avec

g++ -std=c++0x -pthread -O3 -march=native sorter_gaussian_radix.cxx -o sorter_gaussian_radix

Vous pouvez changer la taille de la base en définissant STEP (par exemple, ajoutez -DSTEP = 11). J'ai trouvé le meilleur pour mon ordinateur portable est 8 (par défaut).

Par défaut, il divise le problème en 4 morceaux et l'exécute sur plusieurs threads. Vous pouvez changer cela en passant un paramètre de profondeur à la ligne de commande. Donc, si vous avez deux noyaux, lancez-le en tant que

sorter_gaussian_radix 50000000 1

et si vous avez 16 noyaux

sorter_gaussian_radix 50000000 4

La profondeur maximale actuelle est de 6 (64 fils). Si vous mettez trop de niveaux, vous ralentirez simplement le code.

Une chose que j'ai également essayée était le type de base de la bibliothèque Intel Performance Primitives (IPP). La mise en œuvre d'Alexandru va à l'encontre de l'IPP, qui est environ 30% plus lent. Cette variation est également incluse ici (commentée).

#include <stdlib.h>
#include <string.h>
#include <algorithm>
#include <ctime>
#include <iostream>
#include <thread>
#include <vector>
#include <boost/cstdint.hpp>
// #include "ipps.h"

#ifndef STEP
#define STEP 8
#endif

const int step = STEP;
const int start_step=24;
const int num_steps=(64-start_step+step-1)/step;
int size;
double *dbuf, *copy;

clock_t c1, c2, c3, c4, c5;

const double distrib[]={-2.15387,
                        -1.86273,
                        -1.67594,
                        -1.53412,
                        -1.4178,
                        -1.31801,
                        -1.22986,
                        -1.15035,
                        -1.07752,
                        -1.00999,
                        -0.946782,
                        -0.887147,
                        -0.830511,
                        -0.776422,
                        -0.724514,
                        -0.67449,
                        -0.626099,
                        -0.579132,
                        -0.53341,
                        -0.488776,
                        -0.445096,
                        -0.40225,
                        -0.36013,
                        -0.318639,
                        -0.27769,
                        -0.237202,
                        -0.197099,
                        -0.157311,
                        -0.11777,
                        -0.0784124,
                        -0.0391761,
                        0,
                        0.0391761,
                        0.0784124,
                        0.11777,
                        0.157311,
                        0.197099,
                        0.237202,
                        0.27769,
                        0.318639,
                        0.36013,
                        0.40225,
                        0.445097,
                        0.488776,
                        0.53341,
                        0.579132,
                        0.626099,
                        0.67449,
                        0.724514,
                        0.776422,
                        0.830511,
                        0.887147,
                        0.946782,
                        1.00999,
                        1.07752,
                        1.15035,
                        1.22986,
                        1.31801,
                        1.4178,
                        1.53412,
                        1.67594,
                        1.86273,
                        2.15387};


class Distrib
{
  const int value;
public:
  Distrib(const double &v): value(v) {}

  bool operator()(double a)
  {
    return a<value;
  }
};


void recursive_sort(const int start, const int end,
                    const int index, const int offset,
                    const int depth, const int max_depth)
{
  if(depth<max_depth)
    {
      Distrib dist(distrib[index]);
      const int middle=std::partition(dbuf+start,dbuf+end,dist) - dbuf;

      // const int middle=
      //   std::partition(dbuf+start,dbuf+end,[&](double a)
      //                  {return a<distrib[index];})
      //   - dbuf;

      std::thread lower(recursive_sort,start,middle,index-offset,offset/2,
                        depth+1,max_depth);
      std::thread upper(recursive_sort,middle,end,index+offset,offset/2,
                        depth+1,max_depth);
      lower.join(), upper.join();
    }
  else
    {
  // ippsSortRadixAscend_64f_I(dbuf+start,copy+start,end-start);

      c1=clock();

      double *dbuf_local(dbuf), *copy_local(copy);
      boost::uint64_t mask = (1 << step) - 1;
      int cnt[num_steps][mask+1];

      boost::uint64_t *ibuf = reinterpret_cast<boost::uint64_t *> (dbuf_local);

      for(int i=0;i<num_steps;++i)
        for(uint j=0;j<mask+1;++j)
          cnt[i][j]=0;

      for (int i = start; i < end; i++)
        {
          for (int w = start_step, v = 0; w < 64; w += step, v++)
            {
              int p = (~ibuf[i] >> w) & mask;
              (cnt[v][p])++;
            }
        }

      c2=clock();

      std::vector<int> sum(num_steps,0);
      for (uint i = 0; i <= mask; i++)
        {
          for (int w = start_step, v = 0; w < 64; w += step, v++)
            {
              int tmp = sum[v] + cnt[v][i];
              cnt[v][i] = sum[v];
              sum[v] = tmp;
            }
        }

      c3=clock();

      for (int w = start_step, v = 0; w < 64; w += step, v++)
        {
          ibuf = reinterpret_cast<boost::uint64_t *>(dbuf_local);

          for (int i = start; i < end; i++)
            {
              int p = (~ibuf[i] >> w) & mask;
              copy_local[start+((cnt[v][p])++)] = dbuf_local[i];
            }
          std::swap(copy_local,dbuf_local);
        }

      // Do the last set of reversals
      for (int p = start; p < end; p++)
        if (dbuf_local[p] >= 0.)
          {
            std::reverse(dbuf_local+p, dbuf_local + end);
            break;
          }

      c4=clock();

      // Insertion sort
      for (int i = start+1; i < end; i++) {
        double value = dbuf_local[i];
        if (value < dbuf_local[i - 1]) {
          dbuf_local[i] = dbuf_local[i - 1];
          int p = i - 1;
          for (; p > 0 && value < dbuf_local[p - 1]; p--)
            dbuf_local[p] = dbuf_local[p - 1];
          dbuf_local[p] = value;
        }
      }
      c5=clock();

    }
}


int main(int argc, char **argv) {
  size = atoi(argv[1]);
  copy = new double[size];

  dbuf = new double[size];
  FILE *f = fopen("gaussian.dat", "r");
  fread(dbuf, size, sizeof(double), f);
  fclose(f);

  clock_t c0 = clock();

  const int max_depth= (argc > 2) ? atoi(argv[2]) : 2;

  // ippsSortRadixAscend_64f_I(dbuf,copy,size);

  recursive_sort(0,size,31,16,0,max_depth);

  if(num_steps%2==1)
    std::swap(dbuf,copy);

  // for (int i=0; i<size-1; i++){
  //   if (dbuf[i]>dbuf[i+1])
  //     std::cout << "BAD "
  //               << i << " "
  //               << dbuf[i] << " "
  //               << dbuf[i+1] << " "
  //               << "\n";
  // }

  std::cout << "Finished after "
            << (double) (c1 - c0) / CLOCKS_PER_SEC << " "
            << (double) (c2 - c1) / CLOCKS_PER_SEC << " "
            << (double) (c3 - c2) / CLOCKS_PER_SEC << " "
            << (double) (c4 - c3) / CLOCKS_PER_SEC << " "
            << (double) (c5 - c4) / CLOCKS_PER_SEC << " "
            << "\n";

  // delete [] dbuf;
  // delete [] copy;
  return 0;
}

EDIT : J'ai mis en œuvre les améliorations du cache d'Alexandru, qui ont permis de réduire d'environ 30% le temps passé sur ma machine.

EDIT : Ceci implémente un tri récursif, donc cela devrait bien fonctionner sur la machine à 16 cœurs d’Alexandru. Il utilise également la dernière amélioration d'Alexandru et supprime l'un des inverses. Pour moi, cela a donné une amélioration de 20%.

EDIT : Correction d'un bug de signe qui causait une inefficacité lorsqu'il y a plus de 2 cœurs.

EDIT : Suppression du lambda, donc il sera compilé avec les anciennes versions de gcc. Il inclut la variation de code IPP commentée. J'ai aussi corrigé la documentation pour courir sur 16 cœurs. Autant que je sache, c'est la mise en œuvre la plus rapide.

EDIT : Correction d'un bug quand STEP n'est pas 8. Augmentation du nombre maximum de threads à 64. Ajout de quelques informations de timing.

Je suppose que cela dépend vraiment de ce que vous voulez faire. Si vous voulez trier un groupe de Gaussiens, cela ne vous aidera pas. Mais si vous voulez un groupe de Gaussiens triés, ça ira. Même si cela manque un peu le problème, je pense qu'il sera intéressant de comparer les routines de tri réelles.

Si vous voulez quelque chose pour être rapide, faites-en moins.

Au lieu de générer un groupe d'échantillons aléatoires à partir de la distribution normale, puis de les trier, vous pouvez générer un groupe d'échantillons de la distribution normale dans un ordre trié.

Vous pouvez utiliser la solution ici pour générer n nombres aléatoires uniformes dans un ordre trié. Ensuite, vous pouvez utiliser le cdf inverse (scipy.stats.norm.ppf) de la distribution normale pour transformer les nombres aléatoires uniformes en nombres de la distribution normale via un échantillonnage par transformée inverse .

import scipy.stats
import random

# slightly modified from linked stackoverflow post
def n_random_numbers_increasing(n):
  """Like sorted(random() for i in range(n))),                                
  but faster because we avoid sorting."""
  v = 1.0
  while n:
    v *= random.random() ** (1.0 / n)
    yield 1 - v
    n -= 1

def n_normal_samples_increasing(n):
  return map(scipy.stats.norm.ppf, n_random_numbers_increasing(n))

Si vous voulez avoir les mains plus sales, je suppose que vous pourrez peut-être accélérer les nombreux calculs cdf inverses en utilisant une sorte de méthode itérative et en utilisant le résultat précédent comme première estimation. Puisque les suppositions vont être très proches, une simple itération vous donnera probablement une grande précision.

Essayez cette solution changeante de Guvante avec ce Main (), il commence à trier dès que la lecture 1/4 IO est terminée, le test est plus rapide:

    static void Main(string[] args)
    {
        FileStream filestream = new FileStream(@"..\..\..\gaussian.dat", FileMode.Open, FileAccess.Read);
        doubles = new double[][] { new double[count / 4], new double[count / 4], new double[count / 4], new double[count / 4] };
        Thread[] threads = new Thread[4];

        for (int i = 0; i < 4; i++)
        {
            byte[] bytes = new byte[count * 4];
            filestream.Read(bytes, 0, count * 4);

            for (int j = 0; j < count / 4; j++)
            {
                doubles[i][j] = BitConverter.ToDouble(bytes, i * count/4 + j * 8);
            }

            threads[i] = new Thread(ThreadStart);
            waiting[i] = events[i] = new AutoResetEvent(false);
            threads[i].Start(i);    
        }

        WaitHandle.WaitAll(waiting);
        double[] left = Merge(doubles[0], doubles[1]);
        double[] right = Merge(doubles[2], doubles[3]);
        double[] result = Merge(left, right);
        Console.ReadKey();
    }
}

Puisque vous connaissez la distribution, mon idée serait de créer k compartiments, chacun avec le même nombre d’éléments escompté (puisque vous connaissez la distribution, vous pouvez la calculer). Puis, en temps O (n), balayez le tableau et mettez les éléments dans leurs compartiments.

Ensuite, triez simultanément les seaux. Supposons que vous avez k compartiments et n éléments. Un compartiment prendra (n / k) lg (n / k) pour trier. Supposons maintenant que vous ayez p processeurs utilisables. Étant donné que les compartiments peuvent être triés indépendamment, vous devez traiter un multiplicateur de ceil (k / p). Cela donne un temps d'exécution final de n + ceil (k / p) * (n / k) lg (n / k), ce qui devrait être bien plus rapide que n lg n si vous choisissez bien k.

Une idée d'optimisation de bas niveau consiste à insérer deux doubles dans un registre SSE, afin que chaque thread fonctionne avec deux éléments à la fois. Cela peut être compliqué à faire pour certains algorithmes.

Une autre chose à faire est de trier le tableau en fragments faciles à mettre en cache, puis de fusionner les résultats. Deux niveaux doivent être utilisés: par exemple, 4 Ko pour L1 puis 64 Ko pour L2.

Cela devrait être très convivial pour le cache, car le tri du compartiment ne sortira pas du cache et la fusion finale parcourra la mémoire de manière séquentielle.

De nos jours, le calcul est beaucoup moins cher que les accès en mémoire. Cependant, nous avons un grand nombre d'éléments, il est donc difficile de dire quelle est la taille du tableau lorsque le tri idiote pour le cache est plus lent que pour une version peu complexe sans cache.

Mais je ne fournirai pas d'implémentation de ce qui précède car je le ferais sous Windows (VC ++).

Voici une implémentation du tri du compartiment d'analyse linéaire. Je pense que c'est plus rapide que toutes les implémentations à un seul thread actuelles, à l'exception du tri de base. Le temps d'exécution prévu devrait être linéaire si l'estimation de la cdf est suffisamment précise (j'utilise une interpolation linéaire des valeurs trouvées sur le Web) et je n'ai commis aucune erreur susceptible de provoquer un balayage excessif:

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <algorithm>
#include <ctime>

using std::fill;

const double q[] = {
  0.0,
  9.865E-10,
  2.8665150000000003E-7,
  3.167E-5,
  0.001349898,
  0.022750132,
  0.158655254,
  0.5,
  0.8413447460000001,
  0.9772498679999999,
  0.998650102,
  0.99996833,
  0.9999997133485,
  0.9999999990134999,
  1.0,
};
int main(int argc, char** argv) {
  if (argc <= 1)
    return puts("No argument!");
  unsigned count = atoi(argv[1]);
  unsigned count2 = 3 * count;

  bool *ba = new bool[count2 + 1000];
  fill(ba, ba + count2 + 1000, false);
  double *a = new double[count];
  double *c = new double[count2 + 1000];

  FILE *f = fopen("gaussian.dat", "rb");
  if (fread(a, 8, count, f) != count)
    return puts("fread failed!");
  fclose(f);

  int i;
  int j;
  bool s;
  int t;
  double z;
  double p;
  double d1;
  double d2;
  for (i = 0; i < count; i++) {
    s = a[i] < 0;
    t = a[i];
    if (s) t--;
    z = a[i] - t;
    t += 7;
    if (t < 0) {
      t = 0;
      z = 0;
    } else if (t >= 14) {
      t = 13;
      z = 1;
    }
    p = q[t] * (1 - z) + q[t + 1] * z;
    j = count2 * p;
    while (ba[j] && c[j] < a[i]) {
      j++;
    }
    if (!ba[j]) {
      ba[j] = true;
      c[j] = a[i];
    } else {
      d1 = c[j];
      c[j] = a[i];
      j++;
      while (ba[j]) {
        d2 = c[j];
        c[j] = d1;
        d1 = d2;
        j++;
      }
      c[j] = d1;
      ba[j] = true;
    }
  }
  i = 0;
  int max = count2 + 1000;
  for (j = 0; j < max; j++) {
    if (ba[j]) {
      a[i++] = c[j];
    }
  }
  // for (i = 0; i < count; i += 1) {
  //   printf("here %f\n", a[i]);
  // }
  return 0;
}

Je ne sais pas, pourquoi je ne peux pas éditer mon précédent post, alors voici la nouvelle version, 0,2 seconde plus rapide (mais environ 1,5 s plus rapide en temps CPU (utilisateur)). Cette solution a 2 programmes, calcule d'abord les quantiles pour une distribution normale pour le tri par compartiment et les stocke dans un tableau, t [double * scale] = index du compartiment, où scale est un nombre arbitraire permettant de doubler le moulage. Ensuite, le programme principal peut utiliser ces données pour placer les doublons dans le compartiment approprié. Cela présente un inconvénient: si les données ne sont pas gaussiennes, elles ne fonctionneront pas correctement (et il n’ya presque aucune chance de fonctionner de manière incorrecte pour une distribution normale), mais la modification pour un cas particulier est facile et rapide ::Trier()).

Compilation: g ++ => http://pastebin.com/WG7pZEzH programme d'aide

g ++ -std = c ++ 0x -O3 -march = natif -pthread => http://pastebin.com/T3yzViZP programme de tri principal

Voici une autre solution séquentielle. Celui-ci utilise le fait que les éléments sont distribués normalement, et je pense que l'idée est généralement applicable pour obtenir un tri proche du temps linéaire.

L'algorithme est comme ça:

• CDF approximatif (voir phi()fonction dans l'implémentation)
• Pour tous les éléments, calculez la position approximative dans le tableau trié: size * phi(x)
• Placez les éléments dans un nouveau tableau proche de leur position finale
  • Dans mon implémentation, le tableau de destination comporte des lacunes, de sorte que je n'ai pas à déplacer trop d'éléments lors de l'insertion.
• Utilisez insertsort pour trier les éléments finaux (insertsort est linéaire si la distance par rapport à la position finale est inférieure à une constante).

Malheureusement, la constante cachée est assez grande et cette solution est deux fois plus lente que l'algorithme de tri de base.

Mon jeu personnel favori utilisant les blocs de construction filetés d'Intel a déjà été publié, mais voici une solution parallèle grossière utilisant JDK 7 et sa nouvelle API fork / join:

import java.io.FileInputStream;
import java.nio.channels.FileChannel;
import java.util.Arrays;
import java.util.concurrent.*;
import static java.nio.channels.FileChannel.MapMode.READ_ONLY;
import static java.nio.ByteOrder.LITTLE_ENDIAN;


/**
 * 
 * Original Quicksort: https://github.com/pmbauer/parallel/tree/master/src/main/java/pmbauer/parallel
 *
 */
public class ForkJoinQuicksortTask extends RecursiveAction {

    public static void main(String[] args) throws Exception {

        double[] array = new double[Integer.valueOf(args[0])];

        FileChannel fileChannel = new FileInputStream("gaussian.dat").getChannel();
        fileChannel.map(READ_ONLY, 0, fileChannel.size()).order(LITTLE_ENDIAN).asDoubleBuffer().get(array);

        ForkJoinPool mainPool = new ForkJoinPool();

        System.out.println("Starting parallel computation");

        mainPool.invoke(new ForkJoinQuicksortTask(array));        
    }

    private static final long serialVersionUID = -642903763239072866L;
    private static final int SERIAL_THRESHOLD = 0x1000;

    private final double a[];
    private final int left, right;

    public ForkJoinQuicksortTask(double[] a) {this(a, 0, a.length - 1);}

    private ForkJoinQuicksortTask(double[] a, int left, int right) {
        this.a = a;
        this.left = left;
        this.right = right;
    }

    @Override
    protected void compute() {
        if (right - left < SERIAL_THRESHOLD) {
            Arrays.sort(a, left, right + 1);
        } else {
            int pivotIndex = partition(a, left, right);
            ForkJoinTask<Void> t1 = null;

            if (left < pivotIndex)
                t1 = new ForkJoinQuicksortTask(a, left, pivotIndex).fork();
            if (pivotIndex + 1 < right)
                new ForkJoinQuicksortTask(a, pivotIndex + 1, right).invoke();

            if (t1 != null)
                t1.join();
        }
    }

    public static int partition(double[] a, int left, int right) {
        // chose middle value of range for our pivot
        double pivotValue = a[left + (right - left) / 2];

        --left;
        ++right;

        while (true) {
            do
                ++left;
            while (a[left] < pivotValue);

            do
                --right;
            while (a[right] > pivotValue);

            if (left < right) {
                double tmp = a[left];
                a[left] = a[right];
                a[right] = tmp;
            } else {
                return right;
            }
        }
    }    
}

Avertissement important : J'ai pris l'adaptation de tri rapide pour fork / rejoindre de: https://github.com/pmbauer/parallel/tree/master/src/main/java/pmbauer/parallel

Pour exécuter cela, vous avez besoin d’une version bêta de JDK 7 (http://jdk7.java.net/download.html).

Sur mon Core i7 2.93Ghz Quad (OS X):

Référence python

time python sort.py 50000000
sorting...

real    1m13.885s
user    1m11.942s
sys     0m1.935s

Java JDK 7 fork / join

time java ForkJoinQuicksortTask 50000000
Starting parallel computation

real    0m2.404s
user    0m10.195s
sys     0m0.347s

J'ai aussi essayé de faire des essais de lecture en parallèle et de convertir les octets en doubles, mais je n'y voyais aucune différence.

Mise à jour:

Si quelqu'un veut expérimenter le chargement parallèle des données, la version de chargement parallèle est ci-dessous. En théorie, cela pourrait accélérer encore un peu si votre périphérique IO dispose de suffisamment de capacité en parallèle (les disques SSD le sont généralement). La création de doubles à partir d'octets entraîne également des frais généraux, de sorte que cela pourrait aussi aller plus vite en parallèle. Sur mes systèmes (SSD Ubuntu 10.10 / Nehalem Quad / Intel X25M et SSD OS X 10.6 / i7 Quad / Samsung), je ne voyais aucune différence réelle.

import static java.nio.ByteOrder.LITTLE_ENDIAN;
import static java.nio.channels.FileChannel.MapMode.READ_ONLY;

import java.io.FileInputStream;
import java.nio.DoubleBuffer;
import java.nio.channels.FileChannel;
import java.util.Arrays;
import java.util.concurrent.ForkJoinPool;
import java.util.concurrent.ForkJoinTask;
import java.util.concurrent.RecursiveAction;


/**
 *
 * Original Quicksort: https://github.com/pmbauer/parallel/tree/master/src/main/java/pmbauer/parallel
 *
 */
public class ForkJoinQuicksortTask extends RecursiveAction {

   public static void main(String[] args) throws Exception {

       ForkJoinPool mainPool = new ForkJoinPool();

       double[] array = new double[Integer.valueOf(args[0])];
       FileChannel fileChannel = new FileInputStream("gaussian.dat").getChannel();
       DoubleBuffer buffer = fileChannel.map(READ_ONLY, 0, fileChannel.size()).order(LITTLE_ENDIAN).asDoubleBuffer();

       mainPool.invoke(new ReadAction(buffer, array, 0, array.length));
       mainPool.invoke(new ForkJoinQuicksortTask(array));
   }

   private static final long serialVersionUID = -642903763239072866L;
   private static final int SERIAL_THRESHOLD = 0x1000;

   private final double a[];
   private final int left, right;

   public ForkJoinQuicksortTask(double[] a) {this(a, 0, a.length - 1);}

   private ForkJoinQuicksortTask(double[] a, int left, int right) {
       this.a = a;
       this.left = left;
       this.right = right;
   }

   @Override
   protected void compute() {
       if (right - left < SERIAL_THRESHOLD) {
           Arrays.sort(a, left, right + 1);
       } else {
           int pivotIndex = partition(a, left, right);
           ForkJoinTask<Void> t1 = null;

           if (left < pivotIndex)
               t1 = new ForkJoinQuicksortTask(a, left, pivotIndex).fork();
           if (pivotIndex + 1 < right)
               new ForkJoinQuicksortTask(a, pivotIndex + 1, right).invoke();

           if (t1 != null)
               t1.join();
       }
   }

   public static int partition(double[] a, int left, int right) {
       // chose middle value of range for our pivot
       double pivotValue = a[left + (right - left) / 2];

       --left;
       ++right;

       while (true) {
           do
               ++left;
           while (a[left] < pivotValue);

           do
               --right;
           while (a[right] > pivotValue);

           if (left < right) {
               double tmp = a[left];
               a[left] = a[right];
               a[right] = tmp;
           } else {
               return right;
           }
       }
   }

}

class ReadAction extends RecursiveAction {

   private static final long serialVersionUID = -3498527500076085483L;

   private final DoubleBuffer buffer;
   private final double[] array;
   private final int low, high;

   public ReadAction(DoubleBuffer buffer, double[] array, int low, int high) {
       this.buffer = buffer;
       this.array = array;
       this.low = low;
       this.high = high;
   }

   @Override
   protected void compute() {
       if (high - low < 100000) {
           buffer.position(low);
           buffer.get(array, low, high-low);
       } else {
           int middle = (low + high) >>> 1;

           invokeAll(new ReadAction(buffer.slice(), array, low, middle),  new ReadAction(buffer.slice(), array, middle, high));
       }
   }
}

Update2:

J'ai exécuté le code sur l'une de nos 12 machines de développement avec une légère modification pour définir un nombre fixe de cœurs. Cela a donné les résultats suivants:

Cores  Time
1      7.568s
2      3.903s
3      3.325s
4      2.388s
5      2.227s
6      1.956s
7      1.856s
8      1.827s
9      1.682s
10     1.698s
11     1.620s
12     1.503s

Sur ce système, j'ai également essayé la version Python qui prenait 1m2.994s et la version C ++ de Zjarek qui utilisait 1.925s (pour une raison quelconque, la version C ++ de Zjarek semble fonctionner relativement plus rapidement sur l'ordinateur de static_rtti).

J'ai aussi essayé ce qui se passait si je doublais la taille du fichier à 100 000 000 de fois:

Cores  Time
1      15.056s
2      8.116s
3      5.925s
4      4.802s
5      4.430s
6      3.733s
7      3.540s
8      3.228s
9      3.103s
10     2.827s
11     2.784s
12     2.689s

Dans ce cas, la version C ++ de Zjarek a pris 3.968. Python a juste pris trop de temps ici.

150 000 000 doubles:

Cores  Time
1      23.295s
2      12.391s
3      8.944s
4      6.990s
5      6.216s
6      6.211s
7      5.446s
8      5.155s
9      4.840s
10     4.435s
11     4.248s
12     4.174s

Dans ce cas, la version C ++ de Zjarek était 6.044s. Je n'ai même pas essayé Python.

La version C ++ est très cohérente avec ses résultats, où Java bascule un peu. Tout d'abord, le problème devient un peu plus efficace lorsque le problème s'aggrave, mais il est à nouveau moins efficace.

Une version utilisant des pthreads traditionnels. Code de fusion copié de la réponse de Guvante. Compiler avec g++ -O3 -pthread.

#include <stdio.h>
#include <stdlib.h>
#include <pthread.h>
#include <algorithm>

static unsigned int nthreads = 4;
static unsigned int size = 50000000;

typedef struct {
  double *array;
  int size;
} array_t;


void 
merge(double *left, int leftsize,
      double *right, int rightsize,
      double *result)
{
  int l = 0, r = 0, insertat = 0;
  while (l < leftsize && r < rightsize) {
    if (left[l] < right[r])
      result[insertat++] = left[l++];
    else
      result[insertat++] = right[r++];
  }

  while (l < leftsize) result[insertat++] = left[l++];
  while (r < rightsize) result[insertat++] = right[r++];
}


void *
run_thread(void *input)
{
  array_t numbers = *(array_t *)input;
  std::sort(numbers.array, numbers.array+numbers.size); 
  pthread_exit(NULL);
}

int 
main(int argc, char **argv) 
{
  double *numbers = (double *) malloc(size * sizeof(double));

  FILE *f = fopen("gaussian.dat", "rb");
  if (fread(numbers, sizeof(double), size, f) != size)
    return printf("Reading gaussian.dat failed");
  fclose(f);

  array_t worksets[nthreads];
  int worksetsize = size / nthreads;
  for (int i = 0; i < nthreads; i++) {
    worksets[i].array=numbers+(i*worksetsize);
    worksets[i].size=worksetsize;
  }

  pthread_attr_t attributes;
  pthread_attr_init(&attributes);
  pthread_attr_setdetachstate(&attributes, PTHREAD_CREATE_JOINABLE);

  pthread_t threads[nthreads];
  for (int i = 0; i < nthreads; i++) {
    pthread_create(&threads[i], &attributes, &run_thread, &worksets[i]);
  }

  for (int i = 0; i < nthreads; i++) {
    pthread_join(threads[i], NULL);
  }

  double *tmp = (double *) malloc(size * sizeof(double));
  merge(numbers, worksetsize, numbers+worksetsize, worksetsize, tmp);
  merge(numbers+(worksetsize*2), worksetsize, numbers+(worksetsize*3), worksetsize, tmp+(size/2));
  merge(tmp, worksetsize*2, tmp+(size/2), worksetsize*2, numbers);

  /*
  printf("Verifying result..\n");
  for (int i = 0; i < size - 1; i++) {
    if (numbers[i] > numbers[i+1])
      printf("Result is not correct\n");
  }
  */

  pthread_attr_destroy(&attributes);
  return 0;
}  

Sur mon ordinateur portable, j'obtiens les résultats suivants:

real    0m6.660s
user    0m9.449s
sys     0m1.160s

Voici une implémentation C99 séquentielle qui tente de vraiment utiliser la distribution connue. Il effectue essentiellement un seul tour de tri de seau en utilisant les informations de distribution, puis quelques tours de tri rapide sur chaque seau en supposant une distribution uniforme dans les limites du seau et enfin un tri de sélection modifié pour copier les données dans le tampon d'origine. Le tri rapide mémorise les points de division, de sorte que le tri par sélection ne doit fonctionner que sur des petits cunks. Et malgré (parce que?) De toute cette complexité, ce n'est même pas très rapide.

Pour rendre rapide l'évaluation, les valeurs sont échantillonnées en quelques points et, plus tard, seule une interpolation linéaire est utilisée. En fait, peu importe si Φ est évalué exactement, à condition que l'approximation soit strictement monotone.

Les tailles de bacs sont choisies de telle sorte que le risque de débordement du bac est négligeable. Plus précisément, avec les paramètres actuels, la probabilité qu’un jeu de données de 50000000 éléments provoque un débordement de bac est de 3.65e-09. (Ceci peut être calculé en utilisant la fonction de survie de la distribution de Poisson .)

Pour compiler, veuillez utiliser

gcc -std=c99 -msse3 -O3 -ffinite-math-only

Comme il y a beaucoup plus de calculs que dans les autres solutions, ces indicateurs de compilation sont nécessaires pour le rendre au moins raisonnablement rapide. Sans -msse3les conversions de doublepour intdevenir vraiment lent. Si votre architecture ne prend pas en charge SSE3, ces conversions peuvent également être effectuées à l'aide de la lrint()fonction.

Le code est plutôt moche - je ne sais pas si cela répond à l'exigence d'être "raisonnablement lisible" ...

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

#define N 50000000
#define BINSIZE 720
#define MAXBINSIZE 880
#define BINCOUNT (N / BINSIZE)
#define SPLITS 64
#define PHI_VALS 513

double phi_vals[PHI_VALS];

int bin_index(double x)
{
    double y = (x + 8.0) * ((PHI_VALS - 1) / 16.0);
    int interval = y;
    y -= interval;
    return (1.0 - y) * phi_vals[interval] + y * phi_vals[interval + 1];
}

double bin_value(int bin)
{
    int left = 0;
    int right = PHI_VALS - 1;
    do
    {
        int centre = (left + right) / 2;
        if (bin < phi_vals[centre])
            right = centre;
        else
            left = centre;
    } while (right - left > 1);
    double frac = (bin - phi_vals[left]) / (phi_vals[right] - phi_vals[left]);
    return (left + frac) * (16.0 / (PHI_VALS - 1)) - 8.0;
}

void gaussian_sort(double *restrict a)
{
    double *b = malloc(BINCOUNT * MAXBINSIZE * sizeof(double));
    double **pos = malloc(BINCOUNT * sizeof(double*));
    for (size_t i = 0; i < BINCOUNT; ++i)
        pos[i] = b + MAXBINSIZE * i;
    for (size_t i = 0; i < N; ++i)
        *pos[bin_index(a[i])]++ = a[i];
    double left_val, right_val = bin_value(0);
    for (size_t bin = 0, i = 0; bin < BINCOUNT; ++bin)
    {
        left_val = right_val;
        right_val = bin_value(bin + 1);
        double *splits[SPLITS + 1];
        splits[0] = b + bin * MAXBINSIZE;
        splits[SPLITS] = pos[bin];
        for (int step = SPLITS; step > 1; step >>= 1)
            for (int left_split = 0; left_split < SPLITS; left_split += step)
            {
                double *left = splits[left_split];
                double *right = splits[left_split + step] - 1;
                double frac = (double)(left_split + (step >> 1)) / SPLITS;
                double pivot = (1.0 - frac) * left_val + frac * right_val;
                while (1)
                {
                    while (*left < pivot && left <= right)
                        ++left;
                    while (*right >= pivot && left < right)
                        --right;
                    if (left >= right)
                        break;
                    double tmp = *left;
                    *left = *right;
                    *right = tmp;
                    ++left;
                    --right;
                }
                splits[left_split + (step >> 1)] = left;
            }
        for (int left_split = 0; left_split < SPLITS; ++left_split)
        {
            double *left = splits[left_split];
            double *right = splits[left_split + 1] - 1;
            while (left <= right)
            {
                double *min = left;
                for (double *tmp = left + 1; tmp <= right; ++tmp)
                    if (*tmp < *min)
                        min = tmp;
                a[i++] = *min;
                *min = *right--;
            }
        }
    }
    free(b);
    free(pos);
}

int main()
{
    double *a = malloc(N * sizeof(double));
    FILE *f = fopen("gaussian.dat", "rb");
    assert(fread(a, sizeof(double), N, f) == N);
    fclose(f);
    for (int i = 0; i < PHI_VALS; ++i)
    {
        double x = (i * (16.0 / PHI_VALS) - 8.0) / sqrt(2.0);
        phi_vals[i] =  (erf(x) + 1.0) * 0.5 * BINCOUNT;
    }
    gaussian_sort(a);
    free(a);
}

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <memory.h>
#include <algorithm>

// maps [-inf,+inf] to (0,1)
double normcdf(double x) {
        return 0.5 * (1 + erf(x * M_SQRT1_2));
}

int calcbin(double x, int bins) {
        return (int)floor(normcdf(x) * bins);
}

int *docensus(int bins, int n, double *arr) {
        int *hist = calloc(bins, sizeof(int));
        int i;
        for(i = 0; i < n; i++) {
                hist[calcbin(arr[i], bins)]++;
        }
        return hist;
}

void partition(int bins, int *orig_counts, double *arr) {
        int *counts = malloc(bins * sizeof(int));
        memcpy(counts, orig_counts, bins*sizeof(int));
        int *starts = malloc(bins * sizeof(int));
        int b, i;
        starts[0] = 0;
        for(i = 1; i < bins; i++) {
                starts[i] = starts[i-1] + counts[i-1];
        }
        for(b = 0; b < bins; b++) {
                while (counts[b] > 0) {
                        double v = arr[starts[b]];
                        int correctbin;
                        do {
                                correctbin = calcbin(v, bins);
                                int swappos = starts[correctbin];
                                double tmp = arr[swappos];
                                arr[swappos] = v;
                                v = tmp;
                                starts[correctbin]++;
                                counts[correctbin]--;
                        } while (correctbin != b);
                }
        }
        free(counts);
        free(starts);
}


void sortbins(int bins, int *counts, double *arr) {
        int start = 0;
        int b;
        for(b = 0; b < bins; b++) {
                std::sort(arr + start, arr + start + counts[b]);
                start += counts[b];
        }
}


void checksorted(double *arr, int n) {
        int i;
        for(i = 1; i < n; i++) {
                if (arr[i-1] > arr[i]) {
                        printf("out of order at %d: %lf %lf\n", i, arr[i-1], arr[i]);
                        exit(1);
                }
        }
}


int main(int argc, char *argv[]) {
        if (argc == 1 || argv[1] == NULL) {
                printf("Expected data size as argument\n");
                exit(1);
        }
        int n = atoi(argv[1]);
        const int cachesize = 128 * 1024; // a guess
        int bins = (int) (1.1 * n * sizeof(double) / cachesize);
        if (argc > 2) {
                bins = atoi(argv[2]);
        }
        printf("Using %d bins\n", bins);
        FILE *f = fopen("gaussian.dat", "rb");
        if (f == NULL) {
                printf("Couldn't open gaussian.dat\n");
                exit(1);
        }
        double *arr = malloc(n * sizeof(double));
        fread(arr, sizeof(double), n, f);
        fclose(f);

        int *counts = docensus(bins, n, arr);
        partition(bins, counts, arr);
        sortbins(bins, counts, arr);
        checksorted(arr, n);

        return 0;
}

Cela utilise erf () pour placer chaque élément de manière appropriée dans une corbeille, puis trie chaque corbeille. Il maintient le tableau entièrement en place.

Première passe: docensus () compte le nombre d'éléments dans chaque case.

Deuxième passage: partition () permute le tableau, plaçant chaque élément dans son propre bin

Troisième passage: sortbins () effectue un qsort sur chaque bin.

C'est un peu naïf, et appelle la coûteuse fonction erf () deux fois pour chaque valeur. Les premier et troisième passages sont potentiellement parallélisables. La seconde est hautement série et est probablement ralentie par ses modèles d’accès mémoire très aléatoires. Cela peut également valoir la peine de mettre en cache le numéro de bac de chaque double, en fonction du rapport entre la puissance du processeur et la vitesse de la mémoire.

Ce programme vous permet de choisir le nombre de bacs à utiliser. Ajoutez simplement un second numéro à la ligne de commande. Je l'ai compilé avec gcc -O3, mais ma machine est tellement faible que je ne peux pas vous dire de bons résultats.

Edit: Pouf! Mon programme C s'est transformé comme par magie en un programme C ++ utilisant std :: sort!

Jetez un coup d'œil à l'implémentation du tri de base par Michael Herf ( Radix Tricks ). Sur ma machine, le tri était 5 fois plus rapide que l' std::sortalgorithme de ma première réponse. Le nom de la fonction de tri est RadixSort11.

int main(void)
{
    std::ifstream ifs("C:\\Temp\\gaussian.dat", std::ios::binary | std::ios::in);
    std::vector<float> v;
    v.reserve(50000000);
    double d;
    while (ifs.read(reinterpret_cast<char*>(&d), sizeof(double)))
        v.push_back(static_cast<float>(d));
    std::vector<float> vres(v.size(), 0.0);
    clock_t c0 = clock();
    RadixSort11(&v[0], &vres[0], v.size());
    std::cout << "Finished after: "
              << static_cast<double>(clock() - c0) / CLOCKS_PER_SEC << std::endl;
    return 0;
}