aga-v1.f90

Go to the documentation of this file.

!=======================================================================
!> @file aga-v1.f90
!> @brief The implementation of AGA
!> @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 containing the implementation of the algorithm
!> @details It has to be called from a main program that serves as an 
!! interface if used with a experimental/observational data.
!> @todo Debug OpenMP implementation
module agav1
contains
! ----------------------------------------------------------------------
!> @brief THE algorithm
!> @details This is the implementation of the Asexual Genetic Algorithm
!! (AGA), it requires that the experimental/observational data have been
!! previously loaded (this is done in the main program, and was left out
!! of the main algorithm to improve flexibility), and if MPI is enabled
!! it also needs to be prevously loaded. The algorithm loads the rest of 
!! the runtime parameters and return the fittest individual and an
!! estimate of the errors.
!> @param[in] xobs(ndata) : Vector of data abscissae
!> @param[in] yobs(ndata) : Vector of data ordinates
!> @param[in] ysigma(ndata) : Vector of uncertainties in the ordinate
!> @param[out] q(npar) : Vector of parameters to be adjusted, at the end
!! returns the fittest individul 
!> @param[out] delta(npar) : Window (hyperbox) half-size, at the end
!! the variable returns the uncertainty of q
!> @param[in] ndata : Number of data points
  subroutine Algorithm(xobs,yobs,ysigma,q,delta,npar,ndata)
    use input                   ! located in extfunc.f90
    use func                    ! located in extfunc.f90
    use exfunc                  ! located in extfunc.f90
    use mpi_module              ! located in modules.f90
    use numgaus                 ! located in modules.f90
    implicit none
    integer :: npar, ndata, nsteps, nruns, nparent, nind, nerr,         
               deltaconst, nmig, ngaussdev, igaus, gensize
    real :: beta,reduc,fito,deltafit
    real :: mean, sigma
    !
    real, dimension(:), allocatable :: q0,delta0, q, delta, fit
    real, dimension(:), allocatable :: xobs,yobs,ysigma
    !
    real, dimension(:,:),allocatable :: arrmain,arrpar,arrgen,arrerr,   
                                        qbounds
    !
    integer :: i,j,k,niter,idum,nit,ij,ncrit
    integer :: stk, nproc, omp_get_thread_num, KMP_GET_STACKSIZE_S,nprocs
    integer :: omp_get_num_threads
    !
    integer :: rand_size
    integer, allocatable, dimension(:) :: rand_seed
    character (len=10) :: system_time
    real :: rtime
    !
#ifdef MPIP
    real, dimension(:,:),allocatable :: arrerrtot
#endif
    !
    !     OpenMP initialization
    !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    !$    WRITE (6,*) 'Compiled with OpenMP'
    !$OMP PARALLEL PRIVATE(NPROC)
    !$           stk = KMP_GET_STACKSIZE_S()
    !$           nproc = OMP_GET_THREAD_NUM()
    !$           WRITE(6,*) 'Processor (thread) #:',nproc,' ready'
    !$OMP END PARALLEL
    !$           WRITE(6,*) 'Stack size per thread : ', stk/1048576,' Mb'
    !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 
    !
    !   initializes seed for random number generator
    
    call random_seed(size=rand_size)
    allocate(rand_seed(1:rand_size))
    call date_and_time(time=system_time)
    read(system_time,*) rtime
    rand_seed=int(rtime*1000.)
#ifdef MPIP
    rand_seed=rand_seed*rank
#endif
    call random_seed(put=rand_seed)    
    deallocate(rand_seed)
    !----------------------------------------------
    !
#ifdef MPIP
    do i=0,np-1
                if (rank.eq.i) then
#endif
                        call InitCond(nsteps,nruns,nparent,nind,nerr, deltaconst,   &
                                              nmig, ngaussdev, igaus, npar, gensize, reduc, &
                                          beta, qbounds,q0,delta0,fit,arrmain,arrpar,arrgen )                     
        
                        allocate(arrerr(npar,nerr/np),q(npar),delta(npar))
                        q(:)=q0
                        delta(:)=delta0
#ifdef MPIP
                        allocate(arrerrtot(npar,nerr))
       endif
       call mpi_barrier (mpi_comm_world,err)
    enddo
#endif
    !
    call spawn(npar,nparent,nmig,ngaussdev,delta,q,qbounds,arrpar)
    !   
    ! Main loop
    !  
    do j=1,nerr/np
          Print *,'Nerr=',j,np
       !
        if (nerr .ne. 1) call randdata(ndata,yobs,ysigma)
       !
       nit=1
       !
       do while(nit.le.nruns)
          !
          deltafit=1e20
          fit(1)=1e20
          niter = 1
          ncrit=0
          !
          do while(niter .le. nsteps)
 
            call deltabox(npar,niter,nind,deltaconst,reduc,arrmain,     &
                          delta0,delta)         
             !
             !$omp parallel shared(arrgen,arrmain,arrpar,fit) &
             !$omp           private(nproc,nprocs)
             !$ nproc = OMP_GET_THREAD_NUM()
             !$ nprocs= omp_get_num_threads()
             !$omp do SCHEDULE(static,nparent/nprocs)
             !
             do i=1,nparent
                call spawn(npar,gensize,nmig,ngaussdev,delta,           &
                           arrpar(:,i), qbounds, arrgen )
                !
                arrmain(:,((i-1)*gensize+1):i*gensize)=arrgen(:,:)
             enddo
             !$omp end do
             !
             !$omp do SCHEDULE(static,nind/nprocs)
             do i=1,nind
                !!   this works only for sums of gaussians
                If (igaus .eq. 1) Then
                   !
                   call gaussian(ngaus,i,arrmain)
                   !
                EndIf
                q(:)=arrmain(:,i)
                fit(i)=meritfunc(xobs,yobs,ysigma,arrmain(:,i),         &
                                 npar,ndata)
             enddo
             !$omp end do
             !$omp end parallel
             !
             call pickbest(npar,nind,nparent,arrmain,fit,arrpar)

             niter = niter+1
             deltafit=abs(1-fito/fit(1))
             fito=fit(1)
! CRITERIO NUEVO
!             call criterion(nparn,q,delta,beta,ncrit)
!             print*,nit,niter,ncrit,j
!             if (ncrit .eq. 1) exit
!
!             if (deltafit .lt. beta) exit
              print*,niter,nit
              if (fit(1) .lt. beta) exit
!
          enddo
          nit=nit+1
          !
       enddo
      arrerr(:,j)=arrpar(:,1) !c(:)
    end do
    !
#ifdef MPIP
    call mpi_gather (arrerr,    npar*nerr/np, mpi_real8 ,               &
                     arrerrtot, npar*nerr/np, mpi_real8 ,               &
                      master, mpi_comm_world, err )
    if (rank.eq.master) then
#endif
    !
    open(2,file='coeffit.out')
    !
    do i=1,npar
#ifndef MPIP
       q(i)=med(arrerr(i,:),nerr)
       delta(i)=sigm(arrerr(i,:),nerr)
#else
       q(i)=med(arrerrtot(i,:),nerr)
       delta(i)=sigm(arrerrtot(i,:),nerr)
#endif
    enddo
    if (rank.eq.master) then
       do i=1,npar
          write(2,*)q(i),delta(i)
       enddo
          write(2,*)fit(1)
    endif
    close(2)
#ifdef MPIP
 endif
 dellocate(arrerrtot)
 call mpi_barrier (mpi_comm_world, err)
#endif
 deallocate(fit,arrerr)
 
end subroutine Algorithm
end module agav1
!-----------------------------------------------------------------------