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é.