Discussion :

[MatRem]

Bonjour,

Je lis un fichier binaire qui contient 3601*3601 altitudes codées sur 16 bits au format dted 2 data (environ 25 MO), et j'en fait la moyenne.

En python il faut environ 16 secondes. J'ai codé la même chose en c++ et il me faut 0.5 secondes soit 32 fois plus performant.

Est ce normal, ou est ce mon code python qui est mal écrit ?

python :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
path = "N47E006.hgt"
longitude_line_count = 3601
latitude_point_count = 3601
 
dt1_file = open( path, 'rb' )
 
mean_z = 0
 
for longitude_id in range( longitude_line_count ) :
	for latitude_id in range( latitude_point_count ) :
		value = dt1_file.read(2)
		high_byte = ord( value[0] )
		low_byte = ord( value[1] )
		sign = 1 - ( ( high_byte & 0x80 ) / 0x80 ) * 2
		high_byte = high_byte & 0x7f
		z = sign * ( high_byte * 0x100 + low_byte )
 
		mean_z += z
 
print( "mean_z      : " + str( float(mean_z) / float( longitude_line_count * latitude_point_count ) ) )
 
dt1_file.close()

c++

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
#include <fstream>
#include <iostream>
 
long read( std::ifstream & stream )
{
	char data[2];
	stream.read( data, 2 );
 
	long const sign = 1 - ( ( data[0] & 0x80 ) / 0x80 ) * 2;
	data[0] = data[0] & 0x7f;
 
	return sign * (long)( (data[0] * 0x100) | data[1] );
}
 
int main()
{
	std::ifstream file("N47E006.hgt", std::ios_base::in | std::ios_base::binary );
 
	size_t const longitude_count = 3601;
	size_t const latitude_count = 3601;
 
	long z_sum = 0;
 
	char data[2]; // hight byte, low byte
 
	for( size_t longitude = 0; longitude<longitude_count; ++longitude )
	{
		for( size_t latitude = 0; latitude<latitude_count; ++latitude )
		{
			z_sum += read( file );
		}
	}
 
	std::cout << "mean z      : " << static_cast<float>(z_sum) / static_cast<float>(longitude_count * latitude_count) << std::endl;
 
	file.close();
}

[MatRem]

ps : j'utilise python 2.7 sous ubuntu 11.10, et je ne maîtrise pas le python

[MatRem]

Suite à quelques conseils avisés il s'avère que la lecture du fichier n'a rien à voir avec le problème de performance. Il s'agit juste du très grand nombre d'instructions interprétées lors de l'itération.

En utilisant numpy on arrive au même performance que c++. Ouf ! je commençais à douter de mon utilisation futur de python .

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
import numpy
 
path = "N47E006.hgt"
longitude_line_count = 3601
latitude_point_count = 3601
 
dt1_file = open( path, 'rb' )
 
mean_z = 0
 
data = bytearray( dt1_file.read( longitude_line_count * latitude_point_count * 2 ) );
dt1_file.close()
 
data_array_byte = numpy.array( data, dtype=numpy.int16 )
 
high_byte = data_array_byte[::2]
low_byte = data_array_byte[1::2]
sign = 1 - ( ( high_byte & 0x80 ) / 0x80 ) * 2
high_byte = high_byte & 0x7f
z = sign * ( high_byte * 0x100 + low_byte )
mean_z = z.sum()
 
print( "mean_z      : " + str( float(mean_z) / float( longitude_line_count * latitude_point_count ) ) )

[dividee]

En fait, je pense que la lecture du fichier a beaucoup à voir avec les performances. Les performances du second code ne viennent pas tant de numpy que du fait de lire le fichier en un bloc, au lieu de lire 2 octets à la fois.

Si je lis bien, tu décodes des entiers signés sur 16 bits en complément à 2 ("signed short" en C) à partir des octets. C'est du travail inutile, tu peux directement lire les données dans le bon format avec le module array.

Je n'ai pas testé, mais cela devrait ressembler à ceci:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
import array
 
path = "N47E006.hgt"
longitude_line_count = 3601
latitude_point_count = 3601
 
short_array = array.array('h')      # 'h' = signed short
with open( path, 'rb' ) as dt1_file:
    short_array.fromfile(dt1_file, longitude_line_count * latitude_point_count)
#short_array.byteswap()   # si nécessaire, inverse les octets de poids fort et de poids faible
 
print( "mean_z      : " + str( float(sum(short_array)) / float( longitude_line_count * latitude_point_count ) ) )

[mont29]

Arf*! J’ai trop traîné dans mes tests, me suis fait grillé*!

Mais bon, voici les résultats de mes cogitations*:

D’abords, quel est exactement le codage de ces shorts*? Parce que apparemment, ils sont en big-endian (octet de poids fort d’abord), mais dans ce cas, ta méthode de décodage ne correspond pas au standard pour les entiers négatifs…

Ensuite, je suis d’accord avec dividee, le plus gros problème de performance vient des presque 13 millions d’accès en lecture*!

Voici mon code (compatible py2.7 et py3). r1() utilise le même principe que ton code, mais avec un décodage compatible avec le «*standard*» (i.e. short C). r2() utilise le module struct, qui permet de faire les conversions bytes → valeurs (et réciproquement), en procédant par gros blocs de 3601 éléments. Enfin, r3() fait de même, mais en une seule fois…

r4() est le code de dividee (le plus performant ), r5() est ton code avec numpy, et r6(), le même un poil optimisé

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
# -*- coding: utf-8 -*-
 
import struct
import time
import sys
 
is_py3k = sys.version_info.major == 3
 
def r1():
    num = 3601*3601
    sm = 0
 
    st = time.time()
    with open("test.hgt", "rb") as f:
        if is_py3k:
            for i in range(num):
                val = f.read(2)
                neg = (val[0] // 0x80)
                sm += val[0] * 0x100 + val[1] - (0x10000 if neg else 0)
        else:
            for i in range(num):
                val = f.read(2)
                neg = (ord(val[0]) // 0x80)
                sm += ord(val[0]) * 0x100 + ord(val[1]) - (0x10000 if neg else 0)
 
    print("mean_z      : {}".format(float(sm) / float(num)))
    print("r1 time: {} sec".format(time.time()-st))
 
def r2():
    num = 3601*3601
    col = 3601
    sm = 0
    fmt_half = ">"+"h"*col
 
    st = time.time()
    with open("test.hgt", "rb") as f:
        for i in range(col):
            sm += sum(struct.unpack(fmt_half, f.read(col*2)))
 
    print("mean_z      : {}".format(float(sm) / float(num)))
    print("r2 time: {} sec".format(time.time()-st))
 
def r3():
    num = 3601*3601
    sm = 0
    fmt_full = ">"+"h"*num
 
    st = time.time()
    with open("test.hgt", "rb") as f:
        sm = sum(struct.unpack(fmt_full, f.read(num*2)))
 
    print("mean_z      : {}".format(float(sm) / float(num)))
    print("r3 time: {} sec".format(time.time()-st))
 
def r4():
    import array
    num = 3601*3601
    sm = 0
 
    st = time.time()
    short_array = array.array('h')      # 'h' = signed short
    with open("test.hgt", 'rb') as f:
        short_array.fromfile(f, num)
    short_array.byteswap()   # si nécessaire, inverse les octets de poids fort et de poids faible
    sm = sum(short_array)
 
    print("mean_z      : {}".format(float(sm) / float(num)))
    print("r4 time: {} sec".format(time.time()-st))
 
def r5():
    import numpy
    num = 3601*3601
    sm = 0
 
    st = time.time()
    with open("test.hgt", 'rb') as f:
        data = bytearray(f.read(num*2));
    data_array_byte = numpy.array(data, dtype=numpy.int16)
 
    high_byte = data_array_byte[::2]
    low_byte = data_array_byte[1::2]
    sign = 1 - ((high_byte & 0x80) / 0x80) * 2
    high_byte = high_byte & 0x7f
    z = sign * (high_byte * 0x100 + low_byte)
    sm = z.sum()
 
    print("mean_z      : {}".format(float(sm) / float(num)))
    print("r5 time: {} sec".format(time.time()-st))
 
def r6():
    import numpy
    num = 3601*3601
    sm = 0
 
    st = time.time()
    with open("test.hgt", 'rb') as f:
        data = bytearray(f.read(num*2));
    data_array_byte = numpy.array(data, dtype=numpy.int16)
 
    high_byte = data_array_byte[::2]
    low_byte = data_array_byte[1::2]
    sign = 1 - (high_byte // 0x80) * 2
    high_byte = high_byte & 0x7f
    z = sign * (high_byte * 0x100 + low_byte)
    sm = z.sum()
 
    print("mean_z      : {}".format(float(sm) / float(num)))
    print("r6 time: {} sec".format(time.time()-st))
 
 
r1()
r2()
r3()
r4()
if not is_py3k: # no numpy in my py3k
    r5()
    r6()

Et voici mes résultats (sachant que le code d’origine prenait dans les 16 secondes chez moi aussi)*:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
$ python2.7 ./test2.py 
mean_z      : -4.95617519926
r1 time: 10.6385879517 sec
mean_z      : -4.95617519926
r2 time: 0.332727909088 sec
mean_z      : -4.95617519926
r3 time: 1.2303750515 sec
mean_z      : -4.95617519926
r4 time: 0.193406105042 sec
mean_z      : -6.46038493581
r5 time: 0.471071004868 sec
mean_z      : -6.46038493581
r6 time: 0.364951133728 sec
 
$ python3.2 ./test2.py 
mean_z      : -4.956175199258499
r1 time: 8.517001867294312 sec
mean_z      : -4.956175199258499
r2 time: 0.6926469802856445 sec
mean_z      : -4.956175199258499
r3 time: 1.622196912765503 sec
mean_z      : -4.956175199258499
r4 time: 0.48600101470947266 sec

Donc, si tes nombres sont des shorts codés de façon standard, la r4() de dividee est la meilleure, sinon, c’est la r6() qui gagne…

[dividee]

Apparemment il s'agit de données d'élévations de la surface terrestre, le format est expliqué ici: http://dds.cr.usgs.gov/srtm/version1.../SRTM_Topo.txt

Il est seulement dit que c'est du big-endian (il faut donc effectivement faire le byteswap sur la plupart des architectures), mais pas que c'est en complément à 2. Mais le complément à 2 étant le format le plus classique, c'est probablement celui-là qui est utilisé.

[dividee]

Hum. J'ai récupéré quelques fichiers ici: http://dds.cr.usgs.gov/srtm/version2_1/ mais les seules valeurs négatives que je trouve sont codées 0x80 0x00, qui semble indiquer une absence des données... Si ces "trous" ne sont pas comblés, il faudrait éviter de les compter pour calculer l'élévation moyenne...

Les résultats sont irréalistes avec les versions numpy (r5 et r6):

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
...
mean_z      : 1245.24099426
r4 time: 3.68899989128 sec
...
mean_z      : 4.17507479062
r6 time: 1.47000002861 sec

En visualisant avec un éditeur hexa, les valeurs semblent tourner autour de "0x05 0x00" (1280). numpy tronque sans doute la somme à 32 bits.
Si je remplace

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

    sm = z.sum()

par

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
    sm = 0
    for v in z:
        sm += int(v)

j'obtiens:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
mean_z      : 1329.0461232150253
r6 time: 17.742000102996826 sec

qui semble plus correct. Mais évidemment c'est nettement plus lent... La différence entre 1245 et 1329 est expliquée par le fait que 0x80 0x00 donne -32768 en complément à 2 mais 0 avec l'autre codage...

[MatRem]

Merci d'avoir pris le temps de me lire.

Ces données n'utilisent pas le complément à 2 pour les nombres négatifs seulement un bit de signe. C'est pourquoi je ne peux pas lire les données d'une façon aussi simple que me le propose dividee.

Si ce n'est pas le fait d'itérer qui est coûteux pourquoi cette version tourne aussi en 17 sec :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
path = "N47E006.hgt"
longitude_line_count = 3601
latitude_point_count = 3601
 
dt1_file = open( path, 'rb' )
 
mean_z = 0
 
data = dt1_file.read( 2 * longitude_line_count * latitude_point_count )
 
for longitude_id in range( longitude_line_count ) :
	for latitude_id in range( latitude_point_count ) :
		idx = longitude_id * latitude_point_count + latitude_id
		high_byte = ord( data[idx*2] )
		low_byte = ord( data[idx*2+1] )
		sign = 1 - ( ( high_byte & 0x80 ) / 0x80 ) * 2
		high_byte = high_byte & 0x7f
		z = sign * ( high_byte * 0x100 + low_byte )
 
		mean_z += z
 
print( "mean_z      : " + str( float(mean_z) / float( longitude_line_count * latitude_point_count ) ) )
 
dt1_file.close()

[dividee]

Si ce n'est pas le fait d'itérer qui est coûteux pourquoi cette version tourne aussi en 17 sec

Je me suis trompé, désolé. J'ai sous-estimé la lenteur de python pour l'arithmétique...

Pour me faire pardonner, voici comment résoudre le problème avec numpy:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

sm = z.sum(dtype=numpy.int64)

Pour une solution sans numpy, on peut aussi utiliser le complément à 2 et puis ajuster:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
...
short_array = array.array('h')
with open( path, 'rb' ) as dt1_file:
    short_array.fromfile(dt1_file, longitude_line_count * latitude_point_count)
short_array.byteswap()
 
short_array = [v if v>=0 else -(v+0x8000) for v in short_array] # ajustement
 
mean_z = sum(short_array) / float( longitude_line_count * latitude_point_count)

C'est nettement moins rapide qu'avec numpy, mais tout de même 5 fois plus rapide que le code original.