m_qtlmap_analyse_unitrait
[ Top ] [ ANALYSE ] [ Modules ]
NAME
m_qtlmap_analyse_unitrait
DESCRIPTION
NOTES
optinit,opti_0qtl,opti_1qtl
SEE ALSO
am
[ Top ] [ m_qtlmap_analyse_unitrait ] [ Variables ]
NAME
am
DESCRIPTION
The qtl dam effect found under H1
DIMENSIONS
nm
am2
[ Top ] [ m_qtlmap_analyse_unitrait ] [ Variables ]
NAME
am2
DESCRIPTION
The qtl dam effects found under H2, 1-->first qtl, 2-->2nd sqtl
DIMENSIONS
nm,2
ap
[ Top ] [ m_qtlmap_analyse_unitrait ] [ Variables ]
NAME
ap
DESCRIPTION
The qtl sire effect found under H1
DIMENSIONS
np
ap2
[ Top ] [ m_qtlmap_analyse_unitrait ] [ Variables ]
NAME
ap2
DESCRIPTION
The qtl sire effects found under H2, 1-->first qtl, 2-->2nd sqtl
DIMENSIONS
np,2
carpm
[ Top ] [ m_qtlmap_analyse_unitrait ] [ Variables ]
NAME
carpm
DESCRIPTION
initialization : loop before the minimization of the likelihood : opti_1qtl
DIMENSIONS
carpm1
[ Top ] [ m_qtlmap_analyse_unitrait ] [ Variables ]
NAME
carpm1
DESCRIPTION
initialization : loop before the minimization of the likelihood : opti_1car_2qtl
DIMENSIONS
carpm2
[ Top ] [ m_qtlmap_analyse_unitrait ] [ Variables ]
NAME
carpm2
DESCRIPTION
initialization : loop before the minimization of the likelihood : opti_1car_2qtl
DIMENSIONS
carpp
[ Top ] [ m_qtlmap_analyse_unitrait ] [ Variables ]
NAME
carpp
DESCRIPTION
initialization : loop before the minimization of the likelihood : opti_1qtl
DIMENSIONS
carpp1
[ Top ] [ m_qtlmap_analyse_unitrait ] [ Variables ]
NAME
carpp1
DESCRIPTION
initialization : loop before the minimization of the likelihood : opti_1car_2qtl
DIMENSIONS
carpp2
[ Top ] [ m_qtlmap_analyse_unitrait ] [ Variables ]
NAME
carpp2
DESCRIPTION
initialization : loop before the minimization of the likelihood : opti_1car_2qtl
DIMENSIONS
carppdf
[ Top ] [ m_qtlmap_analyse_unitrait ] [ Variables ]
NAME
carppdf
DESCRIPTION
initialization : loop before the minimization of the likelihood : opti_1qtl
DIMENSIONS
carppdf11
[ Top ] [ m_qtlmap_analyse_unitrait ] [ Variables ]
NAME
carppdf11
DESCRIPTION
initialization : loop before the minimization of the likelihood : opti_1car_2qtl
DIMENSIONS
carppdf12
[ Top ] [ m_qtlmap_analyse_unitrait ] [ Variables ]
NAME
carppdf12
DESCRIPTION
initialization : loop before the minimization of the likelihood : opti_1car_2qtl
DIMENSIONS
carppdf22
[ Top ] [ m_qtlmap_analyse_unitrait ] [ Variables ]
NAME
carppdf22
DESCRIPTION
initialization : loop before the minimization of the likelihood : opti_1car_2qtl
DIMENSIONS
current_chr
[ Top ] [ m_qtlmap_analyse_unitrait ] [ Variables ]
NAME
current_chr
DESCRIPTION
the current chromosome while the likelihood calculus
current_chr2
[ Top ] [ m_qtlmap_analyse_unitrait ] [ Variables ]
NAME
current_chr2
DESCRIPTION
the 2nd chromosome while the likelihood calculus (2QTL)
current_ic
[ Top ] [ m_qtlmap_analyse_unitrait ] [ Variables ]
NAME
current_ic
DESCRIPTION
the current trait
f0
[ Top ] [ m_qtlmap_analyse_unitrait ] [ Variables ]
NAME
f0
DESCRIPTION
The likelihood under H0
f_1
[ Top ] [ m_qtlmap_analyse_unitrait ] [ Variables ]
NAME
f_1
DESCRIPTION
The maximum likelihood founded under H1
fm0
[ Top ] [ m_qtlmap_analyse_unitrait ] [ Variables ]
NAME
fm0
DESCRIPTION
The likelihood by full sib family under H0
DIMENSIONS
nm
fm01
[ Top ] [ m_qtlmap_analyse_unitrait ] [ Variables ]
NAME
fm01
DESCRIPTION
DIMENSIONS
nm
fm1
[ Top ] [ m_qtlmap_analyse_unitrait ] [ Variables ]
NAME
fm1
DESCRIPTION
buffer to get the likelihood founded by full sib family under H1 at the current position
DIMENSIONS
nm
fm_1
[ Top ] [ m_qtlmap_analyse_unitrait ] [ Variables ]
NAME
fm_1
DESCRIPTION
The maximum likelihood founded by full sib family under H1
DIMENSIONS
np
fm_2
[ Top ] [ m_qtlmap_analyse_unitrait ] [ Variables ]
NAME
fm_2
DESCRIPTION
DIMENSIONS
nm
fp0
[ Top ] [ m_qtlmap_analyse_unitrait ] [ Variables ]
NAME
fp0
DESCRIPTION
The likelihood by half sib family under H0
DIMENSIONS
np
fp01
[ Top ] [ m_qtlmap_analyse_unitrait ] [ Variables ]
NAME
fp01
DESCRIPTION
DIMENSIONS
np
fp1
[ Top ] [ m_qtlmap_analyse_unitrait ] [ Variables ]
NAME
fp1
DESCRIPTION
buffer to get the likelihood founded by half sib family under H1 at the current position
DIMENSIONS
np
fp_1
[ Top ] [ m_qtlmap_analyse_unitrait ] [ Variables ]
NAME
fp_1
DESCRIPTION
The maximum likelihood founded by sire family under H1
DIMENSIONS
np
fp_2
[ Top ] [ m_qtlmap_analyse_unitrait ] [ Variables ]
NAME
fp_2
DESCRIPTION
DIMENSIONS
np
sig1
[ Top ] [ m_qtlmap_analyse_unitrait ] [ Variables ]
NAME
sig1
DESCRIPTION
the standart deviation (residual) found under H0
DIMENSIONS
np
sompm
[ Top ] [ m_qtlmap_analyse_unitrait ] [ Variables ]
NAME
sompm
DESCRIPTION
initialization : loop before the minimization of the likelihood : opti_1qtl
DIMENSIONS
sompm1
[ Top ] [ m_qtlmap_analyse_unitrait ] [ Variables ]
NAME
sompm1
DESCRIPTION
initialization : loop before the minimization of the likelihood : opti_1car_2qtl
DIMENSIONS
sompm1m2
[ Top ] [ m_qtlmap_analyse_unitrait ] [ Variables ]
NAME
sompm1m2
DESCRIPTION
initialization : loop before the minimization of the likelihood : opti_1car_2qtl
DIMENSIONS
sompm1y
[ Top ] [ m_qtlmap_analyse_unitrait ] [ Variables ]
NAME
sompm1y
DESCRIPTION
initialization : loop before the minimization of the likelihood : opti_1car_2qtl
DIMENSIONS
sompm2
[ Top ] [ m_qtlmap_analyse_unitrait ] [ Variables ]
NAME
sompm2
DESCRIPTION
initialization : loop before the minimization of the likelihood : opti_1car_2qtl
DIMENSIONS
sompm2y
[ Top ] [ m_qtlmap_analyse_unitrait ] [ Variables ]
NAME
sompm2y
DESCRIPTION
initialization : loop before the minimization of the likelihood : opti_1car_2qtl
DIMENSIONS
sompmy
[ Top ] [ m_qtlmap_analyse_unitrait ] [ Variables ]
NAME
sompmy
DESCRIPTION
initialization : loop before the minimization of the likelihood : opti_1qtl
DIMENSIONS
sompp
[ Top ] [ m_qtlmap_analyse_unitrait ] [ Variables ]
NAME
sompp
DESCRIPTION
initialization : loop before the minimization of the likelihood : opti_1qtl
DIMENSIONS
sompp1
[ Top ] [ m_qtlmap_analyse_unitrait ] [ Variables ]
NAME
sompp1
DESCRIPTION
initialization : loop before the minimization of the likelihood : opti_1car_2qtl
DIMENSIONS
sompp1m1
[ Top ] [ m_qtlmap_analyse_unitrait ] [ Variables ]
NAME
sompp1m1
DESCRIPTION
initialization : loop before the minimization of the likelihood : opti_1car_2qtl
DIMENSIONS
sompp1m2
[ Top ] [ m_qtlmap_analyse_unitrait ] [ Variables ]
NAME
sompp1m2
DESCRIPTION
initialization : loop before the minimization of the likelihood : opti_1car_2qtl
DIMENSIONS
sompp1p2
[ Top ] [ m_qtlmap_analyse_unitrait ] [ Variables ]
NAME
sompp1p2
DESCRIPTION
initialization : loop before the minimization of the likelihood : opti_1car_2qtl
DIMENSIONS
sompp1y
[ Top ] [ m_qtlmap_analyse_unitrait ] [ Variables ]
NAME
sompp1y
DESCRIPTION
initialization : loop before the minimization of the likelihood : opti_1car_2qtl
DIMENSIONS
sompp2
[ Top ] [ m_qtlmap_analyse_unitrait ] [ Variables ]
NAME
sompp2
DESCRIPTION
initialization : loop before the minimization of the likelihood : opti_1car_2qtl
DIMENSIONS
sompp2m1
[ Top ] [ m_qtlmap_analyse_unitrait ] [ Variables ]
NAME
sompp2m1
DESCRIPTION
initialization : loop before the minimization of the likelihood : opti_1car_2qtl
DIMENSIONS
sompp2m2
[ Top ] [ m_qtlmap_analyse_unitrait ] [ Variables ]
NAME
sompp2m2
DESCRIPTION
initialization : loop before the minimization of the likelihood : opti_1car_2qtl
DIMENSIONS
sompp2y
[ Top ] [ m_qtlmap_analyse_unitrait ] [ Variables ]
NAME
sompp2y
DESCRIPTION
initialization : loop before the minimization of the likelihood : opti_1car_2qtl
DIMENSIONS
somppdf
[ Top ] [ m_qtlmap_analyse_unitrait ] [ Variables ]
NAME
somppdf
DESCRIPTION
initialization : loop before the minimization of the likelihood : opti_1qtl
DIMENSIONS
somppdf1
[ Top ] [ m_qtlmap_analyse_unitrait ] [ Variables ]
NAME
somppdf1
DESCRIPTION
initialization : loop before the minimization of the likelihood : opti_1car_2qtl
DIMENSIONS
somppdf2
[ Top ] [ m_qtlmap_analyse_unitrait ] [ Variables ]
NAME
somppdf2
DESCRIPTION
initialization : loop before the minimization of the likelihood : opti_1car_2qtl
DIMENSIONS
somppm
[ Top ] [ m_qtlmap_analyse_unitrait ] [ Variables ]
NAME
somppm
DESCRIPTION
initialization : loop before the minimization of the likelihood : opti_1qtl
DIMENSIONS
somppy
[ Top ] [ m_qtlmap_analyse_unitrait ] [ Variables ]
NAME
somppy
DESCRIPTION
initialization : loop before the minimization of the likelihood : opti_1qtl
DIMENSIONS
somppydf
[ Top ] [ m_qtlmap_analyse_unitrait ] [ Variables ]
NAME
somppydf
DESCRIPTION
initialization : loop before the minimization of the likelihood : opti_1qtl
DIMENSIONS
somppydf1
[ Top ] [ m_qtlmap_analyse_unitrait ] [ Variables ]
NAME
somppydf1
DESCRIPTION
initialization : loop before the minimization of the likelihood : opti_1car_2qtl
DIMENSIONS
somppydf2
[ Top ] [ m_qtlmap_analyse_unitrait ] [ Variables ]
NAME
somppydf2
DESCRIPTION
initialization : loop before the minimization of the likelihood : opti_1car_2qtl
DIMENSIONS
std
[ Top ] [ m_qtlmap_analyse_unitrait ] [ Variables ]
NAME
std
DESCRIPTION
the standart deviation (residual) found under H1
DIMENSIONS
np
std2
[ Top ] [ m_qtlmap_analyse_unitrait ] [ Variables ]
NAME
std2
DESCRIPTION
the standart deviation (residual) found under H2
DIMENSIONS
np
xmoym
[ Top ] [ m_qtlmap_analyse_unitrait ] [ Variables ]
NAME
xmoym
DESCRIPTION
The polygenic dam effect found under H1
DIMENSIONS
nm
xmoym2
[ Top ] [ m_qtlmap_analyse_unitrait ] [ Variables ]
NAME
xmoym2
DESCRIPTION
The polygenic dam effect found under H2
DIMENSIONS
np
xmoyp
[ Top ] [ m_qtlmap_analyse_unitrait ] [ Variables ]
NAME
xmoyp
DESCRIPTION
The polygenic sire effect found under H1
DIMENSIONS
np
xmoyp2
[ Top ] [ m_qtlmap_analyse_unitrait ] [ Variables ]
NAME
xmoyp2
DESCRIPTION
TThe polygenic sire effect found under H2
DIMENSIONS
np
xmu1m
[ Top ] [ m_qtlmap_analyse_unitrait ] [ Variables ]
NAME
xmu1m
DESCRIPTION
The polygenic dam effect found under H0
DIMENSIONS
nm
xmu1p
[ Top ] [ m_qtlmap_analyse_unitrait ] [ Variables ]
NAME
xmu1p
DESCRIPTION
The polygenic sire effect found under H0
DIMENSIONS
np
end_analyse_unitrait_1QTL
[ Top ] [ m_qtlmap_analyse_unitrait ] [ Subroutines ]
NAME
end_analyse_unitrait_1QTL
DESCRIPTION
SOURCE
895 subroutine end_analyse_unitrait_1QTL 896 deallocate (sig1) 897 deallocate (xmu1p) 898 deallocate (xmu1m) 899 deallocate (fp0) 900 deallocate (fm0) 901 deallocate (ap) 902 deallocate (xmoyp) 903 deallocate (std) 904 deallocate (xmoym) 905 deallocate (am) 906 deallocate (fp_1) 907 deallocate (fm_1) 908 end subroutine end_analyse_unitrait_1QTL
end_analyse_unitrait_2qtl
[ Top ] [ m_qtlmap_analyse_unitrait ] [ Subroutines ]
NAME
end_analyse_unitrait_2qtl
DESCRIPTION
SOURCE
918 subroutine end_analyse_unitrait_2qtl 919 deallocate (ap2) 920 deallocate (xmoyp2) 921 deallocate (std2) 922 deallocate (fp_2) 923 deallocate (am2) 924 deallocate (xmoym2) 925 deallocate (fm_2) 926 927 end subroutine end_analyse_unitrait_2qtl
end_sub_1qtl
[ Top ] [ m_qtlmap_analyse_unitrait ] [ Subroutines ]
NAME
end_sub_1qtl
DESCRIPTION
SOURCE
738 subroutine end_sub_1qtl 739 deallocate (sompp) 740 deallocate (sompm) 741 deallocate (sompmy) 742 deallocate (carpp) 743 deallocate (carpm) 744 deallocate (somppy) 745 deallocate (somppm) 746 deallocate (somppdf) 747 deallocate (carppdf) 748 deallocate (somppydf) 749 deallocate (fp1) 750 deallocate (fm1) 751 end subroutine end_sub_1qtl
end_sub_2qtl
[ Top ] [ m_qtlmap_analyse_unitrait ] [ Subroutines ]
NAME
end_sub_2qtl
DESCRIPTION
SOURCE
858 subroutine end_sub_2qtl 859 deallocate (sompp1) 860 deallocate (sompp2) 861 deallocate (sompm1) 862 deallocate (sompm2) 863 deallocate (carpp1) 864 deallocate (carpp2) 865 deallocate (carpm1) 866 deallocate (carpm2) 867 deallocate (sompp1y) 868 deallocate (sompp2y) 869 deallocate (sompm1y) 870 deallocate (sompm2y) 871 deallocate (sompp1p2) 872 deallocate (sompp1m1) 873 deallocate (sompp1m2) 874 deallocate (sompp2m1) 875 deallocate (sompp2m2) 876 deallocate (sompm1m2) 877 deallocate (somppdf1) 878 deallocate (somppdf2) 879 deallocate (carppdf11) 880 deallocate (carppdf12) 881 deallocate (carppdf22) 882 deallocate (somppydf1) 883 deallocate (somppydf2) 884 deallocate (fp01) 885 deallocate (fm01) 886 end subroutine end_sub_2qtl
funct_0qtl
[ Top ] [ m_qtlmap_analyse_unitrait ] [ Subroutines ]
NAME
funct_0qtl
DESCRIPTION
Calcul de la vraisemblance et de ses derivees partielles sous H0
SOURCE
1025 subroutine funct_0qtl(n,x,f,iuser,user) 1026 use m_qtlmap_analyse_gen, only : somy,eff,cary,effdf,somydf,carydf,estmum,somcd,somcddf 1027 1028 implicit none 1029 integer , intent(in) :: n 1030 ! 1031 ! Tableaux dimensionnes selon n, le nombre de parametres a estimer 1032 ! 1033 real (kind=dp) ,dimension(n), intent(in) :: x 1034 real (kind=dp) , intent(inout) :: f 1035 integer , dimension(1), intent(in) :: iuser 1036 real (kind=dp) ,dimension(1), intent(in) :: user 1037 1038 integer :: ip,nm1,nm2,jm 1039 integer :: nestim,nest 1040 1041 real (kind=dp) :: sig,var,vpf,xmup,xmup2 1042 real (kind=dp) :: vdf,somxmu,xmum,xmum2,xmupm 1043 1044 !************nestim=sum(estmum(:ip-1))+nest****************************************************************** 1045 f=0.d0 1046 nestim=0 1047 do ip=1,np 1048 sig=x(ip) 1049 var=sig*sig 1050 xmup=x(ip+np) 1051 xmup2=xmup*xmup 1052 vdf=carydf(ip)+somcddf(ip)*xmup2-2.d0*xmup*somydf(ip) 1053 vdf=0.5d0*vdf/var 1054 fp0(ip)=vdf+(dble(effdf(ip))*dlog(sig)) 1055 f=f+fp0(ip) 1056 ! print *,f 1057 1058 nm1=nmp(ip)+1 1059 nm2=nmp(ip+1) 1060 somxmu=0.d0 1061 nest=0 1062 do jm=nm1,nm2 1063 if(estime(current_ic,jm)) then 1064 nest=nest+1 1065 if(nest.le.estmum(ip)) then 1066 nestim=nestim+1 1067 ! print *,jm,nest,nestim,sum(estmum(:ip-1))+nest,estmum(ip) 1068 xmum=x(2*np+nestim) 1069 somxmu=somxmu+xmum 1070 ! print *,somxmu 1071 else 1072 xmum=-somxmu 1073 end if 1074 xmum2=xmum*xmum 1075 xmupm=2.d0*(xmup+xmum) 1076 vpf=cary(jm)+somcd(jm)*(xmup2+xmum2+2.d0*xmup*xmum)-xmupm*somy(jm) 1077 vpf=0.5d0*vpf/var 1078 fm0(jm)=vpf+(dble(eff(jm))*dlog(sig)) 1079 f=f+fm0(jm) 1080 fp0(ip)=fp0(ip)+fm0(jm) 1081 ! print *,ip,jm,fm0(jm),f 1082 ! stop 1083 else 1084 fm0(jm)=0.d0 1085 end if 1086 end do 1087 end do 1088 1089 return 1090 end subroutine funct_0qtl
funct_0qtl_family
[ Top ] [ m_qtlmap_analyse_unitrait ] [ Subroutines ]
NAME
funct_0qtl_family
DESCRIPTION
Calcul de la vraisemblance et de ses derivees partielles sous H0
SOURCE
1099 subroutine funct_0qtl_family(ip,n,x,f,iuser,user) 1100 implicit none 1101 integer , intent(in) :: ip,n 1102 ! 1103 ! Tableaux dimensionnes selon n, le nombre de parametres a estimer 1104 ! 1105 real (kind=dp) ,dimension(n), intent(in) :: x 1106 real (kind=dp) , intent(inout) :: f 1107 integer , dimension(1), intent(in) :: iuser 1108 real (kind=dp) ,dimension(1), intent(in) :: user 1109 1110 integer :: jm,nestim,nest 1111 1112 real (kind=dp) :: sig,var,vpf,xmup,xmup2 1113 real (kind=dp) :: vdf,somxmu,xmum,xmum2,xmupm 1114 1115 ! nest=count(estime(current_ic,nmp(ip)+1:jm)) 1116 f=0.d0 1117 sig=x(ip) 1118 var=sig*sig 1119 xmup=x(ip+np) 1120 xmup2=xmup*xmup 1121 vdf=carydf(ip)+somcddf(ip)*xmup2-2.d0*xmup*somydf(ip) 1122 vdf=0.5d0*vdf/var 1123 f=vdf+(dble(effdf(ip))*dlog(sig)) 1124 somxmu=0.d0 1125 nest=0 1126 do jm=nmp(ip)+1,nmp(ip+1) 1127 if(estime(current_ic,jm)) then 1128 nest=nest+1 1129 if(nest.le.estmum(ip)) then 1130 if (ip>1) then 1131 nestim=sum(estmum(:ip-1))+nest 1132 else 1133 nestim=nest 1134 end if 1135 ! nestim=nestim+1 1136 ! print *,jm,nest,nestim,sum(estmum(:ip-1))+nest,estmum(ip) 1137 xmum=x(2*np+nestim) 1138 somxmu=somxmu+xmum 1139 ! print *,somxmu 1140 else 1141 xmum=-somxmu 1142 end if 1143 xmum2=xmum*xmum 1144 xmupm=2.d0*(xmup+xmum) 1145 vpf=cary(jm)+somcd(jm)*(xmup2+xmum2+2.d0*xmup*xmum)-xmupm*somy(jm) 1146 vpf=0.5d0*vpf/var 1147 fm0(jm)=vpf+(dble(eff(jm))*dlog(sig)) 1148 f=f+fm0(jm) 1149 else 1150 fm0(jm)=0.d0 1151 end if 1152 end do 1153 1154 end subroutine funct_0qtl_family
funct_1car_2qtl
[ Top ] [ m_qtlmap_analyse_unitrait ] [ Subroutines ]
NAME
funct_1car_2qtl
DESCRIPTION
Calcul de la vraisemblance et de ses derivees partielles sous H1
SOURCE
2067 subroutine funct_1car_2qtl(n,x,f,iuser,user) 2068 integer , intent(in) :: n 2069 ! 2070 ! Tableaux dimensionnes selon n, le nombre de parametres a estimer 2071 ! 2072 real (kind=dp) ,dimension(n), intent(in) :: x 2073 real (kind=dp) , intent(inout) :: f 2074 integer , dimension(1), intent(in) :: iuser 2075 real (kind=dp) ,dimension(1), intent(in) :: user 2076 2077 integer :: ip,nm1,nm2,jm,ngeno1,ngeno2,ig 2078 integer :: nestim,nest,ifem,indam,ngeno21,ngeno22,ig2 2079 2080 real (kind=dp) :: sig,var,vmere,vpf,xmup,xmup2,a1p,a2p,a1p2 2081 real (kind=dp) :: a2p2,vdf,somxmu,xmum,xmum2,xmupm,a1m,a2m 2082 real (kind=dp) :: a1m2,a2m2,z 2083 !**************************************************************************** 2084 f=0.d0 2085 fm01=0.d0 2086 fp01=0.d0 2087 2088 nestim=0 2089 do ip=1,np 2090 sig=x(ip) 2091 var=sig*sig 2092 xmup=x(ip+np) 2093 xmup2=xmup*xmup 2094 a1p=x(2*np+ip) 2095 a2p=x(3*np+ip) 2096 a1p2=a1p*a1p 2097 a2p2=a2p*a2p 2098 vdf=carydf(ip)+dble(effdf(ip))*xmup2-2.d0*xmup*somydf(ip) 2099 vdf=vdf+a1p2*carppdf11(ip)+a2p2*carppdf22(ip) & 2100 -2.d0*a1p*somppydf1(ip)-2.d0*a2p*somppydf2(ip) & 2101 +2.d0*xmup*(a1p*somppdf1(ip)+a2p*somppdf2(ip)) & 2102 +2.d0*a1p*a2p*carppdf12(ip) 2103 vdf=0.5d0*vdf/var 2104 fp01(ip)=vdf+dble(effdf(ip))*dlog(sig) 2105 f=f+fp01(ip) 2106 2107 nm1=nmp(ip)+1 2108 nm2=nmp(ip+1) 2109 somxmu=0.d0 2110 nest=0 2111 do jm=nm1,nm2 2112 if(estime(current_ic,jm)) then 2113 nest=nest+1 2114 if(nest.le.estmum(ip)) then 2115 nestim=nestim+1 2116 xmum=x(4*np+nestim) 2117 somxmu=somxmu+xmum 2118 else 2119 xmum=-somxmu 2120 end if 2121 xmum2=xmum*xmum 2122 xmupm=2.d0*(xmup+xmum) 2123 ifem=repfem(jm) 2124 indam=iam(current_ic,ifem) 2125 a1m=x(4*np+nmumest(current_ic)+indam) 2126 a2m=x(4*np+nmumest(current_ic)+namest(current_ic)+indam) 2127 a1m2=a1m*a1m 2128 a2m2=a2m*a2m 2129 z=cary(jm)+dble(eff(jm))*(xmup2+xmum2+2.d0*xmup*xmum)-xmupm*somy(jm) 2130 vmere=0.d0 2131 ngeno1=ngenom(current_chr,jm)+1 2132 ngeno2=ngenom(current_chr,jm+1) 2133 ngeno21=ngenom(current_chr2,jm)+1 2134 ngeno22=ngenom(current_chr2,jm+1) 2135 2136 if( (ngeno2-ngeno1 == 0) .and. (ngeno22-ngeno21 == 0) )then 2137 ig=ngeno1 2138 ig2=ngeno21 2139 vpf=z +a1p2*carpp1(ig)+a2p2*carpp2(ig)+a1m2*carpm1(ig & 2140 ) +a2m2*carpm2(ig) +2.d0* (a1p*a2p*sompp1p2(ig,ig2)+ & 2141 a1p *a1m*sompp1m1(ig)+a1p*a2m*sompp1m2(ig,ig2) +a2p & 2142 *a1m *sompp2m1(ig,ig2)+a2p*a2m*sompp2m2(ig)+a1m*a2m & 2143 *sompm1m2(ig,ig2)) +xmupm* (a1p*sompp1(ig)+a2p & 2144 *sompp2(ig)+a1m*sompm1(ig)+a2m*sompm2(ig)) -2.d0 & 2145 *(a1p*sompp1y(ig)+ a2p*sompp2y(ig)+ a1m & 2146 *sompm1y(ig)+a2m*sompm2y(ig)) 2147 !print *,"vpf:",vpf 2148 vpf=0.5d0*vpf/var 2149 2150 fm01(jm)=vpf+dble(eff(jm))*dlog(sig) 2151 !print *,"fm01:",fm01(jm) 2152 else 2153 do ig=ngeno1,ngeno2 2154 do ig2=ngeno21,ngeno22 2155 vpf=z +a1p2*carpp1(ig)+a2p2*carpp2(ig)+a1m2*carpm1(ig & 2156 ) +a2m2*carpm2(ig) +2.d0* (a1p*a2p*sompp1p2(ig,ig2)+ & 2157 a1p *a1m*sompp1m1(ig)+a1p*a2m*sompp1m2(ig,ig2) +a2p & 2158 *a1m *sompp2m1(ig,ig2)+a2p*a2m*sompp2m2(ig)+a1m*a2m & 2159 *sompm1m2(ig,ig2)) +xmupm* (a1p*sompp1(ig)+a2p & 2160 *sompp2(ig)+a1m*sompm1(ig)+a2m*sompm2(ig)) -2.d0 & 2161 *(a1p*sompp1y(ig)+ a2p*sompp2y(ig)+ a1m & 2162 *sompm1y(ig)+a2m*sompm2y(ig)) 2163 vpf=-0.5d0*vpf/var 2164 ! note ofi : a voir avec pascale...il y a 2 proba de genotype sur 2 chromosomes.... 2165 if ( current_chr == current_chr2 ) then 2166 vmere=vmere+probg(current_chr,ig)*dexp(vpf) 2167 else 2168 vmere=vmere+probg(current_chr,ig)*probg(current_chr2,ig2)*dexp(vpf) 2169 end if 2170 end do 2171 end do 2172 fm01(jm)=-dlog(vmere)+(dble(eff(jm))*dlog(sig)) 2173 end if 2174 2175 2176 f=f+fm01(jm) 2177 fp01(ip)=fp01(ip)+fm01(jm) 2178 else 2179 fm01(jm)=0.d0 2180 end if 2181 end do 2182 end do 2183 2184 return 2185 end subroutine funct_1car_2qtl
funct_1car_2qtl_family
[ Top ] [ m_qtlmap_analyse_unitrait ] [ Subroutines ]
NAME
funct_1car_2qtl_family
DESCRIPTION
Calcul de la vraisemblance et de ses derivees partielles sous H1
SOURCE
2196 subroutine funct_1car_2qtl_family(ip,n,x,f,iuser,user) 2197 integer , intent(in) :: ip,n 2198 ! 2199 ! Tableaux dimensionnes selon n, le nombre de parametres a estimer 2200 ! 2201 real (kind=dp) ,dimension(n), intent(in) :: x 2202 real (kind=dp) , intent(inout) :: f 2203 integer , dimension(1), intent(in) :: iuser 2204 real (kind=dp) ,dimension(1), intent(in) :: user 2205 2206 integer :: jm,ngeno1,ngeno2,ig 2207 integer :: nestim,nest,ifem,indam,ngeno21,ngeno22,ig2 2208 2209 real (kind=dp) :: sig,var,vmere,vpf,xmup,xmup2,a1p,a2p,a1p2 2210 real (kind=dp) :: a2p2,vdf,somxmu,xmum,xmum2,xmupm,a1m,a2m 2211 real (kind=dp) :: a1m2,a2m2,z 2212 !**************************************************************************** 2213 f=0.d0 2214 fm01=0.d0 2215 sig=x(ip) 2216 var=sig*sig 2217 xmup=x(ip+np) 2218 xmup2=xmup*xmup 2219 a1p=x(2*np+ip) 2220 a2p=x(3*np+ip) 2221 a1p2=a1p*a1p 2222 a2p2=a2p*a2p 2223 vdf=carydf(ip)+dble(effdf(ip))*xmup2-2.d0*xmup*somydf(ip) 2224 vdf=vdf+a1p2*carppdf11(ip)+a2p2*carppdf22(ip) & 2225 -2.d0*a1p*somppydf1(ip)-2.d0*a2p*somppydf2(ip) & 2226 +2.d0*xmup*(a1p*somppdf1(ip)+a2p*somppdf2(ip)) & 2227 +2.d0*a1p*a2p*carppdf12(ip) 2228 vdf=0.5d0*vdf/var 2229 f=vdf+dble(effdf(ip))*dlog(sig) 2230 somxmu=0.d0 2231 nest=0 2232 do jm=nmp(ip)+1,nmp(ip+1) 2233 if(estime(current_ic,jm)) then 2234 nest=nest+1 2235 if(nest.le.estmum(ip)) then 2236 if (ip>1) then 2237 nestim=sum(estmum(:ip-1))+nest 2238 else 2239 nestim=nest 2240 end if 2241 xmum=x(4*np+nestim) 2242 somxmu=somxmu+xmum 2243 else 2244 xmum=-somxmu 2245 end if 2246 xmum2=xmum*xmum 2247 xmupm=2.d0*(xmup+xmum) 2248 ifem=repfem(jm) 2249 indam=iam(current_ic,ifem) 2250 a1m=x(4*np+nmumest(current_ic)+indam) 2251 a2m=x(4*np+nmumest(current_ic)+namest(current_ic)+indam) 2252 a1m2=a1m*a1m 2253 a2m2=a2m*a2m 2254 z=cary(jm)+dble(eff(jm))*(xmup2+xmum2+2.d0*xmup*xmum)-xmupm*somy(jm) 2255 vmere=0.d0 2256 ngeno1=ngenom(current_chr,jm)+1 2257 ngeno2=ngenom(current_chr,jm+1) 2258 ngeno21=ngenom(current_chr2,jm)+1 2259 ngeno22=ngenom(current_chr2,jm+1) 2260 2261 if( (ngeno2-ngeno1 == 0) .and. (ngeno22-ngeno21 == 0) )then 2262 ig=ngeno1 2263 ig2=ngeno21 2264 vpf=z +a1p2*carpp1(ig)+a2p2*carpp2(ig)+a1m2*carpm1(ig & 2265 ) +a2m2*carpm2(ig) +2.d0* (a1p*a2p*sompp1p2(ig,ig2)+ & 2266 a1p *a1m*sompp1m1(ig)+a1p*a2m*sompp1m2(ig,ig2) +a2p & 2267 *a1m *sompp2m1(ig,ig2)+a2p*a2m*sompp2m2(ig)+a1m*a2m & 2268 *sompm1m2(ig,ig2)) +xmupm* (a1p*sompp1(ig)+a2p & 2269 *sompp2(ig)+a1m*sompm1(ig)+a2m*sompm2(ig)) -2.d0 & 2270 *(a1p*sompp1y(ig)+ a2p*sompp2y(ig)+ a1m & 2271 *sompm1y(ig)+a2m*sompm2y(ig)) 2272 !print *,"vpf:",vpf 2273 vpf=0.5d0*vpf/var 2274 2275 fm01(jm)=vpf+dble(eff(jm))*dlog(sig) 2276 !print *,"fm01:",fm01(jm) 2277 else 2278 do ig=ngeno1,ngeno2 2279 do ig2=ngeno21,ngeno22 2280 vpf=z +a1p2*carpp1(ig)+a2p2*carpp2(ig)+a1m2*carpm1(ig & 2281 ) +a2m2*carpm2(ig) +2.d0* (a1p*a2p*sompp1p2(ig,ig2)+ & 2282 a1p *a1m*sompp1m1(ig)+a1p*a2m*sompp1m2(ig,ig2) +a2p & 2283 *a1m *sompp2m1(ig,ig2)+a2p*a2m*sompp2m2(ig)+a1m*a2m & 2284 *sompm1m2(ig,ig2)) +xmupm* (a1p*sompp1(ig)+a2p & 2285 *sompp2(ig)+a1m*sompm1(ig)+a2m*sompm2(ig)) -2.d0 & 2286 *(a1p*sompp1y(ig)+ a2p*sompp2y(ig)+ a1m & 2287 *sompm1y(ig)+a2m*sompm2y(ig)) 2288 vpf=-0.5d0*vpf/var 2289 ! note ofi : a voir avec pascale...il y a 2 proba de genotype sur 2 chromosomes.... 2290 if ( current_chr == current_chr2 ) then 2291 vmere=vmere+probg(current_chr,ig)*dexp(vpf) 2292 else 2293 vmere=vmere+probg(current_chr,ig)*probg(current_chr2,ig2)*dexp(vpf) 2294 end if 2295 end do 2296 end do 2297 fm01(jm)=-dlog(vmere)+(dble(eff(jm))*dlog(sig)) 2298 end if 2299 f=f+fm01(jm) 2300 else 2301 fm01(jm)=0.d0 2302 end if 2303 end do 2304 2305 end subroutine funct_1car_2qtl_family
funct_1qtl
[ Top ] [ m_qtlmap_analyse_unitrait ] [ Subroutines ]
NAME
funct_1qtl
DESCRIPTION
Calcul de la vraisemblance et de ses derivees partielles sous Hypothese 1QTL 1carac
SOURCE
1460 subroutine funct_1qtl(n,x,f,iuser,user) 1461 1462 integer , intent(in) :: n ! number of variables 1463 double precision, dimension(n), intent(in) :: x ! Position 1464 double precision, intent(inout) :: f ! Result value at the position 1465 double precision,dimension(1) :: user 1466 integer ,dimension(1) :: iuser 1467 1468 1469 ! 1470 ! Divers 1471 integer :: ip,nestim,nm1,nm2,nest,jm,ifem,indam,i 1472 integer :: ngeno1,ngeno2,ig 1473 real (kind=dp) :: sig,var,xmup,xmup2,vdf,somxmu,xmum,xmum2 1474 real (kind=dp) :: xmupm,vmere,z,vpf,x_ap,x_ap2,x_am,x_am2 1475 1476 1477 !****************************************************************************** 1478 f=0.d0 1479 nestim=0 1480 do ip=1,np 1481 sig=x(ip) 1482 var=sig*sig 1483 xmup=x(ip+np) 1484 xmup2=xmup*xmup 1485 x_ap=x(2*np+ip) 1486 x_ap2=x_ap*x_ap 1487 vdf=carydf(ip)+somcddf(ip)*xmup2-2.d0*xmup*somydf(ip) 1488 vdf=vdf+x_ap2*carppdf(ip)+2.d0*xmup*x_ap*somppdf(ip)-2.d0*x_ap*somppydf(ip) 1489 vdf=0.5d0*vdf/var 1490 fp1(ip)=vdf+dble(effdf(ip))*dlog(sig) 1491 f=f+fp1(ip) 1492 nm1=nmp(ip)+1 1493 nm2=nmp(ip+1) 1494 somxmu=0.d0 1495 nest=0 1496 do jm=nm1,nm2 1497 if(estime(current_ic,jm)) then 1498 nest=nest+1 1499 if(nest.le.estmum(ip)) then 1500 nestim=nestim+1 1501 xmum=x(3*np+nestim) 1502 somxmu=somxmu+xmum 1503 else 1504 xmum=-somxmu 1505 end if 1506 xmum2=xmum*xmum 1507 xmupm=2.d0*(xmup+xmum) 1508 ifem=repfem(jm) 1509 indam=iam(current_ic,ifem) 1510 x_am=x(3*np+nmumest(current_ic)+indam) 1511 x_am2=x_am*x_am 1512 z=cary(jm)+somcd(jm)*(xmup2+xmum2+2.d0*xmup*xmum)-xmupm*somy(jm) 1513 vmere=0.d0 1514 ngeno1=ngenom(current_chr,jm)+1 1515 ngeno2=ngenom(current_chr,jm+1) 1516 1517 do ig=ngeno1,ngeno2 1518 vpf=z+x_ap2*carpp(ig)+x_am2*carpm(ig) & 1519 +2.d0*x_ap*x_am*somppm(ig) & 1520 +xmupm*(x_ap*sompp(ig)+x_am*sompm(ig)) & 1521 -2.d0*x_ap*somppy(ig)-2.d0*x_am*sompmy(ig) 1522 vpf=-0.5d0*vpf/var 1523 vmere=vmere+probg(current_chr,ig)*dexp(vpf) 1524 1525 end do 1526 1527 if (vmere == 0) then 1528 fm1(jm)=INIFINY_REAL_VALUE 1529 else 1530 fm1(jm)=-dlog(vmere)+(dble(eff(jm))*dlog(sig)) 1531 end if 1532 1533 ! end if 1534 f=f+fm1(jm) 1535 fp1(ip)=fp1(ip)+fm1(jm) 1536 else 1537 fm1(jm)=0.d0 1538 end if 1539 end do 1540 end do 1541 return 1542 end subroutine funct_1qtl
funct_1qtl_family
[ Top ] [ m_qtlmap_analyse_unitrait ] [ Subroutines ]
NAME
funct_1qtl_family
DESCRIPTION
Calcul de la vraisemblance et de ses derivees partielles sous Hypothese 1QTL 1carac
SOURCE
1553 subroutine funct_1qtl_family(ip,n,x,f,iuser,user) 1554 1555 integer , intent(in) :: ip,n ! number of variables 1556 double precision, dimension(n), intent(in) :: x ! Position 1557 double precision, intent(inout) :: f ! Result value at the position 1558 double precision,dimension(1) :: user 1559 integer ,dimension(1) :: iuser 1560 1561 1562 ! 1563 ! Divers 1564 integer :: nestim,nest,jm,ifem,indam,i 1565 integer :: ngeno1,ngeno2,ig 1566 real (kind=dp) :: sig,var,xmup,xmup2,vdf,somxmu,xmum,xmum2 1567 real (kind=dp) :: xmupm,vmere,z,vpf,x_ap,x_ap2,x_am,x_am2 1568 1569 1570 !****************************************************************************** 1571 f=0.d0 1572 nestim=0 1573 sig=x(ip) 1574 var=sig*sig 1575 xmup=x(ip+np) 1576 xmup2=xmup*xmup 1577 x_ap=x(2*np+ip) 1578 x_ap2=x_ap*x_ap 1579 vdf=carydf(ip)+somcddf(ip)*xmup2-2.d0*xmup*somydf(ip) 1580 vdf=vdf+x_ap2*carppdf(ip)+2.d0*xmup*x_ap*somppdf(ip)-2.d0*x_ap*somppydf(ip) 1581 vdf=0.5d0*vdf/var 1582 f=vdf+dble(effdf(ip))*dlog(sig) 1583 somxmu=0.d0 1584 nest=0 1585 do jm=nmp(ip)+1,nmp(ip+1) 1586 if(estime(current_ic,jm)) then 1587 nest=nest+1 1588 if(nest.le.estmum(ip)) then 1589 if (ip>1) then 1590 nestim=sum(estmum(:ip-1))+nest 1591 else 1592 nestim=nest 1593 end if 1594 xmum=x(3*np+nestim) 1595 somxmu=somxmu+xmum 1596 else 1597 xmum=-somxmu 1598 end if 1599 xmum2=xmum*xmum 1600 xmupm=2.d0*(xmup+xmum) 1601 ifem=repfem(jm) 1602 indam=iam(current_ic,ifem) 1603 x_am=x(3*np+nmumest(current_ic)+indam) 1604 x_am2=x_am*x_am 1605 z=cary(jm)+somcd(jm)*(xmup2+xmum2+2.d0*xmup*xmum)-xmupm*somy(jm) 1606 vmere=0.d0 1607 ngeno1=ngenom(current_chr,jm)+1 1608 ngeno2=ngenom(current_chr,jm+1) 1609 1610 do ig=ngeno1,ngeno2 1611 vpf=z+x_ap2*carpp(ig)+x_am2*carpm(ig) & 1612 +2.d0*x_ap*x_am*somppm(ig) & 1613 +xmupm*(x_ap*sompp(ig)+x_am*sompm(ig)) & 1614 -2.d0*x_ap*somppy(ig)-2.d0*x_am*sompmy(ig) 1615 vpf=-0.5d0*vpf/var 1616 vmere=vmere+probg(current_chr,ig)*dexp(vpf) 1617 end do 1618 1619 if (vmere == 0) then 1620 fm1(jm)=INIFINY_REAL_VALUE 1621 else 1622 fm1(jm)=-dlog(vmere)+(dble(eff(jm))*dlog(sig)) 1623 end if 1624 f=f+fm1(jm) 1625 else 1626 fm1(jm)=0.d0 1627 end if 1628 end do 1629 1630 end subroutine funct_1qtl_family
init_analyse_unitrait_1QTL
[ Top ] [ m_qtlmap_analyse_unitrait ] [ Subroutines ]
NAME
init_analyse_unitrait_1QTL
DESCRIPTION
SOURCE
654 subroutine init_analyse_unitrait_1QTL 655 integer :: stat 656 allocate (sig1(np),STAT=stat) 657 call check_allocate(stat,'sig1 [m_qtlmap_analyse_unitrait]') 658 sig1=0.D0 659 allocate (xmu1p(np),STAT=stat) 660 call check_allocate(stat,'xmu1p [m_qtlmap_analyse_unitrait]') 661 xmu1p=0.d0 662 allocate (xmu1m(nm),STAT=stat) 663 call check_allocate(stat,'xmu1m [m_qtlmap_analyse_unitrait]') 664 xmu1m=0.d0 665 allocate (fp0(np),STAT=stat) 666 call check_allocate(stat,'fp0 [m_qtlmap_analyse_unitrait]') 667 668 allocate (fm0(nm),STAT=stat) 669 call check_allocate(stat,'fm0 [m_qtlmap_analyse_unitrait]') 670 671 allocate (ap(np),STAT=stat) 672 call check_allocate(stat,'ap [m_qtlmap_analyse_unitrait]') 673 ap=0.d0 674 allocate (xmoyp(np),STAT=stat) 675 call check_allocate(stat,'xmoyp [m_qtlmap_analyse_unitrait]') 676 xmoyp=0.d0 677 allocate (std(np),STAT=stat) 678 call check_allocate(stat,'std [m_qtlmap_analyse_unitrait]') 679 std=0.d0 680 allocate (am(nm),STAT=stat) 681 call check_allocate(stat,'am [m_qtlmap_analyse_unitrait]') 682 am=0.d0 683 allocate (xmoym(nm),STAT=stat) 684 call check_allocate(stat,'xmoym [m_qtlmap_analyse_unitrait]') 685 xmoym=0.d0 686 allocate (fp_1(np),STAT=stat) 687 call check_allocate(stat,'fp_1 [m_qtlmap_analyse_unitrait]') 688 allocate (fm_1(nm),STAT=stat) 689 call check_allocate(stat,'fm_1 [m_qtlmap_analyse_unitrait]') 690 691 end subroutine init_analyse_unitrait_1QTL
init_analyse_unitrait_2QTL
[ Top ] [ m_qtlmap_analyse_unitrait ] [ Subroutines ]
NAME
init_analyse_unitrait_2QTL
DESCRIPTION
SOURCE
760 subroutine init_analyse_unitrait_2QTL 761 integer :: stat,maxng 762 763 764 allocate (ap2(np,2),STAT=stat) 765 call check_allocate(stat,'ap2 [m_qtlmap_analyse_unitrait_2qtl]') 766 allocate (xmoyp2(np),STAT=stat) 767 call check_allocate(stat,'xmoyp2 [m_qtlmap_analyse_unitrait_2qtl]') 768 allocate (std2(np),STAT=stat) 769 call check_allocate(stat,'std2 [m_qtlmap_analyse_unitrait_2qtl]') 770 allocate (fp_2(np),STAT=stat) 771 call check_allocate(stat,'fp_2 [m_qtlmap_analyse_unitrait_2qtl]') 772 allocate (am2(nm,2),STAT=stat) 773 call check_allocate(stat,'am2 [m_qtlmap_analyse_unitrait_2qtl]') 774 allocate (xmoym2(nm),STAT=stat) 775 call check_allocate(stat,'xmoym2 [m_qtlmap_analyse_unitrait_2qtl]') 776 allocate (fm_2(nm),STAT=stat) 777 call check_allocate(stat,'fm_2 [m_qtlmap_analyse_unitrait_2qtl]') 778 779 end subroutine init_analyse_unitrait_2QTL
init_sub_1qtl
[ Top ] [ m_qtlmap_analyse_unitrait ] [ Subroutines ]
NAME
init_sub_1qtl
DESCRIPTION
SOURCE
700 subroutine init_sub_1qtl 701 integer :: stat,maxng 702 703 maxng = get_maxnbgenotypedam() 704 allocate (sompp(maxng),STAT=stat) 705 call check_allocate(stat,'sompp [m_qtlmap_analyse_unitrait]') 706 allocate (sompm(maxng),STAT=stat) 707 call check_allocate(stat,'sompm [m_qtlmap_analyse_unitrait]') 708 allocate (sompmy(maxng),STAT=stat) 709 call check_allocate(stat,'sompmy [m_qtlmap_analyse_unitrait]') 710 allocate (carpp(maxng),STAT=stat) 711 call check_allocate(stat,'carpp [m_qtlmap_analyse_unitrait]') 712 allocate (carpm(maxng),STAT=stat) 713 call check_allocate(stat,'carpm [m_qtlmap_analyse_unitrait]') 714 allocate (somppy(maxng),STAT=stat) 715 call check_allocate(stat,'somppy [m_qtlmap_analyse_unitrait]') 716 allocate (somppm(maxng),STAT=stat) 717 call check_allocate(stat,'somppm [m_qtlmap_analyse_unitrait]') 718 allocate (somppdf(np),STAT=stat) 719 call check_allocate(stat,'somppdf [m_qtlmap_analyse_unitrait]') 720 allocate (carppdf(np),STAT=stat) 721 call check_allocate(stat,'carppdf [m_qtlmap_analyse_unitrait]') 722 allocate (somppydf(np),STAT=stat) 723 call check_allocate(stat,'somppydf [m_qtlmap_analyse_unitrait]') 724 allocate (fp1(np),STAT=stat) 725 call check_allocate(stat,'fp1 [m_qtlmap_analyse_unitrait]') 726 allocate (fm1(nm),STAT=stat) 727 call check_allocate(stat,'fm1 [m_qtlmap_analyse_unitrait]') 728 729 end subroutine init_sub_1qtl
init_sub_2qtl
[ Top ] [ m_qtlmap_analyse_unitrait ] [ Subroutines ]
NAME
init_sub_2qtl
DESCRIPTION
SOURCE
788 subroutine init_sub_2qtl 789 integer :: stat,maxng 790 791 maxng = get_maxnbgenotypedam() 792 793 allocate (sompp1(maxng),STAT=stat) 794 call check_allocate(stat,'sompp1 [m_qtlmap_analyse_unitrait_2qtl]') 795 allocate (sompp2(maxng),STAT=stat) 796 call check_allocate(stat,'sompp2 [m_qtlmap_analyse_unitrait_2qtl]') 797 allocate (sompm1(maxng),STAT=stat) 798 call check_allocate(stat,'sompm1 [m_qtlmap_analyse_unitrait_2qtl]') 799 allocate (sompm2(maxng),STAT=stat) 800 call check_allocate(stat,'sompm2 [m_qtlmap_analyse_unitrait_2qtl]') 801 allocate (carpp1(maxng),STAT=stat) 802 call check_allocate(stat,'carpp1 [m_qtlmap_analyse_unitrait_2qtl]') 803 allocate (carpp2(maxng),STAT=stat) 804 call check_allocate(stat,'carpp2 [m_qtlmap_analyse_unitrait_2qtl]') 805 allocate (carpm1(maxng),STAT=stat) 806 call check_allocate(stat,'carpm1 [m_qtlmap_analyse_unitrait_2qtl]') 807 allocate (carpm2(maxng),STAT=stat) 808 call check_allocate(stat,'carpm2 [m_qtlmap_analyse_unitrait_2qtl]') 809 allocate (sompp1y(maxng),STAT=stat) 810 call check_allocate(stat,'sompp1y [m_qtlmap_analyse_unitrait_2qtl]') 811 allocate (sompp2y(maxng),STAT=stat) 812 call check_allocate(stat,'sompp2y [m_qtlmap_analyse_unitrait_2qtl]') 813 allocate (sompm1y(maxng),STAT=stat) 814 call check_allocate(stat,'sompm1y [m_qtlmap_analyse_unitrait_2qtl]') 815 allocate (sompm2y(maxng),STAT=stat) 816 call check_allocate(stat,'sompm2y [m_qtlmap_analyse_unitrait_2qtl]') 817 allocate (sompp1p2(maxng,maxng),STAT=stat) 818 call check_allocate(stat,'sompp1p2 [m_qtlmap_analyse_unitrait_2qtl]') 819 allocate (sompp1m1(maxng),STAT=stat) 820 call check_allocate(stat,'sompp1m1 [m_qtlmap_analyse_unitrait_2qtl]') 821 allocate (sompp1m2(maxng,maxng),STAT=stat) 822 call check_allocate(stat,'sompp1m2 [m_qtlmap_analyse_unitrait_2qtl]') 823 allocate (sompp2m1(maxng,maxng),STAT=stat) 824 call check_allocate(stat,'sompp2m1 [m_qtlmap_analyse_unitrait_2qtl]') 825 allocate (sompp2m2(maxng),STAT=stat) 826 call check_allocate(stat,'sompp2m2 [m_qtlmap_analyse_unitrait_2qtl]') 827 allocate (sompm1m2(maxng,maxng),STAT=stat) 828 call check_allocate(stat,'sompm1m2 [m_qtlmap_analyse_unitrait_2qtl]') 829 allocate (somppdf1(np),STAT=stat) 830 call check_allocate(stat,'somppdf1 [m_qtlmap_analyse_unitrait_2qtl]') 831 allocate (somppdf2(np),STAT=stat) 832 call check_allocate(stat,'somppdf2 [m_qtlmap_analyse_unitrait_2qtl]') 833 allocate (carppdf11(np),STAT=stat) 834 call check_allocate(stat,'carppdf11 [m_qtlmap_analyse_unitrait_2qtl]') 835 allocate (carppdf12(np),STAT=stat) 836 call check_allocate(stat,'carppdf12 [m_qtlmap_analyse_unitrait_2qtl]') 837 allocate (carppdf22(np),STAT=stat) 838 call check_allocate(stat,'carppdf22 [m_qtlmap_analyse_unitrait_2qtl]') 839 allocate (somppydf1(np),STAT=stat) 840 call check_allocate(stat,'somppydf1 [m_qtlmap_analyse_unitrait_2qtl]') 841 allocate (somppydf2(np),STAT=stat) 842 call check_allocate(stat,'somppydf2 [m_qtlmap_analyse_unitrait_2qtl]') 843 844 allocate (fp01(np),STAT=stat) 845 call check_allocate(stat,'fp01 [m_qtlmap_analyse_unitrait_2qtl]') 846 allocate (fm01(nm),STAT=stat) 847 call check_allocate(stat,'fm01 [m_qtlmap_analyse_unitrait_2qtl]') 848 849 end subroutine init_sub_2qtl
opti_0qtl
[ Top ] [ m_qtlmap_analyse_unitrait ] [ Subroutines ]
NAME
opti_0qtl
DESCRIPTION
Calcul de la vraisemblance 0 QTL, 1 caractere
SOURCE
936 subroutine opti_0qtl(ic) 937 use m_qtlmap_analyse_gen, only : nmumest,sig0,xmu0p,xmu0m 938 ! 939 ! Divers 940 941 integer , intent(in) :: ic 942 integer iuser(1) 943 double precision ,dimension(:),allocatable :: borni,borns,par 944 real (kind=dp) :: user(1) 945 integer :: npar,ibound,ip,jm,ifail,nest,nestim 946 logical , dimension(:,:),pointer :: filter_inc 947 948 ! 949 current_ic = ic 950 !****************************************************************************** 951 ! Parametres de maximisation 952 npar=(2*np)+nmumest(ic) 953 954 allocate (borni(npar)) 955 allocate (borns(npar)) 956 allocate (par(npar)) 957 allocate (filter_inc(np,npar)) 958 959 ibound=0 960 filter_inc=.false. 961 nestim=0 962 do ip=1,np 963 borni(ip)=SIG_MIN 964 borns(ip)=SIG_MAX 965 filter_inc(ip,ip)=.true. 966 borni(ip+np)=XMU_MIN 967 borns(ip+np)=XMU_MAX 968 filter_inc(ip,ip+np)=.true. 969 970 nest=0 971 do jm=nmp(ip)+1,nmp(ip+1) 972 if(estime(current_ic,jm)) then 973 nest=nest+1 974 if(nest.le.estmum(ip)) then 975 nestim=nestim+1 976 filter_inc(ip,2*np+nestim)=.true. 977 end if 978 end if 979 end do 980 end do 981 982 do jm=1,nmumest(ic) 983 borni(2*np+jm)=XMU_MIN 984 borns(2*np+jm)=XMU_MAX 985 end do 986 987 ! 988 ! Point de depart 989 do ip=1,np 990 par(ip)=sig0(ip) 991 par(ip+np)=xmu0p(ip) 992 end do 993 do jm=1,nmumest(ic) 994 par(2*np+jm)=xmu0m(jm) 995 end do 996 ! 997 ! Optimisation de la vraisemblance 998 ifail=1 999 ! call minimizing_funct(npar,ibound,funct_0qtl,borni,borns,par,f0,iuser,user,ifail) 1000 call minimizing_funct_family_sire(npar,ibound,funct_0qtl_family,filter_inc,fp0,borni,borns,par,f0,iuser,user,ifail) 1001 1002 do ip = 1,np 1003 sig1(ip)=par(ip) 1004 xmu1p(ip)=par(ip+np) 1005 end do 1006 do jm=1,nmumest(ic) 1007 xmu1m(jm)=par(2*np+jm) 1008 end do 1009 ! 1010 deallocate (borni) 1011 deallocate (borns) 1012 deallocate (par) 1013 deallocate (filter_inc) 1014 1015 return 1016 end subroutine opti_0qtl
opti_1car_2qtl
[ Top ] [ m_qtlmap_analyse_unitrait ] [ Subroutines ]
NAME
opti_1car_2qtl
DESCRIPTION
Computing statistical test across the chromosome
SOURCE
1642 subroutine opti_1car_2qtl(ic,lrtsol) 1643 implicit none 1644 1645 integer , intent(in) :: ic 1646 type(TYPE_LRT_SOLUTION) , intent(inout) :: lrtsol 1647 1648 ! 1649 ! Divers 1650 real (kind=dp) :: par(4*np+3*nm),borni(4*np+3*nm) 1651 real (kind=dp) :: borns(4*np+3*nm) 1652 integer :: iuser(1),chr 1653 real (kind=dp) :: user(1) 1654 real (kind=dp) :: pddp(nd*4,2),pddm(nd*4,2) 1655 1656 integer :: ip,ibound,ngeno1,ngeno2,ig,nd1,nd2,kd,kkd,ilong,ntotal,nprime,indexchr(nchr) 1657 integer :: npar,jm,nestim,nm1,nm2,nest,ifem,indam,n,ifail,kd2,nstart 1658 integer :: n1,ix1,ii,iq,npar2,chr2,ngeno21,ngeno22,ilong2,ig1,ig2,nd11,nd12,nd21,nd22 1659 real (kind=dp) :: xmin,xmax,somxmu,f01_t,save_f_1,save_f0 1660 logical , dimension(:,:),pointer :: filter_inc 1661 real (kind=dp) , dimension(:,:,:,:,:),pointer :: maxpar 1662 real (kind=dp) , dimension(:,:,:,:,:),pointer :: fpsave1 1663 real (kind=dp) , dimension(:,:,:,:,:),pointer :: fmsave1 1664 1665 !Trick pour recuperer la valeur de ces tableaux 1666 !La parallelisation sur les caracteres empeche la recopie de ces tableaux (deja instancier NTHREADS fois) 1667 !on recupere le pointeur de ces tableau avant la parallelisation sur le groupe de liaison 1668 !puis on initialise tous les thread avec ces valeurs sauvegarde 1669 real (kind=dp) ,dimension(:),pointer :: save_carydf 1670 real (kind=dp) ,dimension(:),pointer :: save_somcddf 1671 real (kind=dp) ,dimension(:),pointer :: save_somydf 1672 integer ,dimension(:),pointer :: save_effdf 1673 integer ,dimension(:),pointer :: save_estmum 1674 real (kind=dp) ,dimension(:),pointer :: save_somcd 1675 real (kind=dp) ,dimension(:),pointer :: save_somy 1676 real (kind=dp) ,dimension(:),pointer :: save_cary 1677 integer ,dimension(:),pointer :: save_eff 1678 real (kind=dp) ,dimension(:),pointer :: save_fp_1 1679 real (kind=dp) ,dimension(:),pointer :: save_fm_1 1680 1681 ! 1682 !**************************************************************************** 1683 ! 1684 call new(2,lrtsol) 1685 npar2=4*np+3*nm 1686 1687 current_ic = ic 1688 1689 borni=0.d0 1690 borns=0.d0 1691 par=0.d0 1692 lrtsol%lrt0_2=0.d0 1693 lrtsol%lrt1_2=0.d0 1694 1695 save_carydf => carydf 1696 save_somcddf => somcddf 1697 save_somydf => somydf 1698 save_effdf => effdf 1699 save_estmum => estmum 1700 save_somcd => somcd 1701 save_somy => somy 1702 save_cary => cary 1703 save_eff => eff 1704 save_fp_1 => fp_1 1705 save_fm_1 => fm_1 1706 save_f_1 = f_1 1707 save_f0 = f0 1708 ! 1709 ! Calcul de la vraisemblance sous H1 1710 ! Parametres de maximisation 1711 npar=(4*np)+nmumest(ic)+(2*namest(ic)) 1712 ibound=0 1713 allocate (filter_inc(np,npar)) 1714 1715 filter_inc=.false. 1716 nestim=0 1717 do ip=1,np 1718 xmin=-5.d0*std(ip) 1719 xmax=-xmin 1720 borni(ip)=SIG_MIN 1721 borns(ip)=SIG_MAX 1722 filter_inc(ip,ip)=.true. 1723 borni(ip+np)=XMU_MIN 1724 borns(ip+np)=XMU_MAX 1725 borni(2*np+ip)=xmin 1726 borns(2*np+ip)=xmax 1727 filter_inc(ip,2*np+ip)=.true. 1728 borni(3*np+ip)=xmin 1729 borns(3*np+ip)=xmax 1730 filter_inc(ip,3*np+ip)=.true. 1731 1732 nest=0 1733 do jm=nmp(ip)+1,nmp(ip+1) 1734 if(estime(current_ic,jm)) then 1735 nest=nest+1 1736 if(nest.le.estmum(ip)) then 1737 nestim=nestim+1 1738 filter_inc(ip,4*np+nestim)=.true. 1739 end if 1740 ifem=repfem(jm) 1741 indam=iam(ic,ifem) 1742 filter_inc(ip,4*np+nmumest(ic)+indam)=.true. 1743 filter_inc(ip,4*np+nmumest(ic)+namest(ic)+indam)=.true. 1744 end if 1745 end do 1746 1747 end do 1748 do jm=1,nmumest(ic) 1749 borni(4*np+jm)=XMU_MIN 1750 borns(4*np+jm)=XMU_MAX 1751 end do 1752 do jm=1,namest(ic) 1753 borni(4*np+nmumest(ic)+jm)=xmin 1754 borns(4*np+nmumest(ic)+jm)=xmax 1755 borni(4*np+nmumest(ic)+namest(ic)+jm)=xmin 1756 borns(4*np+nmumest(ic)+namest(ic)+jm)=xmax 1757 end do 1758 ! 1759 ! Point de depart 1760 nestim=0 1761 do ip=1,np 1762 par(ip)=std(ip) 1763 par(ip+np)=xmoyp(ip) 1764 par(2*np+ip)=ap(ip)/2.d0 1765 par(3*np+ip)=ap(ip)/2.d0 1766 nm1=nmp(ip)+1 1767 nm2=nmp(ip+1) 1768 nest=0 1769 somxmu=0.d0 1770 do jm=nm1,nm2 1771 if (estime(ic,jm)) then 1772 nest=nest+1 1773 ifem=repfem(jm) 1774 indam=iam(ic,ifem) 1775 if (nest.le.estmum(ip))then 1776 nestim=nestim+1 1777 par(4*np+nestim)=xmoym(jm) 1778 somxmu=somxmu+xmoym(jm) 1779 else 1780 xmoym(jm)=-somxmu 1781 end if 1782 par(4*np+nmumest(ic)+indam)= am(jm)/2.d0 1783 par(4*np+nmumest(ic)+namest(ic)+indam)= am(jm)/2.d0 1784 end if 1785 end do 1786 end do 1787 ! 1788 ! Marche le long du chromosome 1789 ifail=0 1790 lrtsol%lrtmax=-1.d75 1791 1792 ntotal=0 1793 do chr=1,nchr 1794 ntotal=ntotal+get_npo(chr) 1795 indexchr(chr)=ntotal 1796 end do 1797 1798 allocate (maxpar(nchr,nchr,get_maxnpo(),get_maxnpo(),npar)) 1799 allocate (fpsave1(nchr,nchr,get_maxnpo(),get_maxnpo(),np)) 1800 allocate (fmsave1(nchr,nchr,get_maxnpo(),get_maxnpo(),nm)) 1801 1802 !$OMP PARALLEL DEFAULT(SHARED) FIRSTPRIVATE(par) & 1803 !$OMP PRIVATE(ip,jm,chr,chr2,ilong2,n1,n,nstart,ig1,ig2,pddp,pddm,ngeno1,ngeno2,kd,kkd,kd2,ig,ifail,nestim,ii,f01_t) 1804 current_ic = ic 1805 call init_sub_2qtl 1806 carydf=>save_carydf 1807 somcddf => save_somcddf 1808 somydf => save_somydf 1809 effdf => save_effdf 1810 estmum => save_estmum 1811 somcd => save_somcd 1812 somy => save_somy 1813 cary => save_cary 1814 eff => save_eff 1815 fp_1 => save_fp_1 1816 fm_1 => save_fm_1 1817 f_1 = save_f_1 1818 f0 = save_f0 1819 1820 !$OMP DO 1821 do nprime=1,ntotal 1822 ! print *,'start:',nprime 1823 n=nprime 1824 do chr=nchr,1,-1 1825 if ( indexchr(chr) >= nprime) then 1826 current_chr = chr 1827 else 1828 n=n-get_npo(chr) 1829 end if 1830 end do 1831 chr=current_chr 1832 do chr2=chr,nchr 1833 if (chr2 == chr ) then 1834 nstart = n +1 1835 else 1836 nstart = 1 1837 end if 1838 1839 current_chr2 = chr2 1840 do n1=nstart,get_npo(chr2) 1841 ! print *,chr,chr2,n,n1 1842 do ip=1,np 1843 somppdf1(ip)=0.d0 1844 somppdf2(ip)=0.d0 1845 carppdf11(ip)=0.d0 1846 carppdf12(ip)=0.d0 1847 carppdf22(ip)=0.d0 1848 somppydf1(ip)=0.d0 1849 somppydf2(ip)=0.d0 1850 do jm=nmp(ip)+1,nmp(ip+1) 1851 do ig=ngenom(chr,jm)+1,ngenom(chr,jm+1) 1852 sompp1(ig)=0.d0 1853 sompm1(ig)=0.d0 1854 carpp1(ig)=0.d0 1855 carpm1(ig)=0.d0 1856 sompp1y(ig)=0.d0 1857 sompm1y(ig)=0.d0 1858 sompp1m1(ig)=0.d0 1859 do kd=ngend(chr,ig)+1,ngend(chr,ig+1) 1860 kkd=ndesc(chr,kd) 1861 if(presentc(ic,kkd)) then 1862 pddp(kd,1)=-pdd(chr,kd,1,n)-pdd(chr,kd,2,n)+pdd(chr,kd,3,n)+pdd(chr,kd,4,n) 1863 pddm(kd,1)=-pdd(chr,kd,1,n)+pdd(chr,kd,2,n)-pdd(chr,kd,3,n)+pdd(chr,kd,4,n) 1864 ! 1865 if (estime(ic,jm)) then 1866 sompp1(ig)=sompp1(ig)+pddp(kd,1) 1867 carpp1(ig)=carpp1(ig)+pddp(kd,1)*pddp(kd,1) 1868 sompp1y(ig)=sompp1y(ig)+pddp(kd,1)*y(ic,kkd) 1869 sompm1(ig)=sompm1(ig)+pddm(kd,1) 1870 carpm1(ig)=carpm1(ig)+pddm(kd,1)*pddm(kd,1) 1871 sompm1y(ig)=sompm1y(ig)+pddm(kd,1)*y(ic,kkd) 1872 sompp1m1(ig)=sompp1m1(ig)+pddp(kd,1)*pddm(kd,1) 1873 else 1874 somppdf1(ip)=somppdf1(ip)+pddp(kd,1) 1875 carppdf11(ip)=carppdf11(ip)+pddp(kd,1)*pddp(kd,1) 1876 somppydf1(ip)=somppydf1(ip)+pddp(kd,1)*y(ic,kkd) 1877 end if 1878 end if 1879 end do !kd 1880 end do !ig 1881 ! second chromosome 1882 do ig=ngenom(chr2,jm)+1,ngenom(chr2,jm+1) 1883 sompp2(ig)=0.d0 1884 sompm2(ig)=0.d0 1885 carpp2(ig)=0.d0 1886 carpm2(ig)=0.d0 1887 sompp2y(ig)=0.d0 1888 sompm2y(ig)=0.d0 1889 sompp2m2(ig)=0.d0 1890 1891 do kd=ngend(chr2,ig)+1,ngend(chr2,ig+1) 1892 kkd=ndesc(chr2,kd) 1893 if(presentc(ic,kkd)) then 1894 pddp(kd,2)=-pdd(chr2,kd,1,n1)-pdd(chr2,kd,2,n1)+pdd(chr2,kd,3,n1)+pdd(chr2,kd,4,n1) 1895 pddm(kd,2)=-pdd(chr2,kd,1,n1)+pdd(chr2,kd,2,n1)-pdd(chr2,kd,3,n1)+pdd(chr2,kd,4,n1) 1896 if (estime(ic,jm)) then 1897 sompp2(ig)=sompp2(ig)+pddp(kd,2) 1898 carpp2(ig)=carpp2(ig)+pddp(kd,2)*pddp(kd,2) 1899 sompp2y(ig)=sompp2y(ig)+pddp(kd,2)*y(ic,kkd) 1900 sompm2(ig)=sompm2(ig)+pddm(kd,2) 1901 carpm2(ig)=carpm2(ig)+pddm(kd,2)*pddm(kd,2) 1902 sompm2y(ig)=sompm2y(ig)+pddm(kd,2)*y(ic,kkd) 1903 sompp2m2(ig)=sompp2m2(ig)+pddp(kd,2)*pddm(kd,2) 1904 else 1905 somppdf2(ip)=somppdf2(ip)+pddp(kd,2) 1906 carppdf22(ip)=carppdf22(ip)+pddp(kd,2)*pddp(kd,2) 1907 somppydf2(ip)=somppydf2(ip)+pddp(kd,2)*y(ic,kkd) 1908 end if 1909 end if 1910 end do ! kd 1911 end do ! ig
opti_1qtl
[ Top ] [ m_qtlmap_analyse_unitrait ] [ Subroutines ]
NAME
opti_1qtl
DESCRIPTION
Computing statistical test across the chromosome
SOURCE
1166 subroutine opti_1qtl(ic,lrtsol) 1167 integer , intent(in) :: ic 1168 type(TYPE_LRT_SOLUTION) , intent(inout) :: lrtsol 1169 1170 ! 1171 ! Divers 1172 integer :: iuser(1) 1173 double precision ,dimension((3*np)+nmumest(ic)+namest(ic)) :: borni,borns,par 1174 double precision :: user(1) 1175 1176 integer :: chr,ibound,n,ip,jm,ngeno1,ngeno2,ix,ilong,nm1,nm2,ig 1177 integer :: npar,nd1,nd2,kd,kkd,ifail,ii 1178 integer :: nestim,nest,ifem,indam,nprime,ntotal,indexchr(nchr),rang 1179 real (kind=dp) :: pp, pm, somxmu,f1,fsave0 1180 logical , dimension(:,:),pointer :: filter_inc 1181 real (kind=dp) , dimension(:,:,:),pointer :: maxpar 1182 real (kind=dp) , dimension(:,:) ,pointer :: fmax 1183 real (kind=dp) , dimension(:,:,:),pointer :: fpsave1 1184 real (kind=dp) , dimension(:,:,:),pointer :: fmsave1 1185 1186 !Trick pour recuperer la valeur de ces tableaux 1187 !La parallelisation sur les caracteres empeche la recopie de ces tableaux (deja instancier NTHREADS fois) 1188 !on recupere le pointeur de ces tableau avant la parallelisation sur le groupe de liaison 1189 !puis on initialise tous les thread avec ces valeurs sauvegarde 1190 real (kind=dp) ,dimension(:),pointer :: save_carydf 1191 real (kind=dp) ,dimension(:),pointer :: save_somcddf 1192 real (kind=dp) ,dimension(:),pointer :: save_somydf 1193 integer ,dimension(:),pointer :: save_effdf 1194 integer ,dimension(:),pointer :: save_estmum 1195 real (kind=dp) ,dimension(:),pointer :: save_somcd 1196 real (kind=dp) ,dimension(:),pointer :: save_somy 1197 real (kind=dp) ,dimension(:),pointer :: save_cary 1198 integer ,dimension(:),pointer :: save_eff 1199 real (kind=dp) ,dimension(:),pointer :: save_fp0 1200 real (kind=dp) ,dimension(:),pointer :: save_fm0 1201 1202 call new(1,lrtsol) 1203 1204 !****************************************************************************** 1205 ! Calcul de la vraisemblance sous H1 1206 ! Parametres de maximisation 1207 1208 ibound=0 1209 npar=(3*np)+nmumest(ic)+namest(ic) 1210 1211 allocate (filter_inc(np,npar)) 1212 1213 save_carydf => carydf 1214 save_somcddf => somcddf 1215 save_somydf => somydf 1216 save_effdf => effdf 1217 save_estmum => estmum 1218 save_somcd => somcd 1219 save_somy => somy 1220 save_cary => cary 1221 save_eff => eff 1222 fsave0 = f0 1223 save_fp0 => fp0 1224 save_fm0 => fm0 1225 1226 filter_inc=.false. 1227 nestim=0 1228 do ip=1,np 1229 borni(ip)=SIG_MIN 1230 borns(ip)=SIG_MAX 1231 filter_inc(ip,ip)=.true. 1232 borni(ip+np)=XMU_MIN 1233 borns(ip+np)=XMU_MAX 1234 filter_inc(ip,ip+np)=.true. 1235 borni(2*np+ip)=AP_MIN 1236 borns(2*np+ip)=AP_MAX 1237 filter_inc(ip,2*np+ip)=.true. 1238 1239 nest=0 1240 do jm=nmp(ip)+1,nmp(ip+1) 1241 if(estime(current_ic,jm)) then 1242 nest=nest+1 1243 if(nest.le.estmum(ip)) then 1244 nestim=nestim+1 1245 filter_inc(ip,3*np+nestim)=.true. 1246 end if 1247 ifem=repfem(jm) 1248 indam=iam(ic,ifem) 1249 filter_inc(ip,3*np+nmumest(ic)+indam)=.true. 1250 end if 1251 end do 1252 end do 1253 !filter_inc(:,3*np+nmumest(ic):)=.true. 1254 1255 do jm=1,nmumest(ic) 1256 borni(3*np+jm)=XMU_MIN 1257 borns(3*np+jm)=XMU_MAX 1258 end do 1259 do jm=1,namest(ic) 1260 borni(3*np+nmumest(ic)+jm)=AM_MIN 1261 borns(3*np+nmumest(ic)+jm)=AM_MAX 1262 end do 1263 ! 1264 ! Point de depart 1265 do ip=1,np 1266 par(ip)=sig1(ip) 1267 par(ip+np)=xmu1p(ip) 1268 par(2*np+ip)=0.d0 1269 end do 1270 1271 1272 do jm=1,nmumest(ic) 1273 par(3*np+jm)=xmu1m(jm) 1274 end do 1275 do jm=1,namest(ic) 1276 par(3*np+nmumest(ic)+jm)=0.d0 1277 end do 1278 ! 1279 ! Marche le long du chromosome 1280 1281 lrtsol%lrtmax(0)=-1.d75 1282 lrtsol%chrmax(0)=1 1283 1284 ntotal=0 1285 do chr=1,nchr 1286 ntotal=ntotal+get_npo(chr) 1287 indexchr(chr)=ntotal 1288 end do 1289 1290 allocate (maxpar(nchr,get_maxnpo(),npar)) 1291 allocate (fmax(nchr,get_maxnpo())) 1292 allocate (fpsave1(nchr,get_maxnpo(),np)) 1293 allocate (fmsave1(nchr,get_maxnpo(),nm)) 1294 1295 !$OMP PARALLEL DEFAULT(SHARED) FIRSTPRIVATE(par) & 1296 !$OMP PRIVATE(ip,jm,chr,n,nd1,nd2,pp,pm,ngeno1,ngeno2,kd,kkd,ig,ifail,nestim,ii,f1) 1297 current_ic = ic 1298 call init_sub_1qtl 1299 carydf=>save_carydf 1300 somcddf => save_somcddf 1301 somydf => save_somydf 1302 effdf => save_effdf 1303 estmum => save_estmum 1304 somcd => save_somcd 1305 somy => save_somy 1306 cary => save_cary 1307 eff => save_eff 1308 f0 = fsave0 1309 fp0 => save_fp0 1310 fm0 => save_fm0 1311 1312 !$OMP DO 1313 do nprime=1,ntotal 1314 n=nprime 1315 do chr=nchr,1,-1 1316 if ( indexchr(chr) >= nprime) then 1317 current_chr = chr 1318 else 1319 n=n-get_npo(chr) 1320 end if 1321 end do 1322 chr=current_chr 1323 ! print *,chr,n 1324 do ip=1,np 1325 nm1=nmp(ip)+1 1326 nm2=nmp(ip+1) 1327 somppdf(ip)=0.d0 1328 carppdf(ip)=0.d0 1329 somppydf(ip)=0.d0 1330 do jm=nm1,nm2 1331 ngeno1=ngenom(chr,jm)+1 1332 ngeno2=ngenom(chr,jm+1) 1333 do ig=ngeno1,ngeno2 1334 sompp(ig)=0.d0 1335 carpp(ig)=0.d0 1336 somppy(ig)=0.d0 1337 sompm(ig)=0.d0 1338 carpm(ig)=0.d0 1339 sompmy(ig)=0.d0 1340 somppm(ig)=0.d0 1341 nd1=ngend(chr,ig)+1 1342 nd2=ngend(chr,ig+1) 1343 do kd=nd1,nd2 1344 kkd=ndesc(chr,kd) 1345 if(presentc(ic,kkd)) then 1346 pp=-pdd(chr,kd,1,n)-pdd(chr,kd,2,n)+pdd(chr,kd,3,n)+pdd(chr,kd,4,n) 1347 pm=-pdd(chr,kd,1,n)+pdd(chr,kd,2,n)-pdd(chr,kd,3,n)+pdd(chr,kd,4,n) 1348 if(estime(ic,jm)) then 1349 sompp(ig)=sompp(ig)+pp*cd(ic,kkd) 1350 carpp(ig)=carpp(ig)+pp*pp*cd(ic,kkd) 1351 somppy(ig)=somppy(ig)+pp*y(ic,kkd)*cd(ic,kkd) 1352 sompm(ig)=sompm(ig)+pm*cd(ic,kkd) 1353 carpm(ig)=carpm(ig)+pm*pm*cd(ic,kkd) 1354 sompmy(ig)=sompmy(ig)+pm*y(ic,kkd)*cd(ic,kkd) 1355 somppm(ig)=somppm(ig)+pp*pm*cd(ic,kkd) 1356 else 1357 somppdf(ip)=somppdf(ip)+pp*cd(ic,kkd) 1358 carppdf(ip)=carppdf(ip)+pp*pp*cd(ic,kkd) 1359 somppydf(ip)=somppydf(ip)+pp*y(ic,kkd)*cd(ic,kkd) 1360 end if 1361 end if 1362 end do 1363 end do 1364 end do 1365 end do 1366 1367 ! 1368 ! Optimisation de la vraisemblance a la position dx 1369 ifail=1 1370 !call minimizing_funct(npar,ibound,funct_1qtl,borni,borns,par,f1,iuser,user,ifail) 1371 1372 call minimizing_funct_family_sire(npar,ibound,funct_1qtl_family,filter_inc,fp1,borni,borns,par,f1,iuser,user,ifail) 1373 1374 fpsave1(chr,n,:) = fp1 1375 fmsave1(chr,n,:) = fm1 1376 fmax(chr,n) = f1 1377 if ( f1 < INIFINY_REAL_VALUE ) then 1378 lrtsol%lrt1(chr,n)=-2.d0*(f1-f0) 1379 else 1380 lrtsol%lrt1(chr,n)=0 1381 end if 1382 1383 do ii=1,np 1384 lrtsol%xlrp(chr,ii,n)=-2.d0*(fp1(ii)-fp0(ii)) 1385 lrtsol%pater_eff(chr,ii,n)=par(2*np+ii) 1386 end do 1387 1388 do ii=1,nm 1389 lrtsol%xlrm(chr,ii,n)=-2.d0*(fm1(ii)-fm0(ii)) 1390 end do 1391 1392 do ii=1,namest(ic) 1393 lrtsol%mater_eff(chr,ii,n)=par(3*np+nmumest(ic)+ii) 1394 end do 1395 1396 maxpar(chr,n,:)=par 1397 end do 1398 !$OMP END DO 1399 call end_sub_1qtl 1400 !$OMP END PARALLEL 1401 do chr=1,nchr 1402 do n=1,get_npo(chr) 1403 if (lrtsol%lrt1(chr,n)> lrtsol%lrtmax(0)) then 1404 par=maxpar(chr,n,:) 1405 lrtsol%lrtmax(0)=lrtsol%lrt1(chr,n) 1406 lrtsol%nxmax(0)=n 1407 lrtsol%chrmax(0)=chr 1408 f_1=fmax(chr,n) 1409 nestim=0 1410 do ip=1,np 1411 ap(ip)=par(2*np+ip) 1412 xmoyp(ip)=par(np+ip) 1413 std(ip)=par(ip) 1414 fp_1(ip)=fpsave1(chr,n,ip) 1415 nest=0 1416 somxmu=0.d0 1417 do jm=nmp(ip)+1,nmp(ip+1) 1418 fm_1(jm)=fmsave1(chr,n,jm) 1419 if (estime(ic,jm)) then 1420 nest=nest+1 1421 ifem=repfem(jm) 1422 indam=iam(ic,ifem) 1423 am(jm)=par(3*np+nmumest(ic)+indam) 1424 if (nest.le.estmum(ip)) then 1425 nestim=nestim+1 1426 xmoym(jm)=par(3*np+nestim) 1427 somxmu=somxmu+xmoym(jm) 1428 else 1429 xmoym(jm)=-somxmu 1430 end if 1431 else 1432 am(jm)=0.d0 1433 xmoym(jm)=0.d0 1434 end if 1435 end do 1436 end do 1437 end if 1438 end do 1439 end do 1440 1441 1442 deallocate (fpsave1) 1443 deallocate (fmsave1) 1444 deallocate (fmax) 1445 deallocate (maxpar) 1446 deallocate (filter_inc) 1447 1448 return 1449 end subroutine opti_1qtl
set_solution_hypothesis0
[ Top ] [ m_qtlmap_analyse_unitrait ] [ Subroutines ]
NAME
set_solution_hypothesis0
DESCRIPTION
SOURCE
2316 subroutine set_solution_hypothesis0(ic,incsol) 2317 integer ,intent(in) :: ic 2318 type(TYPE_INCIDENCE_SOLUTION) ,intent(inout) :: incsol 2319 2320 integer :: nteff,maxNbPar,jm,ifem,ip,ieff 2321 2322 allocate (incsol%sig(1,np)) 2323 2324 incsol%hypothesis=0 2325 ! Mean family 2326 nteff = 1 2327 if ( count(estime(ic,:)) > 0 ) nteff = 2 2328 2329 maxNbPar = max(np,count(estime(ic,:))) 2330 allocate (incsol%groupeName(nteff)) 2331 allocate (incsol%eqtl_print(nteff)) 2332 allocate (incsol%nbParameterGroup(nteff)) 2333 allocate (incsol%parameterName(nteff,maxNbPar)) 2334 allocate (incsol%paramaterValue(nteff,maxNbPar)) 2335 allocate (incsol%parameterVecsol(nteff,maxNbPar)) 2336 allocate (incsol%parameterPrecis(nteff,maxNbPar)) 2337 incsol%parameterPrecis=0.d0 2338 incsol%parameterVecsol=.true. 2339 incsol%eqtl_print=.true. 2340 do ip=1,np 2341 incsol%sig(1,ip) = sig1(ip)*sigt(ic) 2342 end do 2343 2344 ieff=1 2345 incsol%groupeName(ieff) = 'Mean Sire' 2346 incsol%nbParameterGroup(ieff)=np 2347 2348 do ip=1,np 2349 incsol%parameterName(ieff,ip)='Sire '//trim(pere(ip)) 2350 incsol%paramaterValue(ieff,ip) = xmu1p(ip)*sigt(ic) + xmut(ic) 2351 incsol%parameterVecsol(ieff,ip) = .true. 2352 end do 2353 2354 if ( count(estime(ic,:)) > 0 ) then 2355 ieff = ieff +1 2356 incsol%groupeName(ieff) = 'Mean dam' 2357 incsol%nbParameterGroup(ieff)=count(estime(ic,:)) 2358 ifem=0 2359 2360 do ip=1,np 2361 do jm=nmp(ip)+1,nmp(ip+1) 2362 if ( estime(ic,jm)) then 2363 ifem=ifem+1 2364 incsol%parameterName(ieff,ifem)='Dam '//trim(mere(jm))//" [Sire "//trim(pere(ip))//"]" 2365 incsol%paramaterValue(ieff,ifem) = xmu1m(ifem)*sigt(ic) 2366 incsol%parameterVecsol(ieff,ifem) = .true. 2367 end if 2368 end do 2369 end do 2370 2371 end if 2372 2373 end subroutine set_solution_hypothesis0
set_solution_hypothesis1
[ Top ] [ m_qtlmap_analyse_unitrait ] [ Subroutines ]
NAME
set_solution_hypothesis1
DESCRIPTION
SOURCE
2384 subroutine set_solution_hypothesis1(ic,incsol) 2385 integer , intent(in) :: ic 2386 type(TYPE_INCIDENCE_SOLUTION) ,intent(inout) :: incsol 2387 2388 integer :: nteff,maxNbPar,jm,ifem,ip,ieff 2389 2390 allocate (incsol%sig(1,np)) 2391 2392 incsol%hypothesis=1 2393 ! Mean family , Qtl effect 2394 nteff = 2 2395 if ( count(estime(ic,:)) > 0 ) nteff = 4 2396 2397 maxNbPar = max(np,count(estime(ic,:))) 2398 allocate (incsol%groupeName(nteff)) 2399 allocate (incsol%eqtl_print(nteff)) 2400 allocate (incsol%nbParameterGroup(nteff)) 2401 allocate (incsol%parameterName(nteff,maxNbPar)) 2402 allocate (incsol%paramaterValue(nteff,maxNbPar)) 2403 allocate (incsol%parameterVecsol(nteff,maxNbPar)) 2404 allocate (incsol%parameterPrecis(nteff,maxNbPar)) 2405 allocate (incsol%qtl_groupeName(1,1)) 2406 2407 incsol%parameterPrecis=0.d0 2408 incsol%parameterVecsol=.true. 2409 incsol%eqtl_print=.true. 2410 2411 do ip=1,np 2412 incsol%sig(1,ip) = std(ip)*sigt(ic) 2413 end do 2414 2415 ieff=1 2416 incsol%qtl_groupeName(1,1)=ieff 2417 incsol%groupeName(ieff) = 'Sire Qtl effect' 2418 incsol%nbParameterGroup(ieff)=np 2419 2420 do ip=1,np 2421 incsol%parameterName(ieff,ip)='Sire '//trim(pere(ip)) 2422 incsol%paramaterValue(ieff,ip) = ap(ip)*sigt(ic) 2423 end do 2424 2425 if ( count(estime(ic,:)) > 0 ) then 2426 ieff = ieff +1 2427 incsol%groupeName(ieff) = 'Dam Qtl effect' 2428 incsol%nbParameterGroup(ieff)=namest(ic) 2429 ifem=0 2430 do jm=1,nm 2431 if ( estime(ic,jm)) then 2432 ifem=ifem+1 2433 incsol%parameterName(ieff,ifem)='Dam '//trim(mere(jm)) 2434 incsol%paramaterValue(ieff,ifem) = am(jm)*sigt(ic) 2435 incsol%parameterVecsol(ieff,ifem) = .true. 2436 end if 2437 end do 2438 end if 2439 2440 2441 ieff = ieff +1 2442 incsol%groupeName(ieff) = 'Mean Sire' 2443 incsol%nbParameterGroup(ieff)=np 2444 2445 do ip=1,np 2446 incsol%parameterName(ieff,ip)='Sire '//trim(pere(ip)) 2447 incsol%paramaterValue(ieff,ip) = xmoyp(ip)*sigt(ic) + xmut(ic) 2448 incsol%parameterVecsol(ieff,ip) = .true. 2449 end do 2450 2451 if ( count(estime(ic,:)) > 0 ) then 2452 ieff = ieff +1 2453 incsol%groupeName(ieff) = 'Mean dam' 2454 incsol%nbParameterGroup(ieff)=count(estime(ic,:)) 2455 ifem=0 2456 do ip=1,np 2457 do jm=nmp(ip)+1,nmp(ip+1) 2458 if ( estime(ic,jm)) then 2459 ifem=ifem+1 2460 incsol%parameterName(ieff,ifem)='Dam '//trim(mere(jm))//" [Sire "//trim(pere(ip))//"]" 2461 incsol%paramaterValue(ieff,ifem) = xmoym(jm)*sigt(ic) 2462 incsol%parameterVecsol(ieff,ifem) = .true. 2463 end if 2464 end do 2465 end do 2466 end if 2467 2468 end subroutine set_solution_hypothesis1
set_solution_hypothesis2
[ Top ] [ m_qtlmap_analyse_unitrait ] [ Subroutines ]
NAME
set_solution_hypothesis2
DESCRIPTION
SOURCE
2479 subroutine set_solution_hypothesis2(ic,incsol) 2480 integer , intent(in) :: ic 2481 type(TYPE_INCIDENCE_SOLUTION) ,intent(inout) :: incsol 2482 2483 integer :: nteff,maxNbPar,jm,ifem,ip,ieff,iq 2484 2485 allocate (incsol%sig(1,np)) 2486 2487 incsol%hypothesis=2 2488 ! Mean family , Qtl effect 1, QTL effect2 2489 nteff = 3 2490 if ( count(estime(ic,:)) > 0 ) nteff = 6 2491 2492 maxNbPar = max(np,count(estime(ic,:))) 2493 allocate (incsol%groupeName(nteff)) 2494 allocate (incsol%eqtl_print(nteff)) 2495 allocate (incsol%nbParameterGroup(nteff)) 2496 allocate (incsol%parameterName(nteff,maxNbPar)) 2497 allocate (incsol%paramaterValue(nteff,maxNbPar)) 2498 allocate (incsol%parameterVecsol(nteff,maxNbPar)) 2499 allocate (incsol%parameterPrecis(nteff,maxNbPar)) 2500 allocate (incsol%qtl_groupeName(1,2)) 2501 2502 incsol%parameterPrecis=0.d0 2503 incsol%parameterVecsol=.true. 2504 incsol%nbParameterGroup=0.d0 2505 incsol%eqtl_print=.true. 2506 2507 do ip=1,np 2508 incsol%sig(1,ip) = std2(ip)*sigt(ic) 2509 end do 2510 2511 ieff=0 2512 do iq=1,2 2513 ieff = ieff +1 2514 incsol%qtl_groupeName(1,iq)=ieff 2515 incsol%groupeName(ieff) = 'Sire Qtl effect '//trim(str(iq)) 2516 incsol%nbParameterGroup(ieff)=np 2517 do ip=1,np 2518 incsol%parameterName(ieff,ip)='Sire '//trim(pere(ip)) 2519 incsol%paramaterValue(ieff,ip) = ap2(ip,iq)*sigt(ic) 2520 end do 2521 2522 2523 if ( count(estime(ic,:)) > 0 ) then 2524 ieff = ieff +1 2525 incsol%groupeName(ieff) = 'Dam Qtl effect '//trim(str(iq)) 2526 incsol%nbParameterGroup(ieff)=namest(ic) 2527 ifem=0 2528 do jm=1,nm 2529 if ( estime(ic,jm)) then 2530 ifem=ifem+1 2531 incsol%parameterName(ieff,ifem)='Dam '//trim(mere(jm)) 2532 incsol%paramaterValue(ieff,ifem) = am2(jm,iq)*sigt(ic) 2533 incsol%parameterVecsol(ieff,ifem) = .true. 2534 end if 2535 end do 2536 end if 2537 end do 2538 2539 ieff = ieff +1 2540 incsol%groupeName(ieff) = 'Mean Sire' 2541 incsol%nbParameterGroup(ieff)=np 2542 2543 2544 do ip=1,np 2545 incsol%parameterName(ieff,ip)='Sire '//trim(pere(ip)) 2546 incsol%paramaterValue(ieff,ip) = xmoyp2(ip)*sigt(ic) + xmut(ic) 2547 incsol%parameterVecsol(ieff,ip) = .true. 2548 end do 2549 2550 if ( count(estime(ic,:)) > 0 ) then 2551 ieff = ieff +1 2552 incsol%groupeName(ieff) = 'Mean dam' 2553 incsol%nbParameterGroup(ieff)=count(estime(ic,:)) 2554 ifem=0 2555 do ip=1,np 2556 do jm=nmp(ip)+1,nmp(ip+1) 2557 if ( estime(ic,jm)) then 2558 ifem=ifem+1 2559 incsol%parameterName(ieff,ifem)='Dam '//trim(mere(jm))//" [Sire "//trim(pere(ip))//"]" 2560 incsol%paramaterValue(ieff,ifem) = xmoym2(jm)*sigt(ic) 2561 incsol%parameterVecsol(ieff,ifem) = .true. 2562 end if 2563 end do 2564 end do 2565 2566 end if 2567 2568 end subroutine set_solution_hypothesis2