m_qtlmap_analyse_multitrait_DA
[ Top ] [ ANALYSE ] [ Modules ]
NAME
m_qtlmap_analyse_multitrait_DA
DESCRIPTION
Module analyse multicaractere DA init_analyse_DA_1QTL,perform_DA_1QTL, opti_DA_0qtl, opti_DA_1qtl, end_analyse_DA_1QTL
NOTES
SEE ALSO
am
[ Top ] [ m_qtlmap_analyse_multitrait_DA ] [ Variables ]
NAME
am
DESCRIPTION
The qtl dam effect found der H1
DIMENSIONS
nm
ap
[ Top ] [ m_qtlmap_analyse_multitrait_DA ] [ Variables ]
NAME
ap
DESCRIPTION
The qtl sire effect found der H1
DIMENSIONS
np
carpm
[ Top ] [ m_qtlmap_analyse_multitrait_DA ] [ Variables ]
NAME
carpm
DESCRIPTION
Sum of square of probabilities to receive the qtl dam (for dam with enough dndmin) initialization : loop before the minimization of the likelihood : opti_DA_1qtl
DIMENSIONS
maxval(ngenom)
carpp
[ Top ] [ m_qtlmap_analyse_multitrait_DA ] [ Variables ]
NAME
carpp
DESCRIPTION
Sum of square of probabilities to receive the qtl dam (for dam with enough dndmin) initialization : loop before the minimization of the likelihood : opti_DA_1qtl
DIMENSIONS
maxval(ngenom)
carppdf
[ Top ] [ m_qtlmap_analyse_multitrait_DA ] [ Variables ]
NAME
carppdf
DESCRIPTION
Sum of square of probabilities to receive the qtl sire initialization : loop before the minimization of the likelihood : opti_DA_1qtl
DIMENSIONS
np
current_chr
[ Top ] [ m_qtlmap_analyse_multitrait_DA ] [ Variables ]
NAME
current_chr
DESCRIPTION
The current chromosome while the likelihood calculs
f0
[ Top ] [ m_qtlmap_analyse_multitrait_DA ] [ Variables ]
NAME
f0
DESCRIPTION
The maximum likelihood found under H0
f_1
[ Top ] [ m_qtlmap_analyse_multitrait_DA ] [ Variables ]
NAME
f_1
DESCRIPTION
The maximum likelihood under H1
fm0
[ Top ] [ m_qtlmap_analyse_multitrait_DA ] [ Variables ]
NAME
fm0
DESCRIPTION
The likelihood by full sib family under H0
DIMENSIONS
nm
fm1
[ Top ] [ m_qtlmap_analyse_multitrait_DA ] [ Variables ]
NAME
fm1
DESCRIPTION
The likelihood by full sib family under H1
DIMENSIONS
nm
fm_1
[ Top ] [ m_qtlmap_analyse_multitrait_DA ] [ Variables ]
NAME
fm_1
DESCRIPTION
The maximum likelihood by full sib family
DIMENSIONS
nm
fp0
[ Top ] [ m_qtlmap_analyse_multitrait_DA ] [ Variables ]
NAME
fp0
DESCRIPTION
The likelihood by half sib family under H0
DIMENSIONS
np
fp1
[ Top ] [ m_qtlmap_analyse_multitrait_DA ] [ Variables ]
NAME
fp1
DESCRIPTION
The likelihood by half sib family under H1
DIMENSIONS
np
fp_1
[ Top ] [ m_qtlmap_analyse_multitrait_DA ] [ Variables ]
NAME
fp_1
DESCRIPTION
The maximum likelihood by half sib family
DIMENSIONS
np
sig1
[ Top ] [ m_qtlmap_analyse_multitrait_DA ] [ Variables ]
NAME
sig1
DESCRIPTION
the standart deviation (residual) found under H0
DIMENSIONS
np
sompm
[ Top ] [ m_qtlmap_analyse_multitrait_DA ] [ Variables ]
NAME
sompm
DESCRIPTION
initialization : loop before the minimization of the likelihood : opti_DA_1qtl
DIMENSIONS
maxval(ngenom)
sompmy
[ Top ] [ m_qtlmap_analyse_multitrait_DA ] [ Variables ]
NAME
sompmy
DESCRIPTION
initialization : loop before the minimization of the likelihood : opti_DA_1qtl
DIMENSIONS
maxval(ngenom)
sompp
[ Top ] [ m_qtlmap_analyse_multitrait_DA ] [ Variables ]
NAME
sompp
DESCRIPTION
initialization : loop before the minimization of the likelihood : opti_DA_1qtl
DIMENSIONS
maxval(ngenom)
somppdf
[ Top ] [ m_qtlmap_analyse_multitrait_DA ] [ Variables ]
NAME
somppdf
DESCRIPTION
Sum of probabilities to receive the qtl sire initialization : loop before the minimization of the likelihood : opti_DA_1qtl
DIMENSIONS
somppm
[ Top ] [ m_qtlmap_analyse_multitrait_DA ] [ Variables ]
NAME
somppm
DESCRIPTION
Sum of probabilities to receive the qtl dam
initialization : loop before the minimization of the likelihood : opti_DA_1qtl
DIMENSIONS
maxval(ngenom)
somppy
[ Top ] [ m_qtlmap_analyse_multitrait_DA ] [ Variables ]
NAME
somppy
DESCRIPTION
Sum of probabilities to receive the qtl dam mult with a trait initialization : loop before the minimization of the likelihood : opti_DA_1qtl
DIMENSIONS
maxval(ngenom)
somppydf
[ Top ] [ m_qtlmap_analyse_multitrait_DA ] [ Variables ]
NAME
somppydf
DESCRIPTION
Sum of probabilities to receive the qtl sire mult with a trait initialization : loop before the minimization of the likelihood : opti_DA_1qtl
DIMENSIONS
np
std
[ Top ] [ m_qtlmap_analyse_multitrait_DA ] [ Variables ]
NAME
std
DESCRIPTION
the standart deviation (residual) found under H1
DIMENSIONS
np
xmoym
[ Top ] [ m_qtlmap_analyse_multitrait_DA ] [ Variables ]
NAME
xmoym
DESCRIPTION
The polygenic dam effect found under H1
DIMENSIONS
nm
xmoyp
[ Top ] [ m_qtlmap_analyse_multitrait_DA ] [ Variables ]
NAME
xmoyp
DESCRIPTION
The polygenic sire effect found under H1
DIMENSIONS
np
xmu1m
[ Top ] [ m_qtlmap_analyse_multitrait_DA ] [ Variables ]
NAME
xmu1m
DESCRIPTION
The polygenic dam effect found under H0
DIMENSIONS
nm
xmu1p
[ Top ] [ m_qtlmap_analyse_multitrait_DA ] [ Variables ]
NAME
xmu1p
DESCRIPTION
The polygenic sire effect found under H0
DIMENSIONS
np
end_analyse_DA_1QTL
[ Top ] [ m_qtlmap_analyse_multitrait_DA ] [ Subroutines ]
NAME
end_analyse_DA_1QTL
DESCRIPTION
deallocation of buffer/solution arrays
SOURCE
394 subroutine end_analyse_DA_1QTL 395 deallocate (sig1) 396 deallocate (xmu1p) 397 deallocate (xmu1m) 398 deallocate (fp0) 399 deallocate (fp1) 400 deallocate (fm0) 401 deallocate (fm1) 402 deallocate (sompp) 403 deallocate (sompm) 404 deallocate (sompmy) 405 deallocate (carpp) 406 deallocate (carpm) 407 deallocate (somppy) 408 deallocate (somppm) 409 deallocate (somppdf) 410 deallocate (carppdf) 411 deallocate (somppydf) 412 deallocate (ap) 413 deallocate (xmoyp) 414 deallocate (std) 415 deallocate (xmoym) 416 deallocate (am) 417 deallocate (fp_1) 418 deallocate (fm_1) 419 end subroutine end_analyse_DA_1QTL
funct_0qtl
[ Top ] [ m_qtlmap_analyse_multitrait_DA ] [ Subroutines ]
NAME
funct_0qtl
DESCRIPTION
Calcul de la vraisemblance et de ses derivees partielles sous H0
SOURCE
696 subroutine funct_0qtl(n,x,f,iuser,user) 697 use m_qtlmap_analyse_gen, only : somy,eff,cary,estime,effdf,somydf,carydf,estmum 698 699 implicit none 700 integer , intent(in) :: n 701 ! 702 ! Tableaux dimensionnes selon n, le nombre de parametres a estimer 703 ! 704 real (kind=dp) ,dimension(n), intent(in) :: x 705 real (kind=dp) , intent(inout) :: f 706 integer , dimension(1), intent(in) :: iuser 707 real (kind=dp) ,dimension(1), intent(in) :: user 708 709 integer :: ip,nm1,nm2,jm 710 integer :: nestim,nest 711 712 real (kind=dp) :: sig,var,vpf,xmup,xmup2 713 real (kind=dp) :: vdf,somxmu,xmum,xmum2,xmupm 714 715 !****************************************************************************** 716 f=0.d0 717 nestim=0 718 do ip=1,np 719 sig=x(ip) 720 var=sig*sig 721 xmup=x(ip+np) 722 xmup2=xmup*xmup 723 vdf=carydf(ip)+dble(effdf(ip))*xmup2-2.d0*xmup*somydf(ip) 724 vdf=0.5d0*vdf/var 725 fp0(ip)=vdf+(dble(effdf(ip))*dlog(sig)) 726 f=f+fp0(ip) 727 728 nm1=nmp(ip)+1 729 nm2=nmp(ip+1) 730 somxmu=0.d0 731 nest=0 732 do jm=nm1,nm2 733 !AJOUT OFI. le tableau 'estime' depend d un caractere 734 ! on ne peut donc pas faire de if estime(ic,jm) 735 ! A VALIDEr PAR HELENE : si un caractere est non estimable alors tous les caracteres ne le sont pas... 736 737 if( count(estime(:,jm))==size(carac) ) then 738 nest=nest+1 739 if(nest.le.estmum(ip)) then 740 nestim=nestim+1 741 xmum=x(2*np+nestim) 742 somxmu=somxmu+xmum 743 else 744 xmum=-somxmu 745 end if 746 xmum2=xmum*xmum 747 xmupm=2.d0*(xmup+xmum) 748 vpf=cary(jm)+dble(eff(jm))*(xmup2+xmum2+2.d0*xmup*xmum)-xmupm*somy(jm) 749 vpf=0.5d0*vpf/var 750 fm0(jm)=vpf+(dble(eff(jm))*dlog(sig)) 751 f=f+fm0(jm) 752 fp0(ip)=fp0(ip)+fm0(jm) 753 else 754 fm0(jm)=0.d0 755 end if 756 end do 757 end do 758 return 759 end subroutine funct_0qtl
funct_1qtl
[ Top ] [ m_qtlmap_analyse_multitrait_DA ] [ Subroutines ]
NAME
funct_1qtl
DESCRIPTION
Calcul de la vraisemblance et de ses derivees partielles sous Hypothese 1QTL 1carac
SOURCE
965 subroutine funct_1qtl(n,x,f,iuser,user) 966 967 integer , intent(in) :: n ! number of variables 968 double precision, intent(in) :: x(n) ! Position 969 double precision, intent(out) :: f ! Result value at the position 970 double precision user(1) 971 integer iuser(1) 972 973 974 ! 975 ! Divers 976 integer :: ip,nestim,nm1,nm2,nest,jm,ifem,indam,i 977 integer :: ngeno1,ngeno2,ig 978 real (kind=dp) :: sig,var,xmup,xmup2,vdf,somxmu,xmum,xmum2 979 real (kind=dp) :: xmupm,vmere,z,vpf,x_ap,x_ap2,x_am,x_am2 980 981 982 !****************************************************************************** 983 f=0.d0 984 nestim=0 985 986 do ip=1,np 987 sig=x(ip) 988 var=sig*sig 989 xmup=x(ip+np) 990 xmup2=xmup*xmup 991 x_ap=x(2*np+ip) 992 x_ap2=x_ap*x_ap 993 vdf=carydf(ip)+dble(effdf(ip))*xmup2-2.d0*xmup*somydf(ip) 994 vdf=vdf+x_ap2*carppdf(ip)+2.d0*xmup*x_ap*somppdf(ip)-2.d0*x_ap*somppydf(ip) 995 vdf=0.5d0*vdf/var 996 fp1(ip)=vdf+dble(effdf(ip))*dlog(sig) 997 f=f+fp1(ip) 998 nm1=nmp(ip)+1 999 nm2=nmp(ip+1) 1000 somxmu=0.d0 1001 nest=0 1002 do jm=nm1,nm2 1003 !AJOUT OFI. le tableau 'estime' depend d un caractere 1004 ! on ne peut donc pas faire de if estime(ic,jm) 1005 ! 1006 ! A VALIDEr PAR HELENE : si un caractere est non estimable alors tous les caracteres ne le sont pas... 1007 if( count(estime(:,jm))== size(carac) ) then 1008 nest=nest+1 1009 if(nest.le.estmum(ip)) then 1010 nestim=nestim+1 1011 xmum=x(3*np+nestim) 1012 somxmu=somxmu+xmum 1013 else 1014 xmum=-somxmu 1015 end if 1016 xmum2=xmum*xmum 1017 xmupm=2.d0*(xmup+xmum) 1018 ifem=repfem(jm) 1019 indam=iam(1,ifem) 1020 x_am=x(3*np+nmumest(1)+indam) 1021 x_am2=x_am*x_am 1022 z=cary(jm)+dble(eff(jm))*(xmup2+xmum2+2.d0*xmup*xmum)-xmupm*somy(jm) 1023 vmere=0.d0 1024 ngeno1=ngenom(current_ch,jm)+1 1025 ngeno2=ngenom(current_ch,jm+1) 1026 do ig=ngeno1,ngeno2 1027 vpf=z+x_ap2*carpp(ig)+x_am2*carpm(ig) & 1028 +2.d0*x_ap*x_am*somppm(ig) & 1029 +xmupm*(x_ap*sompp(ig)+x_am*sompm(ig)) & 1030 -2.d0*x_ap*somppy(ig)-2.d0*x_am*sompmy(ig) 1031 vpf=-0.5d0*vpf/var 1032 vmere=vmere+probg(current_ch,ig)*dexp(vpf) 1033 end do 1034 ! if (vmere.lt.1.e-8) then 1035 ! write (nficerr,*) '** WARNING : vmere.eq.0 ', 1036 ! $ 'in funct_1qtl.F **' 1037 ! fm1(jm)=(dble(eff(jm))*dlog(sig)) 1038 ! else 1039 fm1(jm)=-dlog(vmere)+(dble(eff(jm))*dlog(sig)) 1040 ! end if 1041 f=f+fm1(jm) 1042 fp1(ip)=fp1(ip)+fm1(jm) 1043 else 1044 fm1(jm)=0.d0 1045 end if 1046 end do 1047 end do 1048 return 1049 end subroutine funct_1qtl
init_analyse_DA_1QTL
[ Top ] [ m_qtlmap_analyse_multitrait_DA ] [ Subroutines ]
NAME
init_analyse_DA_1QTL
DESCRIPTION
allocation of buffer/solution arrays
SOURCE
311 subroutine init_analyse_DA_1QTL 312 integer :: ic,jm,stat,ng 313 314 ! Get Log Debug Information about estim 315 316 do jm=1,size(estime,2) 317 if( count(estime(:,jm))== size(carac) ) then 318 call log_mess('DAM ['//trim(mere(jm))//"] estime ** OK **",DEBUG_DEF) 319 else 320 call log_mess('DAM ['//trim(mere(jm))//"] estime -- KO --",DEBUG_DEF) 321 do ic=1,size(carac) 322 if ( .not. estime(ic,jm) ) then 323 call log_mess(' --> Trait ['//trim(carac(ic))//'] is not estim ',DEBUG_DEF) 324 end if 325 end do 326 end if 327 end do 328 329 ng = maxval(ngenom) 330 allocate (sig1(np),STAT=stat) 331 call check_allocate(stat,'sig1 [m_qtlmap_analyse_multitrait_DA]') 332 allocate (xmu1p(np),STAT=stat) 333 call check_allocate(stat,'xmu1p [m_qtlmap_analyse_multitrait_DA]') 334 allocate (xmu1m(nm),STAT=stat) 335 call check_allocate(stat,'xmu1m [m_qtlmap_analyse_multitrait_DA]') 336 allocate (fp0(np),STAT=stat) 337 call check_allocate(stat,'fp0 [m_qtlmap_analyse_multitrait_DA]') 338 allocate (fp1(np),STAT=stat) 339 call check_allocate(stat,'fp1 [m_qtlmap_analyse_multitrait_DA]') 340 allocate (fm0(nm),STAT=stat) 341 call check_allocate(stat,'fm0 [m_qtlmap_analyse_multitrait_DA]') 342 allocate (fm1(nm),STAT=stat) 343 call check_allocate(stat,'fm1 [m_qtlmap_analyse_multitrait_DA]') 344 allocate (sompp(ng),STAT=stat) 345 call check_allocate(stat,'sompp [m_qtlmap_analyse_multitrait_DA]') 346 allocate (sompm(ng),STAT=stat) 347 call check_allocate(stat,'sompm [m_qtlmap_analyse_multitrait_DA]') 348 allocate (sompmy(ng),STAT=stat) 349 call check_allocate(stat,'sompmy [m_qtlmap_analyse_multitrait_DA]') 350 allocate (carpp(ng),STAT=stat) 351 call check_allocate(stat,'carpp [m_qtlmap_analyse_multitrait_DA]') 352 allocate (carpm(ng),STAT=stat) 353 call check_allocate(stat,'carpm [m_qtlmap_analyse_multitrait_DA]') 354 allocate (somppy(ng),STAT=stat) 355 call check_allocate(stat,'somppy [m_qtlmap_analyse_multitrait_DA]') 356 allocate (somppm(ng),STAT=stat) 357 call check_allocate(stat,'somppm [m_qtlmap_analyse_multitrait_DA]') 358 allocate (somppdf(np),STAT=stat) 359 call check_allocate(stat,'somppdf [m_qtlmap_analyse_multitrait_DA]') 360 allocate (carppdf(np),STAT=stat) 361 call check_allocate(stat,'carppdf [m_qtlmap_analyse_multitrait_DA]') 362 allocate (somppydf(np),STAT=stat) 363 call check_allocate(stat,'somppydf [m_qtlmap_analyse_multitrait_DA]') 364 allocate (ap(np),STAT=stat) 365 call check_allocate(stat,'ap [m_qtlmap_analyse_multitrait_DA]') 366 ap=0.d0 367 allocate (xmoyp(np),STAT=stat) 368 call check_allocate(stat,'xmoyp [m_qtlmap_analyse_multitrait_DA]') 369 xmoyp=0.d0 370 allocate (std(np),STAT=stat) 371 call check_allocate(stat,'std [m_qtlmap_analyse_multitrait_DA]') 372 std=0.d0 373 allocate (am(nm),STAT=stat) 374 call check_allocate(stat,'am [m_qtlmap_analyse_multitrait_DA]') 375 am=0.d0 376 allocate (xmoym(nm),STAT=stat) 377 call check_allocate(stat,'xmoym [m_qtlmap_analyse_multitrait_DA]') 378 xmoym=0.d0 379 allocate (fp_1(np),STAT=stat) 380 call check_allocate(stat,'fp_1 [m_qtlmap_analyse_multitrait_DA]') 381 allocate (fm_1(nm),STAT=stat) 382 call check_allocate(stat,'fm_1 [m_qtlmap_analyse_multitrait_DA]') 383 384 end subroutine init_analyse_DA_1QTL
opti_DA_0qtl
[ Top ] [ m_qtlmap_analyse_multitrait_DA ] [ Subroutines ]
NAME
opti_DA_0qtl
DESCRIPTION
Calcul de la vraisemblance 0 QTL, 1 caractere
SOURCE
627 subroutine opti_DA_0qtl 628 use m_qtlmap_analyse_gen, only : sig0,xmu0p,xmu0m 629 630 implicit none 631 ! 632 ! Divers 633 634 integer iuser(1) 635 double precision ,dimension(:),allocatable :: borni,borns,par 636 real (kind=dp) :: user(1) 637 integer :: npar,ibound,ip,jm,ifail 638 639 !****************************************************************************** 640 ! Parametres de maximisation 641 npar=(2*np)+nmumest(1) 642 643 allocate (borni(npar)) 644 allocate (borns(npar)) 645 allocate (par(npar)) 646 647 ibound=0 648 do ip=1,np 649 borni(ip)=1.do-6 650 borns(ip)=1.d6 651 borni(ip+np)=-1d6 652 borns(ip+np)=1.d6 653 end do 654 do jm=1,nmumest(1) 655 borni(2*np+jm)=-1.d6 656 borns(2*np+jm)=1.d6 657 end do 658 ! 659 ! Point de depart 660 do ip=1,np 661 par(ip)=sig0(ip) 662 par(ip+np)=xmu0p(ip) 663 end do 664 do jm=1,nmumest(1) 665 par(2*np+jm)=xmu0m(jm) 666 end do 667 ! 668 ! Optimisation de la vraisemblance 669 ifail=0 670 671 call minimizing_funct(npar,ibound,funct_0qtl,borni,borns,par,f0,iuser,user,ifail) 672 673 do ip = 1,np 674 sig1(ip)=par(ip) 675 xmu1p(ip)=par(ip+np) 676 end do 677 do jm=1,nmumest(1) 678 xmu1m(jm)=par(2*np+jm) 679 end do 680 ! 681 deallocate (borni) 682 deallocate (borns) 683 deallocate (par) 684 685 return 686 end subroutine opti_DA_0qtl
opti_DA_1qtl
[ Top ] [ m_qtlmap_analyse_multitrait_DA ] [ Subroutines ]
NAME
opti_DA_1qtl
DESCRIPTION
Computing statistical test across the chromosome
SOURCE
772 subroutine opti_DA_1qtl(ch,ix,n,yda,lrtsol) 773 774 integer , intent(in) :: ch,n,ix 775 type(TYPE_LRT_SOLUTION) , intent(inout) :: lrtsol 776 real (kind=dp) ,intent(in) ,dimension(nd) :: yda 777 ! 778 ! Divers 779 integer :: iuser(1) 780 double precision ,dimension(:),allocatable :: borni,borns,par 781 double precision :: user(1) 782 783 integer :: ibound,ip,jm,ngeno1,ngeno2,ilong,nm1,nm2,ig 784 integer :: npar,nd1,nd2,kd,kkd,ifail,ii 785 integer :: nestim,nest,ifem,indam 786 real (kind=dp) :: pp, pm, somxmu,xlrt_t,f1 787 788 789 !****************************************************************************** 790 ! Calcul de la vraisemblance sous H1 791 ! Parametres de maximisation 792 ibound=0 793 npar=(3*np)+nmumest(1)+namest(1) 794 795 call log_mess("**** OPTI_DA_1QTL *****",DEBUG_DEF) 796 call log_mess("NP:"//str(np),DEBUG_DEF) 797 call log_mess("NPAR:"//str(npar),DEBUG_DEF) 798 799 allocate (borni(npar)) 800 allocate (borns(npar)) 801 allocate (par(npar)) 802 803 current_ch = ch 804 805 do ip=1,np 806 borni(ip)=1.do-6 807 borns(ip)=1.d6 808 borni(ip+np)=-1.d6 809 borns(ip+np)=1.d6 810 borni(2*np+ip)=-1.d6 811 borns(2*np+ip)=1.d6 812 end do 813 do jm=1,nmumest(1) 814 borni(3*np+jm)=-1.d6 815 borns(3*np+jm)=1.d6 816 end do 817 do jm=1,namest(1) 818 borni(3*np+nmumest(1)+jm)=-1.d6 819 borns(3*np+nmumest(1)+jm)=1.d6 820 end do 821 par=0.d0 822 ! 823 ! Point de depart 824 do ip=1,np 825 par(ip)=sig1(ip) 826 par(ip+np)=xmu1p(ip) 827 par(2*np+ip)=0.d0 828 end do 829 830 831 do jm=1,nmumest(1) 832 par(3*np+jm)=xmu1m(jm) 833 end do 834 do jm=1,namest(1) 835 par(3*np+nmumest(1)+jm)=0.d0 836 end do 837 ! 838 !A la position en cours 839 do ip=1,np 840 nm1=nmp(ip)+1 841 nm2=nmp(ip+1) 842 somppdf(ip)=0.d0 843 carppdf(ip)=0.d0 844 somppydf(ip)=0.d0 845 do jm=nm1,nm2 846 ngeno1=ngenom(ch,jm)+1 847 ngeno2=ngenom(ch,jm+1) 848 do ig=ngeno1,ngeno2 849 sompp(ig)=0.d0 850 carpp(ig)=0.d0 851 somppy(ig)=0.d0 852 sompm(ig)=0.d0 853 carpm(ig)=0.d0 854 sompmy(ig)=0.d0 855 somppm(ig)=0.d0 856 nd1=ngend(ch,ig)+1 857 nd2=ngend(ch,ig+1) 858 do kd=nd1,nd2 859 kkd=ndesc(ch,kd) 860 if(presentc(1,kkd)) then 861 pp=-pdd(ch,kd,1,n)-pdd(ch,kd,2,n)+pdd(ch,kd,3,n)+pdd(ch,kd,4,n) 862 pm=-pdd(ch,kd,1,n)+pdd(ch,kd,2,n)-pdd(ch,kd,3,n)+pdd(ch,kd,4,n) 863 !AJOUT OFI. le tableau 'estime' depend d un caractere 864 ! on ne peut donc pas faire de if estime(ic,jm) 865 ! 866 ! A VALIDEr PAR HELENE : si un caractere est non estimable alors tous les caracteres ne le sont pas... 867 868 if( count(estime(:,jm))== size(carac) ) then 869 sompp(ig)=sompp(ig)+pp 870 carpp(ig)=carpp(ig)+pp*pp 871 somppy(ig)=somppy(ig)+pp*yda(kkd) 872 sompm(ig)=sompm(ig)+pm 873 carpm(ig)=carpm(ig)+pm*pm 874 sompmy(ig)=sompmy(ig)+pm*yda(kkd) 875 somppm(ig)=somppm(ig)+pp*pm 876 else 877 somppdf(ip)=somppdf(ip)+pp 878 carppdf(ip)=carppdf(ip)+pp*pp 879 somppydf(ip)=somppydf(ip)+pp*yda(kkd) 880 end if 881 end if 882 end do 883 end do 884 end do 885 end do 886 ! 887 ! Optimisation de la vraisemblance a la position en cours 888 ifail=0 889 890 call minimizing_funct(npar,ibound,funct_1qtl,borni,borns,par,f1,iuser,user,ifail) 891 892 xlrt_t=-2.d0*(f1-f0) 893 call log_mess('N:'//str(n)//" LRT:"//str(xlrt_t)//" F0:"//str(f0)//" F1:"//str(f1),DEBUG_DEF) 894 do ii=1,np 895 lrtsol%xlrp(ch,ii,n)=-2.d0*(fp1(ii)-fp0(ii)) 896 lrtsol%pater_eff(ch,ii,n)=par(2*np+ii) 897 end do 898 899 do ii=1,nm 900 lrtsol%xlrm(ch,ii,n)=-2.d0*(fm1(ii)-fm0(ii)) 901 end do 902 903 do ii=1,namest(1) 904 lrtsol%mater_eff(ch,ii,n)=par(3*np+nmumest(1)+ii) 905 end do 906 907 if(lrtsol%lrtmax(0).le.xlrt_t) then 908 lrtsol%lrtmax(0)= xlrt_t 909 lrtsol%nxmax(0)= n 910 lrtsol%chrmax(0)= ch 911 f_1=f1 ! ic 912 nestim=0 913 do ip=1,np 914 ap(ip)=par(2*np+ip) 915 xmoyp(ip)=par(np+ip) 916 std(ip)=par(ip) 917 fp_1(ip)=fp1(ip) 918 nm1=nmp(ip)+1 919 nm2=nmp(ip+1) 920 nest=0 921 somxmu=0.d0 922 do jm=nm1,nm2 923 fm_1(jm)=fm1(jm) 924 !AJOUT OFI. le tableau 'estime' depend d un caractere 925 ! on ne peut donc pas faire de if estime(ic,jm) 926 ! 927 ! A VALIDEr PAR HELENE : si un caractere est non estimable alors tous les caracteres ne le sont pas... 928 if( count(estime(:,jm))== size(carac) ) then 929 !if (estime(ic,jm)) then 930 nest=nest+1 931 ifem=repfem(jm) 932 indam=iam(1,ifem) 933 am(jm)=par(3*np+nmumest(1)+indam) 934 if (nest.le.estmum(ip)) then 935 nestim=nestim+1 936 xmoym(jm)=par(3*np+nestim) 937 somxmu=somxmu+xmoym(jm) 938 else 939 xmoym(jm)=-somxmu 940 end if 941 else 942 am(jm)=0.d0 943 xmoym(jm)=0.d0 944 end if 945 end do 946 end do 947 end if 948 lrtsol%lrt1(ch,n)=xlrt_t 949 ! 950 deallocate (borni) 951 deallocate (borns) 952 deallocate (par) 953 954 return 955 end subroutine opti_DA_1qtl
perform_DA_1qtl
[ Top ] [ m_qtlmap_analyse_multitrait_DA ] [ Subroutines ]
NAME
perform_DA_1qtl
DESCRIPTION
Calcul des combinaisons lineaires des performances pour la position testee
INPUTS
ch : index of the chromosome tested n : the position tested npo : number of position tested among the chromosome
OUTPUTS
yda : coeff :
SOURCE
439 subroutine perform_DA_1qtl(ch,n,yda,coeff,npo) 440 ! 441 integer ,intent(in) :: ch ! chromosome 442 real (kind=dp) ,intent(out) ,dimension(nd) :: yda 443 integer ,intent(in) :: n,npo 444 real (kind=dp) ,intent(inout),dimension(size(carac),npo) :: coeff 445 446 integer :: ic,k,kkd,nm1,nm2,jm,ip,ngeno1,ngeno2,nd1,nd2,kd,ig 447 integer :: nit,ind,ifail,it,NCV,IRANKX,LDX1,NI,IWKI 448 real(kind=dp) , dimension(:,:),allocatable :: Xy,pb 449 real(kind=dp) , dimension(:),allocatable :: WT, WK 450 integer , dimension(:),allocatable :: NIG, ING, ISX 451 integer , dimension(:,:),allocatable :: NIGP 452 integer :: NGP,LDX,M,NBX,LDCVM,LDE,LDCVX,IWK 453 454 real (kind=dp) , parameter :: TOL=0.01d0 455 real(kind=dp) , dimension(:,:) ,allocatable :: CVM, E, CVX 456 457 character (len=1) :: WEIGHT='W' 458 real (kind=dp) ::s1,s2 459 460 ! Initialisation des parametres 461 NGP=2 462 463 ! Allocation tableau des poids (4*nd est mal ajust�) 464 allocate(pb(maxval(ndesc),NGP)) ! 465 466 do ic=1,ncar 467 ! Initialisation des poids des perf dans chaque groupe genotypique 468 nit=0 469 do ip=1,np 470 nm1=nmp(ip)+1 471 nm2=nmp(ip+1) 472 do jm=nm1,nm2 473 ngeno1=ngenom(ch,jm)+1 474 ngeno2=ngenom(ch,jm+1) 475 do ig=ngeno1,ngeno2 476 nd1=ngend(ch,ig)+1 477 nd2=ngend(ch,ig+1) 478 do kd=nd1,nd2 479 kkd=ndesc(ch,kd) 480 if(count(presentc(:,kkd))==size(carac)) then 481 if (NGP.eq.2) then 482 pb(kd,1)=pdd(ch,kd,1,n)+pdd(ch,kd,2,n) 483 pb(kd,2)=pdd(ch,kd,3,n)+pdd(ch,kd,4,n) 484 else 485 if (NGP.eq.4) then 486 pb(kd,1)=pdd(ch,kd,1,n) 487 pb(kd,2)=pdd(ch,kd,2,n) 488 pb(kd,3)=pdd(ch,kd,3,n) 489 pb(kd,4)=pdd(ch,kd,4,n) 490 else 491 if (NGP.Eq.3) then 492 pb(kd,1)=pdd(ch,kd,1,n) 493 pb(kd,2)=pdd(ch,kd,2,n)+pdd(ch,kd,3,n) 494 pb(kd,3)=pdd(ch,kd,4,n) 495 end if ! NGP 3 496 end if ! NGP 4 497 end if ! NGP 2 498 nit=nit+1 499 end if ! PRESENTC 500 end do 501 end do 502 end do 503 end do 504 end do ! ic 505 ! Initialisation des dimensions du probleme 506 NI=nit*NGP 507 508 ! Initialisation des parametres 509 LDX=NI 510 M=size(carac) 511 NBX=size(carac) 512 LDCVM=NGP 513 if(size(carac)> NGP) then 514 LDE=size(carac) 515 else 516 LDE=NGP 517 end if 518 LDCVX=size(carac) 519 if(NI> NGP-1) then 520 IWKi=NI 521 else 522 IWKi=NGP-1 523 end if 524 IWK=NI*NBX+5*(NBX-1)+(IWKi+1)*NBX+2 525 IFAIL=0 526 527 ! Allocation des vecteurs et tables 528 allocate(Xy(LDX,size(carac))) 529 allocate(WT(LDX)) 530 allocate(WK(IWK)) 531 allocate(NIG(NGP)) 532 allocate(ING(NI)) 533 allocate(ISX(M)) 534 allocate(CVM(LDCVM,NBX)) 535 allocate(E(LDE,6)) 536 allocate(CVX(LDCVX,NGP-1)) 537 538 ! Initialisation des vecteurs et tables 539 Xy=0.d0 540 WT=0.d0 541 WK=0.d0 542 NIG=nit 543 ING=0 544 ISX=1 545 CVM=0.d0 546 E=0.d0 547 CVX=0.d0 548 ! determination de la VC � la position 549 ind=0 550 LDX1=0 551 ISX(:)=1 ! all trait are included in the analysis 552 do while (ind < NGP) 553 do ip=1,np 554 nm1=nmp(ip)+1 555 nm2=nmp(ip+1) 556 do jm=nm1,nm2 557 nd1=ndm(jm)+1 558 nd2=ndm(jm+1) 559 ngeno1=ngenom(ch,jm)+1 560 ngeno2=ngenom(ch,jm+1) 561 do ig=ngeno1,ngeno2 562 nd1=ngend(ch,ig)+1 563 nd2=ngend(ch,ig+1) 564 do kd=nd1,nd2 565 kkd=ndesc(ch,kd) 566 if(count(presentc(:,kkd))==size(carac)) then 567 LDX1=LDX1+1 568 ING(LDX1)=1+ind 569 Xy(LDX1,:)=y(:,kkd) 570 WT(LDX1)=pb(kd,1+ind) 571 end if 572 end do 573 end do 574 end do 575 end do 576 ind=ind+1 577 end do 578 579 s1=0.d0 580 s2=s1 581 do kd=1,nit 582 s1=s1+WT(kd) 583 end do 584 do kd=1+nit,nit*2 585 s2=s2+WT(kd) 586 end do 587 ! print*,'sum', s1,s2 588 call MATH_QTLMAP_G03ACF(WEIGHT,NI,M,Xy,LDX,ISX,NBX,ING,NGP,WT,NIG,CVM,LDCVM,E & 589 ,LDE,NCV,CVX,LDCVX,TOL,IRANKX,IFAIL) 590 if (ifail.ne.0) then 591 do it=1,LDX1 592 ! call log_mess('LDX1 '//trim(str(LDX1))//trim(str(WT(LDX1),DEBUG_DEF))) 593 end do 594 call stop_application('Analyse discriminante sous H1; ifail='//trim(str(ifail))) 595 end if 596 597 do kd=1,nd 598 yda(kd)=0.d0 599 do ic=1,size(carac) 600 yda(kd)=CVX(ic,1)*y(ic,kd)+yda(kd) 601 end do 602 end do 603 do ic=1,size(carac) 604 coeff(ic,n)=CVX(ic,1) 605 end do 606 607 deallocate(pb) 608 deallocate(Xy) 609 deallocate(WT) 610 deallocate(WK) 611 deallocate(NIG) 612 deallocate(ING) 613 deallocate(ISX) 614 deallocate(CVM) 615 deallocate(E) 616 deallocate(CVX) 617 end subroutine perform_DA_1qtl
set_solution_hypothesis1
[ Top ] [ m_qtlmap_analyse_multitrait_DA ] [ Subroutines ]
NAME
set_solution_hypothesis1
DESCRIPTION
SOURCE
1058 subroutine set_solution_hypothesis1(incsol) 1059 type(TYPE_INCIDENCE_SOLUTION) ,intent(inout) :: incsol 1060 1061 integer :: nteff,maxNbPar,jm,ifem,ip,ieff 1062 1063 allocate (incsol%sig(1,np)) 1064 1065 incsol%hypothesis=1 1066 ! Mean family , Qtl effect 1067 nteff = 2 1068 if ( count(estime(1,:)) > 0 ) nteff = 4 1069 1070 maxNbPar = max(np,count(estime(1,:))) 1071 allocate (incsol%groupeName(nteff)) 1072 allocate (incsol%nbParameterGroup(nteff)) 1073 allocate (incsol%parameterName(nteff,maxNbPar)) 1074 allocate (incsol%paramaterValue(nteff,maxNbPar)) 1075 allocate (incsol%parameterVecsol(nteff,maxNbPar)) 1076 allocate (incsol%parameterPrecis(nteff,maxNbPar)) 1077 allocate (incsol%qtl_groupeName(1,1)) 1078 1079 incsol%parameterPrecis=0.d0 1080 incsol%parameterVecsol=.true. 1081 1082 do ip=1,np 1083 incsol%sig(1,ip) = std(ip)*sigt(1) 1084 end do 1085 1086 ieff=1 1087 incsol%qtl_groupeName(1,1)=ieff 1088 incsol%groupeName(ieff) = 'Sire Qtl effect' 1089 incsol%nbParameterGroup(ieff)=np 1090 1091 do ip=1,np 1092 incsol%parameterName(ieff,ip)='Sire '//trim(pere(ip)) 1093 incsol%paramaterValue(ieff,ip) = ap(ip)*sigt(1)! en DA quel sigt prendre ? *sigt(ic) 1094 end do 1095 1096 if ( count(estime(1,:)) > 0 ) then 1097 ieff = ieff +1 1098 incsol%groupeName(ieff) = 'Dam Qtl effect' 1099 incsol%nbParameterGroup(ieff)=namest(1) 1100 ifem=0 1101 do jm=1,nm 1102 if ( estime(1,jm)) then 1103 ifem=ifem+1 1104 incsol%parameterName(ieff,ifem)='Dam '//trim(mere(jm)) 1105 incsol%paramaterValue(ieff,ifem) = am(jm)*sigt(1) 1106 incsol%parameterVecsol(ieff,ifem) = .true. 1107 end if 1108 end do 1109 end if 1110 1111 1112 ieff = ieff +1 1113 incsol%groupeName(ieff) = 'Mean Sire' 1114 incsol%nbParameterGroup(ieff)=np 1115 1116 do ip=1,np 1117 incsol%parameterName(ieff,ip)='Sire '//trim(pere(ip)) 1118 incsol%paramaterValue(ieff,ip) = xmoyp(ip)*sigt(1) + xmut(1) 1119 incsol%parameterVecsol(ieff,ip) = .true. 1120 end do 1121 1122 if ( count(estime(1,:)) > 0 ) then 1123 ieff = ieff +1 1124 incsol%groupeName(ieff) = 'Mean dam' 1125 incsol%nbParameterGroup(ieff)=count(estime(1,:)) 1126 ifem=0 1127 do ip=1,np 1128 do jm=nmp(ip)+1,nmp(ip+1) 1129 if ( estime(1,jm)) then 1130 ifem=ifem+1 1131 incsol%parameterName(ieff,ifem)='Dam '//trim(mere(jm))//" [Sire "//trim(pere(ip))//"]" 1132 incsol%paramaterValue(ieff,ifem) = xmoym(jm)*sigt(1) 1133 incsol%parameterVecsol(ieff,ifem) = .true. 1134 end if 1135 end do 1136 end do 1137 end if 1138 1139 end subroutine set_solution_hypothesis1