test_algen – TOOLIB

Uses

Help

Graph Key

Nodes of different colours represent the following:

Solid arrows point from a submodule to the (sub)module which it is descended from. Dashed arrows point from a module or program unit to modules which it uses.
Where possible, edges connecting nodes are given different colours to make them easier to distinguish in large graphs.

PIKAIA oop example of use

Note

The example chosen is the generation of a surface with given statistical moments. The implementation is part of the analytical developments presented in:
A. Francisco and N. Brunetière, “A hybrid method for fast and efficient rough surface generation”, Proc IMechE Part J:J Engineering Tribology, 2016, Vol. 230(7) 747–768, DOI: 10.1177/1350650115612116

The surface heights are generated thanks the tangent function. This function is indeed very near the common material curves. Thus, providing the right lower and upper bounds makes it possible to match the desired statistical moments.

Note

The main difficulty to cope with is the high non linearity and the convergence only possible near the solution. Classical optimization methods fail in reaching the true solution: the lower and upper limits x(1:2.
The PIKAIA program, which implements a genetic algorithm, is utilized to found a “good” solution, then a very classical optimization routine is used to find the nearly exact solution.

Warning

PIKAIA cost function must be chosen so as to be maximized.
For clarity reasons, the optimization function must be minimized because it is the deviation to the solution. @note Warning

Note

First each thread runs the PIKAIA program to find a rather good starting point x. The parameters (as defined by CTRL) are the same. The only difference is the starting population. Thus, several concurrent populations evolve during CTRL(02) generations.
Then, the optimization subroutine starts from the best population.

Note

make all
./prg
The surface is in ascii mode in the “/out” directory

Variables

Type	Attributes		Name		Initial
real(kind=R8),	parameter	::	SKU	=	+20.0_R8	Imposed surface height kurtosis
real(kind=R8),	parameter	::	SSK	=	-01.0_R8	Imposed surface height skewness
real(kind=R8)		::	f			best cost
real(kind=R8)		::	f_xx			cost function or optimization function at `xx`
integer(kind=I4)		::	i			loop indices
integer(kind=I4)		::	iu			i/o unit
integer(kind=I4)		::	j			loop indices
integer(kind=I4)		::	k			loop indices
type(moment_stat)		::	mom			a statistical moment variable
integer(kind=I4)		::	nb_th			nb threads
integer(kind=I4),	parameter	::	nn	=	512	Surface side dimension
type(pikaia_class)		::	p			PIKAIA class instanciation
integer(kind=I4)		::	status			PIKAIA status
real(kind=R8)		::	tab_sur(1:nn*nn)			output surface heights
real(kind=R8)		::	tend			time variables
real(kind=R8)		::	tstart			time variables
real(kind=R8),	dimension(1:2)	::	xl			lower and upper bonds of xx
real(kind=R8),	dimension(1:2)	::	xu			lower and upper bonds of xx
real(kind=R8),	dimension(1:2)	::	xx			chromosom for PIKAIA

Derived Types

type :: moment_stat

Statistical moment type

Components

Type	Visibility		Name
real(kind=R8),	public	::	ku	kurtosis
real(kind=R8),	public	::	mu	mean
real(kind=R8),	public	::	si	standard deviation
real(kind=R8),	public	::	sk	skewness
real(kind=R8),	public	::	va	variance

Functions

function abs_diff_sk_ku(chrom)

Arguments

Type	Intent	Optional	Attributes		Name
real(kind=R8),	intent(in),		dimension(1:2)	::	chrom

Return Value real(kind=r8)

function add_tang(n, deb, fin, alp, bet, mu, si)

Arguments

Type	Intent		Name
integer(kind=i4),	intent(in)	::	n
integer(kind=i4),	intent(in)	::	deb
integer(kind=i4),	intent(in)	::	fin
real(kind=R8),	intent(in)	::	alp
real(kind=R8),	intent(in)	::	bet
real(kind=R8),	intent(in)	::	mu
real(kind=R8),	intent(in)	::	si

Return Value real(kind=r8)

function cost_func(chrom)

Arguments

Type	Intent	Optional	Attributes		Name
real(kind=R8),	intent(in),		dimension(1:2)	::	chrom

Return Value real(kind=r8)

function tang(xi, n, alp, bet, mu, si)

Arguments

Type	Intent		Name
real(kind=R8),	intent(in)	::	xi
integer(kind=i4),	intent(in)	::	n
real(kind=R8),	intent(in)	::	alp
real(kind=R8),	intent(in)	::	bet
real(kind=R8),	intent(in)	::	mu
real(kind=R8),	intent(in)	::	si

Return Value real(kind=r8)

Subroutines

subroutine calc_moments_1D(tab, mx, nb_mom, lg)

Arguments

Type	Intent	Attributes		Name
real(kind=R8),	intent(in),	dimension(1:lg)	::	tab
type(moment_stat),	intent(out)		::	mx
integer(kind=I4),	intent(in)		::	nb_mom
integer(kind=I4),	intent(in)		::	lg

subroutine cost(me, x, f)

Arguments

Type	Intent	Attributes		Name
class(pikaia_class),	intent(inout)		::	me
real(kind=R8),	intent(in),	dimension(:)	::	x
real(kind=R8),	intent(out)		::	f

subroutine melange(tab, lg)

Arguments

Type	Intent	Optional	Attributes		Name
real(kind=r8),	intent(inout),		dimension(1:lg)	::	tab
integer(kind=i4),	intent(in)			::	lg

subroutine newton_raphson_downhill(x, fvec, eps, ndir, rel)

Arguments

Type	Intent	Attributes		Name
real(kind=R8),	intent(inout),	dimension(1:2)	::	x
real(kind=R8),	intent(out)		::	fvec
real(kind=R8),	intent(in)		::	eps
integer(kind=I4),	intent(in)		::	ndir
real(kind=R8),	intent(in)		::	rel

subroutine profil_theo_trie_1D(tab, lg, x, mx)

Arguments

Type	Intent	Attributes		Name
real(kind=R8),	intent(out),	dimension(1:lg)	::	tab
integer(kind=I4),	intent(in)		::	lg
real(kind=R8),	intent(in),	dimension(1:2)	::	x
type(moment_stat),	intent(out)		::	mx

subroutine sk_ku(xx, sk, ku)

Arguments

Type	Intent	Attributes		Name
real(kind=R8),	intent(in),	dimension(1:2)	::	xx
real(kind=R8),	intent(out)		::	sk
real(kind=R8),	intent(out)		::	ku

Source Code

program test_algen
use omp_lib
use sort_arrays,   only : sort_array2
use data_arch,     only : I4, R8, PI_R8
use miscellaneous, only : get_unit
use pikaia_oop,    only : pikaia_class

implicit none

integer(kind=I4), parameter :: nn  = 512        !! *Surface side dimension*
real(kind=R8),    parameter :: SSK = -01.0_R8   !! *Imposed surface height skewness*
real(kind=R8),    parameter :: SKU = +20.0_R8   !! *Imposed surface height kurtosis*

integer(kind=I4) :: status                      !! **PIKAIA** *status*
integer(kind=I4) :: i, j, k                     !! *loop indices*
integer(kind=I4) :: iu                          !! *i/o unit*
integer(kind=I4) :: nb_th                       !! *nb threads*

real(kind=R8), dimension(1:2) :: xx             !! *chromosom for* **PIKAIA**
real(kind=R8), dimension(1:2) :: xl, xu         !! *lower and upper bonds of xx*

real(kind=R8)    :: f_xx                        !! *cost function or optimization function at* ```xx```
real(kind=R8)    :: tab_sur(1:nn*nn)            !! *output surface heights*

type moment_stat
!! <span style="color:green">Statistical moment type</span>
   real(kind=R8) :: mu !! *mean*
   real(kind=R8) :: va !! *variance*
   real(kind=R8) :: si !! *standard deviation*
   real(kind=R8) :: sk !! *skewness*
   real(kind=R8) :: ku !! *kurtosis*
endtype moment_stat

type(moment_stat) :: mom   !! *a statistical moment variable*

type(pikaia_class)   :: p              !! **PIKAIA** *class instanciation*
real(kind=R8)        :: f              !! *best cost*
real(kind=R8)        :: tstart,tend    !! *time variables*


   !$ real(kind=R8)    :: ostart,oend    !! *openmp time variables*
   !$ integer(kind=I4) :: tid            !! *active thread*

   nb_th = 1

   ! get the number of available threads
   !$OMP PARALLEL PRIVATE(nb_th, TID)
   !$
   !$ tid = omp_get_thread_num()
   !$
   !$ if (tid == 0) then
   !$    nb_th = omp_get_num_threads()
   !$    write(*,'(A)') '--------------'
   !$    write(*,'(A,1X,I5)') 'number of OMP threads: ', nb_th
   !$    write(*,'(A)') '--------------'
   !$ endif
   !$OMP END PARALLEL

   xx(1:2) = 0.0_R8
   xl(1:2) = 0.0_R8
   xu(1:2) = 1.0_R8

   !initialize the class:
   call p%init(           n = 2,                   &  ! IN           ; the parameter space dimension, i.e., the number of adjustable parameters (size of the x vector).
                         xl = xl,                  &  ! IN, DIM(n)   ;  vector of lower bounds for x
                         xu = xu,                  &  ! IN, DIM(n)   ;  vector of upper bounds for x
                          f = cost,                &  !              ; user-supplied scalar function of n variables, which must have the pikaia_func procedure interface.
                     status = status,              &  ! OUT          ; status output flag (0 if there were no errors)
                    !iter_f = report_iteration,    &  !     OPT      ; user-supplied subroutine that will report the best solution for each generation. It must have the iter_func procedure interface.
                         np = 100,                 &  ! IN, OPT      ; number of individuals in a population (default is 100)
                       ngen = 1000,                &  ! IN, OPT      ; maximum number of iterations
                         nd = 9,                   &  ! IN           ; number of significant digits (i.e., number of genes) retained in chromosomal encoding
                     pcross = 0.85_R8,             &  ! IN, OPT      ; crossover probability; must be <= 1.0 (default is 0.85). If crossover takes place, either one or two splicing points are used, with equal probabilities
                     pmutmn = 0.0005_R8,           &  ! IN, OPT      ; minimum mutation rate; must be >= 0.0 (default is 0.0005)
                     pmutmx = 0.25_R8,             &  ! IN, OPT      ; maximum mutation rate; must be <= 1.0 (default is 0.25)
                       pmut = 0.005_R8,            &  ! IN, OPT      ; initial mutation rate; should be small (default is 0.005) (Note: the mutation rate is the probability that any one gene locus will mutate in any one generation.)
                       imut = 2,                   &  ! IN, OPT      ; mutation mode; 1/2/3/4/5 (default is 2).
                                                      !              1=one-point mutation, fixed rate.
                                                      !              2=one-point, adjustable rate based on fitness.
                                                      !              3=one-point, adjustable rate based on distance.
                                                      !              4=one-point+creep, fixed rate.
                                                      !              5=one-point+creep, adjustable rate based on fitness.
                                                      !              6=one-point+creep, adjustable rate based on distance.
                       fdif = 1._R8,               &  ! IN, OPT      ; relative fitness differential; range from 0 (none) to 1 (maximum). (default is 1.0)
                       irep = 3,                   &  ! IN, OPT      ; reproduction plan; 1/2/3=Full generational replacement/Steady-state-replace-random/Steady- state-replace-worst (default is 3)
                     ielite = 0,                   &  ! IN, OPT      ; elitism flag; 0/1=off/on (default is 0) (Applies only to reproduction plans 1 and 2)
                       ivrb = 0,                   &  ! IN, OPT      ; printed output 0/1/2=None/Minimal/Verbose
            convergence_tol = 1.0e-6_R8,           &  ! IN, OPT      ; convergence tolerance; must be > 0.0 (default is 0.0001)
         convergence_window = 200,                 &  ! IN, OPT      ; convergence window; must be >= 0 This is the number of consecutive solutions within the tolerance for convergence to be declared (default is 20)
         initial_guess_frac = 0.1_R8,              &  ! IN, OPT      ; raction of the initial population to set equal to the initial guess. Range from 0 (none) to 1.0 (all). (default is 0.1 or 10%).
                      iseed = 999)                    ! IN, OPT      ; random seed value; must be > 0 (default is 999)

   !Now call pikaia:
   call cpu_time(tstart)
   !$ ostart = omp_get_wtime()

   call p%solve(      x = xx(1:2),     &  ! INOUT, DIM(*) ;
                      f = f,           &  !   OUT         ;
                 status = status) !,   &  !   OUT         ;
!~                     omp = .false. )    ! IN, OPTIONAL

   !$ oend = omp_get_wtime()
   call cpu_time(tend)

   write(*, *) 'Total time spent: ', tend - tstart
   write(*, *) 'Real time spent:  ', oend - ostart

   write(*,*) 'Absolute diff after GA: ', abs_diff_sk_ku(chrom = xx(1:2))

   ! optimization function ran with this "good" solution
   call newton_raphson_downhill( x     = xx(1:2),           & ! starting/ending point
                                 fvec  = f_xx,              & ! best deviation to the wanted solution
                                 eps   = 1.e-6_R8,          & ! dx
                                 ndir  = 32*nb_th,          & ! number of directions to explore
                                 rel   = 0.9_R8)              ! relaxation parameter

   write(*,*) 'Absolute diff after refinement: ', abs_diff_sk_ku(chrom = xx(1:2))


   write(*,*) 'A (nn x nn) SSK, SKU random surface will be now generated'

   ! generation of nn*nn heights with the tangent function with xx(1:2) as bound parameters.
   ! mu = 0.
   ! va = si = 1.
   ! sk = SSK
   ! ku = SKU
   call profil_theo_trie_1D(tab = tab_sur(1:nn*nn),   &  ! resulting vector of heights
                             lg = nn*nn,              &  ! vector length
                              x = xx(1:2),            &  ! best variable couple
                             mx = mom)                   ! resulting moment

   ! height shuffle
   call melange(tab = tab_sur(1:nn*nn),   &  !
                 lg = nn*nn)                 !

   ! random surface output in ascii format
   call get_unit(iunit = iu)
   open(unit = iu, file = 'out/surf.dat')
      k = 1
      do i = 1, nn
      do j = 1, nn
         write(iu, *) i, j, tab_sur(k)
         k = k +1
      enddo
      enddo
   close(unit=iu)

   ! statistical moments check
   write(*,*) 'statistical moments:'
   write(*,*)  mom%mu, mom%si, mom%sk, mom%ku

   ! precision
   write(*,'(a, e8.2)') 'maximum difference (%) ', max( 100*abs((mom%sk -SSK)/SSK), & !
                                                        100*abs((mom%ku -SKU)/SKU) )

stop

contains

   subroutine cost(me, x, f)
   implicit none
   class(pikaia_class), intent(inout)               :: me
   real(kind=R8)      , intent(in   ), dimension(:) :: x
   real(kind=R8)      , intent(  out)               :: f

      f = 1./(1. + abs_diff_sk_ku(chrom = x(1:2)))

   return
   endsubroutine cost


   subroutine newton_raphson_downhill(x, fvec, eps, ndir, rel)
   implicit none
   integer(kind=I4), intent(in   )                 :: ndir
   real(kind=R8),    intent(in   )                 :: eps
   real(kind=R8),    intent(in   )                 :: rel
   real(kind=R8),    intent(inout), dimension(1:2) :: x
   real(kind=R8),    intent(out  )                 :: fvec

      real(kind=R8), dimension(1:2)    :: x0
      real(kind=R8), dimension(1:ndir) :: ct, st, ff, x1, x2

      real(kind=R8)    :: f0, fd1, fd2, ti, dfdx1, dfdx2, dx1, dx2, dt
      integer(kind=I4) :: i

      do i = 1, ndir
         ti = PI_R8*( -1. +(i -1)*2./ndir )
         ct(i) = cos(ti)
         st(i) = sin(ti)
      enddo
      x0 = x

      do
         f0  = abs_diff_sk_ku(chrom = (/x0(1)     , x0(2)     /))
         fd1 = abs_diff_sk_ku(chrom = (/x0(1) +eps, x0(2)     /))
         fd2 = abs_diff_sk_ku(chrom = (/x0(1)     , x0(2) +eps/))

         dfdx1 = (fd1 -f0)/eps
         dfdx2 = (fd2 -f0)/eps

         !$OMP PARALLEL DEFAULT(SHARED) NUM_THREADS(nb_th)
         !$OMP DO PRIVATE(i, dt, dx1, dx2)
            do i = 1, ndir
               dt  = -f0/( eps*(dfdx1*ct(i) +dfdx2*st(i)) )
               dx1 = eps*ct(i)*dt
               dx2 = eps*st(i)*dt

               x1(i) = x0(1) +rel*dx1
               x2(i) = x0(2) +rel*dx2

               ff(i) = abs_diff_sk_ku(chrom = (/x1(i), x2(i)/))
            enddo
         !$OMP END DO
         !$OMP END PARALLEL

         i  = minloc(ff, 1)
         x0 = (/x1(i), x2(i)/)

         if (ff(i)<1.e-8) exit
         write(*,*) x0(1), x0(2), ff(i)
      enddo

      x    = x0
      fvec = ff(i)
   return
   endsubroutine newton_raphson_downhill


   real(kind=R8) function abs_diff_sk_ku(chrom)
   implicit none
   real(kind=R8), intent(in), dimension(1:2):: chrom

      real(kind=R8) :: sk, ku

      call sk_ku(xx = chrom(1:2), sk = sk, ku = ku)
      abs_diff_sk_ku = abs(sk -SSK) +abs(ku -SKU)

   return
   endfunction abs_diff_sk_ku


   real(kind=R8) function cost_func(chrom)
   implicit none
   real(kind=R8), intent(in), dimension(1:2) :: chrom

      cost_func = 1./(1. + abs_diff_sk_ku(chrom(1:2)))

   return
   endfunction cost_func


   subroutine sk_ku(xx, sk, ku)
   implicit none
   real(kind=R8), intent(in ), dimension(1:2) :: xx
   real(kind=R8), intent(out)                 :: sk
   real(kind=R8), intent(out)                 :: ku

      real(kind=R8)    :: xa, xb, mu, si, a, b, un
      real(kind=R8)    :: h, hh, b1, b2, alp, bet
      integer(kind=I4) :: ia, ib, deb, fin, npts, long, i, k

      real(kind=R8), dimension(1:2) :: x
      long  = nn

      do k = 1, 2
        x(k) = max( xx(k), 1.e-4_R8 )
      enddo

      ia = long
      ib = long

      npts = long*long

      deb = 1    +ia
      fin = npts -ib

      a = x(1)
      b = x(2)

      hh = (2._R8 -a -b)/(npts -1)
      h  = (pi_R8/2)*hh

      xa = a +ia*hh
      xb = b +ib*hh

      b1 = -pi_R8/2 *(1.-a)
      b2 = +pi_R8/2 *(1.-b)

      alp = -(b2-npts*b1)/(b2-b1)
      bet =      (npts-1)/(b2-b1)

      un = 1._R8
      !. . . . . . . . . . . . . . . . . . . . . . . . . . . . .
      mu = log(1._R8/sin((pi_R8*xb)/2.)*sin((pi_R8*xa)/2.))
      mu = (un/h)*mu +add_tang(1, deb, fin, alp, bet, mu=0._R8, si=1._R8)
      do i = 1, ia-1
         mu = mu +tang(i*un, 1, alp, bet, mu=0._R8, si=1._R8)
      enddo
      do i = npts, npts -(ib -2), -1
         mu = mu +tang(i*un, 1, alp, bet, mu=0._R8, si=1._R8)
      enddo
      mu = mu/npts
      !. . . . . . . . . . . . . . . . . . . . . . . . . . . . .
      si = (pi_R8*(-2 + xa + xb))/2. - (mu**2*pi_R8*(-2 + xa + xb))/2. + 1._R8/tan((pi_R8*xa)/2.) + 1._R8/tan((pi_R8*xb)/2.) - 2*mu*Log(1._R8/sin((pi_R8*xb)/2.)*Sin((pi_R8*xa)/2.))
      si = (un/h)*si +add_tang(2, deb, fin, alp, bet, mu, si=1._R8)
      do i = 1, ia-1
         si = si+ tang(i*un, 2, alp, bet, mu, si=1._R8)
      enddo
      do i = npts, npts -(ib -2), -1
         si = si +tang(i*un, 2, alp, bet, mu, si=1._R8)
      enddo
      si = si/npts
      si = sqrt(si)
      !. . . . . . . . . . . . . . . . . . . . . . . . . . . . .
      sk = (6*mu*pi_R8 - 2*mu**3*pi_R8 - 3*mu*pi_R8*xa + mu**3*pi_R8*xa - 3*mu*pi_R8*xb + mu**3*pi_R8*xb - 6*mu*1._R8/tan((pi_R8*xa)/2.) - 6*mu*1._R8/tan((pi_R8*xb)/2.) - 1._R8/sin((pi_R8*xa)/2.)**2 + 1._R8/sin((pi_R8*xb)/2.)**2 - 2*Log(Sin((pi_R8*xa)/2.)) + 6*mu**2*Log(Sin((pi_R8*xa)/2.)) + 2*Log(Sin((pi_R8*xb)/2.)) - 6*mu**2*Log(Sin((pi_R8*xb)/2.)))/(2.*si**3)
      sk = (un/h)*sk +add_tang(3, deb, fin, alp, bet, mu, si)
      do i = 1, ia-1
         sk = sk +tang(i*un, 3, alp, bet, mu, si)
      enddo
      do i = npts, npts -(ib -2), -1
         sk = sk +tang(i*un, 3, alp, bet, mu, si)
      enddo
      sk = sk/npts
      !. . . . . . . . . . . . . . . . . . . . . . . . . . . . .
      ku = (6*pi_R8 - 36*mu**2*pi_R8 + 6*mu**4*pi_R8 - 3*pi_R8*xa + 18*mu**2*pi_R8*xa - 3*mu**4*pi_R8*xa - 3*pi_R8*xb + 18*mu**2*pi_R8*xb - 3*mu**4*pi_R8*xb + 4*(-2 + 9*mu**2)*1._R8/tan((pi_R8*xa)/2.) + 4*(-2 + 9*mu**2)*1._R8/tan((pi_R8*xb)/2.) + 12*mu*1._R8/sin((pi_R8*xa)/2.)**2 - 12*mu*1._R8/sin((pi_R8*xb)/2.)**2 + 24*mu*Log(Sin((pi_R8*xa)/2.)) - 24*mu**3*Log(Sin((pi_R8*xa)/2.)) - 24*mu*Log(Sin((pi_R8*xb)/2.)) + 24*mu**3*Log(Sin((pi_R8*xb)/2.)) + 1._R8/sin((pi_R8*xa)/2.)**4*Sin(pi_R8*xa) + 1._R8/sin((pi_R8*xb)/2.)**4*Sin(pi_R8*xb))/(6.*si**4)
      ku = (un/h)*ku +add_tang(4, deb, fin, alp, bet, mu, si)
      do i = 1, ia-1
         ku = ku +tang(i*un, 4, alp, bet, mu, si)
      enddo
      do i = npts, npts -(ib -2), -1
         ku = ku +tang(i*un, 4, alp, bet, mu, si)
      enddo
      ku = ku/npts
      !. . . . . . . . . . . . . . . . . . . . . . . . . . . . .

   return
   endsubroutine sk_ku


   real(kind=R8) function add_tang(n, deb, fin, alp, bet, mu, si)
   implicit none
   real(kind=R8),    intent(in) :: alp, bet, mu, si
   integer(kind=i4), intent(in) :: n, fin, deb

      real(kind=R8) :: xdeb, xfin

      xdeb = deb
      xfin = fin
      add_tang = (1./12)*( +9*(tang(xdeb +0.0_R8, n, alp, bet, mu, si)+tang(xfin -0.0_R8, n, alp, bet, mu, si)) &
                           +1*(tang(xdeb +1.0_R8, n, alp, bet, mu, si)+tang(xfin -1.0_R8, n, alp, bet, mu, si)) &
                           -4*(tang(xdeb +0.5_R8, n, alp, bet, mu, si)+tang(xfin -0.5_R8, n, alp, bet, mu, si)) )
   return
   endfunction add_tang


   real(kind=R8) function tang(xi, n, alp, bet, mu, si)
   implicit none
   real(kind=R8),    intent(in) :: xi, alp, bet, mu, si
   integer(kind=i4), intent(in) :: n

      real(kind=R8) :: tmp

      tmp = (xi +alp)/bet
      tang = ( (tan(tmp) -mu)/si )**n

   return
   endfunction tang


   subroutine calc_moments_1D(tab, mx, nb_mom, lg)
   implicit none
   integer(kind=I4) , intent(in )                  :: lg
   integer(kind=I4) , intent(in )                  :: nb_mom
   real(kind=R8)    , intent(in ), dimension(1:lg) :: tab
   type(moment_stat), intent(out)                  :: mx

      integer(kind=I4) :: i
      real(kind=R8)    :: tmp

      mx%mu = 0
      mx%si = 0
      mx%va = 0
      mx%Sk = 0
      mx%Ku = 0

      do i = 1, lg
         mx%mu = mx%mu +tab(i)/lg
      enddo
      if (nb_mom==1) return

      do i = 1, lg
         mx%va = mx%va +((tab(i) -mx%mu)**2)/lg
      enddo
      mx%si = sqrt( mx%va )
      if (nb_mom==2) return
      if (mx%si < 1.e-15_R8) then
         mx%Sk = 0
         mx%Ku = 0
      else
         do i = 1, lg
            tmp = (tab(i) -mx%mu)/mx%si
            mx%Sk = mx%Sk +(tmp**3)/lg
            mx%Ku = mx%Ku +(tmp**4)/lg
         enddo
      endif

   return
   endsubroutine calc_moments_1D


   subroutine profil_theo_trie_1D(tab, lg, x, mx)
   implicit none
   integer(kind=I4) , intent(in )                  :: lg
   real(kind=R8)    , intent(out), dimension(1:lg) :: tab
   real(kind=R8)    , intent(in ), dimension(1:2)  :: x
   type(moment_stat), intent(out)                  :: mx

      real(kind=R8)    :: b1, b2, alp, bet
      integer(kind=I4) :: i

      b1 = -PI_R8/2 *(1. -x(1))
      b2 = +PI_R8/2 *(1. -x(2))
      alp = -(b2- lg*b1)/(b2 -b1)
      bet =      (lg- 1)/(b2 -b1)
      do i = 1, lg
         tab(i) = tan( (i +alp)/bet )
      enddo
      call calc_moments_1D(tab, mx, nb_mom=4, lg=lg)
      tab(1:lg) = (tab(1:lg) -mx%mu)/mx%si

      mx%mu = 0._R8
      mx%si = 1._R8

   return
   endsubroutine profil_theo_trie_1D


   subroutine melange(tab, lg)
   implicit none
   integer(kind=i4), intent(in   )                  :: lg
   real(kind=r8)   , intent(inout), dimension(1:lg) :: tab

      real(kind=r8), dimension(1:lg) :: tmp
      integer(kind=i4) :: i

      call random_number(harvest=tmp)

      call sort_array2(tab_inout = tmp(1:lg),            &  !
                            tab1 = tab(1:lg), n = lg)       !

   return
   endsubroutine melange

endprogram test_algen

test_algen Program

Contents

Variables

Derived Types

Functions

Subroutines

Source Code

Uses

Calls

Variables

Derived Types

type :: moment_stat

Components

Functions

function abs_diff_sk_ku(chrom)

Arguments

Return Value real(kind=r8)

function add_tang(n, deb, fin, alp, bet, mu, si)

Arguments

Return Value real(kind=r8)

function cost_func(chrom)

Arguments

Return Value real(kind=r8)

function tang(xi, n, alp, bet, mu, si)

Arguments

Return Value real(kind=r8)

Subroutines

subroutine calc_moments_1D(tab, mx, nb_mom, lg)

Arguments

subroutine cost(me, x, f)

Arguments

subroutine melange(tab, lg)

Arguments

subroutine newton_raphson_downhill(x, fvec, eps, ndir, rel)

Arguments

subroutine profil_theo_trie_1D(tab, lg, x, mx)

Arguments

subroutine sk_ku(xx, sk, ku)

Arguments

Source Code