extfunc.f90

Go to the documentation of this file.

!=======================================================================
!> @file extfunc.f90
!> @brief User provided 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 Number of Gaussian components
!> @details Short module to implement the number of Gaussian components
!! used for fit (employed for HH34) as a variable of easy global access.
module numgaus
!> @param ngaus : number of Gaussian components to be fit
  integer :: ngaus
end module numgaus
!-----------------------------------------------------------------------
!> @brief External functions module 
!> @details Here the user needs to provide a
!> target (tarfunc) and a merit function (meritfunc). We recommend to
!! preserve the names and interface in order to avoid changing code 
!! elsewhere.
module exfunc
contains
!---------------------------------------------------------------
!> @brief Target Function
!> @details User provided target function. \n
!! In case some observational or experimental data is being fit,
!! this is the function adjusted.\n
!! For the example implemented here, finding a maxima of a given 
!! function, the evaluation of the level of fitness can be made directly
!! inside the merit function, thus the target function  returns a zero
!! (and it is actually not called). Three examples of target functions 
!! can be found commented below: several Gaussian components, one
!! or two Lorentzians, and a PNe luminosity function.
!> @param[in] x : Independent variable (abscissa)
!> @param[in] npar : Number of independent parameters
!> @param[in] q(npar) : Vector of independent parameters
!> @return Value of the target function
  real function tarfunc(x,npar,q)
    !
    use numgaus
    !
    implicit none
    integer, intent(in) :: npar
    real, dimension(npar), intent(in):: q
    real, intent(in) :: x
!    real, sum0,dintg,intg
!    integer, NumPNs
!
! Funcion Gausianas
!    select case (ngaus)
!    case (1)
!       tarfunc=q(1)*exp(-(x-q(2))**2/q(3))
!    case (2)
!       tarfunc=q(1)*exp(-(x-q(2))**2/q(3))&
!              +q(4)*exp(-(x-q(5))**2/q(6))+q(7)
!    case (3)
!       tarfunc=q(1)*exp(-(x-q(2))**2/q(3))&
!              +q(4)*exp(-(x-q(5))**2/q(6))&
!              +q(7)*exp(-(x-q(8))**2/q(9))
!    case (4)
!       tarfunc=q(1)*exp(-(x-q(2))**2/q(3))&
!             +q(4)*exp(-(x-q(5))**2/q(6))&
!              +q(7)*exp(-(x-q(8))**2/q(9))&
!              +q(10)*exp(-(x-q(11))**2/q(12))
!    case (5)
!       tarfunc=q(1)*exp(-(x-q(2))**2/q(3))&
!              +q(4)*exp(-(x-q(5))**2/q(6))&
!              +q(7)*exp(-(x-q(8))**2/q(9))&
!              +q(10)*exp(-(x-q(11))**2/q(12))&
!              +q(13)*exp(-(x-q(14))**2/q(15))
!    case (6)
!       tarfunc=q(1)*exp(-(x-q(2))**2/q(3))&
!              +q(4)*exp(-(x-q(5))**2/q(6))&
!              +q(7)*exp(-(x-q(8))**2/q(9))&
!              +q(10)*exp(-(x-q(11))**2/q(12))&
!              +q(13)*exp(-(x-q(14))**2/q(15))&
!              +q(16)*exp(-(x-q(17))**2/q(18))
!    end select

! Funcion Lorentzians
!    select case (ngaus)
!       pi=3.141592
!    case (1)
!       tarfunc=q(1)/(pi*q(2)*(1.+((x+q(3)))/q(2)))
!    case (2)
!       tarfunc=q(1)/(pi*q(2)*(1.+((x+q(3)))/q(2)))&
!              +q(4)/(pi*q(5)*(1.+((x+q(6)))/q(5)))
!            
!    end select
!
! Planetary Nebulae Luminosity Function
!    tarfunc=q(1)*exp(0.307*x)*(1.-exp(3.*(q(2)-x)))
!
!    NumPNs=100
!    summ0=0.
!    Do j=0,NumPNs-1 
!       summ0=3.*exp(3*q(1))/(exp(3.*q(1))-exp(3.*A(0,j)))+summ0
!    EndDo
!    dintg=1.114*(exp(3.*q(1)-2.693*ml)-exp(0.307*q(1)))    
!
    tarfunc=0.
    return
  end function tarfunc
!-----------------------------------------------------------------------
!> @brief Merit Function
!> @details User provided merit function. \n
!! This function computes the fitness level of each individual. \n
!! The program assumes that a low value has the highest level of fitness
!! (since its primary application is to fit functions with a merit
!! function that is an error measure, e.g. @f$\chi^2@f$). If one needs
!! to find a maximum of a function, one can use its reciprocal as merit
!! function.\n
!! The interface of the function is designed to adjust a target function
!! to a set of observational/experimental data, there some arguments
!! become dummy if the merit function does not need to be
!! called.
!> @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[in] q(npar) : Vector of parameters to be adjusted
!> @param[in] npar : Number of parameters adjusted
!> @param[in] ndata : Number of data points
!> @return Value of the merit function (lowest value=fittest individual)
    real function meritfunc(xobs,yobs,ysigma,q,npar,ndata)
    implicit none
    real, dimension(ndata), intent(in) :: xobs,yobs,ysigma
    real, dimension(npar), intent(in) :: q
    real :: chi2,yfunc,pi,a,b,sig,x0,y0,r
    integer, intent(in):: npar,ndata
    integer :: i,n
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!                 Chi^2 fitting                                    !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!    chi2=0.
!    do i=1,ndata
!       yfunc=tarfunc(xobs(i),npar,q)
!       chi2=chi2+((yobs(i)-yfunc)/(ysigma(i)))**2
!    enddo
!    meritfunc=chi2
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!                       TEST FUNCTION 1                            !
!       f(x,y)=[16x(1-x)y(1-y)sin(n*pi*x)sin(n*pi*y)]^2            !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    pi=3.141592
    n=9
    meritfunc=16.*q(1)*(1.-q(1))*q(2)*(1.-q(2))*sin(pi*n*q(1))*sin(pi*n*q(2))
    meritfunc=1.-meritfunc
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!                      TEST FUNCTION 2                             !
!                    r=sqrt((x-x0)^2+(y-y0)^2)                     !
!                  f(r)=cos^a(b*pi*r)*exp(-r^2/(2*sig^2))           !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!    a=2.
!    b=3.
!    sig=2.
!    x0=0.25
!    y0=0.25
!    pi=3.141592
!    r=sqrt((q(1)-x0)**2+(q(2)-y0)**2)
!    meritfunc=(cos(b*pi*r))**a*exp(-r**2/(2.*sig))
!    meritfunc=1.-meritfunc
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!     
!                      TEST FUNCTION 3                             !     
!       f(x,y)=(1+cos(12sqrt(x^2+y^2)))/(0.5(x^2+y^2)+2)           ! 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
!    meritfunc=(1+cos(12*sqrt(q(1)**2+q(2)**2)))/(0.5*(q(1)**2+q(2)**2)+2)
!    meritfunc=1.-meritfunc
!
!
    return
  end function meritfunc
!-----------------------------------------------------------------------
end module exfunc
!-----------------------------------------------------------------------
!=======================================================================
!> @brief Input from data and options file.
!> @details contains routines to read the input data and run-time
!! parameters.
module input
contains
!> @brief Observational/experimental data input
!> @details Reads experimental/observational data from file. Data are
!! ordered by columns (see example in ObsData.dat)
!> @param[out] xobs(ndata) : vector of abscissae
!> @param[out] yobs(ndata) : vector of ordinates
!> @param[out] ysigma(ndata) : vector of uncertainties of the ordinates
!> @param[out] ndata : number of experimental points.
  subroutine ReadData(xobs,yobs,ysigma,ndata)
    implicit none
    real, dimension(:), allocatable, intent(out) :: xobs,yobs,ysigma
    integer, intent(out) :: ndata
    !
    real, dimension(:), allocatable :: x1,y1,ys1
    integer :: i, status
    integer, parameter :: n=10000000
    !
    !   read the input file (3 columns, x,y,sigma_y, and determine the
    !   number of lines it contains
    open(1,file='ObsData.dat',status='old',iostat=status)
    allocate(x1(n),y1(n),ys1(n))
    ndata=0
    do i=1,n
       read(1,*,iostat=status) x1(i),y1(i),ys1(i)
       if (status /= 0) exit
       ndata=ndata+1
    enddo
    close(1)
    !
    !   allocate the arrays (size=number of lines in input file)
    allocate(xobs(ndata),yobs(ndata),ysigma(ndata))
    !
    !   fill the arrays with the temporary buffers and free memory
    xobs  (:) = x1 (1:ndata)
    yobs  (:) = y1 (1:ndata)
    ysigma(:) = ys1(1:ndata)
    deallocate(x1,y1,ys1)
    !
    return
  end subroutine ReadData
 !---------------------------------------------------------------------
!> @brief  Reads runtime parameters from file and initializes arrays
!> @details  Reads from file (InitCond.dat) all the runtime input 
!! parameters, initializes the values of the hyperboxes, and allocates
!! memory for all the parents and population.
!> @param[out] nsteps : Number of iterations before the hyperboxes are
!!  restored to their initial size
!> @param[out] nruns : Number of times that the hyperboxes are returned
!! to their initial size
!> @param[out] nparent : Number of parents
!> @param[out] nind : Number of individuals
!> @param[out] nerr : Number of realizations performed to obtain an
!! estimate of the errors
!> @param[out] 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[out] nmig : turns on/off migration. If its value is 0, 
!! migration is disabled, otherwise is enabled
!> @param[out] ngaussdev :  : if value equal to 1, offsprinf is obtained
!! with a Gaussian distribution, otherwise they are obtained with
!! uniform random deviates.
!> @param[out] igaus :  ????? Quiza deberia pasarse al modulo de las 
!! gaussianas
!> @param[out] npar :  number of free parameters
!> @param[out] gensize : size of the new generation from each parent.
!! Since the parent is kept gensize=num of sons + 1
!> @param[out] reduc : reducing factor, if a constant speed is chosen the
!! hyper-box is reduced according to: 
!> @param[out] beta : Tolerance value
!> @param[out] qbounds(2,npar) : limits of the original (user provided)
!! hyperbox
!> @param[out] q0(npar) : Initial guess, computed as the position at the 
!! center of the initial hyperbox
!> @param[out] delta0(npar) : Initial window half-size
!> @param[out] fit(nind) : Array to store the fitness level of the
!! population
!> @param[out] arrmain(npar,nind) : Array for the entire population
!> @param[out] arrpar(npar,nparent) : Array for the parents only
!> @param[out] arrgen(npar,gensize) : Array to store a single generation
!! of each parent
  subroutine InitCond(nsteps,nruns,nparent,nind,nerr, deltaconst, nmig, &
                                          ngaussdev, igaus, npar, gensize, reduc, beta,     &
                                          qbounds,q0,delta0,fit,arrmain,arrpar,arrgen )
    implicit none
    integer, intent(out):: npar,nsteps,nruns,nparent,nind,gensize,      
                           deltaconst,nmig,ngaussdev,nerr,igaus
    real, intent(out)   :: reduc,beta
    real, dimension(:),   allocatable, intent(out):: q0,delta0,fit
    real, dimension(:,:), allocatable, intent(out):: arrmain,arrpar,
                                                     arrgen,qbounds
    real, dimension(:,:), allocatable :: aux
    !
        integer :: i,status
    integer, parameter :: n=10000
    !
    open(1,file='InitCond.dat',status='old',iostat=status)
    allocate(aux(2,1000))
    !
    !
2   continue
    read(1,*,err=2) nsteps,nruns,nparent,nind,nerr
    !
3   continue
    read(1,*,err=3) deltaconst, nmig ,ngaussdev, igaus
    !
4   continue
    read(1,*,err=4) reduc,beta
    !
5      continue
    read(1,*,err=5)aux(1,1),aux(2,1)
    !
    npar=1
    do i=2,n
       read(1,*,iostat=status)aux(1,i),aux(2,i)
       if (status/=0)exit
       npar=npar+1
    enddo
    !
    close(1)
    ! 
    allocate(q0(npar),delta0(npar))
    allocate(qbounds(2,npar))
    qbounds(1,:)=aux(1,1:npar)
    qbounds(2,:)=aux(2,1:npar)
    q0(:)=(qbounds(1,1:npar)+qbounds(2,1:npar))/2.
    delta0(:)=(qbounds(2,1:npar)-qbounds(1,1:npar))/2.
        !
        deallocate(aux)
    !
    gensize=nind/nparent
    allocate(fit(nind), arrmain(npar,nind),arrpar(npar,nparent),         &
             arrgen(npar,gensize))
    !
    return
  end subroutine InitCond
  !--------------------------------------------------------------------
end module input
!----------------------------------------------------------------------