mean and variance of simulated maps

. open several files in batch (in order to...)
. compute mean and variance of simulated maps
. refactoring and fixes

source/sg_griders.f90

@@ -1,9 +1,12 @@
+! subgrider :: Julio Caineta, 2010
+! sg_griders :: modulo de manipulacao de grids
+
 module sg_griders
 
 use sg_gridutils
 implicit none
 
-public::abre,bender,likely,mask,pp,subgrid
+public::abre,bender,likely,mask,pp,subgrid,simmedvar
 
 contains
 
@@ -131,35 +134,35 @@ subroutine likely(res,sims,k,p,mp,id,output,header,timer)
 integer,intent(in)::id
 logical,intent(in)::mp,header
 character(len=256),intent(in)::output
-integer::j,l,a
-real::start,finish,b
+integer::j,l,a,pp,b
+real::start,finish
 call cpu_time(start)
 if (mp) then
     open(id,file=output)
     if (header) call header_skip(id)
 end if
 a=0
-p=0
+p=0.0
+pp=0
 ! open(33,file='likely.dbg',action='write') ! dd
 ! write (33,*) size(res%val), size(sims,1) ! dd
 do j=1,size(res%val)
     if (res%val(j)>=k(1).and.res%val(j)<k(2)) then
         a=a+1
         if (mp) b=0
         do l=1,size(sims,1)
-            if (sims(l,j)>=k(1).and.sims(l,j)<k(2)) then
-                p=p+1
+            if (sims(l,j)>=k(1) .and. sims(l,j)<k(2)) then
+                pp=pp+1
+!               write (33,*) j, res%val(j), sims(l,j)
                 if (mp) b=b+1
             end if
-!            write (33,*) p,b ! dd
         end do
-        if (mp) write (id,*) b/size(sims,1)
+        if (mp) write (id,*) b/real(size(sims,1))
     else
         if (mp) write (id,*) res%nd
     end if
-!    write (33,*) j,a ! dd
 end do
-p=p/(size(sims,1)*a)
+p=pp/real((size(sims,1)*a))
 if (mp) close(id)
 call cpu_time(finish)
 timer=finish-start
@@ -207,4 +210,29 @@ subroutine bender(x,y,z,nd,id,output,header,zef,hor,res,timer)
 timer=finish-start
 end subroutine bender
 
+subroutine simmedvar(do_med,do_var,sims,id,output,header,timer)
+real,dimension(:,:),intent(in)::sims
+integer,intent(in)::id
+character(len=256),intent(in)::output
+logical,intent(in)::header,do_med,do_var
+real,intent(out)::timer
+integer::i,j
+real::start,finish
+!real,dimension(:),allocatable::med,var
+call cpu_time(start)
+!if (do_med) allocate(med(size(sims,2)))
+!if (do_var) allocate(var(size(sims,2)))
+open (id,file=output)
+if (header) call header_skip(id)
+do j=1,size(sims,2)
+!    if (do_med) med(j)=media(sims(:,j))
+!    if (do_var) var(j)=variancia(sims(:,j))
+    if (do_med .and. .not.do_var) write(id,*) media(sims(:,j))
+    if (.not.do_med .and. do_var) write(id,*) variancia(sims(:,j))
+    if (do_med .and. do_var) write(id,*) media(sims(:,j)),variancia(sims(:,j))
+end do
+call cpu_time(finish)
+timer=finish-start
+end subroutine simmedvar
+
 end module

source/sg_gridutils.f90

@@ -1,3 +1,6 @@
+! subgrider :: Julio Caineta, 2010
+! sg_gridutils:: modulo de dependencias para manipulacao de grids (sg_griders)
+
 module sg_gridutils
 
 use qsort_c_module
@@ -61,7 +64,8 @@ function media(data) result(med)
 function variancia(data) result(var)
 real,intent(in),dimension(:)::data
 real::var
-var=media(data**2)-media(data)**2
+var=sum((data-media(data))**2)/size(data)
+!var=media(data**2)-media(data)**2
 end function variancia
 
 ! calcula o percentil p de uma grid

source/sg_main.f90

@@ -1,3 +1,6 @@
+! subgrider :: Julio Caineta, 2010
+! sg_main :: interface
+
 program subgrider
 
 use sg_utils
@@ -7,27 +10,30 @@ program subgrider
 implicit none
 
 ! declaracao de variaveis
-integer::i,nvar,fid,batch,nr_sims,nr_p,zef
+integer::i,nvar,fid,batch,nr_sims(2),nr_p,zef,ask_simmedvar
 integer,dimension(3)::dg,pa,pb,dg_2d
 character(len=256)::ficheiro,output,pocos,base,mp_bend
 character(len=1)::load,mp_p
-logical::mp,header
+logical::mp,header,do_med,do_var,did_load,did_loadsims
 real,allocatable::sims(:,:)
 integer,allocatable::px(:),py(:)
 real::P_x(2),P1,timer,nd,p,q,op,p_sim(2)
 type(grid)::res,hor
 
+did_load=.FALSE.
+did_loadsims=.FALSE.
 !ciclo principal
 do
 batch=-1
 print *,""
-print *,"----|| s u b g r i d e r  v0.2.1 <<jc 2010>> ||----"
+print *,"----|| s u b g r i d e r  v0.2.3 <<jc 2010>> ||----"
 print *,""
 print *,"escolher uma opcao"
 print *,""
 print *,"0 ----- carregar dados"
 print *,"  0.1 - abrir ficheiro"
 print *,"  0.2 - abrir ficheiro anterior"
+print *,"  0.3 - abrir varios ficheiros"
 print *,"1 ----- poco vertical "
 print *,"  1.1 - a partir das coordenadas (x,y)"
 print *,"  1.2 - a partir de um ficheiro com coordenadas (x,y)"
@@ -39,6 +45,8 @@ program subgrider
 print *,"  4.3 - calcular verosimilhanca"
 print *,"5 - planificacao"
 print *,"  5.1 - planificar uma grid a partir de um mapa de planificacao"
+print *,"6 - analise de simulacoes"
+print *,"  6.1 - media e variancia das simulacoes"
 print *,"9 - sair"
 read *,op
 if (op==0.1) then
@@ -60,6 +68,8 @@ program subgrider
     print *,"a carregar o ficheiro ",trim(ficheiro),"..."
     call abre(ficheiro,header,dg,nd,res,fid,timer)
     print *,"ficheiro carregado em ",tempo(timer)
+    did_load=.TRUE.
+    call wait()
 elseif (op==0.2) then
     print *,"abrir ficheiro anterior"
     print *,""
@@ -71,7 +81,25 @@ program subgrider
     print *,"a carregar o ficheiro ",trim(ficheiro),"..."
     call abre(ficheiro,header,dg,nd,res,fid,timer)
     print *,"ficheiro carregado em ",tempo(timer)
-elseif (op==1.1) then
+    did_load=.TRUE.
+    call wait()
+elseif (op==0.3) then
+    print *,"abrir varios ficheiros numerados de A a B"
+    print *,"inserir intervalo (A,B)"
+    read *,nr_sims
+    print *,"base"
+    read *,base
+    call header_ask("r",header,nvar)
+    if (.not. did_load) then
+        print *,"dimensoes das grids: x y z"
+        read (*,*) dg
+    end if
+    print *,"a carregar ficheiros..."
+    call abre_sims(nr_sims,sims,dg,base,'.OUT',header,timer)
+    print *,nr_sims(2)-nr_sims(1)+1,"ficheiros carregados em ",tempo(timer)
+    did_loadsims=.TRUE.
+    call wait()
+elseif (op==1.1 .and. did_load) then
     print *,"extrair poco vertical a partir das coordenadas (x,y)"
     print *,""
     print *,"coordenadas do poco (x,y)"
@@ -81,7 +109,8 @@ program subgrider
     print *,"a furar o poco..."
     call pp(pa(1),pa(2),res,fid,output,header,batch,timer)
     print *,"operacao concluida em ",tempo(timer)
-elseif (op==1.2) then
+    call wait()
+elseif (op==1.2 .and. did_load) then
     print *,"extrair pocos verticais a partir de um ficheiro com coordenadas (x,y)"
     print *,""
     print *,"ficheiro com coordenadas (caminho/nome)"
@@ -97,9 +126,10 @@ program subgrider
     do i=1,size(px)
         call pp(px(i),py(i),res,fid,output,header,batch,timer)
     end do
-    print *,"operacao concluida em ",tempo(timer) !tempo devolvido é do último poço
+    print *,"operacao concluida em ",tempo(timer) !tempo devolvido do ultimo poco
     batch=-1
-elseif (op==2) then
+    call wait()
+elseif (op==2 .and. did_load) then
     print *,"obter uma grid a partir um volume (cubo, paralelipipedo) da grid inicial"
     print *,""
     print *,"coordenadas x (min,max)"
@@ -113,27 +143,18 @@ program subgrider
     print *,"a criar subgrid..."
     call subgrid(pa,pb,res,fid,output,header,timer)
     print *,"operacao concluida em ",tempo(timer)
-elseif (op==3) then
+    call wait()
+elseif (op==3 .and. did_load) then
     print *,"criar uma mascara (data=1, no data=0)"
     print *,""
     call novo(header,nvar,fid,output,batch,timer)
     print *,"a criar mascara..."
     call mask(1,0,header,res,fid,output,timer)
     print *,"operacao concluida em ",tempo(timer)
-elseif (op==4.3) then
+    call wait()
+elseif (op==4.3 .and. did_load .and. did_loadsims) then
     print *,"update bayesiano - calcular verosimilhanca"
     print *,""
-    print *,"numero de simulacoes realizadas (ficheiros *_simX.OUT)"
-    read *,nr_sims
-! <alt>
-    print *,"base"
-    read *,base
-! <\alt>
-!    base=ficheiro(1:len_trim(ficheiro)-4)//'_sim'
-    call header_ask("r",header,nvar)
-    print *,"a carregar grids simuladas..."
-    call abre_sims(nr_sims,sims,dg,base,'.OUT',header,timer)
-    print *,nr_sims,"grids carregadas em ",tempo(timer)
     do
         print *,"produzir grid com mapa de probabilidades? (s/n)"
         read *,mp_p
@@ -160,32 +181,44 @@ program subgrider
     print *,"a probabilidade e' ",trim(rtochar(P1))
     print *,"operacao concluida em ",tempo(timer)
     ! colocado aqui temporariamente, incorporar como batch do percentil
-    print *,""
-    print *,"percentis das simulacoes"
-    call header_ask("w",header,nvar)
-    call novo(header,nvar,fid,output,batch,timer)
-    print *,"inserir intervalo de percentis (pi, pf)"
-    read *,p_sim
-    print *,"a calcular percentis..."
-    call percentil_sims(sims,p_sim,P_x,fid,header,output,timer) !verificar P_x(1)
-    print *,"operacao concluida em ",tempo(timer)
-elseif (op==4.1) then
+    do
+	    print *,"produzir ficheiro com os percentis das simulacoes? (s/n)"
+	    read *,mp_p
+	    if (mp_p=="s" .or. mp_p=="S") then
+		    call header_ask("w",header,nvar)
+		    call novo(header,nvar,fid,output,batch,timer)
+		    print *,"inserir intervalo de percentis (pi, pf)"
+		    read *,p_sim
+		    print *,"a calcular percentis..."
+		    !call percentil_sims(sims,p_sim,P_x,fid,header,output,timer) !verificar P_x(1)
+		    call percentil_sims(sims,p_sim,P_x,fid,header,output,timer) !verificar P_x(1)
+		    print *,"operacao concluida em ",tempo(timer)
+	    elseif (mp_p=="n" .or. mp_p=="N") then
+	        exit
+	    else
+	        print *,"erro - input invalido."
+	    end if
+    end do
+    call wait()
+elseif (op==4.1 .and. did_load) then
     print *,"update bayesiano - calcular percentil"
     print *,""
     print *,"inserir percentil a calcular"
     read *,p
     call percentil(res,p,q,timer)
     print '("o percentil ",f4.2," e ",a)',p,trim(rtochar(q))
     print *,"operacao concluida em ",tempo(timer)
-elseif (op==4.2) then
+    call wait()
+elseif (op==4.2 .and. did_load) then
     print *,"update bayesiano - calcular quantil"
     print *,""
     print *,"inserir valor"
     read *,q
     call percentil_updt(res,q,p,timer)
     print *,"o valor ",trim(rtochar(q))," corresponde ao percentil ",trim(rtochar(p))
     print *,"operacao concluida em ",tempo(timer)
-elseif (op==5.1) then
+    call wait()
+elseif (op==5.1 .and. did_load) then
     print *,"planificar uma grid a partir de um mapa de planificacao"
     print *,"atencao: a grid inicial (alocada) vai ser alterada"
     print *,""
@@ -206,6 +239,36 @@ program subgrider
     print *,"a planificar grid..."
     call bender(dg(1),dg(2),dg(3),-999.25,fid,output,header,zef,hor,res,timer)
     print *,"operacao concluida em ",tempo(timer)
+    call wait()
+elseif (op==6.1 .and. did_loadsims) then
+    print *,"analise de simulacoes"
+    print *,"media e variancia das simulacoes"
+    print *,""
+    print *,"calcular media (1), variancia (2) ou ambas (3)?"
+    read *,ask_simmedvar
+    do
+        if (ask_simmedvar==1) then
+            do_med=.TRUE.
+            do_var=.FALSE.
+            exit
+        elseif (ask_simmedvar==2) then
+            do_med=.FALSE.
+            do_var=.TRUE.
+            exit
+        elseif (ask_simmedvar==3) then
+            do_med=.TRUE.
+            do_var=.TRUE.
+            exit
+        else
+            print *,"erro - input invalido"
+        end if
+    end do
+    call header_ask("w",header,nvar)
+    call novo(header,nvar,fid,output,batch,timer)
+    print *,"a calcular..."
+    call simmedvar(do_med,do_var,sims,fid,output,header,timer)
+    print *,"operacao concluida em ",tempo(timer)
+    call wait()
 elseif (op==9) then
     print *,"programa terminado"
     stop
@@ -214,36 +277,12 @@ program subgrider
 end if
 end do
 
-! funcoes e subrotinas
-
 contains
 
-! devolve o tempo, em string, gasto pela ultima subrotina executada
-function tempo(timer) result(time)
-real,intent(in)::timer
-real::te
-!character(10),allocatable,dimension(:)::time
-character(32)::time
-character(len=6)::format,uni
-if (timer<1) then
-    te=timer*1000
-    format="(f7.3)"
-    uni=" ms"
-elseif (timer>=1 .and. timer<60) then
-    te=timer
-    format="(f6.3)"
-    uni=" s"
-elseif (timer>=60 .and. timer<3600) then
-    te=timer/60
-    format="(f6.3)"
-    uni=" min"
-else
-    te=timer/3600
-    format="(f7.3)"
-    uni=" h"
-end if
-write (time,format) te
-time=time(1:len_trim(time))//trim(uni)
-end function tempo
+subroutine wait()
+print *,""
+print *,"primir [ENTER] para voltar ao menu principal"
+read (*,'()')
+end subroutine wait
 
 end program