m_qtlmap_optimization
NAME
m_qtlmap_optimization
DESCRIPTION
set of functions to find a minimum of a objective function.
this module is an interface to offers :
* a generic interface for all qtlmap analysis
* a specific implementation for full and half sib family structure and likelihood function based on.
- We keep a compatibility with the oldest optimization NAG method e04jyf : http://www.nag.co.uk/numeric/FL/manual/xhtml/E04/e04jyf.xml
The NAG optimization is only available with a NAG compiler (xlf,f95)
source of implementation :
* L-BFGS : http://users.eecs.northwestern.edu/~nocedal/lbfgs.html
* NLOPT : http://ab-initio.mit.edu/wiki/index.php/NLopt
The module offers the following optimization :
1 : * NAG * E04JYF
2 : * L-BFGS-B * Broyden–Fletcher–Goldfarb–Shanno update to approximate the Hessian matrix (L-BFGS stands for (limited memory BFGS))
12 : * NLopt * DIRECT (global, no-derivative)
13 : * NLopt * DIRECT-L (global, no-derivative)
14 : * NLopt * Randomized DIRECT-L (global, no-derivative)
15 : * NLopt * Unscaled DIRECT (global, no-derivative)
16 : * NLopt * Unscaled DIRECT-L (global, no-derivative)
17 : * NLopt * Unscaled Randomized DIRECT-L (global, no-derivative)
18 : * NLopt * Original DIRECT version (global, no-derivative)
19 : * NLopt * Original DIRECT-L version (global, no-derivative)
23 : * NLopt * Limited-memory BFGS (L-BFGS) (local, derivative-based)
24 : * NLopt * Principal-axis, praxis (local, no-derivative)
25 : * NLopt * Limited-memory variable-metric, rank 1 (local, derivative-based)
26 : * NLopt * Limited-memory variable-metric, rank 2 (local, derivative-based)
27 : * NLopt * Truncated Newton (local, derivative-based)
28 : * NLopt * Truncated Newton with restarting (local, derivative-based)
29 : * NLopt * Preconditioned truncated Newton (local, derivative-based)
30 : * NLopt * Preconditioned truncated Newton with restarting (local, derivative-based)
31 : * NLopt * Controlled random search (CRS2) with local mutation (global, no-derivative)
32 : * NLopt * Multi-level single-linkage (MLSL), random (global, no-derivative)
33 : * NLopt * Multi-level single-linkage (MLSL), random (global, derivative)
34 : * NLopt * Multi-level single-linkage (MLSL), quasi-random (global, no-derivative)
35 : * NLopt * Multi-level single-linkage (MLSL), quasi-random (global, derivative)
36 : * NLopt * Method of Moving Asymptotes (MMA) (local, derivative)
37 : * NLopt * COBYLA (Constrained Optimization BY Linear Approximations) (local, no-derivative)
38 : * NLopt * NEWUOA unconstrained optimization via quadratic models (local, no-derivative)
39 : * NLopt * Bound-constrained optimization via NEWUOA-based quadratic models (local, no-derivative)
40 : * NLopt * Nelder-Mead simplex algorithm (local, no-derivative)
41 : * NLopt * Sbplx variant of Nelder-Mead (re-implementation of Rowan's Subplex) (local, no-derivative)
42 : * NLopt * Augmented Lagrangian method (local, no-derivative)
43 : * NLopt * Augmented Lagrangian method (local, derivative)
44 : * NLopt * Augmented Lagrangian method for equality constraints (local, no-derivative)
45 : * NLopt * Augmented Lagrangian method for equality constraints (local, derivative)
46 : * NLopt * BOBYQA bound-constrained optimization via quadratic models (local, no-derivative)
47 : * NLopt * ISRES evolutionary constrained optimization (global, no-derivative)
NOTES
SEE ALSO
LIBNAG_AVAILABLE
[ Top ] [ m_qtlmap_optimization ] [ Constants ]
NAME
LIBNAG_AVAILABLE
DESCRIPTION
flag to know if the nag e04jyf is available
SOURCE
135 #ifdef HAVE_LIBNAG 136 logical , parameter,private :: LIBNAG_AVAILABLE = .true. 137 #else 138 logical , parameter,private :: LIBNAG_AVAILABLE = .false. 139 #endif
OPT_CLASSIC
[ Top ] [ m_qtlmap_optimization ] [ Constants ]
NAME
OPT_CLASSIC
DESCRIPTION
Internal ID to execute a classical optimization
SOURCE
161 integer , parameter , private :: OPT_CLASSIC = 1
OPT_FAM
[ Top ] [ m_qtlmap_optimization ] [ Constants ]
NAME
OPT_FAM
DESCRIPTION
Internal ID to execute a likelihood optimization based on full sib family
SOURCE
170 integer , parameter , private :: OPT_FAM = 2
OPT_FAM_MULTI
[ Top ] [ m_qtlmap_optimization ] [ Constants ]
NAME
OPT_FAM_MULTI
DESCRIPTION
Internal ID to execute a likelihood optimization based on full sib family and multi trait
SOURCE
188 integer , parameter , private :: OPT_FAM_MULTI = 4
OPT_FAM_SIRE
[ Top ] [ m_qtlmap_optimization ] [ Constants ]
NAME
OPT_FAM_SIRE
DESCRIPTION
Internal ID to execute a likelihood optimization based on half sib family
SOURCE
179 integer , parameter , private :: OPT_FAM_SIRE = 3
OPTI_LAST
[ Top ] [ m_qtlmap_optimization ] [ Constants ]
NAME
OPTI_LAST
DESCRIPTION
Last index of optimization listed
SOURCE
104 integer , parameter , public :: OPTI_LAST = OPTI_LBFGS
OPTI_LBFGS
[ Top ] [ m_qtlmap_optimization ] [ Constants ]
NAME
OPTI_LBFGS
DESCRIPTION
constant id L-BFGS (Nocedal) optimization
NOTES
http://users.eecs.northwestern.edu/~nocedal/lbfgs.html
SOURCE
95 integer , parameter , public :: OPTI_LBFGS = 2
OPTI_NAG
[ Top ] [ m_qtlmap_optimization ] [ Constants ]
NAME
OPTI_NAG
DESCRIPTION
constant id e04jyf nag optimization
NOTES
http://www.nag.co.uk/numeric/FL/manual/xhtml/E04/e04jyf.xml
SOURCE
84 integer , parameter , public :: OPTI_NAG = 1
determ
[ Top ] [ m_qtlmap_optimization ] [ Variables ]
NAME
determ
DESCRIPTION
pointer : a determinant vector
DIMENSIONS
np
SOURCE
199 real(kind=dp),dimension(:), private, pointer :: determ
filter_vci
[ Top ] [ m_qtlmap_optimization ] [ Variables ]
NAME
filter_vci
DESCRIPTION
pointer : a boolean vector to active or not the computation of inv_vci
DIMENSIONS
np,nbnivest
inv_vci
[ Top ] [ m_qtlmap_optimization ] [ Variables ]
NAME
inv_vci
DESCRIPTION
pointer : the inverse of the covariance phenotypics trait matrix
DIMENSIONS
np,ncar,ncar
NB_OPTI_NLOPT
[ Top ] [ m_qtlmap_optimization ] [ Variables ]
NAME
NB_OPTI_NLOPT
DESCRIPTION
Number of nplot optimization available
SOURCE
113 integer , private :: NB_OPTI_NLOPT = 0
opti_user
[ Top ] [ m_qtlmap_optimization ] [ Variables ]
NAME
opti_user
DESCRIPTION
default optimization for a classical optimization of a objective function
SOURCE
148 #ifdef HAVE_LIBNAG 149 integer , public :: opti_user = OPTI_NAG 150 #else 151 integer , public :: opti_user = OPTI_LAST+16! OPTI_LBFGS 152 #endif
opti_user_fam
[ Top ] [ m_qtlmap_optimization ] [ Variables ]
NAME
opti_user_fam
DESCRIPTION
default optimization for a specific objective function (likelihood with a family organization)
SOURCE
122 #ifdef HAVE_NLOPT 123 integer , public :: opti_user_fam = OPTI_LAST+16 ! EQUI PNET in NLOPT 124 #else 125 integer , public :: opti_user_fam = OPTI_LBFGS 126 #endif
get_gradient
[ Top ] [ m_qtlmap_optimization ] [ Subroutines ]
NAME
get_gradient
DESCRIPTION
NOTES
SOURCE
1303 subroutine get_gradient(f,n,x,FUNCT1,s1,iuser,s2,user,u,l,g,np,fp,filter_fp,nm,fm,filter_fm) 1304 real(kind=dp) ,intent(in) :: f ! value in X 1305 integer ,intent(in) :: n,s1,s2,np,nm 1306 real(kind=dp) , dimension(n),intent(in) :: x,u,l 1307 real(kind=dp) ,intent(in),dimension(s2) :: user 1308 integer , intent(in),dimension(s1) :: iuser 1309 real(kind=dp) , dimension(n),intent(out):: g ! gradient result 1310 real (kind=dp) ,dimension(np) , intent(inout) :: fp 1311 real (kind=dp) ,dimension(nm) , intent(inout) :: fm 1312 logical ,dimension(np,n) , intent(in) :: filter_fp 1313 logical ,dimension(np,nm,n) , intent(in) :: filter_fm 1314 1315 external :: FUNCT1 1316 1317 !Local variable 1318 integer :: i 1319 double precision , dimension(n) :: xprime 1320 double precision :: fplusH,fmoinH,h2 1321 logical :: bordGaucheOk,bordDroitOk 1322 character(len=1) :: o 1323 #ifdef BENCHMARK_VIEW 1324 count_call_gradient=1+count_call_gradient 1325 #endif 1326 xprime = x 1327 ! do i=1,n 1328 ! print *,'x(',i,')=',x(i) 1329 ! end do 1330 ! stop 1331 ! Precision gradiant df/dx = f(x+h)-f(x-h)/2*h 1332 !calcul dérivé symetrique 1333 h2 = optim_H_PRECISION * 2 1334 g = 0.d0 1335 do i=1,n 1336 bordGaucheOk = l(i)+optim_H_PRECISION < x(i) 1337 bordDroitOk = u(i)-optim_H_PRECISION > x(i) 1338 1339 if ( bordGaucheOk .and. bordDroitOk ) then 1340 xprime(i) = x(i)-optim_H_PRECISION 1341 ! print *,n 1342 ! print *,xprime 1343 ! print *,fmoinH 1344 call FUNCT1(n,xprime,fmoinH,iuser,user) 1345 xprime(i) = x(i)+optim_H_PRECISION 1346 call FUNCT1(n,xprime,fplusH,iuser,user) 1347 g(i)=(fplusH-fmoinH)/h2 1348 xprime(i) = x(i) 1349 #ifdef BENCHMARK_VIEW 1350 count_call_funct=2*nmp(np+1)+count_call_funct 1351 #endif 1352 1353 else if ( bordGaucheOk ) then 1354 ! print *,'-- bord droit --' 1355 xprime(i) = x(i)-optim_H_PRECISION 1356 call FUNCT1(n,xprime,fmoinH,iuser,user) 1357 g(i)=(f-fmoinH)/optim_H_PRECISION 1358 xprime(i) = x(i) 1359 #ifdef BENCHMARK_VIEW 1360 count_call_funct=nmp(np+1)+count_call_funct 1361 #endif 1362 else if ( bordDroitOk ) then 1363 xprime(i) = x(i)+optim_H_PRECISION 1364 call FUNCT1(n,xprime,fplusH,iuser,user) 1365 g(i)=(fplusH-f)/optim_H_PRECISION 1366 xprime(i) = x(i) 1367 #ifdef BENCHMARK_VIEW 1368 count_call_funct=nmp(np+1)+count_call_funct 1369 #endif 1370 end if 1371 end do 1372 1373 end subroutine get_gradient
get_gradient_optim_fam
[ Top ] [ m_qtlmap_optimization ] [ Subroutines ]
NAME
get_gradient_optim_fam
DESCRIPTION
NOTES
SOURCE
1469 subroutine get_gradient_optim_fam(f,n,x,FUNCT_PART,s1,iuser,s2,user,u,l,g,np,fp,filter_fp,nm,fm,filter_fm) 1470 real(kind=dp) ,intent(in) :: f ! value in X 1471 integer ,intent(in) :: n,s1,s2,np,nm 1472 real(kind=dp) , dimension(n),intent(in) :: x,u,l 1473 real(kind=dp) ,intent(inout),dimension(s1) :: user 1474 integer , intent(inout),dimension(s2) :: iuser 1475 real(kind=dp) , dimension(n),intent(out) :: g ! gradient result 1476 real (kind=dp) ,dimension(np) , intent(inout) :: fp 1477 real (kind=dp) ,dimension(nm) , intent(inout) :: fm 1478 logical ,dimension(np,n) , intent(in) :: filter_fp 1479 logical ,dimension(np,nm,n) , intent(in) :: filter_fm 1480 1481 external :: FUNCT_PART 1482 1483 !Local variable 1484 integer :: i,ip,jm 1485 double precision , dimension(n) :: xprime 1486 double precision :: fplusH,fmoinH,h2,ftemp 1487 logical :: bordGaucheOk,bordDroitOk 1488 ! print *,"get_gradient_optim_fam:",size(iuser),size(user) 1489 #ifdef BENCHMARK_VIEW 1490 count_call_gradient=1+count_call_gradient 1491 #endif 1492 xprime = x 1493 h2 = optim_H_PRECISION * 2 1494 g = 0.d0 1495 1496 do i=1,n 1497 bordGaucheOk = l(i)+optim_H_PRECISION < x(i) 1498 bordDroitOk = u(i)-optim_H_PRECISION > x(i) 1499 1500 if ( bordGaucheOk .and. bordDroitOk ) then 1501 xprime(i) = x(i)-optim_H_PRECISION 1502 1503 fmoinH=0 1504 do ip=1,np 1505 do jm=nmp(ip)+1,nmp(ip+1) 1506 ! print *,"**>",ip,jm 1507 if ( filter_fm(ip,jm,i) ) then ! cet element affecte la vraissemblance de la faimme ip,jm 1508 call FUNCT_PART(ip,jm,n,xprime,ftemp,iuser,user) 1509 fmoinH=fmoinH+ftemp 1510 #ifdef BENCHMARK_VIEW 1511 count_call_funct=1+count_call_funct 1512 #endif 1513 else 1514 fmoinH=fmoinH+fm(jm) 1515 #ifdef BENCHMARK_VIEW 1516 count_avoid_funct=1+count_avoid_funct 1517 #endif 1518 end if 1519 end do 1520 end do 1521 1522 xprime(i) = x(i)+optim_H_PRECISION 1523 fplusH=0 1524 do ip=1,np 1525 do jm=nmp(ip)+1,nmp(ip+1) 1526 if ( filter_fm(ip,jm,i) ) then ! cet element affecte la vraissemblance de la faimme ip,jm 1527 call FUNCT_PART(ip,jm,n,xprime,ftemp,iuser,user) 1528 fplusH=fplusH+ftemp 1529 #ifdef BENCHMARK_VIEW 1530 count_call_funct=1+count_call_funct 1531 #endif 1532 else 1533 fplusH=fplusH+fm(jm) 1534 #ifdef BENCHMARK_VIEW 1535 count_avoid_funct=1+count_avoid_funct 1536 #endif 1537 end if 1538 end do 1539 end do 1540 1541 g(i)=(fplusH-fmoinH)/h2 1542 xprime(i) = x(i) 1543 1544 else if ( bordGaucheOk ) then 1545 ! print *,'-- bord droit --' 1546 xprime(i) = x(i)-optim_H_PRECISION 1547 fmoinH=0 1548 do ip=1,np 1549 do jm=nmp(ip)+1,nmp(ip+1) 1550 if ( filter_fm(ip,jm,i) ) then ! cet element affecte la vraissemblance de la famille ip,jm 1551 call FUNCT_PART(ip,jm,n,xprime,ftemp,iuser,user) 1552 fmoinH=fmoinH+ftemp 1553 #ifdef BENCHMARK_VIEW 1554 count_call_funct=1+count_call_funct 1555 #endif 1556 else 1557 fmoinH=fmoinH+fm(jm) 1558 #ifdef BENCHMARK_VIEW 1559 count_avoid_funct=1+count_avoid_funct 1560 #endif 1561 end if 1562 end do 1563 end do 1564 g(i)=(f-fmoinH)/optim_H_PRECISION 1565 xprime(i) = x(i) 1566 else if ( bordDroitOk ) then 1567 xprime(i) = x(i)+optim_H_PRECISION 1568 fplusH=0 1569 do ip=1,np 1570 do jm=nmp(ip)+1,nmp(ip+1) 1571 if ( filter_fm(ip,jm,i) ) then ! cet element affecte la vraissemblance de la famille ip,jm 1572 call FUNCT_PART(ip,jm,n,xprime,ftemp,iuser,user) 1573 fplusH=fplusH+ftemp 1574 #ifdef BENCHMARK_VIEW 1575 count_call_funct=1+count_call_funct 1576 #endif 1577 else 1578 fplusH=fplusH+fm(jm) 1579 #ifdef BENCHMARK_VIEW 1580 count_avoid_funct=1+count_avoid_funct 1581 #endif 1582 end if 1583 end do 1584 end do 1585 1586 g(i)=(fplusH-f)/optim_H_PRECISION 1587 xprime(i) = x(i) 1588 1589 end if 1590 end do 1591 1592 end subroutine get_gradient_optim_fam
get_gradient_optim_fam_multi
[ Top ] [ m_qtlmap_optimization ] [ Subroutines ]
NAME
get_gradient_optim_fam_multi
DESCRIPTION
NOTES
SOURCE
1904 subroutine get_gradient_optim_fam_multi(f,n,x,FUNCT_PART,s1,iuser,s2,user,u,l,g,np,fp,filter_fp,nm,fm,filter_fm) 1905 real(kind=dp) ,intent(in) :: f ! value in X 1906 integer ,intent(in) :: n,s1,s2,np,nm 1907 real(kind=dp) , dimension(n),intent(in) :: x,u,l 1908 real(kind=dp) ,intent(inout),dimension(s1) :: user 1909 integer , intent(inout),dimension(s2) :: iuser 1910 real(kind=dp) , dimension(n),intent(out) :: g ! gradient result 1911 real (kind=dp) ,dimension(np) , intent(in) :: fp 1912 real (kind=dp) ,dimension(nm) , intent(in) :: fm 1913 logical ,dimension(np,n) , intent(in) :: filter_fp 1914 logical ,dimension(np,nm,n) , intent(in) :: filter_fm 1915 logical :: valide 1916 external :: FUNCT_PART 1917 1918 !Local variable 1919 integer :: i,ip,jm,j 1920 real(kind=dp) , dimension(n) :: xprime 1921 real(kind=dp) :: fplusH,fmoinH,h2,ftemp,det 1922 real(kind=dp) , dimension(:,:),pointer :: vc 1923 real(kind=dp) , dimension(ncar,ncar),target :: vct 1924 logical :: bordGaucheOk,bordDroitOk 1925 ! print *,"get_gradient_optim_fam:",size(iuser),size(user) 1926 #ifdef BENCHMARK_VIEW 1927 count_call_gradient=1+count_call_gradient 1928 #endif 1929 xprime = x 1930 h2 = optim_H_PRECISION * 2 1931 g = 0.d0 1932 1933 do i=1,n 1934 bordGaucheOk = l(i)+optim_H_PRECISION < x(i) 1935 bordDroitOk = u(i)-optim_H_PRECISION > x(i) 1936 1937 if ( bordGaucheOk .and. bordDroitOk ) then 1938 xprime(i) = x(i)-optim_H_PRECISION 1939 1940 fmoinH=0 1941 do ip=1,np 1942 valide=.true. 1943 if (filter_vci(ip,i)) then 1944 call get_inv_residual_covariance_matrix(ip,n,xprime,vct,det,valide) 1945 #ifdef BENCHMARK_VIEW 1946 count_call_vci=1+count_call_vci 1947 #endif 1948 vc=>vct 1949 if (.not. valide) then 1950 fmoinH=INIFINY_REAL_VALUE 1951 end if 1952 else 1953 #ifdef BENCHMARK_VIEW 1954 count_avoid_vci=1+count_avoid_vci 1955 #endif 1956 vc => inv_vci(:,:,ip) 1957 det = determ(ip) 1958 if ( det == 0.d0 ) then 1959 fmoinH=INIFINY_REAL_VALUE 1960 valide=.false. 1961 end if 1962 end if 1963 1964 if (valide) then 1965 do jm=nmp(ip)+1,nmp(ip+1) 1966 ! print *,"**>",ip,jm 1967 if ( filter_fm(ip,jm,i) ) then ! cet element affecte la vraissemblance de la faimme ip,jm 1968 call FUNCT_PART(ip,jm,n,xprime,vc,det,ftemp,iuser,user) 1969 fmoinH=fmoinH+ftemp 1970 #ifdef BENCHMARK_VIEW 1971 count_call_funct=1+count_call_funct 1972 #endif 1973 else 1974 fmoinH=fmoinH+fm(jm) 1975 #ifdef BENCHMARK_VIEW 1976 count_avoid_funct=1+count_avoid_funct 1977 #endif 1978 end if 1979 end do 1980 end if 1981 end do 1982 1983 xprime(i) = x(i)+optim_H_PRECISION 1984 fplusH=0 1985 do ip=1,np 1986 valide=.true. 1987 if (filter_vci(ip,i)) then 1988 call get_inv_residual_covariance_matrix(ip,n,xprime,vct,det,valide) 1989 #ifdef BENCHMARK_VIEW 1990 count_call_vci=1+count_call_vci 1991 #endif 1992 vc=>vct 1993 if (.not. valide) then 1994 fmoinH=INIFINY_REAL_VALUE 1995 end if 1996 else 1997 #ifdef BENCHMARK_VIEW 1998 count_avoid_vci=1+count_avoid_vci 1999 #endif 2000 vc => inv_vci(:,:,ip) 2001 det = determ(ip) 2002 if ( det == 0 ) then 2003 fmoinH=INIFINY_REAL_VALUE 2004 valide=.false. 2005 end if 2006 end if 2007 if (valide) then 2008 do jm=nmp(ip)+1,nmp(ip+1) 2009 if ( filter_fm(ip,jm,i) ) then ! cet element affecte la vraissemblance de la faimme ip,jm 2010 call FUNCT_PART(ip,jm,n,xprime,vc,det,ftemp,iuser,user) 2011 fplusH=fplusH+ftemp 2012 #ifdef BENCHMARK_VIEW 2013 count_call_funct=1+count_call_funct 2014 #endif 2015 else 2016 fplusH=fplusH+fm(jm) 2017 #ifdef BENCHMARK_VIEW 2018 count_avoid_funct=1+count_avoid_funct 2019 #endif 2020 end if 2021 end do 2022 end if 2023 end do 2024 2025 g(i)=(fplusH-fmoinH)/h2 2026 xprime(i) = x(i) 2027 ! print *,g(i) 2028 ! stop 2029 2030 else if ( bordGaucheOk ) then 2031 ! print *,'-- bord droit --' 2032 xprime(i) = x(i)-optim_H_PRECISION 2033 fmoinH=0 2034 do ip=1,np 2035 valide=.true. 2036 if (filter_vci(ip,i)) then 2037 call get_inv_residual_covariance_matrix(ip,n,xprime,vct,det,valide) 2038 #ifdef BENCHMARK_VIEW 2039 count_call_vci=1+count_call_vci 2040 #endif 2041 vc=>vct 2042 if (.not. valide) then 2043 fmoinH=INIFINY_REAL_VALUE 2044 end if 2045 else 2046 #ifdef BENCHMARK_VIEW 2047 count_avoid_vci=1+count_avoid_vci 2048 #endif 2049 vc => inv_vci(:,:,ip) 2050 det = determ(ip) 2051 end if 2052 if (valide) then 2053 do jm=nmp(ip)+1,nmp(ip+1) 2054 if ( filter_fm(ip,jm,i) ) then ! cet element affecte la vraissemblance de la famille ip,jm 2055 call FUNCT_PART(ip,jm,n,xprime,vc,det,ftemp,iuser,user) 2056 fmoinH=fmoinH+ftemp 2057 #ifdef BENCHMARK_VIEW 2058 count_call_funct=1+count_call_funct 2059 #endif 2060 else 2061 fmoinH=fmoinH+fm(jm) 2062 #ifdef BENCHMARK_VIEW 2063 count_avoid_funct=1+count_avoid_funct 2064 #endif 2065 end if 2066 end do 2067 end if 2068 end do 2069 g(i)=(f-fmoinH)/optim_H_PRECISION 2070 xprime(i) = x(i) 2071 else if ( bordDroitOk ) then 2072 xprime(i) = x(i)+optim_H_PRECISION 2073 fplusH=0 2074 do ip=1,np 2075 valide=.true. 2076 if (filter_vci(ip,i)) then 2077 call get_inv_residual_covariance_matrix(ip,n,xprime,vct,det,valide) 2078 #ifdef BENCHMARK_VIEW 2079 count_call_vci=1+count_call_vci 2080 #endif 2081 vc=>vct 2082 if (.not. valide) then 2083 fmoinH=INIFINY_REAL_VALUE 2084 end if 2085 else 2086 #ifdef BENCHMARK_VIEW 2087 count_avoid_vci=1+count_avoid_vci 2088 #endif 2089 vc => inv_vci(:,:,ip) 2090 det = determ(ip) 2091 end if 2092 if (valide) then 2093 do jm=nmp(ip)+1,nmp(ip+1) 2094 if ( filter_fm(ip,jm,i) ) then ! cet element affecte la vraissemblance de la famille ip,jm 2095 call FUNCT_PART(ip,jm,n,xprime,vc,det,ftemp,iuser,user) 2096 fplusH=fplusH+ftemp 2097 #ifdef BENCHMARK_VIEW 2098 count_call_funct=1+count_call_funct 2099 #endif 2100 else 2101 fplusH=fplusH+fm(jm) 2102 #ifdef BENCHMARK_VIEW 2103 count_avoid_funct=1+count_avoid_funct 2104 #endif 2105 end if 2106 end do 2107 end if 2108 end do 2109 2110 g(i)=(fplusH-f)/optim_H_PRECISION 2111 xprime(i) = x(i) 2112 2113 end if 2114 end do 2115 2116 2117 end subroutine get_gradient_optim_fam_multi
get_gradient_optim_fam_sire
[ Top ] [ m_qtlmap_optimization ] [ Subroutines ]
NAME
get_gradient_optim_fam_sire
DESCRIPTION
NOTES
SOURCE
1644 subroutine get_gradient_optim_fam_sire(f,n,x,FUNCT_PART,s1,iuser,s2,user,u,l,g,np,fp,filter_fp,nm,fm,filter_fm) 1645 1646 real(kind=dp) ,intent(in) :: f ! value in X 1647 integer ,intent(in) :: n,s1,s2,np,nm 1648 real(kind=dp) , dimension(n),intent(in) :: x,u,l 1649 real(kind=dp) ,intent(inout),dimension(s1) :: user 1650 integer , intent(inout),dimension(s2) :: iuser 1651 real(kind=dp) , dimension(n),intent(out) :: g ! gradient result 1652 real (kind=dp) ,dimension(np) , intent(inout) :: fp 1653 real (kind=dp) ,dimension(nm) , intent(inout) :: fm 1654 logical ,dimension(np,n) , intent(in) :: filter_fp 1655 logical ,dimension(np,nm,n) , intent(in) :: filter_fm 1656 external :: FUNCT_PART 1657 1658 !Local variable 1659 integer :: i,ip 1660 double precision , dimension(n) :: xprime 1661 double precision :: fplusH,fmoinH,h2,ftemp 1662 logical :: bordGaucheOk,bordDroitOk 1663 ! print *,"get_gradient_optim_fam:",size(iuser),size(user) 1664 #ifdef BENCHMARK_VIEW 1665 count_call_gradient=1+count_call_gradient 1666 #endif 1667 xprime = x 1668 h2 = optim_H_PRECISION * 2 1669 g = 0.d0 1670 1671 1672 do i=1,n 1673 bordGaucheOk = l(i)+optim_H_PRECISION < x(i) 1674 bordDroitOk = u(i)-optim_H_PRECISION > x(i) 1675 1676 if ( bordGaucheOk .and. bordDroitOk ) then 1677 xprime(i) = x(i)-optim_H_PRECISION 1678 1679 fmoinH=0 1680 1681 do ip=1,np 1682 if ( filter_fp(ip,i) ) then ! cet element affecte la vraissemblance de la faimme ip 1683 call FUNCT_PART(ip,n,xprime,ftemp,iuser,user) 1684 fmoinH=fmoinH+ftemp 1685 #ifdef BENCHMARK_VIEW 1686 count_call_funct=nmp(ip+1)-nmp(ip)+count_call_funct 1687 #endif 1688 else 1689 fmoinH=fmoinH+fp(ip) 1690 #ifdef BENCHMARK_VIEW 1691 count_avoid_funct=nmp(ip+1)-nmp(ip)+count_avoid_funct 1692 #endif 1693 end if 1694 end do 1695 1696 1697 xprime(i) = x(i)+optim_H_PRECISION 1698 fplusH=0 1699 do ip=1,np 1700 if ( filter_fp(ip,i) ) then ! cet element affecte la vraissemblance de la faimme ip,jm 1701 call FUNCT_PART(ip,n,xprime,ftemp,iuser,user) 1702 fplusH=fplusH+ftemp 1703 #ifdef BENCHMARK_VIEW 1704 count_call_funct=nmp(ip+1)-nmp(ip)+count_call_funct 1705 #endif 1706 else 1707 fplusH=fplusH+fp(ip) 1708 #ifdef BENCHMARK_VIEW 1709 count_avoid_funct=nmp(ip+1)-nmp(ip)+count_avoid_funct 1710 #endif 1711 end if 1712 end do 1713 1714 1715 g(i)=(fplusH-fmoinH)/h2 1716 xprime(i) = x(i) 1717 1718 else if ( bordGaucheOk ) then 1719 ! print *,'-- bord droit --' 1720 xprime(i) = x(i)-optim_H_PRECISION 1721 fmoinH=0 1722 do ip=1,np 1723 if ( filter_fp(ip,i) ) then ! cet element affecte la vraissemblance de la famille ip,jm 1724 call FUNCT_PART(ip,n,xprime,ftemp,iuser,user) 1725 fmoinH=fmoinH+ftemp 1726 #ifdef BENCHMARK_VIEW 1727 count_call_funct=nmp(ip+1)-nmp(ip)+count_call_funct 1728 #endif 1729 else 1730 fmoinH=fmoinH+fp(ip) 1731 #ifdef BENCHMARK_VIEW 1732 count_avoid_funct=nmp(ip+1)-nmp(ip)+count_avoid_funct 1733 #endif 1734 end if 1735 end do 1736 g(i)=(f-fmoinH)/optim_H_PRECISION 1737 xprime(i) = x(i) 1738 else if ( bordDroitOk ) then 1739 xprime(i) = x(i)+optim_H_PRECISION 1740 fplusH=0 1741 do ip=1,np 1742 if ( filter_fp(ip,i) ) then ! cet element affecte la vraissemblance de la famille ip,jm 1743 call FUNCT_PART(ip,n,xprime,ftemp,iuser,user) 1744 fplusH=fplusH+ftemp 1745 #ifdef BENCHMARK_VIEW 1746 count_call_funct=nmp(ip+1)-nmp(ip)+count_call_funct 1747 #endif 1748 else 1749 fplusH=fplusH+fp(ip) 1750 #ifdef BENCHMARK_VIEW 1751 count_avoid_funct=nmp(ip+1)-nmp(ip)+count_avoid_funct 1752 #endif 1753 end if 1754 end do 1755 1756 g(i)=(fplusH-f)/optim_H_PRECISION 1757 xprime(i) = x(i) 1758 1759 end if 1760 end do 1761 1762 1763 1764 end subroutine get_gradient_optim_fam_sire
get_inv_residual_covariance_matrix
[ Top ] [ m_qtlmap_optimization ] [ Subroutines ]
NAME
get_inv_residual_covariance_matrix
DESCRIPTION
NOTES
FUNCTION FOR OPTIMISATION IP-JM-MULTITRAITS
SOURCE
1777 subroutine get_inv_residual_covariance_matrix(ip,n,x,vci,det,ok) 1778 integer , intent(in) :: ip,n 1779 real (kind=dp) ,dimension(n) , intent(in) :: x 1780 real(kind=dp),dimension(ncar,ncar) ,intent(out):: vci 1781 real(kind=dp) ,intent(out):: det 1782 logical ,intent(out) :: ok 1783 1784 real(kind=dp),dimension(ncar+1,ncar) :: VCIInverse 1785 1786 integer :: ic,jc,k,irank 1787 real(kind=dp) :: rh 1788 1789 ok=.true. 1790 1791 !VCi : residual covariance matrix from sire I 1792 ! | SIGi1^2 SIGi1*SIGi2*RHO(1,2) .... SIGi1*SIGip*RHO(1,p) 1793 ! | .. SIGi2^2 1794 ! VCi | .. ... 1795 ! | .. ... 1796 ! | .. 1797 ! | SIGi1*SIGip*RHO(1,p) SIGip^2 1798 k=0 1799 do ic=1,ncar 1800 do jc=ic,ncar 1801 VCIInverse(ic,jc)=x((ic-1)*np+ip)*x((jc-1)*np+ip) 1802 if ( ic /= jc ) then 1803 k=k+1 1804 !VCIInverse(ic,jc)=VCIInverse(ic,jc)*x(ncar*np+k)!rhoc(ic,jc) 1805 ! ** RESCALE OFI *** 1806 rh = x(ncar*np+k) 1807 ! print *,rh 1808 rh = ((dexp(rh)-1.d0)/(dexp(rh)+1.d0)) 1809 VCIInverse(ic,jc)=VCIInverse(ic,jc)*rh 1810 VCIInverse(jc,ic)=VCIInverse(ic,jc) 1811 end if 1812 end do 1813 end do 1814 1815 !Calcul Determinant de la matrice de covariance residuelle du pere ip 1816 !-------------------------------------------------------------------- 1817 irank=0 1818 det=0 1819 !Calcul Inverse de la matrice de covariance residuelle du pere ip 1820 !----------------------------------------------------------------- 1821 CALL MATH_QTLMAP_INVDETMATSYM(ncar,VCIInverse,ncar+1,det,irank) 1822 ! call MATH_QTLMAP_F01ADF(ncar,VCIInverse,ncar+1,irank) 1823 if (irank /= 0) then 1824 ok=.false. 1825 return 1826 end if 1827 1828 do ic=1,ncar 1829 do jc=1,ic 1830 VCI(ic,jc)=VCIInverse(ic+1,jc) 1831 VCI(jc,ic)=VCI(ic,jc) 1832 end do 1833 end do 1834 1835 1836 1837 end subroutine get_inv_residual_covariance_matrix
get_value
[ Top ] [ m_qtlmap_optimization ] [ Subroutines ]
NAME
get_value
DESCRIPTION
NOTES
SOURCE
1387 subroutine get_value(FUNCT1,n,x,f,s1,iuser,s2,user,np,fp,filter_fp,nm,fm,filter_fm) 1388 1389 integer , intent(in) :: n,s1,s2,np,nm 1390 real (kind=dp) ,dimension(n), intent(in) :: x 1391 real (kind=dp) , intent(inout) :: f 1392 integer , dimension(s1), intent(in) :: iuser 1393 real (kind=dp) ,dimension(s2), intent(in) :: user 1394 real (kind=dp) ,dimension(np) , intent(inout) :: fp 1395 real (kind=dp) ,dimension(nm) , intent(inout) :: fm 1396 logical ,dimension(np,n) , intent(in) :: filter_fp 1397 logical ,dimension(np,nm,n) , intent(in) :: filter_fm 1398 external :: FUNCT1 1399 integer :: jm,ip 1400 real (kind=dp) :: ft 1401 1402 ! print *," ** get_value **" 1403 ! stop 1404 call FUNCT1(n,x,f,iuser,user) 1405 1406 #ifdef BENCHMARK_VIEW 1407 count_call_funct=nm+count_call_funct 1408 #endif 1409 end subroutine get_value
get_value_optim_fam
[ Top ] [ m_qtlmap_optimization ] [ Subroutines ]
NAME
get_value_optim_fam
DESCRIPTION
NOTES
SOURCE
1428 subroutine get_value_optim_fam(FUNCT_PART,n,x,f,s1,iuser,s2,user,np,fp,filter_fp,nm,fm,filter_fm) 1429 1430 integer , intent(in) :: n,s1,s2,np,nm 1431 real (kind=dp) ,dimension(n), intent(in) :: x 1432 real (kind=dp) , intent(inout) :: f 1433 integer , dimension(s1), intent(inout) :: iuser 1434 real (kind=dp) ,dimension(s2), intent(inout) :: user 1435 real (kind=dp) ,dimension(np) , intent(inout) :: fp 1436 real (kind=dp) ,dimension(nm) , intent(inout) :: fm 1437 logical ,dimension(np,n) , intent(in) :: filter_fp 1438 logical ,dimension(np,nm,n) , intent(in) :: filter_fm 1439 1440 external :: FUNCT_PART 1441 integer :: jm,ip 1442 1443 f=0 1444 do ip=1,np 1445 do jm=nmp(ip)+1,nmp(ip+1) 1446 call FUNCT_PART(ip,jm,n,x,fm(jm),iuser,user) 1447 f=fm(jm)+f 1448 end do 1449 end do 1450 1451 1452 #ifdef BENCHMARK_VIEW 1453 count_call_funct=nm+count_call_funct 1454 #endif 1455 1456 end subroutine get_value_optim_fam
get_value_optim_fam_multi
[ Top ] [ m_qtlmap_optimization ] [ Subroutines ]
NAME
get_value_optim_fam_multi
DESCRIPTION
NOTES
SOURCE
1850 subroutine get_value_optim_fam_multi(FUNCT_PART,n,x,f,s1,iuser,s2,user,np,fp,filter_fp,nm,fm,filter_fm) 1851 1852 integer , intent(in) :: n,s1,s2,np,nm 1853 real (kind=dp) ,dimension(n), intent(in) :: x 1854 real (kind=dp) , intent(inout) :: f 1855 integer , dimension(s1), intent(inout) :: iuser 1856 real (kind=dp) ,dimension(s2), intent(inout) :: user 1857 real (kind=dp) ,dimension(np) , intent(inout) :: fp 1858 real (kind=dp) ,dimension(nm) , intent(inout) :: fm 1859 logical ,dimension(np,n) , intent(in) :: filter_fp 1860 logical ,dimension(np,nm,n) , intent(in) :: filter_fm 1861 1862 external :: FUNCT_PART 1863 integer :: jm,ip 1864 logical :: valide 1865 ! real (kind=dp) ,dimension(ncar,ncar) :: vct 1866 1867 f=0 1868 do ip=1,np 1869 call get_inv_residual_covariance_matrix(ip,n,x,inv_vci(:,:,ip),determ(ip),valide) 1870 #ifdef BENCHMARK_VIEW 1871 count_call_vci=1+count_call_vci 1872 #endif 1873 if (.not. valide) then 1874 determ(ip)=0.d0 1875 f=INIFINY_REAL_VALUE 1876 ! fp=INIFINY_REAL_VALUE 1877 fm=INIFINY_REAL_VALUE 1878 return 1879 end if 1880 1881 do jm=nmp(ip)+1,nmp(ip+1) 1882 call FUNCT_PART(ip,jm,n,x,inv_vci(:,:,ip),determ(ip),fm(jm),iuser,user) 1883 f=fm(jm)+f 1884 end do 1885 end do 1886 1887 #ifdef BENCHMARK_VIEW 1888 count_call_funct=nm+count_call_funct 1889 #endif 1890 1891 end subroutine get_value_optim_fam_multi
get_value_optim_fam_sire
[ Top ] [ m_qtlmap_optimization ] [ Subroutines ]
NAME
get_value_optim_fam_sire
DESCRIPTION
NOTES
SOURCE
1605 subroutine get_value_optim_fam_sire(FUNCT_PART,n,x,f,s1,iuser,s2,user,np,fp,filter_fp,nm,fm,filter_fm) 1606 1607 integer , intent(in) :: n,s1,s2,np,nm 1608 real (kind=dp) ,dimension(n), intent(in) :: x 1609 real (kind=dp) , intent(inout) :: f 1610 integer , dimension(s1), intent(inout) :: iuser 1611 real (kind=dp) ,dimension(s2), intent(inout) :: user 1612 real (kind=dp) ,dimension(np) , intent(inout) :: fp 1613 real (kind=dp) ,dimension(nm) , intent(inout) :: fm 1614 logical ,dimension(np,n) , intent(in) :: filter_fp 1615 logical ,dimension(np,nm,n) , intent(in) :: filter_fm 1616 external :: FUNCT_PART 1617 integer :: ip,i 1618 logical,dimension(N) :: keep_parameter 1619 1620 f=0 1621 do ip=1,np 1622 call FUNCT_PART(ip,n,x,fp(ip),iuser,user) 1623 f=fp(ip)+f 1624 end do 1625 1626 #ifdef BENCHMARK_VIEW 1627 count_call_funct=nm+count_call_funct 1628 #endif 1629 1630 1631 end subroutine get_value_optim_fam_sire
likelihood_empty
[ Top ] [ m_qtlmap_optimization ] [ Subroutines ]
NAME
likelihood_empty
DESCRIPTION
NOTES
SOURCE
2130 subroutine likelihood_empty(ip,jm,n,x,f,iuser,user) 2131 ! use OMP_LIB 2132 2133 integer , intent(in) :: ip,jm,n 2134 real (kind=dp) ,dimension(n), intent(in) :: x 2135 real (kind=dp) , intent(inout) :: f 2136 integer , dimension(1), intent(in) :: iuser 2137 real (kind=dp) ,dimension(1), intent(in) :: user 2138 2139 2140 end subroutine likelihood_empty
minimizing_e04jyf
[ Top ] [ m_qtlmap_optimization ] [ Subroutines ]
NAME
minimizing_e04jyf
DESCRIPTION
INPUTS
N : the number n of independent variables. N >=1
IBOUND : indicates whether the facility for dealing with bounds of special forms is to be used.
FUNCT1 : SUBROUTINE objective function
interface : FUNCT_PART(n,X,F,IUSER,USER)
BL : lower bounds
BU : upper bound
IUSER : user vector of integer
USER : user vector of real
INPUTS/OUTPUTS
X : start point and the point where rthe minimum was founded
OUTPUTS
F : final point stored in X
fm_user : value of the sub objective function full sib family . F = SUM FM(JM), 1<=JM<=NM
fp_user : value of the sub objective function half sib family . F = SUM FP(IP), 1<=IP<=NP
IFAIL : an error CODE
NOTES
Note: the dimension of the array IUSER and USER must be at least 1.
SOURCE
916 subroutine minimizing_e04jyf(N, & ! the number n of independent variables. N >=1 917 IBOUND, & ! indicates whether the facility for dealing with bounds of special forms is to be used. 918 FUNCT1, & ! SUBROUTINE, supplied by the user. External Procedure 919 BL, & ! lower bounds 920 BU, & ! upper bound 921 X, & ! vector W 922 F, & ! final point stored in X 923 IUSER, & ! Note: the dimension of the array IUSER must be at least 1. 924 USER, & ! Note: the dimension of the array IUSER must be at least 1. 925 FP, & ! Likelihood by Family of Sire 926 FILTER_FP, & ! Filter on The Family Sire 927 FM, & ! Likelihood by Family of SIre-Dam 928 FILTER_FM, & ! Filter on family Sire-Dam 929 IFAIL ) 930 931 932 integer,intent(in) :: N, IBOUND 933 integer,dimension(:) :: IUSER 934 real(kind=dp),dimension(N) :: BL, BU, X 935 real(kind=dp) :: F 936 real(kind=dp),dimension(:) :: USER 937 real(kind=dp),dimension(NP) ,intent(OUT) :: FP ! Likelihood by Family of Sire 938 logical ,dimension(NP,N),intent(IN) :: FILTER_FP ! Filter on The Family Sire 939 real(kind=dp),dimension(NM) ,intent(OUT) :: FM ! Likelihood by Family of SIre-Dam 940 logical ,dimension(NP,NM,N) ,intent(IN):: FILTER_FM ! Filter on family Sire-Dam 941 integer,intent(out) :: IFAIL 942 external :: FUNCT1 943 944 #ifdef HAVE_LIBNAG 945 integer,dimension(:),allocatable :: IW 946 integer :: LIW, LW 947 real(kind=dp),dimension(:),allocatable :: W 948 949 LIW = N+2 950 allocate (IW(LIW)) 951 LW =12*N+N*N+1 952 953 if ( LW < 13 ) LW = 13 954 955 allocate (W(LW)) 956 IFAIL=1 957 958 call e04jyf(N, IBOUND, FUNCT1, BL, BU, X, F, IW, LIW, W, LW,IUSER, USER, IFAIL) 959 !IFAIL = 1 960 ! 961 ! 962 ! On entry, N<1, 963 ! or IBOUND<0, 964 ! or IBOUND>3, 965 ! or IBOUND=0 and BLj>BUj for some j, 966 ! or IBOUND=3 and BL1>BU1, 967 ! or LIW<N+2, 968 ! or 969 ! 970 ! LW 971 ! < 972 ! max13,12×N+N×N-1/2 973 ! . 974 ! 975 !IFAIL = 2 976 ! 977 ! 978 ! There have been 400×n function evaluations, yet the algorithm does not seem to be converging. The calculations can be restarted from the final point held in X. The error may also indicate that Fx has no minimum. 979 ! 980 !IFAIL = 3 981 ! 982 ! 983 ! The conditions for a minimum have not all been met but a lower point could not be found and the algorithm has failed. 984 ! 985 !IFAIL = 4 986 ! 987 ! 988 ! An overflow has occurred during the computation. This is an unlikely failure, but if it occurs you should restart at the latest point given in X. 989 ! 990 !IFAIL = 5 991 ! 992 !IFAIL = 6 993 ! 994 !IFAIL = 7 995 ! 996 !IFAIL = 8 997 ! 998 ! 999 ! There is some doubt about whether the point x found by E04JYF is a minimum. The degree of confidence in the result decreases as IFAIL increases. Thus, when IFAIL =5 it is probable that the final x gives a good estimate of the position of a minimum, but when IFAIL =8 it is very unlikely that the routine has found a minimum. 1000 ! 1001 !IFAIL = 9 1002 ! 1003 ! 1004 ! In the search for a minimum, the modulus of one of the variables has become very large ∼106. This indicates that there is a mistake in FUNCT1, that your problem has no finite solution, or that the problem needs rescaling (see Section 8). 1005 ! 1006 !IFAIL = 10 1007 ! 1008 ! 1009 ! The computed set of forward-difference intervals (stored in W9×N+1,W9×N+2,…, 1010 ! W10×N) is such that Xi+W9×N+i≤Xi for some i. 1011 ! 1012 ! This is an unlikely failure, but if it occurs you should attempt to select another starting point. 1013 ! 1014 if ( IFAIL /= 0 ) then 1015 call log_mess('E04JYF IFAIL:'//str(IFAIL),WARNING_DEF) 1016 end if 1017 1018 deallocate(IW) 1019 deallocate (W) 1020 1021 #else 1022 call stop_application("No NAG environment") 1023 #endif 1024 1025 end subroutine minimizing_e04jyf
minimizing_funct
[ Top ] [ m_qtlmap_optimization ] [ Subroutines ]
NAME
minimizing_funct
DESCRIPTION
apply a classical optimization. The variable opti_user determines which method is used.
INPUTS
N : the number n of independent variables. N >=1
IBOUND : indicates whether the facility for dealing with bounds of special forms is to be used.
FUNCT1 : SUBROUTINE, supplied by the user. External Procedure
interface : FUNCT_PART(ip,jm,n,X,FM,IUSER,USER)
BL : lower bounds
BU : upper bound
IUSER : user vector of integer
USER : user vector of real
INPUTS/OUTPUTS
X : start point and the point where rthe minimum was founded
OUTPUTS
F : final point stored in X
IFAIL : an error CODE
NOTES
Note: the dimension of the array IUSER and USER must be at least 1.
SOURCE
452 subroutine minimizing_funct(N, & ! the number n of independent variables. N >=1 453 IBOUND, & ! indicates whether the facility for dealing with bounds of special forms is to be used. 454 FUNCT1, & ! SUBROUTINE, supplied by the user. External Procedure 455 BL, & ! lower bounds 456 BU, & ! upper bound 457 X, & ! vector W 458 F, & ! final point stored in X 459 IUSER, & ! Note: the dimension of the array IUSER must be at least 1. 460 USER, & ! Note: the dimension of the array IUSER must be at least 1. 461 IFAIL) 462 463 integer,intent(in) :: N, IBOUND 464 integer,dimension(:) :: IUSER 465 real(kind=dp),dimension(N) :: BL, BU, X 466 real(kind=dp) :: F 467 real(kind=dp),dimension(:) :: USER 468 integer,intent(out) :: IFAIL 469 external :: FUNCT1 470 471 real(kind=dp) ,dimension(NP) :: FP 472 real(kind=dp) ,dimension(NM) :: FM 473 logical ,dimension(NP,N) :: FILTER_FP 474 logical ,dimension(NP,NM,N) :: FILTER_FM 475 476 477 478 #ifdef BENCHMARK_VIEW 479 real(kind=dp) :: t1,t2,temp 480 temp=0 481 t1=second(); 482 count_call_funct = 0 483 count_call_gradient = 0 484 count_avoid_funct = 0 485 #endif 486 487 488 ! correction BUG 02/11/2010 : si n=0 aucun parametre a estimer on sort 489 if (N == 0 ) then 490 call log_mess("minimizing_funct is call with a objective function and none parameter",WARNING_DEF) 491 return 492 end if 493 494 select case (opti_user) 495 case (OPTI_NAG) 496 call minimizing_e04jyf(N, IBOUND, FUNCT1, BL, BU, X, F,IUSER, USER,FP,FILTER_FP,FM,FILTER_FM,IFAIL) 497 case (OPTI_LBFGS) 498 call minimizing_lbfgs(N, IBOUND, FUNCT1, BL, BU, X, F,IUSER, USER,FP,FILTER_FP,FM,FILTER_FM,& 499 IFAIL,OPT_CLASSIC,likelihood_empty) 500 case default 501 if ( opti_user > OPTI_LAST .and. opti_user <= (OPTI_LAST+NB_OPTI_NLOPT) ) then 502 call minimizing_nlopt(opti_user-OPTI_LAST,N, IBOUND, FUNCT1, BL, BU, X, F,IUSER, & 503 USER,FP,FILTER_FP,FM,FILTER_FM,IFAIL,OPT_CLASSIC,likelihood_empty) 504 else 505 call stop_application('bad value of opti_user['//str(opti_user)//'].') 506 end if 507 end select 508 509 IF (F .NE. F) THEN 510 F=INIFINY_REAL_VALUE 511 IFAIL=11 512 END IF 513 514 call log_mess(trim(name_optim(opti_user))//' MAX VALUE:'//str(f),VERBOSE_DEF) 515 516 #ifdef BENCHMARK_VIEW 517 t2=second()-t1 518 total_time=total_time+t2 519 if ((count_call_funct+count_avoid_funct)>0) temp=(count_avoid_funct/(count_call_funct+count_avoid_funct))*100 520 print '(a," ",f10.8,1x,f7.3,1x,a,i7,1x,a,i7,1x,a,i7,f4.1,"%",f7.3)',trim(name_optim(opti_user)),t2,f,'call funct:',& 521 count_call_funct/nm,'call gradient:',count_call_gradient,'avoid call:',count_avoid_funct/nm,temp,total_time 522 #endif 523 524 end subroutine minimizing_funct
minimizing_funct_family
[ Top ] [ m_qtlmap_optimization ] [ Subroutines ]
NAME
minimizing_funct_family
DESCRIPTION
apply a likelihood objective function optimization using the full sib family structuration and a flag by family
that determines if a parameter influence the likelihood of the family.
The variable opti_user_fam determines which method is used.
INPUTS
N : the number n of independent variables. N >=1
IBOUND : indicates whether the facility for dealing with bounds of special forms is to be used.
FUNCT_PART : SUBROUTINE of the sub objective function for the family ip,jm
interface : FUNCT_PART(ip,jm,n,X,FM,IUSER,USER)
filter_incid : boolean array that indicates if a parameter influence the likelihood on a specific family.
dimension : np,nm,N
BL : lower bounds
BU : upper bound
IUSER : user vector of integer
USER : user vector of real
INPUTS/OUTPUTS
X : start point and the point where rthe minimum was founded
OUTPUTS
F : final point stored in X
fm_user : value of the sub objective function full sib family . F = SUM FM(JM), 1<=JM<=NM
fp_user : value of the sub objective function half sib family . F = SUM FP(IP), 1<=IP<=NP
IFAIL : an error CODE
NOTES
Note: the dimension of the array IUSER and USER must be at least 1.
SOURCE
559 subroutine minimizing_funct_family(N, & ! the number n of independent variables. N >=1 560 IBOUND, & ! indicates whether the facility for dealing with bounds of special forms is to be used. 561 FUNCT_PART, & ! SUBROUTINE, supplied by the user Likelihood for the famili IP,JM. External Procedure 562 filter_incid,& 563 fm_user, & 564 fp_user, & 565 BL, & ! lower bounds 566 BU, & ! upper bound 567 X, & ! vector W 568 F, & ! final point stored in X 569 IUSER, & ! Note: the dimension of the array IUSER must be at least 1. 570 USER, & ! Note: the dimension of the array IUSER must be at least 1. 571 IFAIL) 572 573 integer,intent(in) :: N, IBOUND 574 integer,dimension(:) :: IUSER 575 real(kind=dp),dimension(N) :: BL, BU, X 576 real(kind=dp) ,dimension(nm) :: fm_user 577 real(kind=dp) ,dimension(np) :: fp_user 578 real(kind=dp) :: F 579 real(kind=dp),dimension(:) :: USER 580 integer,intent(out) :: IFAIL 581 logical, dimension(NP,NM,N) , intent(in) :: filter_incid 582 583 external :: FUNCT_PART 584 585 logical, dimension(NP,N) :: filter_fp 586 integer :: ip 587 588 #ifdef BENCHMARK_VIEW 589 real(kind=dp) :: t1,t2,temp 590 temp=0 591 t1=second(); 592 count_call_funct = 0 593 count_call_gradient = 0 594 count_avoid_funct = 0 595 #endif 596 ! correction BUG 02/11/2010 : si n=0 aucun parametre a estimer on sort 597 if (N == 0 ) then 598 call log_mess("minimizing_funct_family is call with a objective function and none parameter",WARNING_DEF) 599 return 600 end if 601 602 filter_fp = .false. 603 604 select case (opti_user_fam) 605 case (OPTI_NAG) 606 call stop_application('can not use this optimization ['//str(opti_user_fam)//'].') 607 case (OPTI_LBFGS) 608 call minimizing_lbfgs(N, IBOUND, get_value, BL, BU, X, F,IUSER, USER,FP_USER,FILTER_FP,FM_USER,filter_incid,& 609 IFAIL,OPT_FAM,FUNCT_PART) 610 case default 611 if ( opti_user_fam > OPTI_LAST .and. opti_user_fam <= (OPTI_LAST+NB_OPTI_NLOPT) ) then 612 call minimizing_nlopt(opti_user_fam-OPTI_LAST,N, IBOUND, get_value, BL, BU, X, F,IUSER, & 613 USER,FP_USER,FILTER_FP,FM_USER,filter_incid,IFAIL,OPT_FAM,FUNCT_PART) 614 else 615 call stop_application('bad value of opti_user['//str(opti_user_fam)//'].') 616 end if 617 end select 618 619 IF (F .NE. F) THEN 620 F=INIFINY_REAL_VALUE 621 IFAIL=11 622 END IF 623 624 do ip=1,np 625 fp_user(ip) = sum(fm_user(nmp(ip)+1:nmp(ip+1))) 626 end do 627 628 call log_mess(trim(name_optim(opti_user_fam))//' MAX VALUE:'//str(f),VERBOSE_DEF) 629 630 #ifdef BENCHMARK_VIEW 631 t2=second()-t1 632 total_time=total_time+t2 633 !print *,(count_call_funct+count_avoid_funct),count_call_funct,count_avoid_funct 634 if ((count_call_funct+count_avoid_funct)>0) then 635 temp=(real(count_avoid_funct)/real(count_call_funct+count_avoid_funct))*100 636 end if 637 !print *,temp 638 print '(a," ",f10.8,1x,f7.3,1x,a,i7,1x,a,i7,1x,a,i7,1x,f4.1,"%",f7.3)',trim(name_optim(opti_user_fam)),t2,f,'call funct:',& 639 count_call_funct/nm,'call gradient:',count_call_gradient,'avoid call:',count_avoid_funct/nm,temp,total_time 640 #endif 641 642 end subroutine minimizing_funct_family
minimizing_funct_family_multi
[ Top ] [ m_qtlmap_optimization ] [ Subroutines ]
NAME
minimizing_funct_family_multi
DESCRIPTION
apply a likelihood objective function optimization using the half sib family structuration and a flag by sire family
that determines if a parameter influence the likelihood of the sire family.
The variable opti_user_fam determines which method is used.
INPUTS
N : the number n of independent variables. N >=1
IBOUND : indicates whether the facility for dealing with bounds of special forms is to be used.
FUNCT_PART : SUBROUTINE of the sub objective function for the family ip,jm
interface : FUNCT_PART(ip,jm,n,X,vci,determ,FM,IUSER,USER)
filter_incid : boolean array that indicates if a parameter influence the likelihood on a specific sire family.
dimension : np,nm,N
filter_vcinv :
dimension : ncar,ncar,np
BL : lower bounds
BU : upper bound
IUSER : user vector of integer
USER : user vector of real
INPUTS/OUTPUTS
X : start point and the point where rthe minimum was founded
vci_temp : buffer array.dim : ncar,ncar,np
determ_temp : buffer array.dim : np
OUTPUTS
F : final point stored in X
fm_user : value of the sub objective function full sib family . F = SUM FM(JM), 1<=JM<=NM
fp_user : value of the sub objective function half sib family . F = SUM FP(IP), 1<=IP<=NP
IFAIL : an error CODE
NOTES
Note: the dimension of the array IUSER and USER must be at least 1.
SOURCE
790 subroutine minimizing_funct_family_multi(N, & ! the number n of independent variables. N >=1 791 IBOUND, & ! indicates whether the facility for dealing with bounds of special forms is to be used. 792 FUNCT_PART, & ! SUBROUTINE, supplied by the user Likelihood for the famili IP,JM. External Procedure 793 filter_incid,& 794 filter_vcinv,& 795 vci_temp,& 796 determ_temp,& 797 fm_user, & 798 fp_user, & 799 BL, & ! lower bounds 800 BU, & ! upper bound 801 X, & ! vector W 802 F, & ! final point stored in X 803 IUSER, & ! Note: the dimension of the array IUSER must be at least 1. 804 USER, & ! Note: the dimension of the array IUSER must be at least 1. 805 IFAIL) 806 807 integer,intent(in) :: N, IBOUND 808 integer,dimension(:) :: IUSER 809 real(kind=dp),dimension(N) :: BL, BU, X 810 real(kind=dp) ,dimension(NM) :: fm_user 811 real(kind=dp) ,dimension(NP) :: fp_user 812 real(kind=dp) ,dimension(ncar,ncar,np),target :: vci_temp 813 real(kind=dp) ,dimension(np),target :: determ_temp 814 real(kind=dp) :: F 815 real(kind=dp),dimension(:) :: USER 816 integer,intent(out) :: IFAIL 817 logical, dimension(NP,NM,N) , intent(in) :: filter_incid 818 logical, dimension(:,:) ,target , intent(in) ::filter_vcinv 819 820 external :: FUNCT_PART 821 logical, dimension(NP,N) :: filter_fp 822 823 integer :: ip 824 825 #ifdef BENCHMARK_VIEW 826 real(kind=dp) :: t1,t2,temp 827 temp=0 828 t1=second(); 829 count_call_funct = 0 830 count_call_gradient = 0 831 count_avoid_funct = 0 832 count_call_vci = 0 833 count_avoid_vci = 0 834 #endif 835 ! correction BUG 02/11/2010 : si n=0 aucun parametre a estimer on sort 836 if (N == 0 ) then 837 call log_mess("minimizing_funct_family_multi is call with a objective function and none parameter",WARNING_DEF) 838 return 839 end if 840 841 filter_vci=>filter_vcinv 842 determ=>determ_temp 843 determ=0.d0 844 inv_vci=>vci_temp 845 inv_vci=0.d0 846 847 select case (opti_user_fam) 848 case (OPTI_NAG) 849 call stop_application('can not use this optimization ['//str(opti_user_fam)//'].') 850 case (OPTI_LBFGS) 851 call minimizing_lbfgs(N, IBOUND, get_value, BL, BU, X, F,IUSER, USER, fp_user, FILTER_FP,fm_user,filter_incid,& 852 IFAIL,OPT_FAM_MULTI,FUNCT_PART) 853 case default 854 if ( opti_user_fam > OPTI_LAST .and. opti_user_fam <= (OPTI_LAST+NB_OPTI_NLOPT) ) then 855 call minimizing_nlopt(opti_user_fam-OPTI_LAST,N, IBOUND, get_value, BL, BU, X, F,IUSER, & 856 USER,fp_user, FILTER_FP,fm_user,filter_incid,IFAIL,OPT_FAM_MULTI,FUNCT_PART) 857 else 858 call stop_application('bad value of opti_user['//str(opti_user_fam)//'].') 859 end if 860 end select 861 862 863 IF (F .NE. F) THEN 864 F=INIFINY_REAL_VALUE 865 IFAIL=11 866 END IF 867 868 do ip=1,np 869 fp_user(ip) = sum(fm_user(nmp(ip)+1:nmp(ip+1))) 870 end do 871 872 call log_mess(trim(name_optim(opti_user_fam))//' MAX VALUE:'//str(f),VERBOSE_DEF) 873 874 #ifdef BENCHMARK_VIEW 875 t2=second()-t1 876 total_time=total_time+t2 877 !print *,(count_call_funct+count_avoid_funct),count_call_funct,count_avoid_funct 878 if ((count_call_funct+count_avoid_funct)>0) then 879 temp=(real(count_avoid_funct)/real(count_call_funct+count_avoid_funct))*100 880 end if 881 !print *,temp 882 print '(a," ",f10.8,1x,f7.3,1x,a,i10,1x,a,i10,1x,a,i9,1x,f4.1,"%",f7.3,i10,i10)',trim(name_optim(opti_user_fam)),t2,f,& 883 'call funct:',count_call_funct/nm,'call gradient:',count_call_gradient,'avoid call:',count_avoid_funct/nm,temp,total_time,& 884 count_call_vci,count_avoid_vci 885 #endif 886 887 end subroutine minimizing_funct_family_multi
minimizing_funct_family_sire
[ Top ] [ m_qtlmap_optimization ] [ Subroutines ]
NAME
minimizing_funct_family_sire
DESCRIPTION
apply a likelihood objective function optimization using the half sib family structuration and a flag by sire family
that determines if a parameter influence the likelihood of the sire family.
The variable opti_user_fam determines which method is used.
INPUTS
N : the number n of independent variables. N >=1
IBOUND : indicates whether the facility for dealing with bounds of special forms is to be used.
FUNCT_PART : SUBROUTINE of the sub objective function for the family ip,jm
interface : FUNCT_PART(ip,n,X,FM,IUSER,USER)
filter_incid : boolean array that indicates if a parameter influence the likelihood on a specific sire family.
dimension : np,N
BL : lower bounds
BU : upper bound
IUSER : user vector of integer
USER : user vector of real
INPUTS/OUTPUTS
X : start point and the point where rthe minimum was founded
OUTPUTS
F : final point stored in X
fp_user : value of the sub objective function half sib family . F = SUM FP(IP), 1<=IP<=NP
IFAIL : an error CODE
NOTES
Note: the dimension of the array IUSER and USER must be at least 1.
SOURCE
676 subroutine minimizing_funct_family_sire(N, & ! the number n of independent variables. N >=1 677 IBOUND, & ! indicates whether the facility for dealing with bounds of special forms is to be used. 678 FUNCT_PART, & ! SUBROUTINE, supplied by the user Likelihood for the famili IP,JM. External Procedure 679 filter_incid,& 680 fp_user, & 681 BL, & ! lower bounds 682 BU, & ! upper bound 683 X, & ! vector W 684 F, & ! final point stored in X 685 IUSER, & ! Note: the dimension of the array IUSER must be at least 1. 686 USER, & ! Note: the dimension of the array IUSER must be at least 1. 687 IFAIL) 688 689 integer,intent(in) :: N, IBOUND 690 integer,dimension(:) :: IUSER 691 real(kind=dp),dimension(N) :: BL, BU, X 692 real(kind=dp) ,dimension(np):: fp_user 693 real(kind=dp) :: F 694 real(kind=dp),dimension(:) :: USER 695 integer,intent(out) :: IFAIL 696 logical, dimension(NP,N) , intent(in) :: filter_incid 697 698 external :: FUNCT_PART 699 700 integer :: ip 701 real(kind=dp) ,dimension(NM) :: FM 702 logical, dimension(NP,NM,N) :: filter_fm 703 704 #ifdef BENCHMARK_VIEW 705 real(kind=dp) :: t1,t2,temp 706 temp=0 707 t1=second(); 708 count_call_funct = 0 709 count_call_gradient = 0 710 count_avoid_funct = 0 711 #endif 712 ! correction BUG 02/11/2010 : si n=0 aucun parametre a estimer on sort 713 if (N == 0 ) then 714 call log_mess("minimizing_funct_family_sire is call with a objective function and none parameter",WARNING_DEF) 715 return 716 end if 717 718 select case (opti_user_fam) 719 case (OPTI_NAG) 720 call stop_application('can not use this optimization ['//str(opti_user_fam)//'].') 721 case (OPTI_LBFGS) 722 call minimizing_lbfgs(N, IBOUND, get_value, BL, BU, X, F,IUSER, USER, fp_user, filter_incid, FM, filter_fm, & 723 IFAIL,OPT_FAM_SIRE,FUNCT_PART) 724 case default 725 if ( opti_user_fam > OPTI_LAST .and. opti_user_fam <= (OPTI_LAST+NB_OPTI_NLOPT) ) then 726 call minimizing_nlopt(opti_user_fam-OPTI_LAST,N, IBOUND, get_value, BL, BU, X, F,IUSER, USER,& 727 fp_user, filter_incid, FM,filter_fm,IFAIL,OPT_FAM_SIRE,FUNCT_PART) 728 else 729 call stop_application('bad value of opti_user['//str(opti_user_fam)//'].') 730 end if 731 end select 732 733 734 IF (F .NE. F) THEN 735 F=INIFINY_REAL_VALUE 736 IFAIL=11 737 END IF 738 739 call log_mess(trim(name_optim(opti_user_fam))//' MAX VALUE:'//str(f),VERBOSE_DEF) 740 741 #ifdef BENCHMARK_VIEW 742 t2=second()-t1 743 total_time=total_time+t2 744 !print *,(count_call_funct+count_avoid_funct),count_call_funct,count_avoid_funct 745 if ((count_call_funct+count_avoid_funct)>0) then 746 temp=(real(count_avoid_funct)/real(count_call_funct+count_avoid_funct))*100 747 end if 748 !print *,temp 749 print '(a," ",f10.8,1x,f7.3,1x,a,i10,1x,a,i10,1x,a,i9,1x,f4.1,"%",f7.3)',trim(name_optim(opti_user_fam)),t2,f,'call funct:',& 750 count_call_funct/nm,'call gradient:',count_call_gradient,'avoid call:',count_avoid_funct/nm,temp,total_time 751 #endif 752 753 end subroutine minimizing_funct_family_sire
minimizing_lbfgs
[ Top ] [ m_qtlmap_optimization ] [ Subroutines ]
NAME
minimizing_lbfgs
DESCRIPTION
NOTES
Note: the dimension of the array IUSER and USER must be at least 1.
SOURCE
1038 subroutine minimizing_lbfgs(N, & ! the number n of independent variables. N >=1 1039 IBOUND, & ! indicates whether the facility for dealing with bounds of special forms is to be used. 1040 FUNCT1, & ! SUBROUTINE, supplied by the user. External Procedure 1041 BL, & ! lower bounds 1042 BU, & ! upper bound 1043 X, & ! vector W 1044 F, & ! final point stored in X 1045 IUSER, & ! Note: the dimension of the array IUSER must be at least 1. 1046 USER, & ! Note: the dimension of the array IUSER must be at least 1. 1047 FP, & ! Likelihood by Family of Sire 1048 FILTER_FP, & ! Filter on The Family Sire 1049 FM, & ! Likelihood by Family of SIre-Dam 1050 FILTER_FM, & ! Filter on family Sire-Dam 1051 IFAIL, & 1052 OPTIM_FAMILY,& 1053 FUNCT_PART) 1054 1055 1056 integer :: N, IBOUND, IFAIL 1057 integer,dimension(:) :: IUSER 1058 real(kind=dp),dimension(N) :: BL, BU, X 1059 real(kind=dp) :: F 1060 real(kind=dp),dimension(:) :: USER 1061 real(kind=dp),dimension(NP) ,intent(OUT) :: FP ! Likelihood by Family of Sire 1062 logical ,dimension(NP,N),intent(IN) :: FILTER_FP ! Filter on The Family Sire 1063 real(kind=dp),dimension(NM) ,intent(OUT) :: FM ! Likelihood by Family of SIre-Dam 1064 logical ,dimension(NP,NM,N) ,intent(IN):: FILTER_FM ! Filter on family Sire-Dam 1065 integer ,intent(in) :: OPTIM_FAMILY 1066 external :: FUNCT1,FUNCT_PART 1067 1068 ! nmax is the dimension of the largest problem to be solved. 1069 ! mmax is the maximum number of limited memory corrections. 1070 1071 ! Declare the variables needed by the code. 1072 ! A description of all these variables is given at the end of 1073 ! the driver. 1074 1075 ! m is an INTEGER variable that must be set by the user to the 1076 ! number of corrections used in the limited memory matrix. 1077 ! It is not altered by the routine. Values of m < 3 are 1078 ! not recommended, and large values of m can result in excessive 1079 ! computing time. The range 3 <= m <= 20 is recommended. 1080 1081 integer ,parameter :: m = 15 1082 character(len=60) :: task, csave 1083 logical , dimension(4) :: lsave 1084 integer :: iprint,nbd(n), iwa(3*n), isave(44) 1085 double precision :: factr, pgtol 1086 double precision :: l(n), u(n), g(n), dsave(29) 1087 double precision :: wa(2*m*n+4*n+12*m*m+12*m) 1088 integer :: i,c,s1,s2 1089 double precision :: t1,t2 1090 1091 dsave=0.d0 1092 IFAIL = 0 1093 F=0.d0 1094 ! We wish to have output at every iteration. 1095 s1=size(iuser) 1096 s2=size(user) 1097 iprint = -1 1098 1099 ! We specify the tolerances in the stopping criteria. 1100 ! Typical values for factr: 1.d+12 for 1101 ! low accuracy; 1.d+7 for moderate accuracy; 1.d+1 for extremely 1102 ! high accuracy. 1103 factr=1.0d+10 ! par defaut 1.0d+7 1104 pgtol=optim_tolg!5.0d-2 ! par defaut 1.0d-5 1105 ! We specify the dimension n of the sample problem and the number 1106 ! m of limited memory corrections stored. (n and m should not 1107 ! exceed the limits nmax and mmax respectively.) 1108 1109 1110 ! We now provide nbd which defines the bounds on the variables: 1111 ! l specifies the lower bounds, 1112 ! u specifies the upper bounds. 1113 1114 do i=1,n 1115 nbd(i)=2 1116 l(i)=BL(i) 1117 u(i)=BU(i) 1118 end do 1119 1120 1121 ! We start the iteration by initializing task. 1122 ! 1123 task = 'START' 1124 ! ------- the beginning of the loop ---------- 1125 c=0 1126 111 continue 1127 if ( c > optim_maxeval ) then 1128 call log_mess("** WARNING : LBFGS : Maximum iteration was exceeded without convergence **",WARNING_DEF) 1129 goto 112 1130 end if 1131 ! This is the call to the L-BFGS-B code. 1132 ! write (*,fmt='(a1)',advance='no') '-' 1133 call setulb(n,m,x,l,u,nbd,f,g,factr,pgtol,wa,iwa,task,iprint,csave,lsave,isave,dsave) 1134 1135 ! write (*,fmt='(a5)',advance='no') '*' 1136 if (task(1:2) .eq. 'FG') then 1137 ! the minimization routine has returned to request the 1138 ! function f and gradient g values at the current x. 1139 ! Compute function value f for the sample problem. 1140 if ( OPTIM_FAMILY == OPT_FAM) then 1141 call get_value_optim_fam(FUNCT_PART,n,x,f,s1,iuser,s2,user,NP,FP,FILTER_FP,NM,FM,FILTER_FM) 1142 else if ( OPTIM_FAMILY == OPT_FAM_MULTI) then 1143 call get_value_optim_fam_multi(FUNCT_PART,n,x,f,s1,iuser,s2,user,NP,FP,FILTER_FP,NM,FM,FILTER_FM) 1144 else if ( OPTIM_FAMILY == OPT_FAM_SIRE) then 1145 call get_value_optim_fam_sire(FUNCT_PART,n,x,f,s1,iuser,s2,user,NP,FP,FILTER_FP,NM,FM,FILTER_FM) 1146 else if ( OPTIM_FAMILY == OPT_CLASSIC) then 1147 call FUNCT1(n,x,f,iuser,user) 1148 #ifdef BENCHMARK_VIEW 1149 count_call_funct=nmp(np+1)+count_call_funct 1150 #endif 1151 end if 1152 1153 c = c+1 1154 ! Compute gradient g for the sample problem. 1155 if ( OPTIM_FAMILY == OPT_FAM ) then 1156 call get_gradient_optim_fam(f,n,x,FUNCT_PART,s1,iuser,s2,user,u,l,g,NP,FP,FILTER_FP,NM,FM,FILTER_FM) 1157 else if ( OPTIM_FAMILY == OPT_FAM_MULTI ) then 1158 call get_gradient_optim_fam_multi(f,n,x,FUNCT_PART,s1,iuser,s2,user,u,l,g,NP,FP,FILTER_FP,NM,FM,FILTER_FM) 1159 else if ( OPTIM_FAMILY == OPT_FAM_SIRE) then 1160 call get_gradient_optim_fam_sire(f,n,x,FUNCT_PART,s1,iuser,s2,user,u,l,g,NP,FP,FILTER_FP,NM,FM,FILTER_FM) 1161 else if ( OPTIM_FAMILY == OPT_CLASSIC) then 1162 call get_gradient(f,n,x,FUNCT1,s1,iuser,s2,user,u,l,g,NP,FP,FILTER_FP,NM,FM,FILTER_FM) 1163 end if 1164 ! go back to the minimization routine. 1165 goto 111 1166 ! print *,'ENCORE' 1167 endif 1168 1169 ! 1170 ! print *,'test' 1171 if (task(1:5) .eq. 'NEW_X') goto 111 1172 ! call log_mess('LBFGD : task:'//trim(task),DEBUG_DEF) 1173 ! the minimization routine has returned with a new iterate, 1174 ! and we have opted to continue the iteration. 1175 1176 112 continue 1177 1178 ! ---------- the end of the loop ------------- 1179 ! If task is neither FG nor NEW_X we terminate execution. 1180 return 1181 end subroutine minimizing_lbfgs
minimizing_nlopt
[ Top ] [ m_qtlmap_optimization ] [ Subroutines ]
NAME
minimizing_nlopt
DESCRIPTION
NOTES
Note: the dimension of the array IUSER and USER must be at least 1.
SOURCE
1194 subroutine minimizing_nlopt(NLOPT_OPTI,& !optimisateur nlopt 1195 N, & ! the number n of independent variables. N >=1 1196 IBOUND, & ! indicates whether the facility for dealing with bounds of special forms is to be used. 1197 FUNCT1, & ! SUBROUTINE, supplied by the user. External Procedure 1198 BL, & ! lower bounds 1199 BU, & ! upper bound 1200 X, & ! vector W 1201 F, & ! final point stored in X 1202 IUSER, & ! Note: the dimension of the array IUSER must be at least 1. 1203 USER, & ! Note: the dimension of the array USER must be at least 1. 1204 FP, & ! Likelihood by Family of Sire 1205 FILTER_FP, & ! Filter on The Family Sire 1206 FM, & ! Likelihood by Family of SIre-Dam 1207 FILTER_FM, & ! Filter on family Sire-Dam 1208 IFAIL, & 1209 OPTIM_FAMILY,& 1210 FUNCT_PART) 1211 1212 1213 integer,intent(in) :: NLOPT_OPTI,N, IBOUND 1214 integer,dimension(:) :: IUSER 1215 real(kind=dp),dimension(N) :: BL, BU, X 1216 real(kind=dp) :: F 1217 real(kind=dp),dimension(:) :: USER 1218 real(kind=dp),dimension(NP) ,intent(OUT) :: FP ! Likelihood by Family of Sire 1219 logical ,dimension(NP,N),intent(IN) :: FILTER_FP ! Filter on The Family Sire 1220 real(kind=dp),dimension(NM) ,intent(OUT) :: FM ! Likelihood by Family of SIre-Dam 1221 logical ,dimension(NP,NM,N) ,intent(IN):: FILTER_FM ! Filter on family Sire-Dam 1222 integer,intent(out) :: IFAIL 1223 integer ,intent(in) :: OPTIM_FAMILY 1224 external :: FUNCT1,FUNCT_PART 1225 1226 #ifdef HAVE_NLOPT 1227 external :: interface_minimizing_nlopt_partial 1228 #endif 1229 integer :: s1,s2 1230 1231 s1=size(iuser) 1232 s2=size(user) 1233 1234 if ( OPTIM_FAMILY == OPT_FAM ) then 1235 #ifdef HAVE_NLOPT 1236 call interface_minimizing_nlopt_partial(NLOPT_OPTI,get_value_optim_fam,get_gradient_optim_fam,FUNCT_PART,n,x, & 1237 BL,BU,s1,iuser,s2,user,NP,FP,FILTER_FP,NM,FM,FILTER_FM,F,optim_tolx,optim_tolf,optim_maxeval,optim_maxtime,IFAIL) 1238 else if ( OPTIM_FAMILY == OPT_FAM_MULTI ) then 1239 call interface_minimizing_nlopt_partial(NLOPT_OPTI,get_value_optim_fam_multi,get_gradient_optim_fam_multi,& 1240 FUNCT_PART,n,x,BL,BU,s1,iuser,s2,user,NP,FP,FILTER_FP,NM,FM,FILTER_FM,F,& 1241 optim_tolx,optim_tolf,optim_maxeval,optim_maxtime,IFAIL) 1242 else if ( OPTIM_FAMILY == OPT_FAM_SIRE ) then 1243 call interface_minimizing_nlopt_partial(NLOPT_OPTI,get_value_optim_fam_sire,get_gradient_optim_fam_sire,& 1244 FUNCT_PART,n,x,BL,BU,s1,iuser,s2,user,NP,FP,FILTER_FP,NM,FM,FILTER_FM,F,& 1245 optim_tolx,optim_tolf,optim_maxeval,optim_maxtime,IFAIL) 1246 1247 else if ( OPTIM_FAMILY == OPT_CLASSIC ) then 1248 call interface_minimizing_nlopt_partial(NLOPT_OPTI,get_value,get_gradient,& 1249 FUNCT1,n,x,BL,BU,s1,iuser,s2,user,NP,FP,FILTER_FP,NM,FM,FILTER_FM,F,& 1250 optim_tolx,optim_tolf,optim_maxeval,optim_maxtime,IFAIL) 1251 #else 1252 call stop_application("** Can not perform a nlopt optimization ** ") 1253 #endif 1254 end if 1255 1256 if ( IFAIL /= 0 ) then 1257 select case (IFAIL) 1258 case (-4) 1259 call log_mess("NLOPT_FAILURE",ERROR_DEF);!stop; 1260 case (-3) 1261 call log_mess("NLOPT_INVALID_ARGS",ERROR_DEF);!stop; 1262 case (-2) 1263 call log_mess("NLOPT_OUT_OF_MEMORY",ERROR_DEF);!stop; 1264 case (-1) 1265 call log_mess("NLOPT_ROUNDOFF_LIMITED",ERROR_DEF);!stop; 1266 case (0) 1267 call log_mess("NLOPT SUCCESS ("//trim(name_optim(NLOPT_OPTI+OPTI_LAST))//") : "//trim(str(F)),INFO_DEF) 1268 case (1) 1269 call log_mess("NLOPT_MINF_MAX_REACHED ("//trim(name_optim(NLOPT_OPTI+OPTI_LAST))//") : "//trim(str(F)),VERBOSE_DEF) 1270 case (2) 1271 call log_mess("NLOPT_FTOL_REACHED ("//trim(name_optim(NLOPT_OPTI+OPTI_LAST))//") : "//trim(str(F)),VERBOSE_DEF) 1272 case (3) 1273 call log_mess("NLOPT_XTOL_REACHED ("//trim(name_optim(NLOPT_OPTI+OPTI_LAST))//") : "//trim(str(F)),VERBOSE_DEF) 1274 case (4) 1275 call log_mess("NLOPT_MAXEVAL_REACHED ("//trim(name_optim(NLOPT_OPTI+OPTI_LAST))//") : "//trim(str(F)),WARNING_DEF) 1276 case (5) 1277 call log_mess("NLOPT_MAXTIME_REACHED ("//trim(name_optim(NLOPT_OPTI+OPTI_LAST))//") : "//trim(str(F)),WARNING_DEF) 1278 case default 1279 call log_mess("THE METHOD FAILED :"//trim(str(ifail)),WARNING_DEF) 1280 end select 1281 end if 1282 1283 end subroutine minimizing_nlopt
name_optim
[ Top ] [ m_qtlmap_optimization ] [ Subroutines ]
NAME
name_optim
DESCRIPTION
get the optimization name
INPUTS
optim : the optimization id
OUTPUTS
name : the optimization name
SOURCE
343 function name_optim(optim) result(name) 344 integer , intent(in) :: optim 345 character(len=300) :: name 346 SELECT CASE (optim) 347 CASE (OPTI_NAG) 348 name ="E04JYF" 349 CASE (OPTI_LBFGS) 350 name = "LBFGS" 351 CASE (OPTI_LAST+1) 352 name = "NLOPT_GN_DIRECT" 353 CASE (OPTI_LAST+2) 354 name = "NLOPT_GN_DIRECT_L" 355 CASE (OPTI_LAST+3) 356 name = "NLOPT_GN_DIRECT_L_RAND" 357 CASE (OPTI_LAST+4) 358 name = "NLOPT_GN_DIRECT_NOSCAL" 359 CASE (OPTI_LAST+5) 360 name = "NLOPT_GN_DIRECT_L_NOSCAL" 361 CASE (OPTI_LAST+6) 362 name = "NLOPT_GN_DIRECT_L_RAND_NOSCAL" 363 CASE (OPTI_LAST+7) 364 name = "NLOPT_GN_ORIG_DIRECT" 365 CASE (OPTI_LAST+8) 366 name = "NLOPT_GN_ORIG_DIRECT_L" 367 CASE (OPTI_LAST+9) 368 name = "NLOPT_GD_STOGO" 369 CASE (OPTI_LAST+10) 370 name = "NLOPT_GD_STOGO_RAND" 371 CASE (OPTI_LAST+11) 372 name = "NLOPT_LD_LBFGS_NOCEDAL" 373 CASE (OPTI_LAST+12) 374 name = "NLOPT_LD_LBFGS" 375 CASE (OPTI_LAST+13) 376 name = "NLOPT_LN_PRAXIS" 377 CASE (OPTI_LAST+14) 378 name = "NLOPT_LD_VAR1" 379 CASE (OPTI_LAST+15) 380 name = "NLOPT_LD_VAR2" 381 CASE (OPTI_LAST+16) 382 name = "NLOPT_LD_TNEWTON" 383 CASE (OPTI_LAST+17) 384 name = "NLOPT_LD_TNEWTON_RESTART" 385 CASE (OPTI_LAST+18) 386 name = "NLOPT_LD_TNEWTON_PRECOND" 387 CASE (OPTI_LAST+19) 388 name = "NLOPT_LD_TNEWTON_PRECOND_RESTART" 389 CASE (OPTI_LAST+20) 390 name = "NLOPT_GN_CRS2_LM" 391 CASE (OPTI_LAST+21) 392 name = "NLOPT_GN_MLSL" 393 CASE (OPTI_LAST+22) 394 name = "NLOPT_GD_MLSL" 395 CASE (OPTI_LAST+23) 396 name = "NLOPT_GN_MLSL_LDS" 397 CASE (OPTI_LAST+24) 398 name = "NLOPT_GD_MLSL_LDS" 399 CASE (OPTI_LAST+25) 400 name = "NLOPT_LD_MMA" 401 CASE (OPTI_LAST+26) 402 name = "NLOPT_LN_COBYLA" 403 CASE (OPTI_LAST+27) 404 name = "NLOPT_LN_NEWUOA" 405 CASE (OPTI_LAST+28) 406 name = "NLOPT_LN_NEWUOA_BOUND" 407 CASE (OPTI_LAST+29) 408 name = "NLOPT_LN_NELDERMEAD" 409 CASE (OPTI_LAST+30) 410 name = "NLOPT_LN_SBPLX" 411 CASE (OPTI_LAST+31) 412 name = "NLOPT_LN_AUGLAG" 413 CASE (OPTI_LAST+32) 414 name = "NLOPT_LD_AUGLAG" 415 CASE (OPTI_LAST+33) 416 name = "NLOPT_LN_AUGLAG_EQ" 417 CASE (OPTI_LAST+34) 418 name = "NLOPT_LD_AUGLAG_EQ" 419 CASE (OPTI_LAST+35) 420 name = "NLOPT_LN_BOBYQA" 421 CASE (OPTI_LAST+36) 422 name = "NLOPT_GN_ISRES" 423 CASE DEFAULT 424 call stop_application("Unknown name optimize method ["//trim(str(optim))//"]") 425 END SELECT 426 end function
set_optimization
[ Top ] [ m_qtlmap_optimization ] [ Subroutines ]
NAME
set_optimization
DESCRIPTION
set the default optimization
INPUTS
value : a specific optimization id
NOTES
optimization available : see m_qtlmap_optimization
SOURCE
305 subroutine set_optimization(value) 306 integer, intent(in) :: value 307 integer :: i 308 external :: get_nb_opti_nlopt_qtlmap 309 310 #ifdef HAVE_NLOPT 311 call get_nb_opti_nlopt_qtlmap(NB_OPTI_NLOPT) 312 #else 313 NB_OPTI_NLOPT = 0 314 #endif 315 if ( value == OPTI_NAG .and. .not. LIBNAG_AVAILABLE ) then 316 call log_mess('**** Optimization NAG (e04jyf) is not available, qtlmap use :'//str(opti_user)) 317 return 318 end if 319 320 321 if ( (value >= OPTI_NAG) .and. (value <= OPTI_LAST+NB_OPTI_NLOPT) ) then 322 opti_user = value 323 end if 324 325 if ( value /= OPTI_NAG .and. (value <= OPTI_LAST+NB_OPTI_NLOPT)) then 326 opti_user_fam = value 327 end if 328 329 330 end subroutine