modules.f90

Go to the documentation of this file.

!=======================================================================
!> @file modules.f90
!> @brief Global functions and subroutines module
!> @authors Ary Rodriguez-Gonzalez & Alejandro Esquivel
!> @date 23/Jun/2011
!
! Copyright (c) 2011 Ary Rodriguez-Gonzalez & Alejandro Esquivel
!
! This file is part of AGA-V1.
!
! AGA-V1 is free software; you can redistribute it and/or modify
! it under the terms of the GNU General Public License as published by
! the Free Software Foundation; either version 3 of the License, or
! (at your option) any later version.
!
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
! You should have received a copy of the GNU General Public License
! along with this program.  If not, see http://www.gnu.org/licenses/.
!=======================================================================
!> @brief module mpi_module
!> @details This module contains declarations needed to implement the
!! Message Passing Interface (MPI) paralelization.
!-----------------------------------------------------------------------
module mpi_module
        implicit none
#ifdef MPIP
        include "mpif.h"

!> @brief rank of the master node (who is in charge of the I/O)
integer, parameter :: master=0

!> @brief rank of each processor (if MPI is disabled has a value of 0)
integer :: rank
!> @brief total number of processors (if MPI is disabled has a value of 1)
integer :: np
integer :: err                                          !< @brief MPI error variable
integer :: status(MPI_STATUS_SIZE)      !< @brief MPI status variable

#else

integer, parameter :: np=1, master=0, rank=0

#endif

end module mpi_module
!-----------------------------------------------------------------------
!> @brief module func
!! @details This module contains global functions and subroutines used
!! by AGA-V1
module func
contains
!-----------------------------------------------------------------------
!> @brief Updates the size of the hyper-boxes
!> @param[in] npar : number of free parameters (parameter space
!! dimension, denoted by "n" in the paper)
!> @param[in] niter : number of iteration counter
!> @param[in] nind : number of individials in the entire population
!> @param[in] deltaconst : if equal to 1 the reduction of the hyperboxes
!! is at a constant rate (see reduc), otherwise hyperboxes are reduced 
!! with the standard deviation of each parameter.
!> @param[in] reduc : reducing factor, if a constant speed is chosen the
!! hyper-box is reduced according to: 
!! @f$\Delta_{i,niter}=\Delta_{i,niter-1}*reduc^{niter}@f$
!> @param[in] arrmain(npar,nind) : Entire population of "nind" individuals, 
!! each with "npar" parameters.
!> @param[in] delta0(npar) : Original size of the hyperboxes, 
!! array of @f$\Delta_{i,1}@f$.
!> @param[out] delta(npar) : array of @f$\Delta_{i,niter}@f$ (niter 
!! is constant, thus is only a vector w/npar entries). 
!!@f$\Delta_{i,niter}@f$ is the size of the hyper-box corresponding to 
!! the "i-th" parameter
  subroutine deltabox(npar,niter,nind,deltaconst,reduc,arrmain,delta0,delta)
    implicit none
    real,dimension(npar), intent(in)  :: delta0
    real,dimension(npar), intent(out) :: delta
    real,    intent(in),  dimension(npar,nind)       :: arrmain
    integer, intent(in):: npar,niter,deltaconst,nind
    real, intent(in) :: reduc
    real :: mean
    integer :: i
    
    if (deltaconst .eq. 1 .or. niter .eq. 1) then
       do i=1,npar
          delta(i) = delta0(i)*reduc**niter
       enddo
    else
       do i=1,npar
          delta(i)=sigm(arrmain(i,:),nind)
       enddo
    endif

  end subroutine deltabox
!-----------------------------------------------------------------------
!> @brief Computes the mean value
!> @param[in] N : number of data points
!> @param[in] stat(N) : vector from which the mean is to be obtained
!> @return Arithmetic mean value of stat(N)
   real function med(stat,N)
     implicit none
     integer, intent(in)   :: N
     real, dimension(N), intent(in) :: stat
     integer :: i
     !
     med=0.
     do i=1,N
        med=med+stat(i)
     enddo
     !
     med=med/N
     !
    return
    !
  end function med
!-----------------------------------------------------------------------
!> @brief Computes the standard deviation
!> @param[in] N : number of data points
!> @param[in] stat(N) : vector from which the standard deviaiton is to 
!! be obtained
!> @return standard deviation of stat(N)
   real function sigm(stat,N)
     implicit none
     integer, intent(in)   :: N
     real, intent(in), dimension(N) :: stat
     real :: mean
     integer :: i
     !
     mean=med(stat,N)
     !
     sigm=0.
     do i=1,N
        sigm=sigm+(stat(i)-mean)**2
     enddo
     !
     sigm=sqrt(sigm/(N))
     !
     return
     !
   end function sigm
!-----------------------------------------------------------------------
!> @brief Asexual reproduction routine
!> @param[in] npar : number of free parameters
!> @param[in] gensize : size of the new generation from each parent.
!! Since the parent is kept gensize= num of sons + 1
!> @param[in] nmig : turns on/off migration. If its value is 0, 
!! migration is disabled, otherwise is enabled
!> @param[in] ngaussdev : if value equal to 1, offsprinf is obtained
!! with a Gaussian distribution, otherwise they are obtained with
!! uniform random deviates.
!> @param[in] delta(npar) : Size of the parent hyper-box.
!> @param[in] qpar(npar) : parent parameters (its location in the 
!! npar-dimensional parameter space
!> @param[in] qbounds(npar) : bounds of the initial universe
!> @param[out] qfam(npar,gensize) : A new family (progeny+parent)
  subroutine spawn(npar,gensize,nmig,ngaussdev,delta,qpar,qbounds,qfam)
    implicit none
    integer, intent(in) :: npar,gensize,nmig,ngaussdev
    real, dimension(npar), intent(in):: delta, qpar
    real, dimension(2,npar), intent(in) ::qbounds
    real, dimension(npar,gensize) :: qfam
    integer :: i,j,rand_seed
    real :: ran1,ran
    !
    do i=1,npar
       qfam(i,1)=qpar(i)
    enddo

        do j=2,gensize
                i=1
                do while (i.le.npar)
                if (ngaussdev .eq. 1) then
                    ran1=gaussdev(rand_seed)
                else    
             call random_number(ran)
             ran1=2.*ran-1.
                        endif
                        
                        qfam(i,j)=qfam(i,1)+delta(i)*ran1
          
                        if (nmig.ne.0) then
                                i=i+1
                        else if ((qfam(i,j) .ge. qbounds(1,i)).and.                 &
                                 (qfam(i,j) .le. qbounds(2,i)) ) then
                                i=i+1
                        endif
        enddo 
        enddo   

return
end subroutine spawn
!-----------------------------------------------------------------------
!> @brief Evaluates convergence criterion
!> @param[in] npar : Number of free parameters
!> @param[in] q(npar) : Individual data
!> @param[in] delta(npar) : Window size of the individual
!> @param[in] beta : Tolerance value
!> @param[out] ncrit : Return value=1 if criterion met, otherwise is not
!! modified
subroutine criterion(npar,q,delta,beta,ncrit)
  implicit none
  integer, intent(in)  :: npar
  real, intent(in), dimension(npar) :: q, delta
  real, intent(in) :: beta
  logical, dimension(npar) :: crit
  integer, intent(out) :: ncrit
  integer :: i

  crit(:)=.False.
  do i=1,npar
     if (abs(delta(i)/q(i)) .le. beta) crit(i)=.true.
  enddo

  if (all(crit(:)).eq. .true.) ncrit=1

  return
end subroutine criterion
!-----------------------------------------------------------------------
!> @brief Generates new realization from experimental data
!> @param[in] ndata : number of experimental data points
!> @param[in,out] yobs(ndata) : experimental data
!> @param[in] ysigma(ndata) : associated uncertainties of yobs
!> @todo Hay que implementar dispersion Gaussiana y checar por que no se 
!! guarda yobs0
subroutine randdata(ndata, yobs, ysigma)
  implicit none
  integer, intent(in) :: ndata
  real, intent(inout), dimension(ndata) :: yobs
  real, intent(in), dimension(ndata) :: ysigma
  integer :: i
  real :: ran1,ran
  !
  do i=1,ndata
     call random_number(ran)
     ran1 = 2.*ran-1.
     yobs(i)= yobs(i)+ysigma(i)*ran1
  end do
!
return
!
end subroutine randdata
!-----------------------------------------------------------------------
!> @brief Generates pseudo-random (Gaussian distributtion)
!> @details Modified from "Numerical Recipes in Fortran 77, second ed."
!! (Press et al., Cambridge University Press)
!> @param[in] idum : dummy variable (required)
!> @return Pseudo-random number, drawn from a normal distribution
!! (Gaussian with mean=0 and standard deviation=1)
function gaussdev(idum)
  implicit none
  real :: gaussdev
  !  
  integer :: iset,idum
  real :: fac,gset,rsq,v1,v2,ran
  save iset,gset
  data iset/0/
  !
  rsq=0.
  if (idum .lt. 0.) iset=0
  if (iset .eq. 0) then
     do while (rsq .gt. 1.0 .or. rsq .eq. 0.)
        call random_number(ran)
        v1=2.*ran-1.
        call random_number(ran)
        v2=2.*ran-1.
        rsq=v1**2+v2**2
     enddo
     fac=sqrt(-2.*log10(rsq)/rsq)
     gset=v1*fac
     gaussdev=v2*fac
     iset=1
  else
     gaussdev=gset
     iset=0.
  endif
  !
  return
end function gaussdev
!-----------------------------------------------------------------------
!> @brief Select the fittest individuals (nparents number of future
!! parents) from the entire population. Based in the smallest value of
!! the merit function (smaller number=fittest, e.g. @f$\chi^2@f$)
!> @param[in] npar : numbr of free parameters
!> @param[in] nind : number of individuals in the entire population
!> @param[in] nparents: number of parents, total of fittest individuals
!! to survive for the next generation
!> @param[in] arrmain(npar,nind) : entire population
!> @param[in,out] fitness(nind) : fitness level of each individual, 
!! given by the merit function provided b the user in "tarfunc.f90".
!> @param[out] arrpars(npar,nparents) : array of the fittest individuals,
!! which will be parents in the following generation.
subroutine pickbest(npar,nind,nparents,arrmain,fitness,arrpars)
  implicit none
  integer, intent(in) :: npar,nparents,nind
  real,    intent(in),  dimension(npar,nind)       :: arrmain
  real,    intent(inout),  dimension(nind)         :: fitness
  real                  ,  dimension(nind)         :: fit2
  real,    intent(out), dimension(npar,nparents):: arrpars
  integer, dimension(nind) :: index
  integer :: i
  
  call indexx(nind,fitness,index)
  do i=1,nparents
     arrpars(:,i)=arrmain(:,index(i))
     fit2(i)=fitness(index(i))
  end do
  
  fitness(1:nparents)=fit2(1:nparents)
  
contains
!-----------------------------------------------------------------------
!> @brief indexing routine
!> @details Modified from "Numerical Recipes in Fortran 77, second ed."
!! (Press et al., Cambridge University Press)
!> @param[in] n: number of elements in the list
!> @param[in] arr(n) : array to be indexed
!> @param[in] indx(n) : index list, such that arr(indx(j)) is in
!> ascending order for j=1, 2, 3 ...n.  arr and n are not changed
  subroutine indexx(n,arr,indx)
    implicit none
    integer, intent(in)   :: n
    real, intent(in),  dimension(n) :: arr
    integer, intent(out), dimension(n) :: indx
    integer :: i,indxt,ir,itemp,j,jstack,k,l
    real    :: a
    integer, parameter :: m=15, nstack=50
    integer, dimension(nstack) :: istack
    !
    do j=1,n
       indx(j)=j
    end do
    jstack=0
    l=1
    ir=n
    do
       if(ir-l.lt.m)then
          do j=l+1,ir
             indxt=indx(j)
             a=arr(indxt)
             do i=j-1,l,-1
                if(arr(indx(i)).le.a) exit
                indx(i+1)=indx(i)
             end do
             indx(i+1)=indxt
          end do
          if(jstack.eq.0) RETURN
          ir=istack(jstack)
          l=istack(jstack-1)
          jstack=jstack-2
       else
          k=(l+ir)/2
          itemp=indx(k)
          indx(k)=indx(l+1)
          indx(l+1)=itemp
          if (arr(indx(l)).gt.arr(indx(ir))) then
             itemp=indx(l)
             indx(l)=indx(ir)
             indx(ir)=itemp
          end if
          if (arr(indx(l+1)).gt.arr(indx(ir))) then
             itemp=indx(l+1)
             indx(l+1)=indx(ir)
             indx(ir)=itemp
          end if
          if (arr(indx(l)).gt.arr(indx(l+1))) then
             itemp=indx(l)
             indx(l)=indx(l+1)
             indx(l+1)=itemp
          end if
          i=l+1
          j=ir
          indxt=indx(l+1)
          a=arr(indxt)
          do
             do
                i=i+1
                if(arr(indx(i)).ge.a) exit
             end do
             do
                j=j-1
                if(arr(indx(j)).le.a) exit
             end do
             if(j.lt.i) exit   
             itemp=indx(i)
             indx(i)=indx(j)
             indx(j)=itemp
          end do
          indx(l+1)=indx(j)
          indx(j)=indxt
          jstack=jstack+2
          if(jstack.gt.nstack) pause 'NSTACK to small to indexx'
          if(ir-i+1.ge.j-l)then
             istack(jstack)=ir
             istack(jstack-1)=i
             ir=j-1
          else
             istack(jstack)=j-1
             istack(jstack-1)=l
             l=i
          end if
       end if
    end do
    !
  end subroutine indexx

end subroutine pickbest
!-----------------------------------------------------------------------
!> @todo Preguntarle a ary que chingaos hace esta rutina
subroutine gaussian(ngaus,i,arrmain)
implicit none
integer :: ng,ng2
real, dimension(3) :: aux
integer, intent(in) :: i,ngaus
real, intent(inout), dimension(:,:) :: arrmain
!
do ng=ngaus-1,1,-1
   do ng2=0,ng-1
      !
      if(arrmain(3*ng2+2,i) .gt. arrmain(3*ng2+5,i)) then
         aux(1:3)=arrmain(3*ng2+1:3*ng2+3,i)
         arrmain(3*ng2+1:3*ng2+3,i)=arrmain(3*ng2+4:3*ng2+6,i)
         arrmain(3*ng2+4:3*ng2+6,i)=aux(1:3)
      end if
      !
   enddo
enddo
!
end subroutine gaussian
!
end module func
!-----------------------------------------------------------------------