fortran calcul Biais entre 2 fichiers [Débutant(e)]

Invité · 16/04/2013, 13h37

Je viens de débuter un stage dans le domaine de la modélisation ...et je n'ai pas encore toutes les bases en programmation.Je me permets de poster mon message sur ce forum en espérant que vous pourriez m'aider.

Je dois comparer 2 fichiers en calculant le biais:

B = 1/n.Σi (Mi -Oi) Avec Mi la valeur donnée par le model a la date i
Et Oi la valeur donné par le fichier observation à la meme date i.

Mes 2 fichiers:
-Un fichier.txt (Observation.txt) dans lequel se trouve les données observations:
avec la date, l'heure et la concentration d'un polluant à cet heure sur une période comme ce qui suit:
2012081801 47
2012081802 45
2012081804 43
2012081805 40
-Un fichier .txt (Model.txt) dans lequel se trouve les données fournies par un model:
2012081801 67.254
2012081802 62.861
2012081803 62.427
2012081804 63.149
2012081805 63.076

je ne sais pas trop comment m'y prendre, En vous remerciant par avance.

**__dardanos__** · 16/04/2013, 17h10

Bienvenue greendiamond

D'après les extraits des fichiers montrés, je comprends que :

A chaque observation dans le premier fichier correspond une valeur du modèle dans le second,
Les deux fichiers ont le même nombre de lignes,
Les données sont classées dans le même ordre dans les chaque fichier.

Si les trois hypothèses sont vérifiées, le problème est simple.

Connais-tu n (le nombre de valeurs d'observation) ?
Si oui tu peux directement déclarer déclarer des tableaux O et M de dimension n.
Sinon il faut parcourir une première fois l'un des fichier, pour déterminer le déterminer et ensuite allouer dynamiquement O et M.

Il faut après initialiser les tableaux à partir des données lues. Le traitement est le même pour chaque fichier : lire la seconde valeur de la ligne courante et l'affecter à un élément du tableau (les dates ne sont pas utilisées d'après ce que j'ai compris).

Code :

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

1
2
3
4
5
open(unit=11,nom_du_fichier)
do i = 1, n
read(11,*) date, tableau(i)
enddo
close(11)

Le FORTRAN 90 permet la manipulation globale des tableaux. Ainsi d'après ta formule, le biais s'obtient :

Code :

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

B = sum(M-O)/n

On préfère utiliser un forme quadratique pour estimer l'erreur.

Invité · 17/04/2013, 11h37

Salut Dardanos,

Tout d'abord je te remercie pour le temps que tu as consacré pour me répondre...

Oui les données sont bien dans le même ordre. Cependant, le souci est qu'il me manque des données dans le fichier observation,
par exemple dans l'exemple du début de fichier (observation) que j'ai copié, il me manque la concentration correspondant à 3h pour le 18 aout 2012 (2012081803).
Mais comment lui dire de passer à la ligne suivante quand n'y a pas la même date (exple: 2012081803) dans les 2 fichiers?

Donc pour n je devrais certainement prendre le nombre de ligne qu'il y a dans mon premier fichier (le fichier observation.txt).

Oui je commence avec les biais qui me semble le plus simple, puis ensuite j'adapterais pour calculer le RMSE (root mean square error) ainsi que la corrélation de Pearson.

En te remerciant

Invité · 17/04/2013, 11h54

je copie ici, ce que j'ai commencé.....

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
 
 
program Calc_Bias
 
!!   Calcul du Biais
	implicit none
 
!! variables
 
      character (len=10) :: datechO 					! 2012081801 date + heure fichier observation
      character (len=10) :: datechM					! 2012081801 date + heure fichier modele
      real  :: concO							!Concentration d'ozone Observation (ug/m3)
      real  :: ConcM							!Concentration d'ozone Modele (ug/m3)	
      real  :: B                                                            ! biais
      real  :: difference	             		                    				
 
!!    local variables
      integer :: nlines                                     
 
 
 
!!      ouverture du fichier observation
 
		open(unit=120,file="/home/observations.dat",status='old')
		read (120,*) datechO, concO
 
! compter le nombre de lignes dans le fichier observation
		nlines=0
 
		do  ! boucle pour compter le nombre de ligne du fichier
			read(120,*) 
			 nlines = nlines + 1
			print*, 'nombre de lignes', nlines
		enddo
 
!!      ouverture du fichier model
 
		open(unit=121,file="/home/model.dat",status='old')
		read (121,*) datechM, concM
 
 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!               Calcul du Biais                           !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
		if (datechO = datechM)
			do
			difference= concM-concO
 
			else !! passer à la ligne suivante???
 
		endif
 
 
		do	
			B=sum(concM-concO)/nlines
		enddo
 
		write 
 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
			close(120)
			close(121)
end program Calc_Bias

**__dardanos__** · 17/04/2013, 13h32

Salut.
Voici un modèle de base :

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
program main
    implicit none
    integer :: ios_o, ios_m, nb_conc
    integer :: date_o, date_m
    double precision :: conc_o, conc_m, somme
 
    open(unit=11, file='Observation.txt', status='old')
    read(11,*) date_o, conc_o
 
    open(unit=12, file='Model.txt', status='old')
    read(12,*) date_m, conc_m
 
    nb_conc = 0 ; somme = 0
    do
        if (date_m < date_o) then
            read(11,*,iostat=ios_m) date_m, conc_m
            if (ios_m/=0) exit
        endif
        if (date_m == date_o) then
            print '(i10,2e14.6)', date_o, conc_o, conc_m
            somme = somme + conc_o - conc_m
            nb_conc = nb_conc + 1
        else
            print '(a,i10)', 'Observations manquantes : ', date_o
        endif
        read(12,*,iostat=ios_o) date_o, conc_o
        if (ios_o/=0) exit
    enddo
 
    close(11) ; close(12)
 
    print '(a,e14.6)', 'Biais : ', somme / dble(nb_conc)
 
end program main

Stocker les concentrations dans des tableaux est inutile.

Invité · 17/04/2013, 15h36

Je te remercie infiniment dardanos, je le test des demain sur une machine linux et je te tiens au courant.
Bonne soirée à toi.

**__dardanos__** · 17/04/2013, 15h36

Il y avait des erreurs dans ma réponse. Une correction :

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
program main
    implicit none
    integer :: ios_o, ios_m, nb_conc = 0
    integer :: date_o, date_m
    double precision :: conc_o, conc_m, somme = 0
 
    open(unit=11, file='Observation.txt', status='old')
    read(11,*) date_o, conc_o
 
    open(unit=12, file='Model.txt', status='old')
    read(12,*) date_m, conc_m
 
    lect_obs : do
 
        ! Boucler tant que la date lue dans le fichier modele est
        ! inferieure a la date lue dans le fichier observation
        lect_model : do
            if (date_m >= date_o) exit lect_model
            read(12,*,iostat=ios_m) date_m, conc_m
            if (ios_m > 0) exit lect_obs    ! Erreur
            if (ios_m < 0) then             ! Fin de fichier
                print '(a/a/a)', repeat('-',70), &
                    'Toutes les valeurs du fichier modele ont ete traitees', &
                    'Valeurs fichier observation non traitees :'
                do
                    read(11,*,iostat=ios_o) date_o, conc_o
                    if (ios_o/=0) exit lect_obs
                    print '(i10,e14.6)', date_o, conc_o
                enddo
            endif
 
        enddo lect_model
 
        ! Date observation trouvée : ajout contribution
        if (date_m == date_o) then
            print '(i10,2e14.6)', date_o, conc_o, conc_m
            somme = somme + conc_o - conc_m
            nb_conc = nb_conc + 1
        else
            ! Cas : date_m > date_o
            print '(a,i10)', 'Le modele ne permet pas de traiter l''observation a : ', date_o
        endif
 
        ! Lecture d'une nouvelle ligne dans le fichier observation
        read(11,*,iostat=ios_o) date_o, conc_o
        if (ios_o > 0) exit lect_obs    ! Erreur
        if (ios_o < 0) then             ! Fin de fichier
            print '(a/a/a)', repeat('-',70), &
                'Toutes les valeurs du fichier observation ont ete traitees', &
                'Valeurs fichier modele non traitees :'
            do
                read(12,*,iostat=ios_m) date_m, conc_m
                if (ios_m/=0) exit lect_obs
                print '(i10,e14.6)', date_m, conc_m
            enddo
        endif
 
    enddo lect_obs
    close(11) ; close(12)
 
    print '(a/a,e14.6)', repeat('-',70), 'Biais : ', somme / dble(nb_conc)
 
end program main

Invité · 18/04/2013, 14h38

Salut Dardanos,

Hélas, toujours dans mes soucis de Biais...
Les deux programmes précédent, m'affichent la même chose quand je les compile et les exécute:

2012081801 0.470000E+02 0.672540E+02
Biais : -0.202540E+02

En fait j'ai l'impression qu'il me calcule seulement la difference Model-observation correspondant à la première date du fichier (2012081801).

**__dardanos__** · 18/04/2013, 18h45

Salut.
Vérifie bien les fichiers. Avec les 2 fichiers de données que tu as posté, et la dernière version du programme, j'obtiens :

Code :

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

1
2
3
4
5
6
7
8
9
2012081801  0.470000E+02  0.672540E+02
2012081802  0.450000E+02  0.628610E+02
2012081804  0.430000E+02  0.631490E+02
2012081805  0.400000E+02  0.630760E+02
----------------------------------------------------------------------
Toutes les valeurs du fichier observation ont ete traitees
Valeurs fichier modele non traitees :
----------------------------------------------------------------------
Biais :  -0.203350E+02

J'ai vérifié, c'est bien la valeur attendue.

Invité · 24/04/2013, 09h33

Salut Dardanos,
excuse moi pour le retard; J'ai eu quelques soucis de santé; J'ai retesté et de mon coté il m'affiche bien ce que j'avais copié précedement;
Je n'ai pas la même chose que toi...

Invité · 24/04/2013, 09h47

Finalament,

le programme marche à merveille;
Je te remercie infiniment Dardanos; C'est juste qu'un de mes fichiers avait été modifié, il ne restait plus que la première ligne....
Merci mille fois.

Bon courage à toi et bonne continuation

Invité · 24/04/2013, 15h54

salut,
J'ai pu faire un programme calculer le RMSE et maintenant je teste pour calculer la correlation:
j'ai pu faire des programmes séparés (un qui calcul la moyenne, un qui
calcul les eccarts type et un pour la covariance mais en fixant la
valeur des moyennes) mais lorsque je mets tout dans un seul
programme, il ne marche pas. En fait, j'ai l'impression qu'il manque
quelque chose pour lui dire de prendre les valeurs des moyennes qu'il
calcul au début du programme...

Voici ce que j'ai pour l'instant:

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
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
 
program Correlation
    implicit none
    integer :: ios_o, ios_m, nb_conc = 0, i
    integer :: date_o, date_m, ok
    double precision :: conc_o, conc_m, somme_o = 0, somme_m = 0, somme_a= 0, somme_b= 0,  cov, Moy_m , Moy_o
    integer :: ecc_type_model, ecc_type_obs, covariance
 
 
    open(unit=11, file='Observations.dat', status='old')      ! fichier observation
    read(11,*) date_o, conc_o
 
 
    open(unit=12, file='Model.dat', status='old') !fichier model
    read(12,*) date_m, conc_m
 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!                                                                             !!!     
!!!     Moy_o & Moy_m           moyenne Observation & Moyenne Model             !!!
!!!                                                                             !!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
    lect_obs : do
 
        ! Boucler tant que la date lue dans le fichier modele est
        ! inferieure a la date lue dans le fichier observation
        lect_model : do
            if (date_m >= date_o) exit lect_model
            read(12,*,iostat=ios_m) date_m, conc_m
            if (ios_m > 0) exit lect_obs    ! Erreur
            if (ios_m < 0) then             ! Fin de fichier
                print '(a/a/a)', repeat('-',70), &
                    'Toutes les valeurs du fichier modele ont ete traitees', &
                    'Valeurs fichier observation non traitees :'
                do
                    read(11,*,iostat=ios_o) date_o, conc_o
                    if (ios_o/=0) exit lect_obs
                    print '(i10,e14.6)', date_o, conc_o
                enddo
            endif
 
        enddo lect_model
 
        ! Date observation trouvée : ajout contribution
        if (date_m == date_o) then
            print '(i10,2e14.6)', date_o, conc_o, conc_m
	    somme_o = somme_o + conc_o
	    somme_m = somme_m +	conc_m
	    nb_conc = nb_conc + 1
        else
            ! Cas : date_m > date_o
            print '(a,i10)', 'Le modele ne permet pas de traiter l''observation a : ', date_o
        endif
 
        ! Lecture d'une nouvelle ligne dans le fichier observation
        read(11,*,iostat=ios_o) date_o, conc_o
        if (ios_o > 0) exit lect_obs    ! Erreur
        if (ios_o < 0) then             ! Fin de fichier
            print '(a/a/a)', repeat('-',70), &
                'Toutes les valeurs du fichier observation ont ete traitees', &
                'Valeurs fichier modele non traitees :'
            do
                read(12,*,iostat=ios_m) date_m, conc_m
                if (ios_m/=0) exit lect_obs
                print '(i10,e14.6)', date_m, conc_m
            enddo
        endif
 
    enddo lect_obs
 
	Moy_o = somme_o / dble(nb_conc)
	Moy_m = somme_m / dble(nb_conc)
 
    	print '(a/a,e14.6)', repeat('-',70), 'Moyenne Observation: ', Moy_o
	print '(a/a,e14.6)', repeat('-',70), 'Moyenne Model: ', Moy_m
 
 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!                                                                             !!!     
!!!                             //calcul de l'eccart type du model              !!!
!!!                                                                             !!!
!!! ecc_type_model     = racine (1/(nb_conc) * Somme (conc_m - Moy_m)^2         !!!
!!!                                                                             !!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
 
Moy_o = somme_o / dble(nb_conc)
Moy_m = somme_m / dble(nb_conc)
 
       lect_obs_A : do
 
 
        ! Boucler tant que la date lue dans le fichier modele est
 
        ! inferieure a la date lue dans le fichier observation
 
        lect_model_A : do
 
            if (date_m >= date_o) exit lect_model_A
 
            read(12,*,iostat=ios_m) date_m, conc_m
 
            if (ios_m > 0) exit lect_obs_A    ! Erreur
 
            if (ios_m < 0) then             ! Fin de fichier
 
 !             print '(a/a/a)', repeat('-',70), &
 
 !                  'Toutes les valeurs du fichier modele ont ete traitees', &
 
 !                  'Valeurs fichier observation non traitees :'
 
                do
 
                    read(11,*,iostat=ios_o) date_o, conc_o
 
                    if (ios_o/=0) exit lect_obs_A
 
                    print '(i10,e14.6)', date_o, conc_o
 
                enddo
 
            endif
 
 
 
        enddo lect_model_A
 
 
 
        ! Date observation trouvée : ajout contribution
 
        if (date_m == date_o) then
 
            print '(i10,2e14.6)', date_o, conc_m, Moy_m
 
            somme_a = somme_a + (conc_m - Moy_m )**2
 
            nb_conc = nb_conc + 1
 !       else
 
 !          ! Cas : date_m > date_o
 
 !          print '(a,i10)', 'Le modele ne permet pas de traiter l''observation a : ', date_o
 
        endif
 
 
 
        ! Lecture d'une nouvelle ligne dans le fichier observation
 
        read(11,*,iostat=ios_o) date_o, conc_o
 
        if (ios_o > 0) exit lect_obs_A    ! Erreur
 
 !      if (ios_o < 0) then             ! Fin de fichier
 
 !          print '(a/a/a)', repeat('-',70), &
 
 !              'Toutes les valeurs du fichier observation ont ete traitees', &
 
 !              'Valeurs fichier modele non traitees :'
 
            do
 
                read(12,*,iostat=ios_m) date_m, conc_m
 
                if (ios_m/=0) exit lect_obs_A
 
                print '(i10,e14.6)', date_m, conc_m
 
            enddo
 
 !      endif
 
 
 
    enddo lect_obs_A
 
ecc_type_model = sqrt(somme_a/dble(nb_conc))
 
print '(a/a,e14.6)', repeat('-',70), 'eccart type model : ', sqrt(somme_a/dble(nb_conc))
 
 
 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!                                                                             !!!     
!!!                             //calcul de l'eccart type Observation           !!!
!!!                                                                             !!!
!!! ecc_type_obs =      racine (1/(nb_conc) * Somme (conc_o - Moy_o)^2          !!!
!!!                                                                             !!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
 
       lect_obs_B : do
 
 
 
        ! Boucler tant que la date lue dans le fichier modele est
 
        ! inferieure a la date lue dans le fichier observation
 
        lect_model_B : do
 
            if (date_m >= date_o) exit lect_model_B
 
            read(12,*,iostat=ios_m) date_m, conc_m
 
            if (ios_m > 0) exit lect_obs_B    ! Erreur
 
            if (ios_m < 0) then             ! Fin de fichier
 
 !             print '(a/a/a)', repeat('-',70), &
 
 !                  'Toutes les valeurs du fichier modele ont ete traitees', &
 
 !                  'Valeurs fichier observation non traitees :'
 
                do
 
                    read(11,*,iostat=ios_o) date_o, conc_o
 
                    if (ios_o/=0) exit lect_obs_B
 
                    print '(i10,e14.6)', date_o, conc_o
 
                enddo
 
            endif
 
 
 
        enddo lect_model_B
 
 
 
        ! Date observation trouvée : ajout contribution
 
        if (date_m == date_o) then
 
            print '(i10,2e14.6)', date_o, conc_m, Moy_o
 
            somme_b = somme_b + (conc_o - Moy_o )**2
 
            nb_conc = nb_conc + 1
 !       else
 
 !          ! Cas : date_m > date_o
 
 !        print '(a,i10)', 'Le modele ne permet pas de traiter l''observation a : ', date_o
 
        endif
 
 
 
        ! Lecture d'une nouvelle ligne dans le fichier observation
 
        read(11,*,iostat=ios_o) date_o, conc_o
 
        if (ios_o > 0) exit lect_obs_B    ! Erreur
 
 !      if (ios_o < 0) then             ! Fin de fichier
 
 !          print '(a/a/a)', repeat('-',70), &
 
 !              'Toutes les valeurs du fichier observation ont ete traitees', &
 
 !              'Valeurs fichier modele non traitees :'
 
            do
 
                read(12,*,iostat=ios_m) date_m, conc_m
 
                if (ios_m/=0) exit lect_obs_B
 
                print '(i10,e14.6)', date_m, conc_m
 
            enddo
 
 !      endif
 
 
 
    enddo lect_obs_B
 
ecc_type_obs = sqrt(somme_b/dble(nb_conc))
 
print '(a/a,e14.6)', repeat('-',70), 'eccart type Observation : ', sqrt(somme_b/dble(nb_conc))
 
 
 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!                                                                             !!!     
!!!                             //calcul de la covariance cov                   !!!
!!!                                                                             !!!
!!!covariance = 1/nb_conc* somme[(conc_m-Moy_m)*(conc_o-Moy_o)]                 !!!
!!!                                                                             !!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
   lect_obs_C : do
 
 
        ! Boucler tant que la date lue dans le fichier modele est
 
        ! inferieure a la date lue dans le fichier observation
 
        lect_model_C : do
 
            if (date_m >= date_o) exit lect_model_C
 
            read(12,*,iostat=ios_m) date_m, conc_m
 
            if (ios_m > 0) exit lect_obs_C    ! Erreur
 
  !         if (ios_m < 0) then             ! Fin de fichier
 
  !            print '(a/a/a)', repeat('-',70), &
 
  !                 'Toutes les valeurs du fichier modele ont ete traitees', &
 
  !                 'Valeurs fichier observation non traitees :'
 
                do
 
                    read(11,*,iostat=ios_o) date_o, conc_o
 
                    if (ios_o/=0) exit lect_obs_C
 
                    print '(i10,e14.6)', date_o, conc_o
 
                enddo
 
!           endif
 
 
 
        enddo lect_model_C
 
 
 
        ! Date observation trouvée : ajout contribution
 
        if (date_m == date_o) then
 
            print '(i10,2e14.6)', date_o, conc_m, Moy_o
 
            cov = cov + (conc_o - Moy_o )* (conc_m - Moy_m )
 
            nb_conc = nb_conc + 1
  !      else
 
  !         ! Cas : date_m > date_o
 
  !         print '(a,i10)', 'Le modele ne permet pas de traiter l''observation a : ', date_o
 
        endif
 
 
 
        ! Lecture d'une nouvelle ligne dans le fichier observation
 
        read(11,*,iostat=ios_o) date_o, conc_o
 
        if (ios_o > 0) exit lect_obs_C    ! Erreur
 
  !     if (ios_o < 0) then             ! Fin de fichier
 
  !         print '(a/a/a)', repeat('-',70), &
 
  !             'Toutes les valeurs du fichier observation ont ete traitees', &
 
  !             'Valeurs fichier modele non traitees :'
 
            do
 
                read(12,*,iostat=ios_m) date_m, conc_m
 
                if (ios_m/=0) exit lect_obs_C
 
                print '(i10,e14.6)', date_m, conc_m
 
            enddo
 
 !      endif
 
 
 
    enddo lect_obs_C
 
covariance = cov/dble(nb_conc)
 
print '(a/a,e14.6)', repeat('-',70), 'Covariance : ', cov/dble(nb_conc)
 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
  		close(11) ; close(12)
 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!                                                                             !!!     
!!!                             //calcul de la correlation                      !!!
!!!                                                                             !!!
!!!correlation =covariance/(ecc_type_obs*ecc_type_model)                        !!!
!!!                                                                             !!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
print '(a/a,e14.6)', repeat('-',70), 'Correlation :',  covariance/(ecc_type_obs*ecc_type_model)
 
end program Correlation

Merci par avance

**__dardanos__** · 24/04/2013, 17h28

Salut,
Ça n'est pas la bonne méthode. Il vaut mieux lire les données une seule fois, et les stocker dans un tableau. Le traitement est ensuite simple avec les fonctions intrinsèques.

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
program main
    implicit none
    type enreg
        integer :: date
        double precision :: conc_o
        double precision :: conc_m
    end type enreg
 
    integer :: i, nb_conc
    character(len=20) :: file_o = 'Observation.txt'
    character(len=20) :: file_m = 'Model.txt'
    type(enreg), dimension(:), allocatable :: tab
    double precision :: moy_o, moy_m, ecc_type_m, ecc_type_o, cov
 
    call extraire_donnees(file_o,file_m)
 
    ! Affichage
    do i = 1, size(tab)
        print *, tab(i)
    enddo
    print '(a)', repeat('-',70)
 
    ! Biais
    print '(a,e14.6)', 'Biais : ', sum( tab%conc_o - tab%conc_m ) / nb_conc
 
    ! Moyennes
    moy_o = sum(tab%conc_o) / nb_conc
    moy_m = sum(tab%conc_m) / nb_conc
    print '(a,e14.6)', 'Moyenne observations : ', moy_o
    print '(a,e14.6)', 'Moyenne modele : '      , moy_m
 
    ! Ecarts types
    ecc_type_m = sqrt( sum( (tab%conc_m - moy_m)**2 ) / nb_conc )
    ecc_type_o = sqrt( sum( (tab%conc_o - moy_o)**2 ) / nb_conc )
    print '(a,e14.6)', 'Ecart type du modele : ',   ecc_type_m
    print '(a,e14.6)', 'Ecart type observation : ', ecc_type_o
 
    ! Covariance
    cov = sum( (tab%conc_m - moy_m)*(tab%conc_o - moy_o) ) / nb_conc
    print '(a,e14.6)', 'Covariance : ', cov
 
    ! Correlation
    print '(a,e14.6)', 'Correlation : ', cov / (ecc_type_o * ecc_type_m)
 
    deallocate(tab)
 
contains
 
    subroutine extraire_donnees (file_o,file_m)
        character(len=*) :: file_o, file_m
        integer :: ios_o, ios_m
        integer :: date_o, date_m
        double precision :: conc_o, conc_m
        logical :: premier = .true.
 
        do
            open(unit=11, file= file_o, status='old')
            read(11,*) date_o, conc_o
 
            open(unit=12, file= file_m, status='old')
            read(12,*) date_m, conc_m
 
            nb_conc = 0
            lect_obs : do
 
                ! Boucler tant que la date lue dans le fichier modele est
                ! inferieure a la date lue dans le fichier observation
                lect_model : do
                    if (date_m >= date_o) exit lect_model
                    read(12,*,iostat=ios_m) date_m, conc_m
                    if (ios_m /= 0) exit lect_obs
                enddo lect_model
 
                ! Date observation trouvée : ajout contribution
                if (date_m == date_o) then
                    nb_conc = nb_conc + 1
                    if (.not. premier) then
                        tab(nb_conc) = enreg( date_o, conc_o, conc_m)
                    endif
                endif
 
                ! Lecture d'une nouvelle ligne dans le fichier observation
                read(11,*,iostat=ios_o) date_o, conc_o
                if (ios_o /= 0) exit lect_obs
 
            enddo lect_obs
            close(11) ; close(12)
 
            if (premier) then
                premier = .false.
                allocate(tab(nb_conc))
            else
                return
            endif
 
        enddo
 
    end subroutine extraire_donnees
 
end program main

Invité · 25/04/2013, 10h21

Salut Dardanos,

merci infiniment, tout marche très bien, c'est parfait.

J'en demande peut être un peu trop mais je tente

Maintenant que j'ai ce programme qui me calcul les différents indicateurs pour une station donnée ( a un endroit donnée), je dois en créer un qui me fait la même chose mais pour toutes mes stations. (J'en ai 1330 au total).

Je m'explique:
j'ai plusieurs stations que l'on distingue par des codes( de 7 characteres): AT60190, ES1675A, IT1539A, ...etc

J'ai mes fichiers observations dont le nom est sous la forme:
AD0942A_hor_EUR_O3.dat
AD0944A_hor_EUR_O3.dat
AT0ENK1_hor_EUR_O3.dat
AT0ILL1_hor_EUR_O3.dat
AT0PIL1_hor_EUR_O3.dat
AT0SON1_hor_EUR_O3.dat
AT0VOR1_hor_EUR_O3.dat
AT0ZIL1_hor_EUR_O3.dat
AT0ZOE2_hor_EUR_O3.dat
(les 7 premiers caractères de chaque fichier correspondent au code de la station)

J'ai mes fichiers modèle dont le nom est sous la fomre suivante:
AT0ENK12012_CHIMERE_REF.dat
AT0ILL12012_CHIMERE_REF.dat
AT0PIL12012_CHIMERE_REF.dat
AT0SON12012_CHIMERE_REF.dat
AT0VOR12012_CHIMERE_REF.dat
AT0ZOE22012_CHIMERE_REF.dat
AT100022012_CHIMERE_REF.dat
AT2F1012012_CHIMERE_REF.dat
(les 7 premiers caractères de chaque fichier correspondent au code de la station)

Ce que je dois réussir à faire c'est une boucle sur les stations: c'est à dire à chaque station "model" il me faut faire la comparaison (calcul du biais, rmse, correlation) avec la station correspondant aux données d'observations
Cependant à certaines stations, les fichiers observations n'existent pas.

Le programme devrait m'ecrire dans un nouveau fichier sous forme de tableau:
le code de la station Biais RMSE ecc_type_model ecc_type_obs correlation:
par exemple pour une premiere ligne:

Code :

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

AT80503  0.352190E+02  0.436179E+02  0.174788E+02  0.229360E+02  0.211310E+00

Mon souci est donc comment faire les boucles sur les stations afin d'arriver à ce que je souhaite obtenir?
En te remerciant

**__dardanos__** · 25/04/2013, 15h03

Salut,
J'ai mis ce qui avait été fait dans le module comparaison_fichiers (à compiler d'abord) :

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
module comparaison_fichiers
    implicit none
    type enreg
        integer :: date
        double precision :: conc_o
        double precision :: conc_m
    end type enreg
 
    integer, parameter :: nb_res = 4
    integer :: nb_conc
    type(enreg), dimension(:), allocatable :: tab
 
contains
 
    function comparaison (file_o, file_m)
        character(len=*), intent(in) :: file_o, file_m
        double precision :: moy_o, moy_m, ecc_type_m, ecc_type_o, cov
        double precision, dimension(nb_res) :: comparaison
 
        call extraire_donnees(file_o, file_m)
 
        ! Biais
        comparaison(1) = sum( tab%conc_o - tab%conc_m ) / nb_conc
 
        ! Moyennes
        moy_o = sum(tab%conc_o) / nb_conc
        moy_m = sum(tab%conc_m) / nb_conc
 
        ! Ecarts types
        ecc_type_m = sqrt( sum( (tab%conc_m - moy_m)**2 ) / nb_conc )
        ecc_type_o = sqrt( sum( (tab%conc_o - moy_o)**2 ) / nb_conc )
        comparaison(2) = ecc_type_m ; comparaison(3) = ecc_type_o
 
        ! Covariance
        cov = sum( (tab%conc_m - moy_m)*(tab%conc_o - moy_o) ) / nb_conc
 
        ! Correlation
        comparaison(4) = cov / (ecc_type_o * ecc_type_m)
 
    end function comparaison
 
    subroutine extraire_donnees (file_o, file_m)
        character(len=*) :: file_o, file_m
        integer :: ios_o, ios_m
        integer :: date_o, date_m
        double precision :: conc_o, conc_m
        logical :: premier = .true.
 
        do
            open(unit=11, file= file_o, status='old')
            read(11,*) date_o, conc_o
 
            open(unit=12, file= file_m, status='old')
            read(12,*) date_m, conc_m
 
            nb_conc = 0
            lect_obs : do
 
                ! Boucler tant que la date lue dans le fichier modele est
                ! inferieure a la date lue dans le fichier observation
                lect_model : do
                    if (date_m >= date_o) exit lect_model
                    read(12,*,iostat=ios_m) date_m, conc_m
                    if (ios_m /= 0) exit lect_obs
                enddo lect_model
 
                ! Date observation trouvée : ajout contribution
                if (date_m == date_o) then
                    nb_conc = nb_conc + 1
                    if (.not. premier) then
                        tab(nb_conc) = enreg( date_o, conc_o, conc_m)
                    endif
                endif
 
                ! Lecture d'une nouvelle ligne dans le fichier observation
                read(11,*,iostat=ios_o) date_o, conc_o
                if (ios_o /= 0) exit lect_obs
 
            enddo lect_obs
            close(11) ; close(12)
 
            if (premier) then
                premier = .false.
                allocate(tab(nb_conc))
            else
                return
            endif
 
        enddo
 
    end subroutine extraire_donnees
 
end module comparaison_fichiers

Le programme principal doit faire ce que tu as demandé, à condition de placer les fichiers dans le répertoire courant :

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
program main
    ! ---------------------------------------------------------------
    ! Les fichiers modele et d'observation doivent etre places dans
    ! le repertoire courant.
    ! ---------------------------------------------------------------
    use comparaison_fichiers
    implicit none
    character(len=20), parameter :: suffixe_o = '_hor_EUR_O3.dat'
 
    integer :: ios, idx
    logical :: file_exists
    character(len=20) :: rad
    character(len=50) :: fichier_m, fichier_o
 
    open( unit=21, file='liste_fichiers_modele.txt',       status='old')
    open( unit=22, file='liste_fichiers_observations.txt', status='old')
    open( unit=23, file='comparaisons.txt',                status='unknown')
    open( unit=24, file='fichiers_modele_non_traites.txt', status='unknown')
 
    write(23,'(10a15)') 'Station', 'Biais', 'ecc_type_model', 'ecc_type_obs', 'correlation'
 
    do
        ! Lecture du nom du fichier modele a traiter
        read(21,*,iostat=ios) fichier_m
        if (ios/=0) exit
 
        ! Nom du fichier d'observation
        idx = index(fichier_m,'_')
        rad = fichier_m(:idx-1)
        fichier_o = trim(rad) // trim(suffixe_o)
 
        ! Test de la presence d'un fichier d'observation
        inquire(file= fichier_o, exist=file_exists)
        if (file_exists) then
            write(23,'(a15,10e15.6)') rad, comparaison( fichier_o, fichier_m)
        else
            write(24,'(a)') trim(fichier_m)
        endif
    enddo
 
    close(21) ; close(22) ; close(23) ; close (24)
 
end program main

Invité · 25/04/2013, 15h37

Lorsque je compile comparaison_fichiers, voila ce qu'il m'indique:

Code :

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

1
2
3
4
 
/usr/lib/gcc/x86_64-redhat-linux/4.1.2/libgfortranbegin.a(fmain.o): In function `main':
(.text+0xa): undefined reference to `MAIN__'
collect2: ld a retourné 1 code d'état d'exécution

Invité · 25/04/2013, 15h48

Que veux tu dire par repertoire courant? Est-ce que cela voudrait dire que je dois avoir les programme comparaison_fichier et main dans le même repertoire que mes fichiers de données et de model?

Je commence avec trois staions, j'ai créer un dossier dans lequel j'y ai mis les programmes "comparaison_fichiers.f90", "main.f90, ainsi que mes fichiers model (ES1038A2012_CHIMERE_REF.dat, FR240132012_CHIMERE_REF.dat, IT1247A2012_CHIMERE_REF.dat) et mes fichiers observations (ES1038A_hor_EUR_O3.dat, FR24013_hor_EUR_O3.dat, IT1247A_hor_EUR_O3.dat).
Puis j'ai rajouté deux fichiers:

-

liste_station_model.dat dans lequel on trouve:

Code :

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

1
2
3
ES1038A2012_CHIMERE_REF.dat
FR240132012_CHIMERE_REF.dat
IT1247A2012_CHIMERE_REF.dat

-

liste_station_observation.dat dans lequel on trouve:

Code :

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

1
2
3
ES1038A_hor_EUR_O3.dat
FR24013_hor_EUR_O3.dat
IT1247A_hor_EUR_O3.dat

**__dardanos__** · 25/04/2013, 18h02

Pour compiler à la main il faut faire :

Code :

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

1
2
3
gfortran -c comparaison_fichiers.f90
gfortran -c main.f90
gfortran -o main main.o comparaison_fichiers.o

Que veux tu dire par repertoire courant?

L'exécutable main, les fichiers liste_fichiers_modele.txt et liste_fichiers_observations.txt (et pas liste_station_model.dat et liste_station_observation.dat), et l'ensemble des fichiers d'observation et du modèle doivent être dans la même répertoire.

Invité · 25/04/2013, 18h34

En tout cas je te remercie Dardanos, je teste ca dès demain matin sur linux et je te tiens au courant.
J'avais fait gfortran comparaison_fichiers.f90 comme c'est ce que j'ai l'habitude de faire pour compiler; donc il me manquait surement -c.
Je te souhaite une bonne soirée

Invité · 26/04/2013, 09h49

Salut,

Ce que tu m'a dit de faire marche bien, cependant il ne traite aucun de mes fichiers.
J'ai bien un fichier comparaison.dat mais dans lequel il est écrit seulement:

Code :

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

        Station          Biais ecc_type_model   ecc_type_obs    correlation

quand a l'autre fichier; fichiers_modele_non_traites.dat il m'indique les noms des 3 fichiers correspondant aux stations que je voulais étudier (pourtant les fichiers observations existent bien pour ces stations):

Code :

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

1
2
3
ES1038A2012_CHIMERE_REF.dat
FR240132012_CHIMERE_REF.dat
IT1247A2012_CHIMERE_REF.dat

J'ai mis ici ce que j'avais dans mon repertoire, peut être cela pourrait aider

Pièce jointe 117047

fortran calcul Biais entre 2 fichiers [Débutant(e)]

Fortran

Discussions similaires

Partager

Partager