m_qtlmap_analyse_multitrait
[ Top ] [ ANALYSE ] [ Modules ]
NAME
m_qtlmap_analyse_multitrait
DESCRIPTION
NOTES
SEE ALSO
am
[ Top ] [ m_qtlmap_analyse_multitrait ] [ Variables ]
NAME
am
DESCRIPTION
The qtl dam effect found der H1
DIMENSIONS
ncar,nm
ap
[ Top ] [ m_qtlmap_analyse_multitrait ] [ Variables ]
NAME
ap
DESCRIPTION
The qtl sire effect found der H1
DIMENSIONS
ncar,np
carpm11
[ Top ] [ m_qtlmap_analyse_multitrait ] [ Variables ]
NAME
carpm11
DESCRIPTION
initialization : loop before the minimization of the likelihood : opti_mcar_1qtl
DIMENSIONS
carpp11
[ Top ] [ m_qtlmap_analyse_multitrait ] [ Variables ]
NAME
carpp11
DESCRIPTION
initialization : loop before the minimization of the likelihood : opti_mcar_1qtl
DIMENSIONS
carppdf11
[ Top ] [ m_qtlmap_analyse_multitrait ] [ Variables ]
NAME
carppdf11
DESCRIPTION
initialization : loop before the minimization of the likelihood : opti_mcar_1qtl
DIMENSIONS
cary11
[ Top ] [ m_qtlmap_analyse_multitrait ] [ Variables ]
NAME
cary11
DESCRIPTION
DIMENSIONS
carydf01
[ Top ] [ m_qtlmap_analyse_multitrait ] [ Variables ]
NAME
carydf01
DESCRIPTION
initialization : optinit_mcar
DIMENSIONS
corcd
[ Top ] [ m_qtlmap_analyse_multitrait ] [ Variables ]
NAME
corcd
DESCRIPTION
DIMENSIONS
current_chr
[ Top ] [ m_qtlmap_analyse_multitrait ] [ Variables ]
NAME
current_chr
DESCRIPTION
DIMENSIONS
eff
[ Top ] [ m_qtlmap_analyse_multitrait ] [ Variables ]
NAME
eff
DESCRIPTION
DIMENSIONS
effdf
[ Top ] [ m_qtlmap_analyse_multitrait ] [ Variables ]
NAME
effdf
DESCRIPTION
DIMENSIONS
effp
[ Top ] [ m_qtlmap_analyse_multitrait ] [ Variables ]
NAME
effp
DESCRIPTION
DIMENSIONS
estmum
[ Top ] [ m_qtlmap_analyse_multitrait ] [ Variables ]
NAME
estmum
DESCRIPTION
DIMENSIONS
f01
[ Top ] [ m_qtlmap_analyse_multitrait ] [ Variables ]
NAME
f01
DESCRIPTION
DIMENSIONS
fm2
[ Top ] [ m_qtlmap_analyse_multitrait ] [ Variables ]
NAME
fm2
DESCRIPTION
the likelihood by full sib family
DIMENSIONS
nm
fm21
[ Top ] [ m_qtlmap_analyse_multitrait ] [ Variables ]
NAME
fm21
DESCRIPTION
maximum value of likelihood by full sib family under H0
DIMENSIONS
nm
fp2
[ Top ] [ m_qtlmap_analyse_multitrait ] [ Variables ]
NAME
fp2
DESCRIPTION
the likelihood by half sib family
DIMENSIONS
np
fp21
[ Top ] [ m_qtlmap_analyse_multitrait ] [ Variables ]
NAME
fp21
DESCRIPTION
maximum value of likelihood by half sib family under H0
DIMENSIONS
np
rhoi
[ Top ] [ m_qtlmap_analyse_multitrait ] [ Variables ]
NAME
rhoi
DESCRIPTION
DIMENSIONS
rhoi2
[ Top ] [ m_qtlmap_analyse_multitrait ] [ Variables ]
NAME
rhoi2
DESCRIPTION
DIMENSIONS
rhoi2_1
[ Top ] [ m_qtlmap_analyse_multitrait ] [ Variables ]
NAME
rhoi2_1
DESCRIPTION
DIMENSIONS
rhoi2_2
[ Top ] [ m_qtlmap_analyse_multitrait ] [ Variables ]
NAME
rhoi2_2
DESCRIPTION
DIMENSIONS
sig01
[ Top ] [ m_qtlmap_analyse_multitrait ] [ Variables ]
NAME
sig01
DESCRIPTION
DIMENSIONS
sig2
[ Top ] [ m_qtlmap_analyse_multitrait ] [ Variables ]
NAME
sig2
DESCRIPTION
DIMENSIONS
sompm11
[ Top ] [ m_qtlmap_analyse_multitrait ] [ Variables ]
NAME
sompm11
DESCRIPTION
initialization : loop before the minimization of the likelihood : opti_mcar_1qtl
DIMENSIONS
sompmy11
[ Top ] [ m_qtlmap_analyse_multitrait ] [ Variables ]
NAME
sompmy11
DESCRIPTION
initialization : loop before the minimization of the likelihood : opti_mcar_1qtl
DIMENSIONS
sompp11
[ Top ] [ m_qtlmap_analyse_multitrait ] [ Variables ]
NAME
sompp11
DESCRIPTION
initialization : loop before the minimization of the likelihood : opti_mcar_1qtl
DIMENSIONS
somppdf11
[ Top ] [ m_qtlmap_analyse_multitrait ] [ Variables ]
NAME
somppdf11
DESCRIPTION
initialization : loop before the minimization of the likelihood : opti_mcar_1qtl
DIMENSIONS
somppm11
[ Top ] [ m_qtlmap_analyse_multitrait ] [ Variables ]
NAME
somppm11
DESCRIPTION
initialization : loop before the minimization of the likelihood : opti_mcar_1qtl
DIMENSIONS
somppy11
[ Top ] [ m_qtlmap_analyse_multitrait ] [ Variables ]
NAME
somppy11
DESCRIPTION
initialization : loop before the minimization of the likelihood : opti_mcar_1qtl
DIMENSIONS
somppydf11
[ Top ] [ m_qtlmap_analyse_multitrait ] [ Variables ]
NAME
somppydf11
DESCRIPTION
initialization : loop before the minimization of the likelihood : opti_mcar_1qtl
DIMENSIONS
somy11
[ Top ] [ m_qtlmap_analyse_multitrait ] [ Variables ]
NAME
somy11
DESCRIPTION
DIMENSIONS
somydf01
[ Top ] [ m_qtlmap_analyse_multitrait ] [ Variables ]
NAME
somydf01
DESCRIPTION
initialization : optinit_mcar
DIMENSIONS
somyy11
[ Top ] [ m_qtlmap_analyse_multitrait ] [ Variables ]
NAME
somyy11
DESCRIPTION
DIMENSIONS
somyydf01
[ Top ] [ m_qtlmap_analyse_multitrait ] [ Variables ]
NAME
somyydf01
DESCRIPTION
initialization : optinit_mcar
DIMENSIONS
std
[ Top ] [ m_qtlmap_analyse_multitrait ] [ Variables ]
NAME
std
DESCRIPTION
the standart deviation (residual) found under H1
DIMENSIONS
ncar,np
xmoym
[ Top ] [ m_qtlmap_analyse_multitrait ] [ Variables ]
NAME
xmoym
DESCRIPTION
The polygenic dam effect found under H1
DIMENSIONS
nm
xmoyp
[ Top ] [ m_qtlmap_analyse_multitrait ] [ Variables ]
NAME
xmoyp
DESCRIPTION
The polygenic sire effect found under H1
DIMENSIONS
ncar,np
xmu01m
[ Top ] [ m_qtlmap_analyse_multitrait ] [ Variables ]
NAME
xmu01m
DESCRIPTION
DIMENSIONS
xmu01p
[ Top ] [ m_qtlmap_analyse_multitrait ] [ Variables ]
NAME
xmu01p
DESCRIPTION
DIMENSIONS
xmu2m
[ Top ] [ m_qtlmap_analyse_multitrait ] [ Variables ]
NAME
xmu2m
DESCRIPTION
DIMENSIONS
xmu2p
[ Top ] [ m_qtlmap_analyse_multitrait ] [ Variables ]
NAME
xmu2p
DESCRIPTION
DIMENSIONS
end_analyse_multitrait
[ Top ] [ m_qtlmap_analyse_multitrait ] [ Subroutines ]
NAME
end_analyse_multitrait
DESCRIPTION
deallocation of buffer/solution arrays
SOURCE
572 subroutine end_analyse_multitrait 573 574 deallocate (somydf01) 575 deallocate (carydf01) 576 deallocate (somyydf01) 577 deallocate (sig01) 578 deallocate (xmu01p) 579 deallocate (xmu01m) 580 deallocate (rhoi) 581 deallocate (fp21) 582 deallocate (fm21) 583 deallocate (somy11) 584 deallocate (cary11) 585 deallocate (somyy11) 586 deallocate (xmu2p) 587 deallocate (sig2) 588 deallocate (xmu2m) 589 deallocate (rhoi2) 590 deallocate (std) 591 deallocate (ap) 592 deallocate (xmoyp) 593 deallocate (am) 594 deallocate (xmoym) 595 deallocate (rhoi2_2) 596 deallocate (rhoi2_1) 597 598 deallocate (effdf) 599 deallocate (estmum) 600 deallocate (eff) 601 deallocate (effp) 602 603 deallocate (corcd) 604 605 end subroutine end_analyse_multitrait
end_sub_1qtl
[ Top ] [ m_qtlmap_analyse_multitrait ] [ Subroutines ]
NAME
end_sub_1qtl
DESCRIPTION
deallocation of buffer/solution arrays
SOURCE
549 subroutine end_sub_1qtl 550 deallocate (sompp11) 551 deallocate (sompm11) 552 deallocate (carpp11) 553 deallocate (carpm11) 554 deallocate (somppm11) 555 deallocate (sompmy11) 556 deallocate (somppy11) 557 deallocate (somppdf11) 558 deallocate (carppdf11) 559 deallocate (somppydf11) 560 deallocate (fp2) 561 deallocate (fm2) 562 563 end subroutine end_sub_1qtl
funct_mcar_0qtl
[ Top ] [ m_qtlmap_analyse_multitrait ] [ Subroutines ]
NAME
funct_mcar_0qtl
DESCRIPTION
Calcul de la statistique de test multicaractere le long du chromosome
SOURCE
919 subroutine funct_mcar_0qtl(n,x,f,iuser,user) 920 integer , intent(in) :: n 921 real (kind=dp) ,dimension(n), intent(inout) :: x 922 real (kind=dp) , intent(inout) :: f 923 integer , dimension(1), intent(in) :: iuser 924 real (kind=dp) ,dimension(1), intent(in) :: user 925 integer :: ip 926 927 do ip=1,np 928 call funct_mcar_0qtl_family(ip,n,x,fp21(ip) ,iuser,user) 929 end do 930 931 f = sum(fp21) 932 return 933 934 end subroutine funct_mcar_0qtl
funct_mcar_0qtl_family
[ Top ] [ m_qtlmap_analyse_multitrait ] [ Subroutines ]
NAME
funct_mcar_0qtl_family
DESCRIPTION
Calcul de la statistique de test multicaractere le long du chromosome
SOURCE
943 subroutine funct_mcar_0qtl_family(ip,n,x,f,iuser,user) 944 integer , intent(in) :: ip,n 945 real (kind=dp) ,dimension(n), intent(inout) :: x 946 real (kind=dp) , intent(inout) :: f 947 integer , dimension(1), intent(in) :: iuser 948 real (kind=dp) ,dimension(1), intent(in) :: user 949 ! 950 ! Tableaux dimensionnes selon np, le nombre de peres 951 real (kind=dp) ,dimension(np) :: det 952 953 ! 954 ! Tableaux dimensionnes selon le nombre nc de caracteres 955 real (kind=dp) ,dimension(ncar,np) :: sigc,xmupc,xmup2c 956 real (kind=dp) ,dimension(ncar,nm) :: xmumc,xmum2c,xmupmc 957 real (kind=dp) ,dimension(ncar,ncar,np) ::covc,vc 958 real (kind=dp) ,dimension(ncar+1,ncar) :: A 959 real (kind=dp) ,dimension(ncar,ncar) :: rhoc,Ab 960 real (kind=dp) ,dimension(ncar) :: wrkspce,somxmu 961 ! 962 ! Divers 963 integer ::i,k,j,nt,nestim,ifail,ic,jc,nest 964 integer :: jm 965 real (kind=dp) :: rh,zdf,vdf,zpf,vpf,determ 966 !*********************************************************************************** 967 f=0.d0 968 ! R�cup�ration des corr�lations 969 k=0 970 do i=2,ncar 971 do j=1,i-1 972 k=k+1 973 nt=2*ncar*np+ncar*nmumest(1)+k 974 rh=x(nt) 975 ! print *,"rh:",rh 976 rhoc(j,i)=((dexp(rh)-1.d0)/(dexp(rh)+1.d0)) 977 rhoc(i,j)=rhoc(j,i) 978 ! print *,"rhoc ",i,j,rhoc(i,j) 979 end do 980 end do 981 ! 982 ! Calcul de l'inverse et du determinant de la matrice de var-cov 983 nestim=0 984 ! 985 ! Initialisation des parametres matrice de variance cov 986 do i=1,ncar 987 sigc(i,ip)=x((i-1)*np+ip) 988 covc(i,i,ip)=sigc(i,ip)*sigc(i,ip) 989 if(i.gt.1) then 990 do j=1,i-1 991 covc(i,j,ip)=rhoc(i,j)*sigc(i,ip)*sigc(j,ip) 992 covc(j,i,ip)=covc(i,j,ip) 993 end do 994 end if 995 end do 996 ! 997 ! Si 2 caracteres, calcul a la main 998 if(ncar.eq.2) then 999 det(ip)=covc(1,1,ip)*covc(2,2,ip)-covc(1,2,ip)*covc(2,1,ip) 1000 1001 vc(1,1,ip)=covc(2,2,ip)/det(ip) 1002 vc(2,2,ip)=covc(1,1,ip)/det(ip) 1003 vc(1,2,ip)=-covc(1,2,ip)/det(ip) 1004 vc(2,1,ip)=-covc(2,1,ip)/det(ip) 1005 else 1006 ! sinon appelle nag 1007 do i=1,ncar 1008 do j=1,ncar 1009 A(i,j)=covc(i,j,ip) 1010 Ab(i,j)=covc(i,j,ip) 1011 end do 1012 end do 1013 ifail=1 1014 ! call MATH_QTLMAP_F03ABF(Ab,ncar,ncar,determ,ifail) 1015 ! if (ifail.ne.0) then 1016 ! MODIFICATION COMPORTEMENT 1017 ! OFI:Si le determinant n'est pas calculable, on considere que les parametre de la fonction 1018 ! ne sont pas viable, on met une valeur tres haute au LRT 1019 1020 ! fp21 = INIFINY_REAL_VALUE 1021 ! f = INIFINY_REAL_VALUE 1022 ! fm21 = INIFINY_REAL_VALUE 1023 ! return 1024 ! end if 1025 CALL MATH_QTLMAP_INVDETMATSYM(ncar,A,ncar+1,determ,ifail) 1026 det(ip)=determ 1027 1028 !ifail=1 1029 !call MATH_QTLMAP_F01ADF(ncar,A,ncar+1,ifail) 1030 if (ifail.ne.0 .or. det(ip)< 1.do-9) then 1031 !print*,'ifail=',ifail,' ; inversion matrice0' 1032 !stop 1033 ! MODIFICATION COMPORTEMENT 1034 ! OFI:Si le determinant n'est pas calculable, on considere que les parametre de la fonction 1035 ! ne sont pas viable, on met une valeur tres haute au LRT 1036 fp21 = INIFINY_REAL_VALUE 1037 f = INIFINY_REAL_VALUE 1038 fm21 = INIFINY_REAL_VALUE 1039 return 1040 end if 1041 1042 do i=1,ncar 1043 do j=1,i 1044 vc(i,j,ip)=A(i+1,j) 1045 vc(j,i,ip)=vc(i,j,ip) 1046 end do 1047 end do 1048 end if 1049 ! 1050 ! Initialisation polyg�nique pere 1051 do i=1,ncar 1052 xmupc(i,ip)=x(ncar*np+(i-1)*np+ip) 1053 xmup2c(i,ip)=xmupc(i,ip)*xmupc(i,ip) 1054 end do 1055 ! 1056 ! Calcul de la vraisemblance 1057 ! 1058 ! Vraisemblance demi-fr�res 1059 zdf=0.d0 1060 vdf=0.d0 1061 do ic=1,ncar 1062 vdf=carydf01(ic,ip)+dble(effdf(ip))*xmup2c(ic,ip)-2.d0*xmupc(ic,ip)*somydf01(ic,ip) 1063 zdf=zdf + vc(ic,ic,ip)*vdf 1064 do jc=1,ic-1 1065 vdf=somyydf01(ic,jc,ip) - xmupc(ic,ip)*somydf01(jc,ip) & 1066 - xmupc(jc,ip)*somydf01(ic,ip)+ dble(effdf(ip))*xmupc(ic,ip)*xmupc(jc,ip) 1067 zdf = zdf + 2.d0*vdf*vc(ic,jc,ip) 1068 end do 1069 end do 1070 f= 0.5d0 * zdf + 0.5*dble(effdf(ip))*dlog(det(ip)) 1071 ! 1072 ! Vraisemblance plein-fr�res 1073 nest=0 1074 somxmu=0.d0 1075 do jm=nmp(ip)+1,nmp(ip+1) 1076 zpf=0.d0 1077 ! WARNING******************************** 1078 ! On estime par rapport a quel caractere 1079 if(estime(ncar,jm)) then 1080 nest=nest+1 1081 if (ip>1) then 1082 nestim=sum(estmum(:ip-1))+nest 1083 else 1084 nestim=nest 1085 end if 1086 1087 do ic=1,ncar 1088 ! print *,"ip:",ip,nest,estmum(ip) 1089 if(nest.le.estmum(ip)) then 1090 ! if(ic.eq.1) nestim=nestim+1 1091 ! print *,"polygenique mere:",x(2*ncar*np+(ic-1)*nmumest(1)+nestim),' ic:',ic 1092 xmumc(ic,jm)=x(2*ncar*np+(ic-1)*nmumest(1)+nestim) 1093 somxmu(ic)=somxmu(ic)+xmumc(ic,jm) 1094 else 1095 xmumc(ic,jm)=-somxmu(ic) 1096 end if 1097 xmum2c(ic,jm)=xmumc(ic,jm)*xmumc(ic,jm) 1098 xmupmc(ic,jm)=xmupc(ic,ip)+xmumc(ic,jm) 1099 1100 vpf=dble(eff(jm))*(xmup2c(ic,ip)+xmum2c(ic,jm) & 1101 +2.d0*xmupc(ic,ip)*xmumc(ic,jm)) & 1102 -2.d0*xmupmc(ic,jm)*somy11(ic,jm)+cary11(ic,jm) 1103 zpf=zpf + vpf*vc(ic,ic,ip) 1104 if (ic.gt.1) then 1105 do jc=1,ic-1 1106 vpf=dble(eff(jm))*xmupmc(ic,jm)*xmupmc(jc,jm) & 1107 +somyy11(ic,jc,jm)-xmupmc(ic,jm)*somy11(jc,jm)-xmupmc(jc,jm)*somy11(ic,jm) 1108 zpf = zpf + 2.d0*vpf*vc(ic,jc,ip) 1109 end do 1110 end if 1111 end do 1112 fm21(jm)= 0.5d0*zpf+0.5d0*dble(eff(jm))*dlog(det(ip)) 1113 else 1114 fm21(jm)=0.d0 1115 end if 1116 f=f+fm21(jm) 1117 end do 1118 return 1119 end subroutine funct_mcar_0qtl_family
funct_mcar_1qtl
[ Top ] [ m_qtlmap_analyse_multitrait ] [ Subroutines ]
NAME
funct_mcar_1qtl
DESCRIPTION
Calcul de la vraisemblance et de ses derivees partielles sous H1
SOURCE
1406 subroutine funct_mcar_1qtl(n,x,f,iuser,user) 1407 1408 integer , intent(in) :: n 1409 real (kind=dp) ,dimension(n), intent(in) :: x 1410 real (kind=dp) , intent(inout) :: f 1411 integer , dimension(1) ,intent(in) :: iuser 1412 real (kind=dp) ,dimension(1),intent(in) :: user 1413 1414 integer :: ip !fp2 1415 1416 do ip=1,np 1417 call funct_mcar_1qtl_family(ip,n,x,fp2(ip),iuser,user) 1418 end do 1419 f = sum(fp2) 1420 1421 end subroutine funct_mcar_1qtl
funct_mcar_1qtl_family
[ Top ] [ m_qtlmap_analyse_multitrait ] [ Subroutines ]
NAME
funct_mcar_1qtl_family
DESCRIPTION
Calcul de la vraisemblance et de ses derivees partielles sous H1
SOURCE
1430 subroutine funct_mcar_1qtl_family(ip,n,x,f,iuser,user) 1431 1432 integer , intent(in) :: ip,n 1433 real (kind=dp) ,dimension(n), intent(in) :: x 1434 real (kind=dp) , intent(inout) :: f 1435 integer , dimension(1) ,intent(in) :: iuser 1436 real (kind=dp) ,dimension(1),intent(in) :: user 1437 1438 real (kind=dp) ,dimension(np) :: det 1439 ! 1440 ! Tableaux dimensionnes selon le nombre nc de caracteres 1441 real (kind=dp),dimension(ncar,np):: sigc,xmupc,xmup2c,apc,ap2c 1442 real (kind=dp),dimension(ncar) :: wrkspce,somxmu 1443 real (kind=dp),dimension(ncar,nm):: xmumc,xmum2c,xmupmc,amc,am2c 1444 real (kind=dp),dimension(ncar,ncar,np) ::covc,vc 1445 real (kind=dp),dimension(ncar+1,ncar) :: A 1446 real (kind=dp),dimension(ncar,ncar) :: rhoc,Ab 1447 1448 real (kind=dp)::rh,zdf,vdf,vmere,zpf,vpf,determ,savezpf 1449 integer :: i,j,k,nt,nestim,ifail,it,nest,jm 1450 integer :: ifem,indam,ic,ngeno1,ngeno2,ig,ipar,itindi,jc 1451 1452 !*********************************************************************************** 1453 f=0.d0 1454 k=0 1455 rhoc=0.d0 1456 ! Recuperation des correlations 1457 do i=2,ncar 1458 do j=1,i-1 1459 k=k+1 1460 nt=3*ncar*np+ncar*nmumest(1)+ncar*namest(1)+k 1461 rh=x(nt) 1462 rhoc(j,i)=((dexp(rh)-1.d0)/(dexp(rh)+1.d0)) 1463 rhoc(i,j)=rhoc(j,i) 1464 end do 1465 end do 1466 ! 1467 ! Calcul de l'inverse et du determinant de la matrice de var-cov 1468 nestim=0 1469 ! 1470 ! Initialisation des parametres 1471 do i=1,ncar 1472 sigc(i,ip)=dabs(x((i-1)*np+ip)) 1473 covc(i,i,ip)=sigc(i,ip)*sigc(i,ip) 1474 if(i.gt.1) then 1475 do j=1,i-1 1476 covc(i,j,ip)=rhoc(i,j)*sigc(i,ip)*sigc(j,ip) 1477 covc(j,i,ip)=covc(i,j,ip) 1478 end do 1479 end if 1480 xmupc(i,ip)=x(ncar*np+(i-1)*np+ip) 1481 xmup2c(i,ip)=xmupc(i,ip)*xmupc(i,ip) 1482 apc(i,ip)=x(ip+2*ncar*np+(i-1)*np) 1483 ap2c(i,ip)=apc(i,ip)*apc(i,ip) 1484 end do !!i 1485 ! Si 2 caracteres, calcul a la main 1486 if(ncar.eq.2) then 1487 det(ip)=covc(1,1,ip)*covc(2,2,ip)-covc(1,2,ip)*covc(2,1,ip) 1488 vc(1,1,ip)=covc(2,2,ip)/det(ip) 1489 vc(2,2,ip)=covc(1,1,ip)/det(ip) 1490 vc(1,2,ip)=-covc(1,2,ip)/det(ip) 1491 vc(2,1,ip)=-covc(2,1,ip)/det(ip) 1492 else 1493 ! sinon appelle nag 1494 do i=1,ncar 1495 do j=1,ncar 1496 A(i,j)=covc(i,j,ip) 1497 Ab(i,j)=covc(i,j,ip) 1498 end do 1499 end do 1500 ifail=1 1501 ! call MATH_QTLMAP_F03ABF(Ab,ncar,ncar,determ,ifail) 1502 1503 ! if (ifail.ne.0) then 1504 ! MODIFICATION COMPORTEMENT 1505 ! OFI:Si le determinant n'est pas calculable, on considere que les parametre de la fonction 1506 ! ne sont pas viable, on met une valeur tres haute au LRT 1507 ! fp21 = INIFINY_REAL_VALUE 1508 ! f = INIFINY_REAL_VALUE 1509 ! fm21 = INIFINY_REAL_VALUE 1510 ! return 1511 ! end if 1512 CALL MATH_QTLMAP_INVDETMATSYM(ncar,A,ncar+1,determ,ifail) 1513 det(ip)=determ 1514 ! ifail=1 1515 ! call MATH_QTLMAP_F01ADF(ncar,A,ncar+1,ifail) 1516 if (ifail.ne.0 .or. det(ip)<=0 ) then 1517 !call stop_application('ifail='//trim(str(ifail))//' ; inversion matrice') 1518 ! MODIFICATION COMPORTEMENT 1519 ! OFI:Si le determinant n'est pas calculable, on considere que les parametre de la fonction 1520 ! ne sont pas viable, on met une valeur tres haute au LRT 1521 fp21 = INIFINY_REAL_VALUE 1522 f = INIFINY_REAL_VALUE 1523 fm21 = INIFINY_REAL_VALUE 1524 return 1525 end if 1526 do i=1,ncar 1527 do j=1,i 1528 vc(i,j,ip)=A(i+1,j) 1529 vc(j,i,ip)=vc(i,j,ip) 1530 ! print *,vc(i,j,ip) 1531 end do 1532 end do 1533 end if 1534 ! Initialisation parametres des m�res 1535 nest=0 1536 somxmu=0.d0 1537 do jm=nmp(ip)+1,nmp(ip+1) 1538 !WARNING -->utilisation de estime 1539 if(count(estime(:,jm)) == ncar) then 1540 nest=nest+1 1541 if (ip>1) then 1542 nestim=sum(estmum(:ip-1))+nest 1543 else 1544 nestim=nest 1545 end if 1546 1547 do i=1,ncar 1548 if(nest.le.estmum(ip)) then 1549 ! if (i.eq.1) nestim=nestim+1 1550 xmumc(i,jm)=x(3*ncar*np+(i-1)*nmumest(1)+nestim) 1551 somxmu(i)=somxmu(i)+xmumc(i,jm) 1552 else 1553 xmumc(i,jm)=-somxmu(i) 1554 end if 1555 xmum2c(i,jm)=xmumc(i,jm)*xmumc(i,jm) 1556 xmupmc(i,jm)=xmupc(i,ip)+xmumc(i,jm) 1557 ifem=repfem(jm) 1558 indam=iam(1,ifem) 1559 amc(i,jm)=x(3*ncar*np+ncar*nmumest(1)+(i-1)*namest(1)+indam) 1560 am2c(i,jm)=amc(i,jm)*amc(i,jm) 1561 end do 1562 end if 1563 end do 1564 ! end do !!ip 1565 ! 1566 ! Calcul de la vraisemblance 1567 ! do ip=1,np 1568 ! Vraisemblance demi-fr�res 1569 zdf=0.d0 1570 vdf=0.d0 1571 do ic=1,ncar 1572 vdf=carydf01(ic,ip)+dble(effdf(ip))*xmup2c(ic,ip)-2.d0*xmupc(ic,ip)*somydf01(ic,ip) 1573 vdf = vdf + ap2c(ic,ip)*carppdf11(ip) & 1574 +2.d0*xmupc(ic,ip)*apc(ic,ip)*somppdf11(ip) & 1575 -2.d0*apc(ic,ip)*somppydf11(ic,ip) 1576 zdf=zdf + vc(ic,ic,ip)*vdf 1577 do jc=1,ic-1 1578 vdf=somyydf01(ic,jc,ip) - xmupc(ic,ip)*somydf01(jc,ip) & 1579 - xmupc(jc,ip)*somydf01(ic,ip) & 1580 + dble(effdf(ip))*xmupc(ic,ip)*xmupc(jc,ip) 1581 vdf=vdf - apc(ic,ip)*somppydf11(jc,ip) & 1582 - apc(jc,ip)*somppydf11(ic,ip) & 1583 + apc(ic,ip)*apc(jc,ip)*carppdf11(ip) & 1584 + xmupc(ic,ip)*apc(jc,ip)*somppdf11(ip) & 1585 + xmupc(jc,ip)*apc(ic,ip)*somppdf11(ip) 1586 zdf = zdf + 2.d0*vdf*vc(ic,jc,ip) 1587 end do 1588 end do 1589 1590 f= 0.5d0 * zdf + 0.5d0*dble(effdf(ip))*dlog(det(ip)) 1591 1592 ! 1593 ! Vraisemblance plein-fr�res 1594 do jm=nmp(ip)+1,nmp(ip+1) 1595 itindi=0 1596 zpf=0.d0 1597 ! WARNING : utilisation de estime 1598 if(estime(ncar,jm)) then 1599 do ic=1,ncar 1600 vpf=dble(eff(jm))*(xmup2c(ic,ip)+xmum2c(ic,jm) & 1601 +2.d0*xmupc(ic,ip)*xmumc(ic,jm)) & 1602 -2.d0*xmupmc(ic,jm)*somy11(ic,jm)+cary11(ic,jm) 1603 zpf=zpf + vpf*vc(ic,ic,ip) 1604 if (ic.gt.1) then 1605 do jc=1,ic-1 1606 vpf=dble(eff(jm))*xmupmc(ic,jm)*xmupmc(jc,jm) & 1607 +somyy11(ic,jc,jm)-xmupmc(ic,jm)*somy11(jc,jm) & 1608 -xmupmc(jc,jm)*somy11(ic,jm) 1609 zpf = zpf + 2.d0*vpf*vc(ic,jc,ip) 1610 end do 1611 end if 1612 end do 1613 ngeno1=ngenom(current_chr,jm)+1 1614 ngeno2=ngenom(current_chr,jm+1) 1615 ! si un seul genotype possible pour la mere, ne passe pas par les exp(vpf) 1616 if ((ngeno2-ngeno1).eq.0) then 1617 ! if (.true.) then 1618 ig=ngeno1 1619 do ic=1,ncar 1620 vpf= ap2c(ic,ip)*carpp11(ig) & 1621 +am2c(ic,jm)*carpm11(ig) & 1622 +2.d0*apc(ic,ip)*amc(ic,jm)*somppm11(ig) & 1623 +2.d0*xmupmc(ic,jm)*(apc(ic,ip)*sompp11(ig) & 1624 +amc(ic,jm)*sompm11(ig)) & 1625 -2.d0*apc(ic,ip)*somppy11(ic,ig) & 1626 -2.d0*amc(ic,jm)*sompmy11(ic,ig) 1627 zpf= zpf + vpf*vc(ic,ic,ip) 1628 ! if ((ngeno2-ngeno1)>0) print *,'ic:',ic,zpf,vpf,vc(ic,ic,ip) 1629 if (ic.gt.1) then 1630 do jc=1,ic-1 1631 vpf= xmupmc(ic,jm)*(apc(jc,ip)*sompp11(ig) & 1632 +amc(jc,jm)*sompm11(ig)) & 1633 +xmupmc(jc,jm)*(apc(ic,ip)*sompp11(ig) & 1634 +amc(ic,jm)*sompm11(ig)) & 1635 -apc(jc,ip)*somppy11(ic,ig) & 1636 -amc(jc,jm)*sompmy11(ic,ig) & 1637 -apc(ic,ip)*somppy11(jc,ig) & 1638 -amc(ic,jm)*sompmy11(jc,ig) & 1639 +apc(ic,ip)*apc(jc,ip)*carpp11(ig) & 1640 +amc(ic,jm)*amc(jc,jm)*carpm11(ig) & 1641 +(apc(ic,ip)*amc(jc,jm)+amc(ic,jm)*apc(jc,ip)) & 1642 *somppm11(ig) 1643 zpf= zpf + 2.d0*vpf*vc(ic,jc,ip) 1644 end do !! jc 1645 end if 1646 end do !! ic 1647 1648 fm2(jm)= 0.5d0*zpf+0.5d0*dble(eff(jm))*dlog(det(ip)) 1649 else 1650 ! si plusieurs genotypes possibles pour la mere, 1651 vmere=0.d0 1652 savezpf=zpf 1653 do ig=ngeno1,ngeno2 !attention NGENO2 ******************************************************************************* 1654 zpf=savezpf 1655 do ic=1,ncar 1656 vpf= ap2c(ic,ip)*carpp11(ig) & 1657 +am2c(ic,jm)*carpm11(ig) & 1658 +2.d0*apc(ic,ip)*amc(ic,jm)*somppm11(ig) & 1659 +2.d0*xmupmc(ic,jm)*(apc(ic,ip)*sompp11(ig) & 1660 +amc(ic,jm)*sompm11(ig)) & 1661 -2.d0*apc(ic,ip)*somppy11(ic,ig) & 1662 -2.d0*amc(ic,jm)*sompmy11(ic,ig) 1663 zpf= zpf + vpf*vc(ic,ic,ip) 1664 1665 if (ic.gt.1) then 1666 do jc=1,ic-1 1667 vpf= xmupmc(ic,jm)*(apc(jc,ip)*sompp11(ig) & 1668 +amc(jc,jm)*sompm11(ig)) & 1669 +xmupmc(jc,jm)*(apc(ic,ip)*sompp11(ig) & 1670 +amc(ic,jm)*sompm11(ig)) & 1671 -apc(jc,ip)*somppy11(ic,ig) & 1672 -amc(jc,jm)*sompmy11(ic,ig) & 1673 -apc(ic,ip)*somppy11(jc,ig) & 1674 -amc(ic,jm)*sompmy11(jc,ig) & 1675 +apc(ic,ip)*apc(jc,ip)*carpp11(ig) & 1676 +amc(ic,jm)*amc(jc,jm)*carpm11(ig) & 1677 +(apc(ic,ip)*amc(jc,jm)+amc(ic,jm)*apc(jc,ip)) & 1678 *somppm11(ig) 1679 zpf= zpf + 2.d0*vpf*vc(ic,jc,ip) 1680 ! print *,'zpf:',zpf,ig,ic,jc 1681 end do !! jc 1682 end if 1683 end do !! ic 1684 ! print *,'=>',ip,jm,zpf 1685 zpf=-0.5d0*zpf 1686 ! print *,'-0.5*zpf:',zpf 1687 zpf=dexp(zpf) 1688 ! print *,'dexp(zpf):',zpf 1689 if ( zpf == 0 ) zpf = 1/huge(zpf) 1690 1691 if ((zpf .ne. zpf) ) then 1692 fm2=INIFINY_REAL_VALUE 1693 fp2=INIFINY_REAL_VALUE 1694 f=INIFINY_REAL_VALUE 1695 ! print *,"FIN AVANT..." 1696 ! stop 1697 return 1698 end if 1699 1700 1701 vmere=vmere+probg(current_chr,ig)*zpf 1702 ! print *,'vmere:',vmere 1703 1704 end do!! ig 1705 1706 1707 if ( vmere /= 0) then 1708 fm2(jm)=-dlog(vmere)+0.5d0*dble(eff(jm))*dlog(det(ip)) 1709 ! print *,'fm2:',fm2(jm),-dlog(vmere) 1710 ! stop 1711 else 1712 fm2=INIFINY_REAL_VALUE 1713 fp2=INIFINY_REAL_VALUE 1714 f=INIFINY_REAL_VALUE 1715 return 1716 end if 1717 1718 end if 1719 ! 1720 else 1721 fm2(jm)=0.d0 1722 end if 1723 f=f+fm2(jm) 1724 end do 1725 1726 end subroutine funct_mcar_1qtl_family
init_analyse_multitrait
[ Top ] [ m_qtlmap_analyse_multitrait ] [ Subroutines ]
NAME
init_analyse_multitrait
DESCRIPTION
allocation of buffer/solution arrays
SOURCE
439 subroutine init_analyse_multitrait 440 integer :: stat 441 442 allocate (std(ncar,np),STAT=stat) 443 call check_allocate(stat,'xmoyp [m_qtlmap_analyse_multitrait]') 444 allocate (ap(ncar,np),STAT=stat) 445 call check_allocate(stat,'xmoyp [m_qtlmap_analyse_multitrait]') 446 allocate (xmoyp(ncar,np),STAT=stat) 447 call check_allocate(stat,'xmoyp [m_qtlmap_analyse_multitrait]') 448 allocate (am(ncar,nm),STAT=stat) 449 call check_allocate(stat,'xmoyp [m_qtlmap_analyse_multitrait]') 450 allocate (xmoym(ncar,nm),STAT=stat) 451 call check_allocate(stat,'xmoyp [m_qtlmap_analyse_multitrait]') 452 allocate (somydf01(ncar,np),STAT=stat) 453 call check_allocate(stat,'somydf01 [m_qtlmap_analyse_multitrait]') 454 allocate (carydf01(ncar,np),STAT=stat) 455 call check_allocate(stat,'carydf01 [m_qtlmap_analyse_multitrait]') 456 allocate (somyydf01(ncar,ncar,np),STAT=stat) 457 call check_allocate(stat,'somyydf01 [m_qtlmap_analyse_multitrait]') 458 allocate (sig01(ncar,np),STAT=stat) 459 call check_allocate(stat,'sig01 [m_qtlmap_analyse_multitrait]') 460 allocate (xmu01p(ncar,np),STAT=stat) 461 call check_allocate(stat,'xmu01p [m_qtlmap_analyse_multitrait]') 462 allocate (xmu01m(ncar,nm),STAT=stat) 463 call check_allocate(stat,'xmu01m [m_qtlmap_analyse_multitrait]') 464 allocate (rhoi(ncar,ncar),STAT=stat) 465 call check_allocate(stat,'rhoi [m_qtlmap_analyse_multitrait]') 466 467 allocate (fp21(np),STAT=stat) 468 call check_allocate(stat,'fp21 [m_qtlmap_analyse_multitrait]') 469 allocate (fm21(nm),STAT=stat) 470 call check_allocate(stat,'fm21 [m_qtlmap_analyse_multitrait]') 471 allocate (somy11(ncar,nm),STAT=stat) 472 call check_allocate(stat,'somy11 [m_qtlmap_analyse_multitrait]') 473 allocate (cary11(ncar,nm),STAT=stat) 474 call check_allocate(stat,'cary11 [m_qtlmap_analyse_multitrait]') 475 allocate (somyy11(ncar,ncar,nm),STAT=stat) 476 call check_allocate(stat,'somyy11 [m_qtlmap_analyse_multitrait]') 477 allocate (xmu2p(ncar,np),STAT=stat) 478 call check_allocate(stat,'xmu2p [m_qtlmap_analyse_multitrait]') 479 allocate (sig2(ncar,np),STAT=stat) 480 call check_allocate(stat,'sig2 [m_qtlmap_analyse_multitrait]') 481 allocate (xmu2m(ncar,nm),STAT=stat) 482 call check_allocate(stat,'xmu2m [m_qtlmap_analyse_multitrait]') 483 xmu2m=0.d0 484 allocate (rhoi2(ncar,ncar),STAT=stat) 485 call check_allocate(stat,'rhoi2 [m_qtlmap_analyse_multitrait]') 486 rhoi2=0.d0 487 allocate (rhoi2_2(ncar,ncar)) 488 rhoi2_2=0.d0 489 allocate (rhoi2_1(ncar,ncar)) 490 rhoi2_1=0.d0 491 492 allocate (effdf(np)) 493 allocate (estmum(np)) 494 allocate (eff(nm)) 495 allocate (effp(np)) 496 497 allocate (corcd(ncar,ncar,nd)) 498 corcd=1.d0 499 500 end subroutine init_analyse_multitrait
init_sub_1qtl
[ Top ] [ m_qtlmap_analyse_multitrait ] [ Subroutines ]
NAME
init_sub_1qtl
DESCRIPTION
allocation of buffer/solution arrays
SOURCE
509 subroutine init_sub_1qtl 510 integer :: stat,ng 511 512 ng = get_maxnbgenotypedam() 513 514 allocate (sompp11(ng),STAT=stat) 515 call check_allocate(stat,'sompp11 [m_qtlmap_analyse_multitrait]') 516 allocate (sompm11(ng),STAT=stat) 517 call check_allocate(stat,'sompm11 [m_qtlmap_analyse_multitrait]') 518 allocate (carpp11(ng),STAT=stat) 519 call check_allocate(stat,'carpp11 [m_qtlmap_analyse_multitrait]') 520 allocate (carpm11(ng),STAT=stat) 521 call check_allocate(stat,'carpm11 [m_qtlmap_analyse_multitrait]') 522 allocate (somppm11(ng),STAT=stat) 523 call check_allocate(stat,'somppm11 [m_qtlmap_analyse_multitrait]') 524 allocate (sompmy11(ncar,ng),STAT=stat) 525 call check_allocate(stat,'sompmy11 [m_qtlmap_analyse_multitrait]') 526 allocate (somppy11(ncar,ng),STAT=stat) 527 call check_allocate(stat,'somppy11 [m_qtlmap_analyse_multitrait]') 528 allocate (somppdf11(ng),STAT=stat) 529 call check_allocate(stat,'somppdf11 [m_qtlmap_analyse_multitrait]') 530 allocate (carppdf11(ng),STAT=stat) 531 call check_allocate(stat,'carppdf11 [m_qtlmap_analyse_multitrait]') 532 allocate (somppydf11(ncar,ng),STAT=stat) 533 call check_allocate(stat,'somppydf11 [m_qtlmap_analyse_multitrait]') 534 535 allocate (fp2(np),STAT=stat) 536 call check_allocate(stat,'fp2 [m_qtlmap_analyse_multitrait]') 537 allocate (fm2(nm),STAT=stat) 538 call check_allocate(stat,'fm2 [m_qtlmap_analyse_multitrait]') 539 540 end subroutine init_sub_1qtl
opti_mcar_0qtl
[ Top ] [ m_qtlmap_analyse_multitrait ] [ Subroutines ]
NAME
opti_mcar_0qtl
DESCRIPTION
Calcul de la statistique de test multicaractere le long du chromosome
SOURCE
786 subroutine opti_mcar_0qtl 787 ! Divers 788 real (kind=dp) , dimension(:), allocatable :: par,borni,borns 789 integer iuser(1) 790 real (kind=dp)user(1) 791 792 integer :: ibound,npar,nrho,ip,jm,j,k,nest,nm1,nm2 793 integer :: r1,imumest,i,nestim,ifail,ic 794 real (kind=dp) ::somxmu 795 logical , dimension(:,:),pointer :: filter_inc 796 797 ! 798 !**************************************************************************** 799 ! 800 ! Parametres de maximisation 801 nrho=ncar*(ncar-1)/2 802 npar=(2*ncar*np)+(ncar*nmumest(1))+nrho 803 allocate (par(npar)) 804 allocate (borni(npar)) 805 allocate (borns(npar)) 806 allocate (filter_inc(np,npar)) 807 808 ibound=0 809 filter_inc=.false. 810 nestim=0 811 do ip=1,np 812 do j=1,ncar 813 borni((j-1)*np+ip)=SIG_MIN 814 borns((j-1)*np+ip)=SIG_MAX 815 filter_inc(ip,(j-1)*np+ip)=.true. 816 borni(ip+ncar*np+(j-1)*np)=XMU_MIN 817 borns(ip+ncar*np+(j-1)*np)=XMU_MAX 818 filter_inc(ip,ip+ncar*np+(j-1)*np)=.true. 819 end do 820 821 nest=0 822 do jm=nmp(ip)+1,nmp(ip+1) 823 if(estime(1,jm)) then 824 nest=nest+1 825 if(nest.le.estmum(ip)) then 826 nestim=nestim+1 827 do j=1,ncar 828 filter_inc(ip,2*ncar*np+(j-1)*nmumest(1)+nestim)=.true. 829 end do 830 end if 831 end if 832 end do 833 834 end do 835 do jm=1,nmumest(1) 836 do j=1,ncar 837 borni(2*ncar*np+(j-1)*nmumest(1)+jm)=XMU_MIN 838 borns(2*ncar*np+(j-1)*nmumest(1)+jm)=XMU_MAX 839 end do 840 end do 841 do j=1,nrho 842 borni(2*ncar*np+ncar*nmumest(1)+j)=DEFAULT_PARAM_MIN 843 borns(2*ncar*np+ncar*nmumest(1)+j)=DEFAULT_PARAM_MAX 844 filter_inc(:,2*ncar*np+ncar*nmumest(1)+j)=.true. 845 end do 846 ! 847 ! Point de depart 848 k=0 849 imumest=0 850 do i=1,ncar 851 nestim=0 852 do ip=1,np 853 par((i-1)*np+ip)=sig01(i,ip) 854 par(ncar*np+(i-1)*np+ip)=xmu01p(i,ip) 855 nm1=nmp(ip)+1 856 nm2=nmp(ip+1) 857 nest=0 858 somxmu=0.d0 859 do jm=nm1,nm2 860 if (estime(i,jm)) then 861 nest=nest+1 862 if (nest.le.estmum(ip)) then 863 nestim=nestim+1 864 par(2*ncar*np+(i-1)*nmumest+nestim)=xmu01m(i,nestim) 865 end if 866 end if 867 end do 868 end do 869 if (i.gt.1) then 870 do j=1,i-1 871 k=k+1 872 r1=rhoi(j,i) 873 par(2*ncar*np+ncar*nmumest(1)+k)=dlog((1.d0+r1)/(1.d0-r1)) 874 end do 875 end if 876 end do 877 ! 878 ! Optimisation de la vraisemblance 879 ifail=1 880 ! call minimizing_funct(npar,ibound,funct_mcar_0qtl,borni,borns,par,f01,iuser,user,ifail) 881 call minimizing_funct_family_sire(npar,ibound,funct_mcar_0qtl_family,filter_inc,fp21,borni,borns,par,f01,iuser,user,ifail) 882 883 k=0 884 do i=2,ncar 885 do j=1,i-1 886 k=k+1 887 rhoi2(j,i)=par(2*ncar*np+ncar*nmumest(1)+k) 888 !**ajout ofi pour affichage sous H0 889 rhoi2_1(j,i)=(dexp(rhoi2(j,i))-1.d0)/(dexp(rhoi2(j,i))+1.d0) 890 rhoi2_1(i,j)=rhoi2_1(j,i) 891 !**fin ajout 892 rhoi2(i,j)=rhoi2(j,i) 893 end do 894 end do 895 896 do i=1,ncar 897 do ip=1,np 898 sig2(i,ip)=par((i-1)*np+ip) 899 xmu2p(i,ip)=par(ncar*np+(i-1)*np+ip) 900 end do 901 do jm=1,nmumest(1) 902 xmu2m(i,jm)=par(2*ncar*np+(i-1)*nmumest(1)+jm) 903 end do 904 end do 905 906 deallocate (par) 907 deallocate (borni) 908 deallocate (borns) 909 deallocate (filter_inc) 910 end subroutine opti_mcar_0qtl
opti_mcar_1qtl
[ Top ] [ m_qtlmap_analyse_multitrait ] [ Subroutines ]
NAME
opti_mcar_1qtl
DESCRIPTION
Calcul de la statistique de test multicaractere le long du chromosome
SOURCE
1128 subroutine opti_mcar_1qtl(lrtsol) 1129 1130 type(TYPE_LRT_SOLUTION) , intent(inout) :: lrtsol 1131 1132 real (kind=dp) :: r(ncar,ncar),rhof(ncar,ncar) 1133 1134 1135 ! 1136 ! Tableaux dimensionnes selon nm, le nombre de meres 1137 1138 ! Divers 1139 real (kind=dp), dimension((3*ncar*np)+ncar*nmumest(1)+ncar*namest(1)+ncar*(ncar-1)/2) :: par,borni,borns 1140 1141 integer iuser(1) 1142 double precision user(1) 1143 integer :: nrho,npar,ifail,ibound,i,k,ip,chr,ntotal,indexchr(nchr),nprime 1144 integer :: jm,nm1,nm2,nd1,nd2,kd,kkd,ifem,j,ilong,n,ix 1145 integer :: ngeno1,ngeno2,ig,nestim,nest,indam 1146 1147 real (kind=dp) , dimension(:,:,:),pointer :: listepar 1148 real (kind=dp) :: somxmu,pp,pm,xlrtc_t,f2,f2_2 1149 logical , dimension(:,:),pointer :: filter_inc 1150 1151 call new(1,lrtsol) 1152 1153 !**************************************************************************** 1154 ! Calcul de la vraisemblance sous H2 1155 ! Parametres de maximisation 1156 nrho=ncar*(ncar-1)/2 1157 npar=(3*ncar*np)+ncar*nmumest(1)+ncar*namest(1)+nrho 1158 1159 allocate (filter_inc(np,npar)) 1160 1161 par=0.d0 1162 filter_inc=.false. 1163 1164 ibound=0 1165 rhof = 0.d0 1166 k=0 1167 do i=1,ncar 1168 nestim=0 1169 do ip=1,np 1170 borni((i-1)*np+ip)=SIG_MIN ! variance de la famille ip pour le caractere i 1171 borns((i-1)*np+ip)=SIG_MAX 1172 filter_inc(ip,(i-1)*np+ip)=.true. 1173 borni(ip+ncar*np+(i-1)*np)=DEFAULT_PARAM_MIN ! moyenne de la famille ip pour le caractere i 1174 borns(ip+ncar*np+(i-1)*np)=DEFAULT_PARAM_MAX 1175 filter_inc(ip,ip+ncar*np+(i-1)*np)=.true. 1176 borni(ip+2*ncar*np+(i-1)*np)=AM_MIN ! effet qtl de la famille ip pour le caractere i 1177 borns(ip+2*ncar*np+(i-1)*np)= AM_MAX 1178 filter_inc(ip,ip+2*ncar*np+(i-1)*np)=.true. 1179 1180 nest=0 1181 do jm=nmp(ip)+1,nmp(ip+1) 1182 if(estime(1,jm)) then 1183 nest=nest+1 1184 if(nest.le.estmum(ip)) then 1185 nestim=nestim+1 1186 filter_inc(ip,3*ncar*np+(i-1)*nmumest(1)+nestim)=.true. 1187 end if 1188 ifem=repfem(jm) 1189 filter_inc(ip,3*ncar*np+ncar*nmumest(1)+(i-1)*namest(1)+iam(1,ifem))=.true. 1190 end if 1191 end do 1192 end do ! ip 1193 do jm=1,nmumest(1) 1194 borni(3*ncar*np+(i-1)*nmumest(1)+jm)=DEFAULT_PARAM_MIN ! moyenne de la famille ip/jm pour le caractere i 1195 borns(3*ncar*np+(i-1)*nmumest(1)+jm)=DEFAULT_PARAM_MAX 1196 end do 1197 do jm=1,namest(1) 1198 borni(3*ncar*np+ncar*nmumest(1)+(i-1)*namest(1)+jm)=AM_MIN ! effet qtl de la famille ip/ifem pour le caractere i 1199 borns(3*ncar*np+ncar*nmumest(1)+(i-1)*namest(1)+jm)=AM_MAX 1200 end do 1201 if (i.gt.1) then 1202 do j=1,i-1 1203 k=k+1 1204 borni(3*ncar*np+ncar*nmumest(1)+ncar*namest(1)+k)=DEFAULT_PARAM_MIN 1205 borns(3*ncar*np+ncar*nmumest(1)+ncar*namest(1)+k)=DEFAULT_PARAM_MAX 1206 end do 1207 end if 1208 end do 1209 1210 filter_inc(:,3*ncar*np+ncar*nmumest(1)+ncar*namest(1)+1:)=.true. 1211 ! 1212 ! 1213 ! Point de depart 1214 k=0 1215 do i=1,ncar 1216 do ip=1,np 1217 par((i-1)*np+ip)=sig2(i,ip) 1218 par(ip+ncar*np+(i-1)*np)=xmu2p(i,ip) 1219 par(ip+2*ncar*np+(i-1)*np)=0.d0 1220 end do 1221 do jm=1,nmumest(1) 1222 par(3*ncar*np+(i-1)*nmumest+jm)=xmu2m(i,jm) 1223 end do 1224 do jm=1,namest(1) 1225 par(3*ncar*np+ncar*nmumest+(i-1)*namest+jm)= 0.d0 1226 end do 1227 if (i.gt.1) then 1228 do j=1,i-1 1229 k=k+1 1230 par(3*ncar*np+ncar*nmumest+ncar*namest+k)=rhoi2(j,i) 1231 end do 1232 end if 1233 end do 1234 1235 1236 allocate (listepar(nchr,get_maxnpo(),npar)) 1237 ibound=0 1238 ! Marche le long du chromosome 1239 lrtsol%lrtmax=-1.d75 1240 1241 ntotal=0 1242 do chr=1,nchr 1243 ntotal=ntotal+get_npo(chr) 1244 indexchr(chr)=ntotal 1245 end do 1246 1247 !$OMP PARALLEL DEFAULT(SHARED) FIRSTPRIVATE(par) & 1248 !$OMP PRIVATE(ip,jm,chr,i,n,nd1,nd2,pp,pm,kd,kkd,ig,ifail,nestim,f2) 1249 call init_sub_1qtl 1250 !$OMP DO 1251 do nprime=1,ntotal 1252 n=nprime 1253 do chr=nchr,1,-1 1254 if ( indexchr(chr) >= nprime) then 1255 current_chr = chr 1256 else 1257 n=n-get_npo(chr) 1258 end if 1259 end do 1260 chr=current_chr 1261 1262 do ip=1,np 1263 somppdf11(ip)=0.d0 1264 carppdf11(ip)=0.d0 1265 do i=1,ncar 1266 somppydf11(i,ip)=0.d0 1267 end do 1268 do jm=nmp(ip)+1,nmp(ip+1) 1269 do ig=ngenom(current_chr,jm)+1,ngenom(current_chr,jm+1) 1270 sompp11(ig)=0.d0 1271 carpp11(ig)=0.d0 1272 do i=1,ncar 1273 somppy11(i,ig)=0.d0 1274 sompmy11(i,ig)=0.d0 1275 end do 1276 sompm11(ig)=0.d0 1277 carpm11(ig)=0.d0 1278 somppm11(ig)=0.d0 1279 nd1=ngend(current_chr,ig)+1 1280 nd2=ngend(current_chr,ig+1) 1281 do kd=nd1,nd2 1282 kkd=ndesc(current_chr,kd) 1283 if(count(presentc(:,kkd)) == ncar )then 1284 pp=-pdd(current_chr,kd,1,n)-pdd(current_chr,kd,2,n)+pdd(current_chr,kd,3,n)+pdd(current_chr,kd,4,n) 1285 pm=-pdd(current_chr,kd,1,n)+pdd(current_chr,kd,2,n)-pdd(current_chr,kd,3,n)+pdd(current_chr,kd,4,n) 1286 if(estime(ncar,jm)) then 1287 sompp11(ig)=sompp11(ig)+pp 1288 carpp11(ig)=carpp11(ig)+pp*pp 1289 do i=1,ncar 1290 somppy11(i,ig)=somppy11(i,ig)+pp*y(i,kkd) 1291 sompmy11(i,ig)=sompmy11(i,ig)+pm*y(i,kkd) 1292 end do 1293 sompm11(ig)=sompm11(ig)+pm 1294 carpm11(ig)=carpm11(ig)+pm*pm 1295 somppm11(ig)=somppm11(ig)+pp*pm 1296 else 1297 somppdf11(ip)=somppdf11(ip)+pp 1298 carppdf11(ip)=carppdf11(ip)+pp*pp 1299 do i=1,ncar 1300 somppydf11(i,ip)=somppydf11(i,ip)+pp*y(i,kkd) 1301 end do 1302 end if 1303 end if 1304 end do !! fin kd 1305 end do !! fin ig 1306 end do !! fin jm 1307 end do !! fin ip 1308 1309 ! Optimisation de la vraisemblance a la position dx 1310 ifail=1 1311 ! call minimizing_funct(npar,ibound,funct_mcar_1qtl,borni,borns,par,f2,iuser,user,ifail) 1312 call minimizing_funct_family_sire(npar,ibound,funct_mcar_1qtl_family,filter_inc,fp2,borni,borns,par,f2,iuser,user,ifail) 1313 1314 listepar(chr,n,:npar) = par(:npar) 1315 if ( f2 < INIFINY_REAL_VALUE ) then 1316 lrtsol%lrt1(chr,n)=-2.d0*(f2-f01) 1317 else 1318 lrtsol%lrt1(chr,n)=0 1319 end if 1320 1321 ! seulement le dernier effet qtl (du dernier caracteres)........ 1322 ! A corriger... 1323 do i=1,ncar 1324 do ip=1,np 1325 lrtsol%xlrp(chr,ip,n)=-2.d0*(fp2(ip)-fp21(ip)) 1326 lrtsol%pater_eff(chr,ip,n)=par(ip+2*ncar*np+(i-1)*np) 1327 end do 1328 do jm=1,nm 1329 lrtsol%xlrm(chr,jm,n)=-2.d0*(fm2(jm)-fm21(jm)) 1330 end do 1331 do jm=1,namest(1) 1332 lrtsol%mater_eff(chr,jm,n)=par(3*ncar*np+ncar*nmumest(1)+(i-1)*namest(1)+jm) 1333 end do 1334 end do 1335 end do 1336 !$OMP END DO 1337 call end_sub_1qtl 1338 !$OMP END PARALLEL 1339 1340 !recherche du maximum 1341 do chr=1,nchr 1342 do n=1,get_npo(chr) 1343 if(lrtsol%lrtmax(0) < lrtsol%lrt1(chr,n)) then 1344 lrtsol%lrtmax(0)=lrtsol%lrt1(chr,n) 1345 lrtsol%nxmax(0)=n 1346 lrtsol%chrmax(0)=chr 1347 end if 1348 end do 1349 end do 1350 !initilisation des valeurs par rapport au maximum trouve 1351 par(:npar)=listepar(lrtsol%chrmax(0),lrtsol%nxmax(0),:npar) 1352 k=0 1353 rhof=0.d0 1354 do i=2,ncar 1355 do j=1,i-1 1356 k=k+1 1357 r(j,i)=par(3*ncar*np+ncar*nmumest(1)+ncar*namest(1)+k) 1358 rhof(j,i)=(dexp(r(j,i))-1.d0)/(dexp(r(j,i))+1.d0) 1359 rhof(i,j)=rhof(j,i) 1360 end do 1361 end do 1362 do i=1,ncar 1363 do j=1,ncar 1364 rhoi2_2(i,j)=rhof(i,j) 1365 end do 1366 nestim=0 1367 do ip=1,np 1368 std(i,ip)=par(ip+(i-1)*np) 1369 xmoyp(i,ip)=par(ip+ncar*np+(i-1)*np) 1370 ap(i,ip)=par(ip+2*ncar*np+(i-1)*np) 1371 nest=0 1372 somxmu = 0.d0 1373 do jm=nmp(ip)+1,nmp(ip+1) 1374 if (count(estime(:,jm))==ncar) then 1375 nest=nest+1 1376 ifem=repfem(jm) 1377 indam=iam(1,ifem) 1378 am(i,jm)=par(3*ncar*np+ncar*nmumest(1)+(i-1)*namest(1)+indam) 1379 if (nest.le.estmum(ip)) then 1380 nestim=nestim+1 1381 xmoym(i,jm)=par(3*ncar*np+(i-1)*nmumest(1)+nestim) 1382 somxmu=somxmu+xmoym(i,jm) 1383 else 1384 xmoym(i,jm)=-somxmu 1385 end if 1386 else 1387 am(i,jm)=0.d0 1388 xmoym(i,jm)=0.d0 1389 end if 1390 end do !jm 1391 end do !ip 1392 end do !icar 1393 1394 deallocate (listepar) 1395 deallocate (filter_inc) 1396 1397 end subroutine opti_mcar_1qtl
optinit_mcar
[ Top ] [ m_qtlmap_analyse_multitrait ] [ Subroutines ]
NAME
optinit_mcar
DESCRIPTION
Initialise les ecart types et moyennes intra-famille
SOURCE
614 subroutine optinit_mcar 615 integer :: efft(ncar),i 616 integer :: ic,nm1,nm2,jm,nest,ifem,imumest,jc,ip,nd1,nd2,kd 617 real (kind=dp) :: somt(ncar),cov(ncar,ncar),xmut(ncar),sigt(ncar) 618 619 real (kind=dp) :: somy(nm),cary(nm),somyp 620 621 cov=0.d0 622 somyydf01=0.d0 623 somyy11=0.d0 !! que moitie inf 624 ! 625 do ic=1,ncar 626 627 628 effp=0 629 imumest=0 630 estfem=.false. 631 do ip=1,np 632 somyp=0.d0 633 sig01(ic,ip)=0.d0 634 effdf(ip)=0 635 somydf01(ic,ip)=0.d0 636 carydf01(ic,ip)=0.d0 637 estmum(ip)=0 638 nm1=nmp(ip)+1 639 nm2=nmp(ip+1) 640 do jm=nm1,nm2 641 eff(jm)=0 642 somy(jm)=0.d0 643 cary(jm)=0.d0 644 nd1=ndm(jm)+1 645 nd2=ndm(jm+1) 646 ifem=repfem(jm) 647 do kd=nd1,nd2 648 if(count(presentc(:,kd)) == ncar ) then 649 eff(jm)=eff(jm)+1 650 somy(jm)=somy(jm)+y(ic,kd)*cd(ic,kd) 651 cary(jm)=cary(jm)+(y(ic,kd)*y(ic,kd))*cd(ic,kd) 652 end if 653 end do 654 655 effp(ip)=effp(ip)+eff(jm) 656 somyp=somyp+somy(jm) 657 if( estime(ic,jm) ) then 658 imumest=imumest+1 659 estmum(ip)=estmum(ip)+1 660 xmu01m(ic,imumest)=somy(jm)/dble(eff(jm)) 661 estfem(ifem)=.true. 662 else 663 effdf(ip)=effdf(ip)+eff(jm) 664 somydf01(ic,ip)=somydf01(ic,ip)+somy(jm) 665 carydf01(ic,ip)=carydf01(ic,ip)+cary(jm) 666 end if 667 end do 668 if (effp(ip) == 0.d0) then 669 call stop_application('Father ['//trim(pere(ip))//'] has got no child with trait value') 670 else 671 xmu01p(ic,ip)=somyp/dble(effp(ip)) 672 end if 673 do jm=nm1,nm2 674 nd1=ndm(jm)+1 675 nd2=ndm(jm+1) 676 do kd=nd1,nd2 677 if(count(presentc(:,kd)) == ncar ) then 678 sig01(ic,ip)=sig01(ic,ip)+(y(ic,kd)-xmu01p(ic,ip))*(y(ic,kd)-xmu01p(ic,ip)) 679 end if 680 end do 681 end do 682 if ( effp(ip) == 1) then 683 call stop_application('Father ['//trim(pere(ip))//'] has got only one child with trait value') 684 else 685 sig01(ic,ip)=dsqrt(sig01(ic,ip)/dble(effp(ip)-1)) 686 end if 687 688 if(estmum(ip).gt.0) then 689 estmum(ip)=estmum(ip)-1 690 nmumest(ic)=nmumest(ic)-1 691 end if 692 end do 693 694 ! Initialisation des param�tres dans des vecteurs d�pendants du nombre de caracteres 695 efft(ic)=0 696 somt(ic)=0.d0 697 xmut(ic)=0.d0 698 sigt(ic)=0.d0 699 imumest=0 700 do ip=1,np 701 nm1=nmp(ip)+1 702 nm2=nmp(ip+1) 703 nest=0 704 do jm=nm1,nm2 705 706 efft(ic)=efft(ic)+eff(jm) 707 somt(ic)=somt(ic)+somy(jm) 708 709 somy11(ic,jm)=somy(jm) 710 cary11(ic,jm)=cary(jm) 711 712 if(estime(ic,jm)) then 713 nest=nest+1 714 if (nest.le.estmum(ip)) then 715 imumest=imumest+1 716 end if 717 end if 718 end do !! fin jm 719 end do !! fin ip 720 721 xmut(ic)=somt(ic)/dble(efft(ic)) 722 723 do ip=1,np 724 nm1=nmp(ip)+1 725 nm2=nmp(ip+1) 726 do jm=nm1,nm2 727 nd1=ndm(jm)+1 728 nd2=ndm(jm+1) 729 do kd=nd1,nd2 730 if(count(presentc(:,kd)) == ncar ) then 731 sigt(ic)=sigt(ic)+(y(ic,kd)-xmut(ic))*(y(ic,kd)-xmut(ic)) 732 if (ic.gt.1) then 733 do jc=1,ic-1 734 if(count(presentc(:,kd)) == ncar ) then 735 somyy11(jc,ic,jm)=somyy11(jc,ic,jm)+(y(jc,kd)*y(ic,kd)) !! que moitie inf 736 somyy11(ic,jc,jm)=somyy11(jc,ic,jm) 737 cov(jc,ic)=cov(jc,ic)+(y(jc,kd)-xmut(jc))*(y(ic,kd)-xmut(ic)) 738 end if 739 end do 740 end if 741 end if 742 end do 743 if(.not.estime(ic,jm)) then 744 if (ic.gt.1) then 745 do jc=1,ic-1 746 somyydf01(jc,ic,ip)=somyydf01(jc,ic,ip)+somyy11(jc,ic,jm) 747 somyydf01(ic,jc,ip)=somyydf01(jc,ic,ip) 748 end do 749 end if 750 end if 751 end do 752 end do 753 sigt(ic)=sigt(ic)/dble(efft(ic)) 754 sigt(ic)=dsqrt(sigt(ic)) 755 756 if(ic.gt.1) then 757 do jc=1,ic-1 758 if (efft(jc).ne.efft(ic)) then 759 do i=1,nd 760 if (count(presentc(:,i))<ncar) then 761 call log_mess("*Remove animal :"//trim(animal(i))) 762 end if 763 end do 764 call log_mess('Missing trait data for traits '//trim(carac(ic))// & 765 ' and '//trim(carac(jc))//', not possible to perform multiple trait analysis',ERROR_DEF) 766 call stop_application(" ** ") 767 end if 768 cov(jc,ic)=cov(jc,ic)/dble(efft(ic)) 769 cov(ic,jc)=cov(jc,ic) 770 rhoi(jc,ic)=cov(jc,ic)/(sigt(jc)*sigt(ic)) 771 rhoi(ic,jc)=rhoi(jc,ic) 772 end do 773 end if 774 end do 775 !! fin ic 776 return 777 end subroutine optinit_mcar
set_solution_hypothesis0
[ Top ] [ m_qtlmap_analyse_multitrait ] [ Subroutines ]
NAME
set_solution_hypothesis0
DESCRIPTION
SOURCE
1735 subroutine set_solution_hypothesis0(incsol) 1736 type(TYPE_INCIDENCE_SOLUTION) ,intent(inout) :: incsol 1737 1738 integer :: nteff,maxNbPar,jm,ifem,ip,ieff,ic 1739 1740 incsol%hypothesis=0 1741 allocate (incsol%sig(ncar,np)) 1742 ! Mean family 1743 nteff = ncar 1744 maxNbPar = np 1745 do ic=1,ncar 1746 if ( count(estime(ic,:)) > 0 ) then 1747 nteff = nteff + 1 1748 maxNbPar = max(maxNbPar,count(estime(ic,:))) 1749 end if 1750 end do 1751 1752 allocate (incsol%groupeName(nteff)) 1753 allocate (incsol%nbParameterGroup(nteff)) 1754 allocate (incsol%parameterName(nteff,maxNbPar)) 1755 allocate (incsol%paramaterValue(nteff,maxNbPar)) 1756 allocate (incsol%parameterVecsol(nteff,maxNbPar)) 1757 allocate (incsol%parameterPrecis(nteff,maxNbPar)) 1758 1759 incsol%parameterVecsol=.true. 1760 incsol%parameterPrecis=0.d0 1761 1762 ieff=0 1763 do ic=1,ncar 1764 do ip=1,np 1765 incsol%sig(ic,ip) = sig2(ic,ip)*sigt(ic) 1766 end do 1767 1768 ieff=ieff+1 1769 incsol%groupeName(ieff) = 'Mean Sire ['//trim(carac(ic))//"]" 1770 incsol%nbParameterGroup(ieff)=np 1771 1772 do ip=1,np 1773 incsol%parameterName(ieff,ip)='Sire '//trim(pere(ip)) 1774 incsol%paramaterValue(ieff,ip) = xmu2p(ic,ip)*sigt(ic) + xmut(ic) 1775 end do 1776 1777 if ( count(estime(ic,:)) > 0 ) then 1778 ieff = ieff +1 1779 incsol%groupeName(ieff) = 'Mean dam ['//trim(carac(ic))//"]" 1780 incsol%nbParameterGroup(ieff)= count(estime(ic,:)) 1781 ifem=0 1782 do ip=1,np 1783 do jm=nmp(ip)+1,nmp(ip+1) 1784 if ( estime(ic,jm)) then 1785 ifem=ifem+1 1786 incsol%parameterName(ieff,ifem)='Dam '//trim(mere(jm))//" [Sire "//trim(pere(ip))//"]" 1787 incsol%paramaterValue(ieff,ifem) = xmu2m(ic,ifem)*sigt(ic) 1788 end if 1789 end do 1790 end do 1791 end if 1792 1793 end do 1794 1795 allocate (incsol%rhoi(ncar,ncar)) 1796 incsol%rhoi=rhoi2_1 1797 1798 end subroutine set_solution_hypothesis0
set_solution_hypothesis1
[ Top ] [ m_qtlmap_analyse_multitrait ] [ Subroutines ]
NAME
set_solution_hypothesis1
DESCRIPTION
SOURCE
1809 subroutine set_solution_hypothesis1(incsol) 1810 type(TYPE_INCIDENCE_SOLUTION) ,intent(inout) :: incsol 1811 1812 integer :: nteff,maxNbPar,jm,ifem,ip,ieff,ic 1813 1814 1815 incsol%hypothesis=1 1816 allocate (incsol%sig(ncar,np)) 1817 ! Mean family , Qtl effect 1818 nteff = 2 * ncar 1819 maxNbPar = np 1820 do ic=1,ncar 1821 if ( count(estime(ic,:)) > 0 ) then 1822 nteff = nteff + 2 1823 maxNbPar = max(maxNbPar,count(estime(ic,:))) 1824 end if 1825 end do 1826 1827 allocate (incsol%groupeName(nteff)) 1828 allocate (incsol%nbParameterGroup(nteff)) 1829 allocate (incsol%parameterName(nteff,maxNbPar)) 1830 allocate (incsol%paramaterValue(nteff,maxNbPar)) 1831 allocate (incsol%parameterVecsol(nteff,maxNbPar)) 1832 allocate (incsol%parameterPrecis(nteff,maxNbPar)) 1833 allocate (incsol%qtl_groupeName(ncar,1)) ! 1 qtl 1834 1835 incsol%parameterVecsol=.true. 1836 incsol%parameterPrecis=0.d0 1837 1838 ieff=0 1839 do ic=1,ncar 1840 1841 do ip=1,np 1842 incsol%sig(ic,ip) = std(ic,ip)*sigt(ic) 1843 end do 1844 1845 ieff=ieff+1 1846 incsol%groupeName(ieff) = 'Mean Sire ['//trim(carac(ic))//"]" 1847 incsol%nbParameterGroup(ieff)=np 1848 1849 do ip=1,np 1850 incsol%parameterName(ieff,ip)='Sire '//trim(pere(ip)) 1851 incsol%paramaterValue(ieff,ip) = xmoyp(ic,ip)*sigt(ic) + xmut(ic) 1852 end do 1853 1854 if ( count(estime(ic,:)) > 0 ) then 1855 ieff = ieff +1 1856 incsol%groupeName(ieff) = 'Mean dam ['//trim(carac(ic))//"]" 1857 incsol%nbParameterGroup(ieff)= count(estime(ic,:)) 1858 ifem=0 1859 do ip=1,np 1860 do jm=nmp(ip)+1,nmp(ip+1) 1861 if ( estime(ic,jm)) then 1862 ifem=ifem+1 1863 incsol%parameterName(ieff,ifem)='Dam '//trim(mere(jm)) 1864 incsol%paramaterValue(ieff,ifem) = xmoym(ic,jm)*sigt(ic)+ xmut(ic) 1865 end if 1866 end do 1867 end do 1868 end if 1869 1870 ieff = ieff +1 1871 incsol%qtl_groupeName(ic,1) = ieff 1872 incsol%groupeName(ieff) = 'Sire Qtl effect ['//trim(carac(ic))//"]" 1873 incsol%nbParameterGroup(ieff)=np 1874 do ip=1,np 1875 incsol%parameterName(ieff,ip)='Sire '//trim(pere(ip)) 1876 incsol%paramaterValue(ieff,ip) = ap(ic,ip)*sigt(ic) 1877 end do 1878 1879 if ( count(estime(ic,:)) > 0 ) then 1880 ieff = ieff +1 1881 incsol%groupeName(ieff) = 'Dam Qtl effect ['//trim(carac(ic))//"]" 1882 incsol%nbParameterGroup(ieff)=count(estime(ic,:)) 1883 ifem=0 1884 do jm=1,nm 1885 if ( estime(ic,jm)) then 1886 ifem=ifem+1 1887 incsol%parameterName(ieff,ifem)='Dam '//trim(mere(jm)) 1888 incsol%paramaterValue(ieff,ifem) = am(ic,jm)*sigt(ic) 1889 incsol%parameterVecsol(ieff,ifem) = .true. 1890 end if 1891 end do 1892 end if 1893 1894 end do 1895 1896 allocate (incsol%rhoi(ncar,ncar)) 1897 incsol%rhoi=rhoi2_2 1898 1899 end subroutine set_solution_hypothesis1