m_qtlmap_incidence_multi
[ Top ] [ ANALYSE ] [ Modules ]
NAME
m_qtlmap_incidence_multi
DESCRIPTION
NOTES
SEE ALSO
current_chr
[ Top ] [ m_qtlmap_incidence_multi ] [ Variables ]
NAME
current_chr
DESCRIPTION
current linkage group tested
NOTES
dataset
[ Top ] [ m_qtlmap_incidence_multi ] [ Variables ]
NAME
dataset
DESCRIPTION
The current dataset using bi the likelihood function
NOTES
estfem_multi
[ Top ] [ m_qtlmap_incidence_multi ] [ Variables ]
NAME
estfem_multi
DESCRIPTION
Estimability per progeny genotyped for all trait combined
NOTES
estime_multi
[ Top ] [ m_qtlmap_incidence_multi ] [ Variables ]
NAME
estime_multi
DESCRIPTION
Estimability per progeny phenotyped for all trait combined
NOTES
iam_multi
[ Top ] [ m_qtlmap_incidence_multi ] [ Variables ]
NAME
iam_multi
DESCRIPTION
Indexe of the femal for all trait combined
NOTES
my_listdesc
[ Top ] [ m_qtlmap_incidence_multi ] [ Variables ]
NAME
my_listdesc
DESCRIPTION
description of each contingence matrix
NOTES
my_xincreduitmul
[ Top ] [ m_qtlmap_incidence_multi ] [ Variables ]
NAME
my_xincreduitmul
DESCRIPTION
contingence matrix for each trait
NOTES
namest_multi
[ Top ] [ m_qtlmap_incidence_multi ] [ Variables ]
NAME
namest_multi
DESCRIPTION
Number of female with estimability for all trait combined
NOTES
ntnivmaxtotal
[ Top ] [ m_qtlmap_incidence_multi ] [ Variables ]
NAME
ntnivmaxtotal
DESCRIPTION
maximum level finded from all contingence matrix
NOTES
add_effcov_multi
[ Top ] [ m_qtlmap_incidence_multi ] [ Subroutines ]
NAME
add_effcov_multi
DESCRIPTION
NOTE
SOURCE
916 subroutine add_effcov_multi(xinc,listDesc,meff,mcov) 917 real (kind=dp) , dimension(ncar,nd,ntnivmaxtotal) , intent(inout) :: xinc 918 type(INCIDENCE_TYPE) , dimension(ncar) , intent(inout) :: listDesc 919 integer , intent(in) :: meff ! indice on effet fixed list to remove (<=0 otherwise) 920 integer , intent(in) :: mcov ! indice on covariate list to remove (<=0 otherwise) 921 integer :: ic,ntniv,ip,jm,kd,ntnifix,nbt,nbef,nbco,nteff,jef,jcov,ieff,nbniv,nivd,i,j,nbef2,nbco2,nkd 922 923 call log_mess("incidence : ** Add eff and cov ** ",DEBUG_DEF) 924 925 do ic=1,ncar 926 927 nbef=modele(ic,1) 928 nbco=modele(ic,2) 929 930 nbef2=nbef 931 if ( meff /= 0 ) nbef2=nbef-1 932 nbco2=nbco 933 if ( mcov /=0 ) nbco2=nbco-1 934 935 936 if (nbef2+nbco2 <=0)return 937 938 nbniv=0 939 ieff=listDesc(ic)%nteff 940 941 if ( nbef2 > 0 ) then 942 ieff=ieff+1 943 listDesc(ic)%eqtl_print(ieff)=.false. 944 do i=1,nbef 945 if ( i == meff ) cycle 946 nbniv=nlev(modele(ic,3+i))+nbniv 947 end do 948 listDesc(ic)%desc(ieff)%name="Fixed effects" 949 listDesc(ic)%desc(ieff)%start=listDesc(ic)%ntniv+1 950 listDesc(ic)%desc(ieff)%end=listDesc(ic)%ntniv+nbniv 951 listDesc(ic)%desc(ieff)%haveSubDesc=.true. 952 953 !init of xinc 954 xinc(ic,:,listDesc(ic)%desc(ieff)%start:listDesc(ic)%desc(ieff)%end)=0.d0 955 956 allocate ( listDesc(ic)%desc(ieff)%listSubDesc(nbniv) ) 957 958 j=0 959 do jef=1,nbef 960 if ( jef == meff ) cycle 961 do i=1,nlev(modele(ic,3+jef)) 962 j=j+1 963 listDesc(ic)%desc(ieff)%listSubDesc(j)%name=trim(namefix(modele(ic,3+jef)))//' level '//trim(str(i))!&//"["//trim(carac(ic))//"]" 964 listDesc(ic)%desc(ieff)%listSubDesc(j)%start=listDesc(ic)%ntniv+j 965 listDesc(ic)%desc(ieff)%listSubDesc(j)%end=listDesc(ic)%ntniv+j 966 end do 967 end do 968 listDesc(ic)%nteff = listDesc(ic)%nteff+1 969 end if 970 971 if ( nbco2 > 0 ) then 972 ieff=ieff+1 973 listDesc(ic)%eqtl_print(ieff)=.false. 974 listDesc(ic)%desc(ieff)%name="Covariates" 975 listDesc(ic)%desc(ieff)%start=listDesc(ic)%ntniv+nbniv+1 976 listDesc(ic)%desc(ieff)%end=listDesc(ic)%ntniv+nbniv+nbco2 977 listDesc(ic)%desc(ieff)%haveSubDesc=.true. 978 979 !init of xinc 980 xinc(ic,:,listDesc(ic)%desc(ieff)%start:listDesc(ic)%desc(ieff)%end)=0.d0 981 982 allocate ( listDesc(ic)%desc(ieff)%listSubDesc(nbco2) ) 983 984 j=0 985 do jcov=1,nbco 986 if ( jcov == mcov ) cycle 987 j=j+1 988 listDesc(ic)%desc(ieff)%listSubDesc(j)%name=trim(namecov(modele(ic,3+modele(ic,1)+jcov)))!//"["//trim(carac(ic))//"]" 989 listDesc(ic)%desc(ieff)%listSubDesc(j)%start=listDesc(ic)%ntniv+nbniv+j 990 listDesc(ic)%desc(ieff)%listSubDesc(j)%end=listDesc(ic)%ntniv+nbniv+j 991 end do 992 listDesc(ic)%nteff = listDesc(ic)%nteff+1 993 end if 994 995 nkd = 0 996 nbniv=0 997 do ip=1,np 998 do jm=nmp(ip)+1,nmp(ip+1) 999 do kd=ndm(jm)+1,ndm(jm+1) 1000 1001 if ( count(presentc(:,kd))==ncar ) then 1002 nkd=nkd+1 1003 ntniv=listDesc(ic)%ntniv 1004 1005 ! 1006 ! autres effets fixes 1007 ! 1008 nbniv=0 1009 do jef=1,nbef 1010 if ( jef == meff ) cycle ! this effet fixed is not adding 1011 ! listDesc(ic)%nivdir(kd,nteff+jef)=ntniv+nbniv+nivx(kd,modele(ic,3+jef)) 1012 nivd=ntniv+nbniv+nivx(kd,modele(ic,3+jef)) 1013 xinc(ic,nkd,nivd)=1 1014 nbniv=nbniv+nlev(modele(ic,3+jef)) 1015 end do 1016 1017 ntnifix=ntniv+nbniv 1018 nbt=3+nbef 1019 !covariate 1020 j=0 1021 do jcov=1,nbco 1022 if ( jcov == mcov ) cycle ! this covariate is not adding 1023 j=j+1 1024 ! listDesc(ic)%covdir(kd,jcov)=covar(kd,modele(ic,nbt+jcov)) 1025 xinc(ic,nkd,ntnifix+j)=covar(kd,modele(ic,nbt+jcov)) 1026 end do 1027 !ntniv=ntnifix+nbco 1028 !nbt=nbt+nbco 1029 end if 1030 end do 1031 end do 1032 end do 1033 1034 listDesc(ic)%ntniv=listDesc(ic)%ntniv+nbniv+nbco2 1035 end do !ic 1036 1037 call log_mess("FIN ** Add eff and cov ** ",DEBUG_DEF) 1038 1039 end subroutine add_effcov_multi
add_general_mean_multi
[ Top ] [ m_qtlmap_incidence_multi ] [ Subroutines ]
NAME
add_general_mean_multi
DESCRIPTION
NOTE
SOURCE
560 subroutine add_general_mean_multi(xinc,listDesc) 561 real (kind=dp) , dimension(ncar,nd,ntnivmaxtotal) , intent(inout) :: xinc 562 type(INCIDENCE_TYPE) , dimension(ncar) , intent(inout) :: listDesc 563 integer :: ic,ntniv 564 call log_mess("multi incidence : ** Add General mean ** ",DEBUG_DEF) 565 566 do ic=1,ncar 567 listDesc(ic)%ntniv = listDesc(ic)%ntniv+1 568 ntniv = listDesc(ic)%ntniv 569 listDesc(ic)%nteff = listDesc(ic)%nteff+1 570 listDesc(ic)%desc(listDesc(ic)%nteff)%name="General Mean" 571 listDesc(ic)%desc(listDesc(ic)%nteff)%start=ntniv 572 listDesc(ic)%desc(listDesc(ic)%nteff)%end=ntniv 573 listDesc(ic)%desc(listDesc(ic)%nteff)%haveSubDesc=.true. 574 allocate (listDesc(ic)%desc(listDesc(ic)%nteff)%listSubDesc(1)) 575 listDesc(ic)%desc(listDesc(ic)%nteff)%listSubDesc(1)%name="General Mean"!//"["//trim(carac(ic))//"]" 576 listDesc(ic)%desc(listDesc(ic)%nteff)%listSubDesc(1)%start=ntniv 577 listDesc(ic)%desc(listDesc(ic)%nteff)%listSubDesc(1)%end=ntniv 578 xinc(ic,:,ntniv) = 1 579 end do 580 !incidenceDesc%nivdir(:,incidenceDesc%nteff) = 1 581 end subroutine add_general_mean_multi
add_polygenic_effect_multi
[ Top ] [ m_qtlmap_incidence_multi ] [ Subroutines ]
NAME
add_polygenic_effect_multi
DESCRIPTION
NOTE
SOURCE
592 subroutine add_polygenic_effect_multi(xinc,listDesc) 593 real (kind=dp) , dimension(ncar,nd,ntnivmaxtotal) , intent(inout) :: xinc 594 type(INCIDENCE_TYPE) , dimension(ncar) , intent(inout) :: listDesc 595 integer :: nkd 596 integer :: ip,jm,kd,nt,j,fem,nteff,ic 597 call log_mess("incidence : ** Multi : Add Polygenic effect to estimate **",DEBUG_DEF) 598 599 do ic=1,ncar 600 nt=listDesc(ic)%ntniv 601 nteff = listDesc(ic)%nteff 602 listDesc(ic)%desc(listDesc(ic)%nteff+1)%name="Sire polygenic effects" 603 listDesc(ic)%desc(listDesc(ic)%nteff+1)%start=nt+1 604 listDesc(ic)%desc(listDesc(ic)%nteff+1)%end=nt+np 605 listDesc(ic)%desc(listDesc(ic)%nteff+1)%haveSubDesc=.true. 606 allocate (listDesc(ic)%desc(listDesc(ic)%nteff+1)%listSubDesc(np)) 607 do ip=1,np 608 listDesc(ic)%desc(listDesc(ic)%nteff+1)%listSubDesc(ip)%name="Sire "//trim(pere(ip))!//"["//trim(carac(ic))//"]" 609 listDesc(ic)%desc(listDesc(ic)%nteff+1)%listSubDesc(ip)%start=nt+ip 610 listDesc(ic)%desc(listDesc(ic)%nteff+1)%listSubDesc(ip)%end=nt+ip 611 end do 612 listDesc(ic)%nteff = listDesc(ic)%nteff+1 613 listDesc(ic)%ntniv = listDesc(ic)%ntniv+np 614 615 if ( count(estime_multi(:))>0) then 616 listDesc(ic)%desc(listDesc(ic)%nteff+1)%name="Dam polygenic effects" 617 listDesc(ic)%desc(listDesc(ic)%nteff+1)%start=listDesc(ic)%desc(listDesc(ic)%nteff)%end+1 618 listDesc(ic)%desc(listDesc(ic)%nteff+1)%end=nt+np+count(estime_multi(:))!nbfemLoc(ic) 619 listDesc(ic)%desc(listDesc(ic)%nteff+1)%haveSubDesc=.true. 620 allocate (listDesc(ic)%desc(listDesc(ic)%nteff+1)%listSubDesc(count(estime_multi(:)))) 621 fem=0 622 do ip=1,np 623 do jm=nmp(ip)+1,nmp(ip+1) 624 if ( estime_multi(jm) ) then 625 fem=fem+1 626 listDesc(ic)%desc(listDesc(ic)%nteff+1)%listSubDesc(fem)%name="Dam "//trim(mere(jm))//" [Sire "//trim(pere(ip))//"]"!&//"["//trim(carac(ic))//"]" 627 listDesc(ic)%desc(listDesc(ic)%nteff+1)%listSubDesc(fem)%start=nt+np+fem 628 listDesc(ic)%desc(listDesc(ic)%nteff+1)%listSubDesc(fem)%end=nt+np+fem 629 end if 630 end do 631 end do 632 listDesc(ic)%nteff = listDesc(ic)%nteff+1 633 listDesc(ic)%ntniv = listDesc(ic)%ntniv+count(estime_multi(:))!nbfemLoc(ic) 634 else 635 636 end if 637 j=0 638 639 !initialisation 640 xinc(ic,:,nt+1:nt+np+count(estime_multi(:nm))) = 0 641 642 nkd=0 643 do ip=1,np 644 xinc(ic,:,nt+ip) = 0 645 do jm=nmp(ip)+1,nmp(ip+1) 646 do kd=ndm(jm)+1,ndm(jm+1) 647 if ( count(presentc(:,kd))==ncar ) then 648 nkd =nkd +1 649 xinc(ic,nkd,nt+ip) = 1 650 ! listDesc(ic)%nivdir(kd,nteff+1)= nt+ip 651 if (estime_multi(jm)) then 652 fem=count(estime_multi(:jm)) ! on compte le nombre estimable dans la famille de pere ip 653 xinc(ic,nkd,nt+np+fem) = 1 654 ! listDesc(ic)%nivdir(kd,nteff+2)= nt+np+fem 655 end if 656 end if 657 end do 658 end do 659 end do 660 end do !ic 661 end subroutine add_polygenic_effect_multi
add_qtleffect_multi
[ Top ] [ m_qtlmap_incidence_multi ] [ Subroutines ]
NAME
add_qtleffect_multi
DESCRIPTION
NOTE
SOURCE
672 subroutine add_qtleffect_multi(nteff,numqtl,chr,n,xinc,listDesc,mint,linear) 673 integer ,intent(in) :: nteff,numqtl,chr,n 674 real (kind=dp) , dimension(ncar,nd,ntnivmaxtotal) ,intent(inout) :: xinc 675 type(INCIDENCE_TYPE) , dimension(ncar) , intent(inout) :: listDesc 676 integer , intent(in) :: mint ! index of interaction to remove from incidence 677 logical , intent(in) :: linear 678 679 integer :: ic,ntniv,nteffstart,ip,fem,jm,l,s,u,nbint,nbt,jef 680 logical ,dimension(nfem) :: femInit 681 682 call log_mess("multi incidence : ** Add Qtl effect to estimate **",DEBUG_DEF) 683 684 do ic=1,ncar 685 686 nteffstart = listDesc(ic)%nteff 687 listDesc(ic)%nteff = listDesc(ic)%nteff+1 688 if ( count(estime_multi(:)) > 0 ) listDesc(ic)%nteff = listDesc(ic)%nteff+1 689 690 nbint=modele(ic,3) 691 nbt=3+modele(ic,1)+modele(ic,2) 692 listDesc(ic)%ntlev=1 693 do jef=1,nbint 694 if ( mint == jef) cycle ! this interaction is not adding 695 listDesc(ic)%ntlev=listDesc(ic)%ntlev*nlev(modele(ic,nbt+jef)) 696 end do 697 698 listDesc(ic)%desc(nteffstart+1)%name="Sire QTL effects ["//trim(str(numqtl))//"]" 699 listDesc(ic)%desc(nteffstart+1)%start=listDesc(ic)%ntniv+1 700 listDesc(ic)%ntniv_qtlsires(numqtl)=listDesc(ic)%desc(nteffstart+1)%start 701 listDesc(ic)%desc(nteffstart+1)%end=listDesc(ic)%ntniv+listDesc(ic)%ntlev*np 702 703 listDesc(ic)%desc(nteffstart+1)%haveSubDesc=.true. 704 allocate (listDesc(ic)%desc(nteffstart+1)%listSubDesc(np*listDesc(ic)%ntlev)) 705 u=0 706 do ip=1,np 707 l=0 708 do s=listDesc(ic)%ntniv+(ip-1)*listDesc(ic)%ntlev+1,listDesc(ic)%ntniv+ip*listDesc(ic)%ntlev 709 l=l+1 710 u=u+1 711 listDesc(ic)%desc(nteffstart+1)%listSubDesc(u)%name="Sire "//trim(pere(ip))//" Level "//trim(str(l))!&//"["//trim(carac(ic))//"]" 712 listDesc(ic)%desc(nteffstart+1)%listSubDesc(u)%start=s 713 listDesc(ic)%desc(nteffstart+1)%listSubDesc(u)%end=s 714 end do 715 end do 716 717 if ( count(estime_multi(:)) > 0 ) then 718 listDesc(ic)%desc(nteffstart+2)%name="Dam QTL effects ["//trim(str(numqtl))//"]" 719 listDesc(ic)%desc(nteffstart+2)%start=listDesc(ic)%desc(nteffstart+1)%end+1 720 listDesc(ic)%ntniv_qtldams(numqtl)=listDesc(ic)%desc(nteffstart+2)%start 721 femInit=.false. 722 l=0 723 do jm=1,nm 724 if (estime_multi(jm) ) then 725 fem=iam_multi(repfem(jm)) 726 if (femInit(fem)) then 727 cycle 728 end if 729 femInit(fem)=.true. 730 l=l+1 731 end if 732 end do 733 734 allocate (listDesc(ic)%desc(nteffstart+2)%listSubDesc(l*listDesc(ic)%ntlev)) 735 listDesc(ic)%desc(nteffstart+2)%end=listDesc(ic)%ntniv+listDesc(ic)%ntlev*(np+l) 736 listDesc(ic)%desc(nteffstart+2)%haveSubDesc=.true. 737 !fem=0 738 femInit=.false. 739 u=0 740 do jm=1,nm 741 if (estime_multi(jm) ) then 742 fem=iam_multi(repfem(jm)) 743 if (femInit(fem)) then 744 cycle 745 end if 746 femInit(fem)=.true. 747 l=0 748 do s=listDesc(ic)%desc(nteffstart+2)%start+(fem-1)*listDesc(ic)%ntlev,& 749 listDesc(ic)%desc(nteffstart+2)%start-1+fem*listDesc(ic)%ntlev 750 l=l+1 751 u=u+1 752 !fem=fem+1 753 listDesc(ic)%desc(nteffstart+2)%listSubDesc(u)%name="Dam "//trim(mere(jm))//" Level "//trim(str(l))!&//"["//trim(carac(ic))//"]" 754 listDesc(ic)%desc(nteffstart+2)%listSubDesc(u)%start=s 755 listDesc(ic)%desc(nteffstart+2)%listSubDesc(u)%end=s 756 end do 757 end if 758 end do 759 760 listDesc(ic)%ntniv = listDesc(ic)%desc(nteffstart+2)%end 761 else 762 listDesc(ic)%ntniv = listDesc(ic)%desc(nteffstart+1)%end 763 end if 764 end do !ic 765 766 call change_qtleffect_multi(nteff,numqtl,chr,n,xinc,listDesc,mint) 767 768 end subroutine add_qtleffect_multi
change_qtleffect_multi
[ Top ] [ m_qtlmap_incidence_multi ] [ Subroutines ]
NAME
change_qtleffect_multi
DESCRIPTION
NOTE
SOURCE
779 subroutine change_qtleffect_multi(nteff,iq,chr,n,xinc,listDesc,mint) 780 integer ,intent(in) :: nteff,iq,chr,n 781 real (kind=dp) , dimension(ncar,nd,ntnivmaxtotal) ,intent(inout) :: xinc 782 type(INCIDENCE_TYPE) , dimension(ncar) , intent(inout) :: listDesc 783 integer ,intent(in) :: mint ! index of interaction to remove from incidence 784 785 !local 786 integer :: nbef,nbco,nbint,nbt,jef,nivint,ntlev,ip,jm,kd,ntt,km,ntniv,nkd,ic 787 real (kind=dp) ,dimension(:,:),allocatable:: pp,pm 788 789 allocate (pp(nd,ncar)) 790 allocate (pm(nd,ncar)) 791 792 nkd=0 793 794 do ip=1,np 795 do jm=nmp(ip)+1,nmp(ip+1) 796 !if ( linear ) then 797 !call getprob_lin(chr,n,ic,ip,jm,ndm(jm)+1,ndm(jm+1),pp,pm) 798 !else 799 do ic=1,ncar 800 call getglobalprob(chr,n,ic,ip,jm,ndm(jm)+1,ndm(jm+1),pp(:,ic),pm(:,ic)) 801 end do 802 !end if 803 do kd=ndm(jm)+1,ndm(jm+1) 804 if ( count(presentc(:,kd)) == ncar ) then 805 nkd=nkd+1 806 807 ! si des effets fixes sont en interaction avec l'allele paternel au qtl 808 ! il faut estimer les effets qtl pour chacun des niveaux concernes 809 ! 810 ! on commence donc par creer un effet composites regroupant l'ensemble de 811 ! ces effets 812 ! 813 do ic=1,ncar 814 ntniv = listDesc(ic)%desc(nteff)%start -1 815 nivint=1 816 ntlev=1 817 nbt=3+modele(ic,1)+modele(ic,2) 818 !nt=ntniv 819 do jef=1,modele(ic,3) 820 if ( mint == jef) cycle ! this interaction is not adding 821 nivint=nivint+ntlev*(nivx(kd,modele(ic,nbt+jef))-1) 822 ntlev=ntlev*nlev(modele(ic,nbt+jef)) 823 end do 824 ntt=ntniv+nivint+ntlev*(ip-1) 825 xinc(ic,nkd,ntt)=pp(kd,ic) 826 827 ! incidenceDesc%nivdir(kd,nteff+1)=ntt 828 ! print *,'nt:',nivdir(kd,nteff+1),' vp:',pp(kd),' nteff:',nteff,' kd:',kd 829 830 if(estime_multi(jm)) then 831 km=iam_multi(repfem(jm)) 832 ntt= ntniv+ntlev*np+nivint+ntlev*(km-1) 833 xinc(ic,nkd,ntt)=pm(kd,ic) 834 ! incidenceDesc%nivdir(kd,nteff+2)=ntt 835 ! print *,'nt:',nivdir(kd,nteff+2),' vp:',pm(kd) 836 end if 837 end do !ic 838 end if !presentc 839 end do 840 end do 841 end do 842 843 844 deallocate (pp) 845 deallocate (pm) 846 ! a modifier 847 ! --> on devrait integerer l intialisation des pdds dans la boucle precedente 848 ! attention les interaction qtl ne sont pas encore pris en compte 849 ! if (.not. linear) then 850 call pdd_at_position_multi(listDesc(1)%dataset,iq,chr,n) 851 ! end if 852 end subroutine change_qtleffect_multi
confusion_multi_type1
[ Top ] [ m_qtlmap_incidence_multi ] [ Subroutines ]
NAME
confusion_multi_type1
DESCRIPTION
NOTE
SOURCE
1050 subroutine confusion_multi_type1(listDesc,xxx,workstruct) 1051 type(INCIDENCE_TYPE) , dimension(ncar) , intent(inout) :: listDesc 1052 real (kind=dp) ,dimension(ncar,ntnivmaxtotal,ntnivmaxtotal) ,intent(in) :: xxx 1053 type(INCIDENCE_GEN_STRUCT) ,intent(inout) :: workstruct 1054 1055 real (kind=dp) ,dimension(:,:,:),allocatable :: xcor 1056 type(CORR_ALERT_TYPE) ,dimension(:),allocatable :: alerts 1057 type(CORR_ALERT_TYPE) ,dimension(:),allocatable :: as 1058 integer :: ic,ialert,sizealert,eff,i,j,starteff 1059 1060 allocate (xcor(ncar,ntnivmaxtotal,ntnivmaxtotal)) 1061 1062 sizealert=0 1063 do ic=1,ncar 1064 call get_correlations(listDesc(ic),xxx(ic,:,:),xcor(ic,:,:)) 1065 sizealert=sizealert+np*namest_multi*listDesc(ic)%ntlev*listDesc(ic)%nbnivest 1066 end do 1067 1068 workstruct%corrmax=0.d0 1069 1070 1071 1072 allocate (alerts(sizealert)) 1073 allocate (as(sizealert*workstruct%nqtl)) 1074 1075 workstruct%nalert=0 1076 starteff = workstruct%listnteff(workstruct%nqtl)+1 1077 if (namest_multi>0) starteff = starteff+1 1078 1079 do ic=1,ncar 1080 do eff=starteff,listDesc(ic)%nteff 1081 call log_mess("computing confusion beetween qtls and:"//listDesc(ic)%desc(eff)%name,VERBOSE_DEF) 1082 call confusion_between_qtl_and_effect(listDesc(ic),xcor(ic,:,:),& 1083 eff,sizealert,alerts,ialert,workstruct%corrmax) 1084 if (ialert>0) then 1085 ! print *,"ialert:",alerts(1)%corr,alerts(2)%corr,alerts(3)%corr 1086 as(workstruct%nalert+1:workstruct%nalert+ialert)=alerts(:ialert) 1087 workstruct%nalert=workstruct%nalert+ialert 1088 end if 1089 end do 1090 end do 1091 1092 allocate (workstruct%alertQtl(workstruct%nalert)) 1093 workstruct%alertQtl=as(:workstruct%nalert) 1094 1095 deallocate (alerts) 1096 deallocate (as) 1097 deallocate (xcor) 1098 1099 1100 end subroutine confusion_multi_type1
end_incidence_multi
[ Top ] [ m_qtlmap_incidence_multi ] [ Subroutines ]
NAME
end_incidence_multi
DESCRIPTION
NOTE
SOURCE
295 subroutine end_incidence_multi(listDesc) 296 type(INCIDENCE_TYPE) ,dimension(ncar), intent(inout) :: listDesc 297 integer :: ic 298 299 deallocate (estime_multi) 300 deallocate (iam_multi) 301 deallocate (estfem_multi) 302 303 call end_dataset(listDesc(1)%dataset) 304 deallocate (listDesc(1)%dataset) 305 do ic=2,ncar 306 listDesc(ic)%dataset=>null() 307 listDesc(ic)%borni=>null() 308 listDesc(ic)%borns=>null() 309 listDesc(ic)%par=>null() 310 end do 311 312 do ic=1,ncar 313 call end_incidence(listDesc(ic)) 314 end do 315 316 end subroutine end_incidence_multi
gen_loop_opti_nqtl_multi
[ Top ] [ m_qtlmap_incidence_multi ] [ Subroutines ]
NAME
gen_loop_opti_nqtl_multi
DESCRIPTION
NOTE
SOURCE
3380 recursive subroutine gen_loop_opti_nqtl_multi(iqtl, & 3381 nqtl, & 3382 curPos, & 3383 workstruct,& 3384 xinc, & 3385 listDesc , & 3386 sigsquareEstimeMax,& 3387 rhoiestimMax, & 3388 BestimMax, & 3389 lrtsol, & 3390 printpercent, & 3391 FUNCT_MODEL) 3392 integer , intent(in) :: nqtl ! under hypothesis nqtl 3393 integer , intent(in) :: iqtl ! iteration of inside loop 3394 type(POSITION_LRT_INCIDENCE) ,intent(inout) :: curPos 3395 type(INCIDENCE_GEN_STRUCT) ,intent(inout) :: workstruct 3396 real (kind=dp) , dimension(ncar,nd,ntnivmaxtotal) , intent(inout) :: xinc ! incidence matrix 3397 type(INCIDENCE_TYPE) , dimension(ncar) , intent(inout) :: listDesc ! description of the incidence matrix 3398 real (kind=dp) , dimension(ncar,ntnivmaxtotal) , intent(inout) :: BestimMax ! estimation of parameter on the maximum finded 3399 real (kind=dp) , dimension(ncar,ncar) ,intent(out) :: rhoiestimMax 3400 real (kind=dp) , dimension(ncar,np) , intent(inout) :: sigsquareEstimeMax ! estimation of parameter on the maximum finded 3401 type(TYPE_LRT_SOLUTION) ,intent(inout) :: lrtsol 3402 logical , intent(in) :: printpercent 3403 external :: FUNCT_MODEL 3404 3405 3406 real (kind=dp) , dimension(ncar,ntnivmaxtotal) :: Bestim 3407 real (kind=dp) , dimension(ncar,np) :: sigsquareEstime 3408 real (kind=dp) , dimension(ncar,ncar) :: rhoiestim 3409 integer :: n,hypothesis,chr,chr_start,iip,ifem,ii,lnpo,chr1,n1,ip 3410 logical hypotheseOne, hypotheseTwo 3411 real (kind=dp) ,dimension(:,:,:),pointer :: xxx 3412 3413 allocate (xxx(ncar,ntnivmaxtotal,ntnivmaxtotal)) 3414 ! print *," ************** loop_opti_nqtl ***********************" 3415 ! incidenceDescSave = incidenceDesc 3416 if ( iqtl == nqtl ) then 3417 hypotheseOne = ( lrtsol%nqtl == 1 ) 3418 hypotheseTwo = ( lrtsol%nqtl == 2 ) 3419 if (nqtl >= 2 ) then 3420 chr_start = curPos%listChr(nqtl-1) 3421 else 3422 chr_start = 1 3423 end if 3424 3425 do chr=chr_start,nchr 3426 ! Step above the chromosome 3427 curPos%listChr(nqtl)=chr 3428 if ( (nqtl >= 2) .and. (chr == chr_start ) ) then 3429 n = curPos%listN(nqtl-1) 3430 else 3431 n=0 3432 end if 3433 3434 do while (n < get_npo(chr)) 3435 n=n+1 3436 curPos%listN(nqtl)=n 3437 ! if (nqtl == 3)print *,"chr:",chr," n:",n 3438 ! 3439 !add qtl effect/interaction at position n to estim 3440 call change_qtleffect_multi(workstruct%listnteff(iqtl),iqtl,chr,n,xinc,listDesc,0) 3441 !call debug_write_incidence(xinc,incidenceDesc) 3442 ! do ii=1,ncar 3443 ! call debug_write_incidence(xinc(ii,:,:),listDesc(ii)) 3444 ! end do 3445 !call model 3446 call FUNCT_MODEL(xinc,listDesc,curPos,workstruct,sigsquareEstime,rhoiestim,Bestim,.false.,.false.,xxx,.false.) 3447 3448 if ( hypotheseOne ) then 3449 lrtsol%lrt1(chr,n)=curPos%lrt(1) 3450 iip=1 3451 do ip=1,np 3452 if ( listDesc(1)%vecsol(1+ip)) then 3453 iip=iip+1 3454 ! ATTENTION VALIDE DANS LE CAS SANS INTERACTION AVEC QTL 3455 lrtsol%pater_eff(chr,ip,n)=Bestim(1,iip) 3456 end if 3457 end do 3458 ifem=0 3459 do ii=1,nm 3460 if ( estime_multi(ii) ) then 3461 ifem=ifem+1 3462 if ( listDesc(1)%vecsol(np+1+ifem) ) then 3463 iip=iip+1 3464 ! ATTENTION VALIDE DANS LE CAS SANS INTERACTION AVEC QTL 3465 lrtsol%mater_eff(chr,ifem,n)=Bestim(1,iip) 3466 end if 3467 end if 3468 end do 3469 end if 3470 if ( hypotheseTwo ) then 3471 ! print *,n,size(lrtsol%lrt1_2,1),listChr(nqtl-1),listChr(nqtl-2),& 3472 ! size(lrtsol%lrt1_2,2),size(lrtsol%lrt1_2,3),size(lrtsol%lrt1_2,4) 3473 chr1=curPos%listChr(1) 3474 n1 = curPos%listN(1) 3475 lrtsol%lrt1_2(chr1,chr,n1,n)=curPos%lrt(2) 3476 lrtsol%lrt0_2(chr1,chr,n1,n)=curPos%lrt(1) 3477 3478 end if 3479 !print *,curPos%listDx,curPos%lrt(nqtl) 3480 !keep maximum 3481 if(lrtsol%lrtmax(nqtl-1) < curPos%lrt(nqtl)) then 3482 lrtsol%nxmax=curPos%listN 3483 lrtsol%chrmax=curPos%listChr 3484 BestimMax = Bestim 3485 sigsquareEstimeMax=sigsquareEstime 3486 rhoiestimMax=rhoiestim 3487 lrtsol%lrtmax=curPos%lrt 3488 !print *,"max sigsquareestime:",workstruct%sigsquareEstime 3489 end if 3490 end do 3491 end do 3492 else 3493 3494 if ( iqtl-2 >= 0 ) then 3495 chr_start=curPos%listChr(iqtl-1) 3496 else 3497 chr_start=1 3498 end if 3499 3500 do chr=chr_start,nchr 3501 ! Step above the chromosome 3502 3503 if ( iqtl-2 >= 0 .and. chr == chr_start) then 3504 n=curPos%listN(iqtl-1) 3505 else 3506 n=0 3507 end if 3508 3509 curPos%listChr(iqtl)=chr 3510 lnpo=get_npo(chr) 3511 do while (n < lnpo) !ix=startx,ilong,pas 3512 if (printpercent) call log_mess( trim(str((float(n)/float(lnpo))*100.d0))//"%", VERBOSE_DEF ) 3513 n=n+1 3514 curPos%listN(iqtl)=n 3515 !add qtl effect/interaction at position n to estim 3516 call change_qtleffect_multi(workstruct%listnteff(iqtl),iqtl,chr,n,xinc,listDesc,0) 3517 3518 ! ****recursivite **** 3519 call gen_loop_opti_nqtl_multi(iqtl+1,nqtl,curPos,workstruct,& 3520 xinc,listDesc,sigsquareEstimeMax,rhoiestimMax,BestimMax,lrtsol,.false.,FUNCT_MODEL) 3521 3522 end do 3523 end do 3524 end if 3525 deallocate (xxx) 3526 3527 3528 end subroutine gen_loop_opti_nqtl_multi
gen_opti_nqtl_multi
[ Top ] [ m_qtlmap_incidence_multi ] [ Subroutines ]
NAME
gen_opti_nqtl_multi
DESCRIPTION
NOTE
SOURCE
3539 subroutine gen_opti_nqtl_multi( nqtl, & 3540 curPosMax, & 3541 workstruct,& 3542 xinc, & 3543 listDescMax , & 3544 sigsquareEstimeMax,& 3545 rhoiestimMax, & 3546 BestimMax, & 3547 lrtsol, & 3548 FUNCT_MODEL) 3549 integer , intent(in) :: nqtl ! under hypothesis nqtl 3550 type(POSITION_LRT_INCIDENCE) ,intent(inout) :: curPosMax 3551 type(INCIDENCE_GEN_STRUCT) ,intent(inout) :: workstruct 3552 real (kind=dp) , dimension(ncar,nd,ntnivmaxtotal) , intent(inout) :: xinc ! incidence matrix 3553 type(INCIDENCE_TYPE) , dimension(ncar) , intent(inout) :: listDescMax ! description of the incidence matrix 3554 real (kind=dp) , dimension(ncar,ntnivmaxtotal) , intent(inout) :: BestimMax ! estimation of parameter on the maximum finded 3555 real (kind=dp) , dimension(ncar,ncar) ,intent(out) :: rhoiestimMax 3556 real (kind=dp) , dimension(ncar,np) , intent(inout) :: sigsquareEstimeMax ! estimation of parameter on the maximum finded 3557 type(TYPE_LRT_SOLUTION) ,intent(inout) :: lrtsol 3558 external :: FUNCT_MODEL 3559 3560 3561 real (kind=dp) , dimension(ncar,ntnivmaxtotal) :: Bestim 3562 real (kind=dp) , dimension(ncar,np) :: sigsquareEstime 3563 real (kind=dp) , dimension(ncar,ncar) :: rhoiestim 3564 3565 real (kind=dp) ,dimension(:,:,:),pointer :: xxx 3566 real (kind=dp) , dimension(:,:,:),allocatable :: xinc2 3567 3568 logical :: hypotheseOne,hypotheseTwo,ok 3569 integer :: i,chr,nGL,nGLTotal,nGLFact,InGL,nlinMax,nlin,nmax(nqtl),chmax(nqtl) 3570 integer :: nlinValide,iqtl,ic,ii,ifem,iip,ip 3571 integer , dimension(:,:,:) ,allocatable :: bestim_save,sigsq_save,rhoi_save 3572 integer , dimension(:,:) ,allocatable :: lchr,lnpos 3573 real(kind=dp) , dimension(:,:) ,allocatable :: lrt_save 3574 type(POSITION_LRT_INCIDENCE) :: curPos 3575 type(INCIDENCE_TYPE) , dimension(ncar) :: listDesc 3576 3577 hypotheseOne = ( lrtsol%nqtl == 1 ) 3578 hypotheseTwo = ( lrtsol%nqtl == 2 ) 3579 3580 nGL = 0 3581 3582 !Nombre de point sur le groupe de liaison 3583 !on cherche le nombre de point dependant de l hyp nqtl du numbre de point sur le groupe de liason 3584 3585 do chr=1,nchr 3586 nGL=nGL + get_npo(chr) 3587 end do 3588 3589 nGLTotal=1 3590 do iqtl=1,nqtl 3591 nGLTotal=nGLTotal*nGL 3592 end do 3593 3594 curPosMax%listN(nqtl)=0 3595 curPosMax%listN(1:nqtl-1)=1 3596 curPosMax%listChr=1 3597 3598 allocate (lchr(nGLTotal,nqtl)) 3599 allocate (lnpos(nGLTotal,nqtl)) 3600 3601 nlinValide=0 3602 do nlin=1,nGLTotal 3603 i=nqtl 3604 ok=.false. 3605 do while ( (.not. ok) .and. (i>=1)) 3606 curPosMax%listN(i) = curPosMax%listN(i)+1 3607 if ( curPosMax%listN(i) > get_npo(curPosMax%listChr(i)) ) then 3608 if ( nchr > curPosMax%listChr(i) ) then 3609 curPosMax%listN(i)=1 3610 curPosMax%listChr(i)=curPosMax%listChr(i)+1 3611 ok=.true. 3612 else 3613 curPosMax%listN(i)=1 3614 curPosMax%listChr(i)=1 3615 i=i-1 3616 end if 3617 else 3618 ok = .true. 3619 end if 3620 end do 3621 ok=.true. 3622 do iqtl=2,nqtl 3623 ok = ok .and. ( curPosMax%listN(iqtl-1)<=curPosMax%listN(iqtl) ) 3624 end do 3625 3626 if ( ok ) then 3627 nlinValide=nlinValide+1 3628 lchr(nlinValide,:)=curPosMax%listChr 3629 lnpos(nlinValide,:)=curPosMax%listN 3630 end if 3631 end do 3632 nGLTotal=nlinValide 3633 3634 allocate (bestim_save(nGLTotal,ncar,ntnivmaxtotal)) 3635 allocate (sigsq_save(nGLTotal,ncar,np)) 3636 allocate (rhoi_save(nGLTotal,ncar,ncar)) 3637 allocate (lrt_save(nGLTotal,nqtl)) 3638 lrt_save=0.d0 3639 !$OMP PARALLEL DEFAULT(SHARED) & 3640 !$OMP PRIVATE(curPos,i,ic,iqtl,ok,rhoiestim,sigsquareEstime,Bestim,xxx,ip,iip,ifem,listDesc,xinc2) 3641 call init_position (nqtl,curPos) 3642 allocate (xxx(ncar,ntnivmaxtotal,ntnivmaxtotal)) 3643 allocate (xinc2(ncar,nd,ntnivmaxtotal)) 3644 xinc2=xinc 3645 !initialisation of incidence matrix 3646 do ic=1,ncar 3647 call copy_incidence_desc (listDescMax(ic),listDesc(ic)) 3648 end do 3649 3650 curPos%lrt(nqtl)=0 3651 curPos%listN(nqtl)=0 3652 curPos%listN(1:nqtl-1)=1 3653 curPos%listChr=1 3654 3655 !$OMP DO 3656 do nlin=1,nGLTotal 3657 call log_mess( trim(str((float(nlin)/float(nGLTotal))*100.d0))//"%", INFO_DEF ) 3658 3659 curPos%listN=lnpos(nlin,:) 3660 curPos%listChr=lchr(nlin,:) 3661 3662 do iqtl=1,nqtl 3663 call change_qtleffect_multi(workstruct%listnteff(iqtl),iqtl,curPos%listChr(iqtl),curPos%listN(iqtl),xinc2,listDesc,0) 3664 end do 3665 ! do ic=1,ncar 3666 ! call debug_write_incidence(xinc2(ic,:,:),listDesc(ic)) 3667 ! end do 3668 3669 call FUNCT_MODEL(xinc2,listDesc,curPos,workstruct,sigsquareEstime,rhoiestim,Bestim,.false.,.false.,xxx,.false.) 3670 3671 3672 3673 !save value 3674 bestim_save(nlin,:,:)=Bestim 3675 sigsq_save(nlin,:,:)=sigsquareEstime 3676 rhoi_save(nlin,:,:)=rhoiestim 3677 lrt_save(nlin,:)=curPos%lrt 3678 3679 if ( hypotheseOne ) then 3680 lrtsol%lrt1(curPos%listChr(1),curPos%listN(1))=curPos%lrt(1) 3681 iip=1 3682 do ip=1,np 3683 if ( listDesc(1)%vecsol(1+ip)) then 3684 iip=iip+1 3685 ! ATTENTION VALIDE DANS LE CAS SANS INTERACTION AVEC QTL 3686 lrtsol%pater_eff(curPos%listChr(1),ip,curPos%listN(1))=Bestim(1,iip) 3687 end if 3688 end do 3689 ifem=0 3690 do ii=1,nm 3691 if ( estime_multi(ii)) then 3692 ifem=ifem+1 3693 if ( listDesc(1)%vecsol(np+1+ifem) ) then 3694 iip=iip+1 3695 ! ATTENTION VALIDE DANS LE CAS SANS INTERACTION AVEC QTL 3696 lrtsol%mater_eff(curPos%listChr(1),ifem,curPos%listN(1))=Bestim(1,iip) 3697 end if 3698 end if 3699 end do 3700 end if 3701 if ( hypotheseTwo ) then 3702 lrtsol%lrt1_2(curPos%listChr(1),curPos%listChr(2),curPos%listN(1),curPos%listN(2))=curPos%lrt(2) 3703 lrtsol%lrt0_2(curPos%listChr(1),curPos%listChr(2),curPos%listN(1),curPos%listN(2))=curPos%lrt(1) 3704 end if 3705 ! end if !ok 3706 end do ! nlin 3707 !$OMP END DO 3708 call end_position (curPos) 3709 deallocate (xxx) 3710 deallocate (xinc2) 3711 do ic=1,ncar 3712 call release_copy_incidence_desc(listDesc(ic)) 3713 end do 3714 !$OMP END PARALLEL 3715 3716 nlinMax=1 3717 lrtsol%lrtmax(0)=lrt_save(nlinMax,nqtl) 3718 nmax=1 3719 chmax=1 3720 3721 call init_position (nqtl,curPos) 3722 3723 curPos%listN(nqtl)=0 3724 curPos%listN(1:nqtl-1)=1 3725 curPos%listChr=1 3726 3727 do nlin=1,nGLTotal 3728 curPos%listN=lnpos(nlin,:) 3729 curPos%listChr=lchr(nlin,:) 3730 3731 if(lrtsol%lrtmax(nqtl-1) < lrt_save(nlin,nqtl)) then 3732 nmax=curPos%listN 3733 chmax=curPos%listChr 3734 nlinMax=nlin 3735 lrtsol%lrtmax=lrt_save(nlin,:) 3736 end if 3737 end do 3738 3739 do iqtl=1,nqtl 3740 lrtsol%nxmax(iqtl-1)=nmax(iqtl) 3741 lrtsol%chrmax(iqtl-1)=chmax(iqtl) 3742 BestimMax = bestim_save(nlinMax,:,:) 3743 sigsquareEstimeMax=sigsq_save(nlinMax,:,:) 3744 rhoiestimMax=rhoi_save(nlinMax,:,:) 3745 curPosMax%listN=curPos%listN 3746 curPosMax%listChr=curPos%listChr 3747 end do 3748 3749 deallocate (lchr) 3750 deallocate (lnpos) 3751 3752 call end_position (curPos) 3753 3754 deallocate (bestim_save) 3755 deallocate (sigsq_save) 3756 deallocate (rhoi_save) 3757 deallocate (lrt_save) 3758 3759 end subroutine gen_opti_nqtl_multi
get_inv_residual_covariance_matrix
[ Top ] [ m_qtlmap_incidence_multi ] [ Subroutines ]
NAME
get_inv_residual_covariance_matrix
DESCRIPTION
NOTE
SOURCE
1658 subroutine get_inv_residual_covariance_matrix(ip,n,x,vci,determ,ok) 1659 integer , intent(in) :: ip,n 1660 real (kind=dp) ,dimension(n) , intent(in) :: x 1661 real(kind=dp),dimension(ncar,ncar),intent(out) :: VCI 1662 real(kind=dp), intent(out) :: determ 1663 logical ,intent(out) :: ok 1664 real(kind=dp),dimension(ncar,ncar) :: rhoc 1665 real(kind=dp),dimension(ncar+1,ncar) :: VCIInverse 1666 1667 real(kind=dp) :: rh 1668 integer :: ic,jc,k,irank 1669 1670 ok=.true. 1671 !RHO(I,J) avec i/=j => x((i-1)*(ncar-i+1)) 1672 k=0 1673 do ic=2,ncar 1674 do jc=1,ic-1 1675 k=k+1 1676 ! *** RESCALE OFI *** 1677 !rhoc(jc,ic)=x(ncar*np+k) 1678 rh = x(ncar*np+k) 1679 rhoc(jc,ic) = ((dexp(rh)-1.d0)/(dexp(rh)+1.d0)) 1680 rhoc(ic,jc)=rhoc(jc,ic) 1681 end do 1682 end do 1683 1684 !VCi : residual covariance matrix from sire I 1685 ! | SIGi1^2 SIGi1*SIGi2*RHO(1,2) .... SIGi1*SIGip*RHO(1,p) 1686 ! | .. SIGi2^2 1687 ! VCi | .. ... 1688 ! | .. ... 1689 ! | .. 1690 ! | SIGi1*SIGip*RHO(1,p) SIGip^2 1691 1692 do ic=1,ncar 1693 do jc=ic,ncar 1694 VCI(ic,jc)=x((ic-1)*np+ip)*x((jc-1)*np+ip) 1695 if ( ic /= jc ) then 1696 VCI(ic,jc)=VCI(ic,jc)*rhoc(ic,jc) 1697 VCI(jc,ic)=VCI(ic,jc) 1698 end if 1699 end do 1700 end do 1701 1702 !Calcul Determinant de la matrice de covariance residuelle du pere ip 1703 !-------------------------------------------------------------------- 1704 irank=0 1705 determ=0 1706 ! call MATH_QTLMAP_F03ABF(VCI,ncar,ncar,determ,irank) 1707 ! if (irank /= 0) then 1708 ! print *,"mauvais determinant matrice" 1709 ! ok=.false. 1710 ! return 1711 ! end if 1712 !Calcul Inverse de la matrice de covariance residuelle du pere ip 1713 !----------------------------------------------------------------- 1714 VCIInverse(:ncar,:ncar)=VCI 1715 CALL MATH_QTLMAP_INVDETMATSYM(ncar,VCIInverse,ncar+1,determ,irank) 1716 ! call MATH_QTLMAP_F01ADF(ncar,VCIInverse,ncar+1,irank) 1717 if (irank /= 0) then 1718 print *,"mauvaise inversion matrice" 1719 ok=.false. 1720 return 1721 end if 1722 1723 do ic=1,ncar 1724 do jc=1,ic 1725 VCI(ic,jc)=VCIInverse(ic+1,jc) 1726 VCI(jc,ic)=VCI(ic,jc) 1727 end do 1728 end do 1729 1730 1731 1732 end subroutine get_inv_residual_covariance_matrix
get_inv_residual_covariance_matrix_cd
[ Top ] [ m_qtlmap_incidence_multi ] [ Subroutines ]
NAME
get_inv_residual_covariance_matrix_cd
DESCRIPTION
NOTE
SOURCE
2006 subroutine get_inv_residual_covariance_matrix_cd(ip,kd,n,x,vci,determ,ok) 2007 integer , intent(in) :: ip,kd,n 2008 real (kind=dp) ,dimension(n) , intent(in) :: x 2009 real(kind=dp),dimension(ncar,ncar),intent(out) :: VCI 2010 real(kind=dp), intent(out) :: determ 2011 logical ,intent(out) :: ok 2012 real(kind=dp),dimension(ncar,ncar) :: rhoc 2013 real(kind=dp),dimension(ncar+1,ncar) :: VCIInverse 2014 2015 integer , parameter :: MAX_COUNT = 200 2016 2017 real(kind=dp) :: rh 2018 integer :: ic,jc,k,irank,countt 2019 2020 ok=.true. 2021 !RHO(I,J) avec i/=j => x((i-1)*(ncar-i+1)) 2022 k=0 2023 do ic=2,ncar 2024 do jc=1,ic-1 2025 k=k+1 2026 ! *** RESCALE OFI *** 2027 !rhoc(jc,ic)=x(ncar*np+k) 2028 rh = x(ncar*np+k) 2029 rhoc(jc,ic) = ((dexp(rh)-1.d0)/(dexp(rh)+1.d0)) 2030 rhoc(ic,jc)=rhoc(jc,ic) 2031 end do 2032 end do 2033 2034 !VCi : residual covariance matrix from sire I 2035 ! | SIGi1^2 SIGi1*SIGi2*RHO(1,2) .... SIGi1*SIGip*RHO(1,p) 2036 ! | .. SIGi2^2 2037 ! VCi | .. ... 2038 ! | .. ... 2039 ! | .. 2040 ! | SIGi1*SIGip*RHO(1,p) SIGip^2 2041 !print *,'kd:',kd,cd(1,kd),my_corcd(1,2,kd) 2042 do ic=1,ncar 2043 do jc=ic,ncar 2044 VCI(ic,jc)=x((ic-1)*np+ip)*x((jc-1)*np+ip) 2045 if ( ic /= jc ) then 2046 VCI(ic,jc)=VCI(ic,jc)*rhoc(ic,jc)*datasetUser%corcd(ic,jc,kd) 2047 VCI(jc,ic)=VCI(ic,jc) 2048 else 2049 VCI(ic,jc)=VCI(ic,jc)*cd(ic,kd) 2050 end if 2051 end do 2052 end do 2053 2054 !Calcul Determinant de la matrice de covariance residuelle du pere ip 2055 !-------------------------------------------------------------------- 2056 irank=0 2057 determ=0 2058 ! call MATH_QTLMAP_F03ABF(VCI,ncar,ncar,determ,irank) 2059 ! if (irank /= 0) then 2060 ! print *,"mauvais determinant matrice" 2061 ! ok=.false. 2062 ! return 2063 ! end if 2064 !Calcul Inverse de la matrice de covariance residuelle du pere ip 2065 !----------------------------------------------------------------- 2066 VCIInverse(:ncar,:ncar)=VCI 2067 irank = 1 2068 countt=0 2069 2070 ! do while (irank /= 0 .and. countt <= MAX_COUNT) 2071 do while (irank /= 0) 2072 CALL MATH_QTLMAP_INVDETMATSYM(ncar,VCIInverse,ncar+1,determ,irank) 2073 ! call MATH_QTLMAP_F01ADF(ncar,VCIInverse,ncar+1,irank) 2074 if (irank /= 0) then 2075 countt=countt+1 2076 ! print *,"mauvaise inversion matrice:",irank,' determ:',determ 2077 ! ok=.false. 2078 ! print *," *************** " 2079 ! print *," ** START VCI ** " 2080 ! do ic=1,ncar 2081 ! do jc=ic,ncar 2082 ! print *,x((ic-1)*np+ip)*x((jc-1)*np+ip) 2083 ! end do 2084 ! end do 2085 ! print *," ** RHOC ** " 2086 ! do ic=1,ncar 2087 ! print *,rhoc(ic,:) 2088 ! end do 2089 ! print *," ** VCI **" 2090 ! do ic=1,ncar 2091 ! print *,VCI(ic,:) 2092 ! end do 2093 ! print *," ** CORCD pour KD:",kd," **" 2094 ! do ic=1,ncar 2095 ! do jc=ic,ncar 2096 ! print *,ic,jc,datasetUser%corcd(ic,jc,kd) 2097 ! end do 2098 ! end do 2099 ! print *," ** CD pour KD:",kd," **" 2100 ! do ic=1,ncar 2101 ! print *,ic,cd(ic,kd) 2102 ! end do 2103 VCIInverse(:ncar,:ncar)=VCI 2104 do ic=1,ncar 2105 VCIInverse(ic,ic)=VCIInverse(ic,ic)+0.1d0*real(countt) 2106 end do 2107 ! stop 2108 ! return 2109 end if 2110 end do 2111 2112 ! if ( countt > MAX_COUNT ) then 2113 ! print *,"mauvaise inversion matrice countt:",countt 2114 ! return 2115 ! end if 2116 2117 do ic=1,ncar 2118 do jc=1,ic 2119 VCI(ic,jc)=VCIInverse(ic+1,jc) 2120 VCI(jc,ic)=VCI(ic,jc) 2121 end do 2122 end do 2123 2124 2125 2126 end subroutine get_inv_residual_covariance_matrix_cd
get_inv_residual_covariance_matrix_LU
[ Top ] [ m_qtlmap_incidence_multi ] [ Subroutines ]
NAME
get_inv_residual_covariance_matrix_LU
DESCRIPTION
NOTE
SOURCE
2546 subroutine get_inv_residual_covariance_matrix_LU(ip,n,x,vci,determ,ok) 2547 integer , intent(in) :: ip,n 2548 real (kind=dp) ,dimension(n) , intent(in) :: x 2549 real(kind=dp),dimension(ncar,ncar),intent(out) :: VCI 2550 real(kind=dp) , intent(out) :: determ 2551 logical ,intent(out) :: ok 2552 real(kind=dp),dimension(ncar,ncar) :: mat_L,V 2553 real(kind=dp),dimension(ncar) :: mat_H 2554 real(kind=dp),dimension(ncar,ncar) :: VCITemp 2555 real(kind=dp),dimension(ncar+1,ncar) :: VCIInverse 2556 2557 integer :: ic,jc,k,irank 2558 2559 ok=.true. 2560 !RHO(I,J) avec i/=j => x((i-1)*(ncar-i+1)) 2561 k=0 2562 mat_L=0.d0 2563 do ic=1,ncar 2564 do jc=1,ic 2565 k=k+1 2566 mat_L(ic,jc) = x(ncar*(np-1)+k) 2567 end do 2568 end do 2569 2570 V = matmul(mat_L,transpose(mat_L)) 2571 mat_H=0.d0 2572 ! on fixe pour ip=1 H 2573 if ( ip == 1) then 2574 mat_H(:)=1.d0 2575 else 2576 ! do ip=2,np 2577 do jc=1,ncar 2578 mat_H(jc)=x((ip-2)*ncar+jc) 2579 end do 2580 ! end do 2581 end if 2582 2583 determ=0 2584 2585 do ic=1,ncar 2586 do jc=1,ncar 2587 VCIInverse(ic,jc) = V(ic,jc)*mat_H(ic)*mat_H(jc) 2588 end do 2589 end do 2590 2591 !Calcul Determinant de la matrice de covariance residuelle du pere ip 2592 !-------------------------------------------------------------------- 2593 irank=0 2594 2595 ! call MATH_QTLMAP_F03ABF(VCI,ncar,ncar,determ,irank) 2596 ! if (irank /= 0) then 2597 ! print *,"mauvais determinant matrice" 2598 ! ok=.false. 2599 ! return 2600 ! end if 2601 !Calcul Inverse de la matrice de covariance residuelle du pere ip 2602 !----------------------------------------------------------------- 2603 !VCIInverse(:ncar,:ncar)=VCI 2604 CALL MATH_QTLMAP_INVDETMATSYM(ncar,VCIInverse,ncar+1,determ,irank) 2605 ! call MATH_QTLMAP_F01ADF(ncar,VCIInverse,ncar+1,irank) 2606 if (irank /= 0) then 2607 print *,"mat_L" 2608 do ic=1,ncar 2609 print *,VCIInverse(ic,:ncar) 2610 end do 2611 print *,"mauvaise inversion matrice" 2612 ok=.false. 2613 stop 2614 return 2615 end if 2616 2617 do ic=1,ncar 2618 do jc=1,ic 2619 VCI(ic,jc)=VCIInverse(ic+1,jc) 2620 VCI(jc,ic)=VCI(ic,jc) 2621 end do 2622 end do 2623 2624 2625 end subroutine get_inv_residual_covariance_matrix_LU
get_inv_residual_covariance_matrix_LU_cd
[ Top ] [ m_qtlmap_incidence_multi ] [ Subroutines ]
NAME
get_inv_residual_covariance_matrix_LU_cd
DESCRIPTION
NOTE
SOURCE
2636 subroutine get_inv_residual_covariance_matrix_LU_cd(ip,kd,n,x,vci,determ,ok) 2637 integer , intent(in) :: ip,kd,n 2638 real (kind=dp) ,dimension(n) , intent(in) :: x 2639 real(kind=dp),dimension(ncar,ncar),intent(out) :: VCI 2640 real(kind=dp), intent(out) :: determ 2641 logical ,intent(out) :: ok 2642 2643 real(kind=dp),dimension(ncar+1,ncar) :: VCIInverse 2644 2645 real(kind=dp),dimension(ncar,ncar) :: mat_L,V 2646 real(kind=dp),dimension(ncar) :: mat_H 2647 2648 real(kind=dp) :: rh 2649 integer :: ic,jc,k,irank 2650 2651 ok=.true. 2652 !RHO(I,J) avec i/=j => x((i-1)*(ncar-i+1)) 2653 k=0 2654 mat_L=0.d0 2655 do ic=1,ncar 2656 do jc=1,ic 2657 k=k+1 2658 mat_L(ic,jc) = x(ncar*(np-1)+k) 2659 end do 2660 end do 2661 2662 V = matmul(mat_L,transpose(mat_L)) 2663 2664 mat_H=0.d0 2665 ! on fixe pour ip=1 H 2666 if ( ip == 1) then 2667 mat_H(:)=1.d0 2668 else 2669 ! do ip=2,np 2670 do jc=1,ncar 2671 mat_H(jc)=x((ip-2)*ncar+jc) 2672 end do 2673 ! end do 2674 end if 2675 2676 2677 determ=0 2678 VCIInverse=0.d0 2679 do ic=1,ncar 2680 do jc=1,ncar 2681 VCIInverse(ic,jc) = V(ic,jc)*mat_H(ic)*mat_H(jc) 2682 if ( ic /= jc ) then 2683 VCIInverse(ic,jc) = VCIInverse(ic,jc)*datasetUser%corcd(ic,jc,kd) 2684 else 2685 VCIInverse(ic,jc) = VCIInverse(ic,jc)*cd(ic,kd) 2686 end if 2687 ! VCIInverse(jc,ic) = VCIInverse(ic,jc) 2688 end do 2689 end do 2690 2691 !Calcul Determinant de la matrice de covariance residuelle du pere ip 2692 !-------------------------------------------------------------------- 2693 irank=0 2694 2695 ! call MATH_QTLMAP_F03ABF(VCI,ncar,ncar,determ,irank) 2696 ! if (irank /= 0) then 2697 ! print *,"mauvais determinant matrice" 2698 ! ok=.false. 2699 ! return 2700 ! end if 2701 !Calcul Inverse de la matrice de covariance residuelle du pere ip 2702 !----------------------------------------------------------------- 2703 ! VCIInverse(:ncar,:ncar)=VCI 2704 2705 ! CALL MATH_QTLMAP_INVDETMAT(ncar,VCIInverse,determ,irank) 2706 ! 2707 ! ! call MATH_QTLMAP_F01ADF(ncar,VCIInverse,ncar+1,irank) 2708 ! if (irank /= 0) then 2709 ! do ic=1,ncar 2710 ! print *,VCIInverse(ic,:ncar) 2711 ! end do 2712 ! print *,"mauvaise inversion matrice" 2713 ! ok=.false. 2714 ! stop 2715 ! return 2716 ! end if 2717 ! 2718 ! do ic=1,ncar 2719 ! do jc=ic,ncar 2720 ! VCI(ic,jc)=VCIInverse(ic,jc) 2721 ! VCI(jc,ic)=VCI(ic,jc) 2722 ! end do 2723 ! end do 2724 2725 CALL MATH_QTLMAP_INVDETMATSYM(ncar,VCIInverse,ncar+1,determ,irank) 2726 ! call MATH_QTLMAP_F01ADF(ncar,VCIInverse,ncar+1,irank) 2727 if (irank /= 0) then 2728 print *,"mauvaise inversion matrice" 2729 ok=.false. 2730 return 2731 end if 2732 2733 do ic=1,ncar 2734 do jc=1,ic 2735 VCI(ic,jc)=VCIInverse(ic+1,jc) 2736 VCI(jc,ic)=VCI(ic,jc) 2737 end do 2738 end do 2739 2740 end subroutine get_inv_residual_covariance_matrix_LU_cd
init_dataset_ncar
[ Top ] [ m_qtlmap_incidence_multi ] [ Subroutines ]
NAME
init_dataset_ncar
DESCRIPTION
NOTE constition des famille de plein et demi frere via des indexes pas de notion d estimable...on comptabilise seulement les effectifs et on repere les lignes de la matrice de contingence
SOURCE
329 subroutine init_dataset_ncar(nqtl,dataset) 330 type (DATASET_TYPE) , intent(inout) :: dataset 331 integer , intent(in) :: nqtl 332 logical ,dimension(nd) :: pres 333 integer :: ip,jm,kd,ifem,kdd,s1,s2,chr,iq,ic,kkd,effp,efft,jc,eff 334 real(kind=dp) :: somyp(ncar),somy(ncar) 335 336 call log_mess("build_family for incidence multitrait analysis...",DEBUG_DEF) 337 338 estime_multi=.false. 339 allocate(dataset%lSires(np)) 340 kdd=0 341 do ip=1,np 342 allocate(dataset%lSires(ip)%full_sib(nmp(ip+1)-nmp(ip))) 343 ifem=0 344 do jm=nmp(ip)+1,nmp(ip+1) 345 ifem=ifem+1 346 kd=ndm(jm)+1 347 do while ( kd<=ndm(jm+1) .and. .not. (count(presentc(:,kd))==ncar) ) 348 kd=kd+1 349 end do 350 if (kd>ndm(jm+1)) then 351 dataset%lSires(ip)%full_sib(ifem)%firstKd=-1 352 dataset%lSires(ip)%full_sib(ifem)%lastKd=-2 353 cycle 354 end if 355 356 kdd=kdd+1 357 ! on a trouve le premier kd valide du pere ip dans la matrice d incidence 358 dataset%lSires(ip)%full_sib(ifem)%firstKd=kdd 359 do kd=kd+1,ndm(jm+1) 360 if ( (count(presentc(:,kd))==ncar) ) then 361 kdd=kdd + 1 362 end if 363 end do 364 dataset%lSires(ip)%full_sib(ifem)%lastKd=kdd 365 366 !NOUVEAU ESTIME.... 367 estime_multi(jm)=(dataset%lSires(ip)%full_sib(ifem)%lastKd-dataset%lSires(ip)%full_sib(ifem)%firstKd)>=NDMIN 368 369 call log_mess("Full Sib ["//trim(pere(ip))//"-"//trim(mere(jm))//"] :"//& 370 trim(str(dataset%lSires(ip)%full_sib(ifem)%firstKd))//"->"& 371 //trim(str(dataset%lSires(ip)%full_sib(ifem)%lastKd)),VERBOSE_DEF) 372 s1=0 373 do chr=1,nchr 374 s1=max(s1,ngenom(chr,jm+1)-ngenom(chr,jm)) 375 end do 376 s2=dataset%lSires(ip)%full_sib(ifem)%lastKd-dataset%lSires(ip)%full_sib(ifem)%firstKd+1 377 if ( nqtl >= 1) then 378 allocate(dataset%lSires(ip)%full_sib(ifem)%ppt(nqtl,s1,s2)) 379 allocate(dataset%lSires(ip)%full_sib(ifem)%pmt(nqtl,s1,s2)) 380 end if 381 end do 382 ifem=0 383 do jm=nmp(ip)+1,nmp(ip+1) 384 ifem=ifem+1 385 if (dataset%lSires(ip)%full_sib(ifem)%firstKd<0) cycle 386 dataset%lSires(ip)%half_sib%firstKd= dataset%lSires(ip)%full_sib(ifem)%firstKd 387 exit 388 end do 389 ifem=nmp(ip+1)-nmp(ip)+1 390 do jm=nmp(ip+1),nmp(ip)+1,-1 391 ifem=ifem-1 392 if (dataset%lSires(ip)%full_sib(ifem)%lastKd<0) cycle 393 dataset%lSires(ip)%half_sib%lastKd=dataset%lSires(ip)%full_sib(ifem)%lastKd 394 exit 395 end do 396 397 call log_mess("Half Sib ["//trim(pere(ip))//"] :"//& 398 trim(str(dataset%lSires(ip)%half_sib%firstKd))//"->"//trim(str(dataset%lSires(ip)%half_sib%lastKd)),& 399 VERBOSE_DEF) 400 !s1=lSires(ip)%half_sib%lastKd-lSires(ip)%half_sib%firstKd+1 401 if ( nqtl >= 1) then 402 allocate(dataset%lSires(ip)%ppt(nqtl,& 403 dataset%lSires(ip)%half_sib%firstKd:dataset%lSires(ip)%half_sib%lastKd)) 404 end if 405 end do ! ip 406 407 do ic=2,ncar 408 dataset%lSires = dataset%lSires 409 end do 410 411 dataset%nkd=0 412 dataset%nkd_max_by_fam=0 413 414 do ip=1,np 415 do jm=nmp(ip)+1,nmp(ip+1) 416 eff=0 417 do kd=ndm(jm)+1,ndm(jm+1) 418 if (count(presentc(:,kd))==ncar) then 419 dataset%nkd=dataset%nkd+1 420 eff=eff+1 421 end if 422 end do !kd 423 dataset%nkd_max_by_fam=max(dataset%nkd_max_by_fam,eff) 424 end do ! jm 425 end do ! ip 426 427 allocate (dataset%Y(ncar,dataset%nkd)) 428 allocate (dataset%CD(ncar,dataset%nkd)) 429 430 do ic=1,ncar 431 kkd=0 432 do ip=1,np 433 do jm=nmp(ip)+1,nmp(ip+1) 434 do kd=ndm(jm)+1,ndm(jm+1) 435 if (count(presentc(:,kd))==ncar) then 436 kkd=kkd+1 437 if ( natureY(ic) /= 'i' ) then 438 dataset%Y(ic,kkd)=y(ic,kd) 439 dataset%CD(ic,kkd)=cd(ic,kd) 440 else 441 call stop_application("Type of trait 'i' are not managed in multi-trait analysis!") 442 dataset%YDISCRETORD(ic,kkd)=ydiscretord(ic,kd) 443 end if 444 end if 445 end do !kd 446 end do ! jm 447 end do ! ip 448 end do !ic 449 450 somy=0 451 efft=0 452 do ip=1,np 453 allocate(dataset%lSires(ip)%sig0(ncar)) 454 allocate(dataset%lSires(ip)%xmu0(ncar)) 455 effp=0 456 somyp=0 457 do jm=nmp(ip)+1,nmp(ip+1) 458 do kd=ndm(jm)+1,ndm(jm+1) 459 if (count(presentc(:,kd))==ncar) then 460 effp=effp+1 461 do ic=1,ncar 462 somyp(ic)=somyp(ic)+y(ic,kd)!*cd(ic,kd) 463 end do 464 end if 465 end do 466 end do 467 468 !somme total sur le caractere 469 do ic=1,ncar 470 somy(ic)=somy(ic)+somyp(ic) 471 end do 472 !effectif total pour le caractere 473 efft=efft+effp 474 475 if (effp == 0.d0) then 476 call stop_application('Father ['//trim(pere(ip))//'] has got no child with trait value ['//trim(carac(ic))//"]") 477 end if 478 479 !moyenne du caractere dans la famille du pere ip 480 do ic=1,ncar 481 dataset%lSires(ip)%xmu0(ic)=somyp(ic)/dble(effp) 482 call log_mess("Mean for Sire/trait ["//trim(pere(ip))//","//trim(carac(ic))& 483 //"]="//str(dataset%lSires(ip)%xmu0(ic)),INFO_DEF) 484 end do 485 !ecart type du caractere dans la famille du pere ip 486 dataset%lSires(ip)%sig0(:)=0.d0 487 do jm=nmp(ip)+1,nmp(ip+1) 488 do kd=ndm(jm)+1,ndm(jm+1) 489 if (count(presentc(:,kd))==ncar) then 490 do ic=1,ncar 491 dataset%lSires(ip)%sig0(ic)=dataset%lSires(ip)%sig0(ic)+(y(ic,kd)-& 492 dataset%lSires(ip)%xmu0(ic))*(y(ic,kd)-dataset%lSires(ip)%xmu0(ic)) 493 end do 494 end if 495 end do !kd 496 end do !jm 497 498 do ic=1,ncar 499 dataset%lSires(ip)%sig0(ic)=sqrt(dataset%lSires(ip)%sig0(ic)/dble(effp-1)) 500 call log_mess("Standart deviation for Sire/trait ["//trim(pere(ip))//","//trim(carac(ic))& 501 //"]="//str(dataset%lSires(ip)%sig0(ic)),INFO_DEF) 502 end do 503 end do !ip 504 505 !moyenne du caractere sur l ensemble de la population 506 allocate(dataset%xmu(ncar)) 507 do ic=1,ncar 508 dataset%xmu(ic)=somy(ic)/dble(efft) 509 call log_mess("Mean for trait ["//trim(carac(ic))//"]="//str(dataset%xmu(ic)),INFO_DEF) 510 end do 511 512 allocate(dataset%sig(ncar)) 513 allocate(dataset%cov(ncar,ncar)) 514 allocate(dataset%rhoi(ncar,ncar)) 515 dataset%cov=0.d0 516 dataset%sig=0.d0 517 518 do ip=1,np 519 do jm=nmp(ip)+1,nmp(ip+1) 520 do kd=ndm(jm)+1,ndm(jm+1) 521 if (count(presentc(:,kd))==ncar) then 522 do ic=1,ncar 523 dataset%sig(ic)=dataset%sig(ic)+(y(ic,kd)-dataset%xmu(ic))*(y(ic,kd)-dataset%xmu(ic)) 524 do jc=1,ic-1 525 dataset%cov(jc,ic)=dataset%cov(jc,ic)+(y(jc,kd)-dataset%xmu(jc))*(y(ic,kd)-dataset%xmu(ic)) 526 end do 527 end do 528 end if 529 end do 530 end do 531 end do 532 533 !variance sur l ensemble de la population 534 do ic=1,ncar 535 dataset%sig(ic)=sqrt(dataset%sig(ic)/dble(efft-1)) 536 call log_mess("Standart deviation for trait ["//trim(carac(ic))//"]="//str(dataset%sig(ic)),INFO_DEF) 537 end do 538 !calcul des covariances 539 do ic=1,ncar 540 do jc=1,ic-1 541 dataset%cov(jc,ic)=dataset%cov(jc,ic)/dble(efft) 542 dataset%cov(ic,jc)=dataset%cov(jc,ic) 543 dataset%rhoi(jc,ic)=dataset%cov(jc,ic)/(dataset%sig(ic)*dataset%sig(jc)) 544 dataset%rhoi(ic,jc)=dataset%rhoi(jc,ic) 545 end do 546 end do 547 548 ! stop 549 end subroutine init_dataset_ncar
init_incidence_multi
[ Top ] [ m_qtlmap_incidence_multi ] [ Subroutines ]
NAME
init_incidence_multi
DESCRIPTION
NOTE
SOURCE
127 subroutine init_incidence_multi(nqtl,xinc,listDesc,type_incidence,resolution_LU) 128 integer , intent(in) :: nqtl,type_incidence 129 real (kind=dp) , dimension(:,:,:), pointer, intent(inout) :: xinc 130 type(INCIDENCE_TYPE) ,dimension(ncar), intent(inout) :: listDesc 131 logical, intent(in) :: resolution_LU 132 133 integer :: ip,jm,kd,ifem,eff,nteffmaxtotal,ic,nparmax,k,jc,s 134 135 ! General array to get index information about femal estimable 136 ! *** 137 allocate (estime_multi(nm)) 138 allocate (estfem_multi(nfem)) 139 allocate (iam_multi(nfem)) 140 estime_multi=.false. 141 estfem_multi=.false. 142 namest_multi=0 143 iam_multi=0 144 145 do ip=1,np 146 do jm=nmp(ip)+1,nmp(ip+1) 147 eff=0 148 do kd=ndm(jm)+1,ndm(jm+1) 149 if(count(presentc(:,kd))==ncar) then 150 eff=eff+1 151 end if 152 end do 153 ifem=repfem(jm) 154 if( eff>=NDMIN ) then 155 estime_multi(jm)=.true. 156 estfem_multi(ifem)=.true. 157 end if 158 end do !jm 159 end do !ip 160 161 do ifem=1,nfem 162 if(estfem_multi(ifem)) then 163 namest_multi=namest_multi+1 164 iam_multi(ifem)=namest_multi 165 end if 166 end do 167 168 ntnivmaxtotal=0 169 nteffmaxtotal=0 170 do ic=1,ncar 171 listDesc(ic)%ic=ic 172 listDesc(ic)%nqtl=nqtl 173 call set_ntnivmax(ic,nqtl,listDesc(ic)%ntnivmax,listDesc(ic)%nteffmax,listDesc(ic)%ntlev,listDesc(ic)%nbniv) 174 ntnivmaxtotal=max(ntnivmaxtotal,listDesc(ic)%ntnivmax) 175 nteffmaxtotal=max(nteffmaxtotal,listDesc(ic)%nteffmax) 176 177 if ( type_incidence==INCIDENCE_TRAIT_COV ) then 178 listDesc(ic)%ntnivmax=listDesc(ic)%ntnivmax+ncar 179 listDesc(ic)%nteffmax=listDesc(ic)%nteffmax+1 180 end if 181 182 183 listDesc(ic)%ntniv = 0 184 listDesc(ic)%nbnivest = 0 185 listDesc(ic)%nteff = 0 186 187 call log_mess("********************* INIT CARAC="// trim(str(ic))//" *****************************",DEBUG_DEF) 188 call log_mess("NTNIVMAX:"//str(listDesc(ic)%ntnivmax),DEBUG_DEF) 189 call log_mess("NTEFFMAX:"//str(listDesc(ic)%nteffmax),DEBUG_DEF) 190 call log_mess("NTLEVQTL / Reproductor:"//str(listDesc(ic)%ntlev),DEBUG_DEF) 191 call log_mess("NBLEV FIXED EFFECT:"//str(listDesc(ic)%nbniv),DEBUG_DEF) 192 call log_mess("*********************************************************",DEBUG_DEF) 193 194 allocate(listDesc(ic)%desc(listDesc(ic)%nteffmax)) 195 allocate(listDesc(ic)%vecsol(listDesc(ic)%ntnivmax)) 196 listDesc(ic)%vecsol=.false. 197 ! allocate(listDesc(ic)%corniv(ntnivmax)) 198 allocate(listDesc(ic)%precis(listDesc(ic)%ntnivmax)) 199 listDesc(ic)%precis=0.d0 200 201 if ( nqtl > 0) then 202 allocate (listDesc(ic)%ntniv_qtlsires(nqtl)) 203 allocate (listDesc(ic)%ntniv_qtldams(nqtl)) 204 end if 205 allocate (listDesc(ic)%corr_niv_nivb(listDesc(ic)%ntnivmax)) 206 end do 207 208 allocate (xinc(ncar,nd,ntnivmaxtotal)) 209 xinc=0.d0 210 211 allocate (listDesc(1)%dataset) 212 call init_dataset_ncar(nqtl,listDesc(1)%dataset) 213 214 do ic=1,ncar 215 allocate (listDesc(ic)%nonull(listDesc(1)%dataset%nkd,ntnivmaxtotal)) 216 allocate (listDesc(ic)%nbnonull(listDesc(1)%dataset%nkd)) 217 end do 218 219 nparmax=ncar*np+ncar*(ncar-1)/2+ncar*ntnivmaxtotal 220 221 allocate (listDesc(1)%borni(nparmax)) 222 allocate (listDesc(1)%borns(nparmax)) 223 allocate (listDesc(1)%par(nparmax)) 224 225 if ( resolution_LU ) then 226 listDesc(1)%par(1:(np-1)*ncar)=1.d0 ! matH 227 listDesc(1)%par((np-1)*ncar+1:(np-1)*ncar+ncar*(ncar+1)/2)=0.d0 ! matL 228 229 s=(np-1)*ncar+1 230 listDesc(1)%par(s)=1.d0 231 232 do ic=2, ncar 233 listDesc(1)%par(s+ic)=1.d0 234 s = s + ic 235 end do 236 237 listDesc(1)%borni(:)=DEFAULT_PARAM_MIN 238 listDesc(1)%borns(:)=DEFAULT_PARAM_MAX 239 240 listDesc(1)%borni(1:(np-1)*ncar)=0.1d0 241 242 listDesc(1)%par(np*ncar+ncar*(ncar-1)/2+1:)=0.d0 243 244 ! print *,"**********>PAR:",listDesc(1)%par(:(np-1)*ncar) 245 else 246 247 listDesc(1)%borni(1:ncar*np)=SIG_MIN 248 listDesc(1)%borns(1:ncar*np)=SIG_MAX 249 250 k=ncar*np 251 !bornes des correlations residuelles 252 do ic=2,ncar 253 do jc=1,ic-1 254 k=k+1 255 ! ** MODIF RESCALE ** OFI => 256 !listDesc(1)%borni(k)=-1.d0 257 !listDesc(1)%borns(k)=1.d0 258 listDesc(1)%borni(k)=DEFAULT_PARAM_MIN 259 listDesc(1)%borns(k)=DEFAULT_PARAM_MAX 260 end do 261 end do 262 263 listDesc(1)%borni(k+1:)=DEFAULT_PARAM_MIN 264 listDesc(1)%borns(k+1:)=DEFAULT_PARAM_MAX 265 266 !initialisation des variances 267 do ic=1,ncar 268 do ip=1,np 269 listDesc(1)%par((ic-1)*np+ip)=listDesc(1)%dataset%lSires(ip)%sig0(ic) 270 end do 271 end do 272 273 listDesc(1)%par(ncar*np+1:)=0.d0 274 275 end if 276 277 do ic=2,ncar 278 listDesc(ic)%dataset=>listDesc(1)%dataset 279 listDesc(ic)%borni=>listDesc(1)%borni 280 listDesc(ic)%borns=>listDesc(1)%borns 281 listDesc(ic)%par=>listDesc(1)%par 282 end do 283 284 end subroutine init_incidence_multi
init_startpoint
[ Top ] [ m_qtlmap_incidence_multi ] [ Subroutines ]
NAME
init_startpoint
DESCRIPTION
NOTE initialise le point de depart pour la minimisation
SOURCE
1186 subroutine init_startpoint(workStruct,incDesc) 1187 type(INCIDENCE_TYPE) ,intent(inout) :: incDesc 1188 type(INCIDENCE_GEN_STRUCT) ,intent(in) :: workstruct 1189 1190 real(kind=dp) :: r1,tempLin(ncar*ncar) 1191 integer :: k,ic,jc,ip,i 1192 1193 !initialisation d'un point de depart pour les variances 1194 k=0 1195 do ic=1,ncar 1196 do ip=1,np 1197 k=k+1 1198 incDesc%par(k)=sqrt(workstruct%sigsquare(workstruct%nqtl,ip,ic)) 1199 end do 1200 end do 1201 1202 i=0 1203 do jc=1,ncar-1 1204 do ic=jc+1,ncar 1205 i = i + 1 1206 tempLin(i) = workstruct%rhoi(workstruct%nqtl,ic,jc) 1207 end do 1208 end do 1209 1210 i=0 1211 !initialisation d'un point de depart pour les correlations residuelles 1212 do ic=2,ncar 1213 do jc=1,ic-1 1214 k=k+1 1215 i=i+1 1216 ! *** RESCALE OFI *** 1217 r1=tempLin(i) 1218 incDesc%par(k)=dlog((1.d0+r1)/(1.d0-r1)) 1219 !incDesc%par(k)=workstruct%rhoi(workstruct%nqtl,jc,ic) 1220 end do 1221 end do 1222 1223 1224 1225 1226 1227 incDesc%par(ncar*np+1+ncar*(ncar-1)/2:)=0.d0 1228 1229 end subroutine init_startpoint
likelihood_ncar_h0_family
[ Top ] [ m_qtlmap_incidence_multi ] [ Subroutines ]
NAME
likelihood_ncar_h0_family
DESCRIPTION
NOTE
SOURCE
1743 subroutine likelihood_ncar_h0_family(ip,jm,n,x,vci,determ,f,iuser,user) 1744 integer , intent(in) :: ip,jm,n 1745 real (kind=dp) ,dimension(n), intent(in) :: x 1746 real (kind=dp) , intent(inout) :: f 1747 integer , dimension(1), intent(inout) :: iuser 1748 real (kind=dp) ,dimension(1), intent(inout) :: user 1749 real (kind=dp) ,intent(in) :: determ 1750 real(kind=dp),dimension(ncar,ncar) ,intent(in) :: VCI 1751 real (kind=dp) :: vpf 1752 integer :: ifem 1753 real(kind=dp),dimension(dataset%nkd_max_by_fam,ncar) :: V 1754 1755 ! 1756 integer :: ic,jc,s,nbnivest,kd1,kd2,kd,na 1757 logical :: ok,valide 1758 1759 ! print *,x 1760 f = 0.d0 1761 V = 0.d0 1762 s=np*ncar+ncar*(ncar-1)/2 1763 1764 ifem=jm-nmp(ip) 1765 kd1=dataset%lSires(ip)%full_sib(ifem)%firstkd 1766 kd2=dataset%lSires(ip)%full_sib(ifem)%lastkd 1767 na = kd2-kd1+1 1768 if (kd2<kd1) then 1769 f=0 1770 return 1771 end if 1772 1773 ! Dans la famille du pere ip, on calcul le Y (performance auquel on soustrait l'ensemble des effets) 1774 ! pour chaque caractere 1775 do ic=1,ncar 1776 nbnivest=my_listdesc(ic)%nbnivest 1777 call matmul_incidence(kd1,kd2,my_listdesc(ic),nd,ntnivmaxtotal,& 1778 My_XINCREDUITMUL(:,:,ic),X(s+1:s+nbnivest),dataset%nkd_max_by_fam,1,V(:,ic)) 1779 V(:na,ic) = my_listdesc(ic)%dataset%Y(ic,kd1:kd2)-V(:na,ic) 1780 s=s+my_listdesc(ic)%nbnivest 1781 end do 1782 1783 1784 vpf=0 1785 !calcul de la vraissemblance de plein frere 1786 do kd=1,na 1787 vpf=vpf+dot_product(matmul(V(kd,:),VCI),V(kd,:)) 1788 end do 1789 !finallement, la vraissemblance de la famille ip/jm: 1790 f=+0.5*vpf+dble(kd2-kd1+1)*0.5*log(determ) 1791 end subroutine likelihood_ncar_h0_family
likelihood_ncar_h0_family_LU
[ Top ] [ m_qtlmap_incidence_multi ] [ Subroutines ]
NAME
likelihood_ncar_h0_family_LU
DESCRIPTION
NOTE
SOURCE
2752 subroutine likelihood_ncar_h0_family_LU(ip,jm,n,x,f,iuser,user) 2753 integer , intent(in) :: ip,jm,n 2754 real (kind=dp) ,dimension(n), intent(in) :: x 2755 real (kind=dp) , intent(inout) :: f 2756 integer , dimension(1), intent(inout) :: iuser 2757 real (kind=dp) ,dimension(1), intent(inout) :: user 2758 real (kind=dp) :: determ 2759 real(kind=dp) ,dimension(ncar,ncar) :: VCI 2760 ! real(kind=dp),dimension(ncar,ncar) :: VCITemp 2761 real (kind=dp) :: vpf 2762 integer :: ifem 2763 real(kind=dp),dimension(dataset%nkd_max_by_fam,ncar) :: V 2764 2765 ! 2766 integer :: ic,jc,s,nbnivest,kd1,kd2,kd,na 2767 logical :: ok,valide 2768 2769 call get_inv_residual_covariance_matrix_LU(ip,n,x,vci,determ,ok) 2770 f = 0.d0 2771 2772 ! do ip=1,np 2773 ! VCITemp = VCI(:,:) 2774 ! do jm=nmp(ip)+1,nmp(ip+1) 2775 V = 0.d0 2776 s=np*ncar+ncar*(ncar-1)/2 2777 2778 ifem=jm-nmp(ip) 2779 kd1=dataset%lSires(ip)%full_sib(ifem)%firstkd 2780 kd2=dataset%lSires(ip)%full_sib(ifem)%lastkd 2781 na = kd2-kd1+1 2782 if (kd2<kd1) then 2783 f=0 2784 return 2785 end if 2786 2787 ! Dans la famille du pere ip, on calcul le Y (performance auquel on soustrait l'ensemble des effets) 2788 ! pour chaque caractere 2789 do ic=1,ncar 2790 nbnivest=my_listdesc(ic)%nbnivest 2791 call matmul_incidence(kd1,kd2,my_listdesc(ic),nd,ntnivmaxtotal,& 2792 My_XINCREDUITMUL(:,:,ic),X(s+1:s+nbnivest),dataset%nkd_max_by_fam,1,V(:,ic)) 2793 V(:na,ic) = my_listdesc(ic)%dataset%Y(ic,kd1:kd2)-V(:na,ic) 2794 s=s+my_listdesc(ic)%nbnivest 2795 end do 2796 2797 2798 vpf=0 2799 !calcul de la vraissemblance de plein frere 2800 do kd=1,na 2801 vpf=vpf+dot_product(matmul(V(kd,:),VCI),V(kd,:)) 2802 end do 2803 !finallement, la vraissemblance de la famille ip/jm: 2804 f=0.5*vpf+dble(kd2-kd1+1)*0.5*log(determ) 2805 ! print *,f 2806 ! end do 2807 ! end do 2808 !stop 2809 2810 end subroutine likelihood_ncar_h0_family_LU
likelihood_ncar_h0_family_withcd
[ Top ] [ m_qtlmap_incidence_multi ] [ Subroutines ]
NAME
likelihood_ncar_h0_family_withcd
DESCRIPTION
NOTE
multivariate analysis
(indicecar-1)*np +ip : SIG ip,indice_caracetere
ncar*np + (ic-1)*ncar+jc : RHO (ic,jc) residual correlation between two traits....
ncar*np + ncar*ncar + nbnivest
SOURCE
1929 subroutine likelihood_ncar_h0_family_withcd(ip,jm,n,x,f,iuser,user) 1930 integer , intent(in) :: ip,jm,n 1931 real (kind=dp) ,dimension(n), intent(in) :: x 1932 real (kind=dp) , intent(inout) :: f 1933 integer , dimension(1), intent(inout) :: iuser 1934 real (kind=dp) ,dimension(1), intent(inout) :: user 1935 1936 1937 real (kind=dp),dimension(dataset%nkd_max_by_fam) :: determ 1938 real(kind=dp),dimension(ncar,ncar,dataset%nkd_max_by_fam) :: VCI 1939 real (kind=dp) :: vpf 1940 integer :: ifem,kkk 1941 real(kind=dp),dimension(dataset%nkd_max_by_fam,ncar) :: V 1942 1943 ! 1944 integer :: ic,jc,s,nbnivest,kd1,kd2,kd,na 1945 logical :: ok,valide 1946 1947 determ=0.d0 1948 kd=1 1949 do kkk=ndm(jm)+1,ndm(jm+1) 1950 if ( count(presentc(:,kkk)) == ncar ) then 1951 call get_inv_residual_covariance_matrix_cd(ip,kkk,n,x,vci(:,:,kd),determ(kd),ok) 1952 if ( .not. ok ) then 1953 f=INIFINY_REAL_VALUE 1954 return 1955 end if 1956 kd = kd + 1 1957 end if 1958 end do 1959 1960 f = 0.d0 1961 V = 0.d0 1962 s=np*ncar+ncar*(ncar-1)/2 1963 1964 ifem=jm-nmp(ip) 1965 1966 kd1=dataset%lSires(ip)%full_sib(ifem)%firstkd 1967 kd2=dataset%lSires(ip)%full_sib(ifem)%lastkd 1968 na = kd2-kd1+1 1969 1970 if (kd2<kd1) then 1971 f=0 1972 return 1973 end if 1974 1975 ! Dans la famille du pere ip, on calcul le Y (performance auquel on soustrait l'ensemble des effets) 1976 ! pour chaque caractere 1977 do ic=1,ncar 1978 nbnivest=my_listdesc(ic)%nbnivest 1979 call matmul_incidence(kd1,kd2,my_listdesc(ic),nd,ntnivmaxtotal,& 1980 My_XINCREDUITMUL(:,:,ic),X(s+1:s+nbnivest),dataset%nkd_max_by_fam,1,V(:,ic)) 1981 V(:na,ic) = my_listdesc(ic)%dataset%Y(ic,kd1:kd2)-V(:na,ic) 1982 s=s+my_listdesc(ic)%nbnivest 1983 end do 1984 1985 vpf=0 1986 !calcul de la vraissemblance de plein frere 1987 do kd=1,na 1988 vpf=vpf+dot_product(matmul(V(kd,:),VCI(:,:,kd)),V(kd,:)) 1989 end do 1990 1991 !finallement, la vraissemblance de la famille ip/jm: 1992 f=+0.5*vpf+0.5*sum(log(determ(:na))) 1993 !print *,f 1994 1995 end subroutine likelihood_ncar_h0_family_withcd
likelihood_ncar_h0_LU
[ Top ] [ m_qtlmap_incidence_multi ] [ Subroutines ]
NAME
likelihood_ncar_h0_LU
DESCRIPTION
NOTE
SOURCE
2821 subroutine likelihood_ncar_h0_LU(n,x,f,iuser,user) 2822 integer , intent(in) :: n 2823 real (kind=dp) ,dimension(n), intent(in) :: x 2824 real (kind=dp) , intent(inout) :: f 2825 integer , dimension(1), intent(inout) :: iuser 2826 real (kind=dp) ,dimension(1), intent(inout) :: user 2827 real (kind=dp) :: determ 2828 real(kind=dp) ,dimension(ncar,ncar) :: VCI 2829 ! real(kind=dp),dimension(ncar,ncar) :: VCITemp 2830 real (kind=dp) :: vpf 2831 integer :: ifem,ip,jm 2832 real(kind=dp),dimension(dataset%nkd_max_by_fam,ncar) :: V 2833 2834 ! 2835 integer :: ic,jc,s,nbnivest,kd1,kd2,kd,na 2836 logical :: ok,valide 2837 2838 f = 0.d0 2839 2840 do ip=1,np 2841 call get_inv_residual_covariance_matrix_LU(ip,n,x,vci,determ,ok) 2842 ! VCITemp = VCI(:,:) 2843 do jm=nmp(ip)+1,nmp(ip+1) 2844 V = 0.d0 2845 s=np*ncar+ncar*(ncar-1)/2 2846 2847 ifem=jm-nmp(ip) 2848 kd1=dataset%lSires(ip)%full_sib(ifem)%firstkd 2849 kd2=dataset%lSires(ip)%full_sib(ifem)%lastkd 2850 na = kd2-kd1+1 2851 if (kd2<kd1) then 2852 cycle 2853 end if 2854 2855 ! Dans la famille du pere ip, on calcul le Y (performance auquel on soustrait l'ensemble des effets) 2856 ! pour chaque caractere 2857 do ic=1,ncar 2858 nbnivest=my_listdesc(ic)%nbnivest 2859 call matmul_incidence(kd1,kd2,my_listdesc(ic),nd,ntnivmaxtotal,& 2860 My_XINCREDUITMUL(:,:,ic),X(s+1:s+nbnivest),dataset%nkd_max_by_fam,1,V(:,ic)) 2861 V(:na,ic) = my_listdesc(ic)%dataset%Y(ic,kd1:kd2)-V(:na,ic) 2862 s=s+my_listdesc(ic)%nbnivest 2863 end do 2864 2865 2866 vpf=0 2867 !calcul de la vraissemblance de plein frere 2868 do kd=1,na 2869 vpf=vpf+dot_product(matmul(V(kd,:),VCI),V(kd,:)) 2870 end do 2871 !finallement, la vraissemblance de la famille ip/jm: 2872 f=f+0.5*vpf+dble(kd2-kd1+1)*0.5*log(determ) 2873 ! print *,f 2874 end do 2875 end do 2876 2877 ! print *,f 2878 !stop 2879 2880 end subroutine likelihood_ncar_h0_LU
likelihood_ncar_h0_LU_family_withcd
[ Top ] [ m_qtlmap_incidence_multi ] [ Subroutines ]
NAME
likelihood_ncar_h0_LU_family_withcd
DESCRIPTION
NOTE
multivariate analysis
(indicecar-1)*np +ip : SIG ip,indice_caracetere
ncar*np + (ic-1)*ncar+jc : RHO (ic,jc) residual correlation between two traits....
ncar*np + ncar*ncar + nbnivest
SOURCE
2895 subroutine likelihood_ncar_h0_LU_family_withcd(ip,jm,n,x,f,iuser,user) 2896 integer , intent(in) :: ip,jm,n 2897 real (kind=dp) ,dimension(n), intent(in) :: x 2898 real (kind=dp) , intent(inout) :: f 2899 integer , dimension(1), intent(inout) :: iuser 2900 real (kind=dp) ,dimension(1), intent(inout) :: user 2901 2902 2903 real (kind=dp),dimension(dataset%nkd_max_by_fam) :: determ 2904 real(kind=dp),dimension(ncar,ncar,dataset%nkd_max_by_fam) :: VCI 2905 real (kind=dp) :: vpf 2906 integer :: ifem,kkk 2907 real(kind=dp),dimension(dataset%nkd_max_by_fam,ncar) :: V 2908 2909 ! 2910 integer :: ic,jc,s,nbnivest,kd1,kd2,kd,na 2911 logical :: ok,valide 2912 2913 determ=0.d0 2914 kd=1 2915 do kkk=ndm(jm)+1,ndm(jm+1) 2916 if ( count(presentc(:,kkk)) == ncar ) then 2917 call get_inv_residual_covariance_matrix_LU_cd(ip,kkk,n,x,vci(:,:,kd),determ(kd),ok) 2918 if ( .not. ok ) then 2919 f=INIFINY_REAL_VALUE 2920 return 2921 end if 2922 kd = kd + 1 2923 end if 2924 end do 2925 2926 f = 0.d0 2927 V = 0.d0 2928 s=np*ncar+ncar*(ncar-1)/2 2929 2930 ifem=jm-nmp(ip) 2931 2932 kd1=dataset%lSires(ip)%full_sib(ifem)%firstkd 2933 kd2=dataset%lSires(ip)%full_sib(ifem)%lastkd 2934 na = kd2-kd1+1 2935 2936 if (kd2<kd1) then 2937 f=0 2938 return 2939 end if 2940 2941 ! Dans la famille du pere ip, on calcul le Y (performance auquel on soustrait l'ensemble des effets) 2942 ! pour chaque caractere 2943 do ic=1,ncar 2944 nbnivest=my_listdesc(ic)%nbnivest 2945 call matmul_incidence(kd1,kd2,my_listdesc(ic),nd,ntnivmaxtotal,& 2946 My_XINCREDUITMUL(:,:,ic),X(s+1:s+nbnivest),dataset%nkd_max_by_fam,1,V(:,ic)) 2947 V(:na,ic) = my_listdesc(ic)%dataset%Y(ic,kd1:kd2)-V(:na,ic) 2948 s=s+my_listdesc(ic)%nbnivest 2949 end do 2950 2951 vpf=0 2952 !calcul de la vraissemblance de plein frere 2953 do kd=1,na 2954 vpf=vpf+dot_product(matmul(V(kd,:),VCI(:,:,kd)),V(kd,:)) 2955 end do 2956 2957 !finallement, la vraissemblance de la famille ip/jm: 2958 f=+0.5*vpf+0.5*sum(log(determ(:na))) 2959 !print *,f 2960 2961 end subroutine likelihood_ncar_h0_LU_family_withcd
likelihood_ncar_hn_family
[ Top ] [ m_qtlmap_incidence_multi ] [ Subroutines ]
NAME
likelihood_ncar_hn_family
DESCRIPTION
NOTE
multivariate analysis
(indicecar-1)*np +ip : SIG ip,indice_caracetere
ncar*np + (ic-1)*ncar+jc : RHO (ic,jc) residual correlation between two traits....
ncar*np + ncar*ncar + nbnivest
SOURCE
1805 subroutine likelihood_ncar_hn_family(ip,jm,n,x,VCI,determ,f,iuser,user) 1806 integer , intent(in) :: ip,jm,n 1807 real (kind=dp) ,dimension(n), intent(in) :: x 1808 real (kind=dp) , intent(inout) :: f 1809 integer , dimension(1), intent(inout) :: iuser 1810 real (kind=dp) ,dimension(1), intent(inout) :: user 1811 real(kind=dp),dimension(ncar,ncar) ,intent(in) :: VCI 1812 real(kind=dp) ,intent(in) :: determ 1813 1814 real(kind=dp),dimension(dataset%nkd_max_by_fam,ncar) :: V 1815 1816 1817 ! 1818 integer :: jjm,ig(my_listdesc(1)%nqtl),ifem,indf,indm,iq,ngg,iig,z,irank,ic,jc 1819 real(kind=dp) :: effm,pbr,vmere,vpf 1820 integer :: kd1,kd2,jj,nqtl,nbnivest,s,kd,na 1821 logical :: ok,valide 1822 1823 f=0.d0 1824 nqtl=my_listdesc(1)%nqtl 1825 1826 jjm=jm-nmp(ip) 1827 1828 if (dataset%lSires(ip)%full_sib(jjm)%lastKd>0) then 1829 kd1=dataset%lSires(ip)%full_sib(jjm)%firstKd 1830 kd2=dataset%lSires(ip)%full_sib(jjm)%lastKd 1831 na = kd2-kd1+1 1832 effm=dble(na) 1833 else 1834 return 1835 end if 1836 1837 vmere=0.d0 1838 if ( estime_multi(jm) ) ifem=iam_multi(repfem(jm)) 1839 !ngg : le nombre de genotype possible sur tous les qtls... 1840 ngg=1 1841 do iq=1,nqtl 1842 ig(iq)=ngenom(current_chr(iq),jm)+1 1843 ngg=ngg*(ngenom(current_chr(iq),jm+1)-ngenom(current_chr(iq),jm)) 1844 end do 1845 !pour toutes les combinaisons possibles des genotypes 1846 do iig=1,ngg 1847 pbr=1 1848 !on modifie la matrice d incidence pour les n qtl 1849 do iq=1,nqtl 1850 do ic=1,ncar 1851 indm=my_listdesc(ic)%ntniv_qtlsires(iq)+ip-1 1852 if ( estime_multi(jm) ) then 1853 !Si la mere est estimable, on place dans la matrice les 1854 !pdds dam et pdds sires (si ces effets sont estimables) 1855 indf=my_listdesc(ic)%ntniv_qtldams(iq)+ifem-1 1856 if ( my_listdesc(ic)%vecsol(indf)) then 1857 indf=my_listdesc(ic)%corr_niv_nivb(indf) ! vrai correspondance dans le tableau si estimable 1858 My_XINCREDUITMUL(kd1:kd2,indf,ic)=dataset%lSires(ip)%full_sib(jjm)%pmt(iq,(ig(iq)-ngenom(current_chr(iq),jm)),:) 1859 end if 1860 if ( my_listdesc(ic)%vecsol(indm) ) then 1861 indm=my_listdesc(ic)%corr_niv_nivb(indm) ! vrai correspondance dans le tableau si estimable 1862 My_XINCREDUITMUL(kd1:kd2,indm,ic)=dataset%lSires(ip)%full_sib(jjm)%ppt(iq,(ig(iq)-ngenom(current_chr(iq),jm)),:) 1863 end if 1864 else 1865 !la mere n est pas estimable, on place seulement les pdds males 1866 if ( my_listdesc(ic)%vecsol(indm) ) then 1867 indm=my_listdesc(ic)%corr_niv_nivb(indm) ! vrai correspondance dans le tableau si estimable 1868 My_XINCREDUITMUL(kd1:kd2,indm,ic)=dataset%lSires(ip)%ppt(iq,kd1:kd2) 1869 end if 1870 end if 1871 end do ! ic 1872 !print *,ig(iq)-ngenom(current_chr(iq),jm),ig(iq) 1873 pbr=pbr*probg(current_chr(iq),ig(iq)) 1874 end do ! iq 1875 ! Dans la famille du pere ip, on calcul le Y (performance auquel on soustrait l'ensemble des effets) 1876 ! pour chaque caractere 1877 s=np*ncar+ncar*(ncar-1)/2 1878 do ic=1,ncar 1879 nbnivest=my_listdesc(ic)%nbnivest 1880 call matmul_incidence(kd1,kd2,my_listdesc(ic),nd,ntnivmaxtotal,& 1881 My_XINCREDUITMUL(:,:,ic),X(s+1:s+nbnivest),dataset%nkd_max_by_fam,1,V(:,ic)) 1882 V(:na,ic) = my_listdesc(ic)%dataset%Y(ic,kd1:kd2)-V(:na,ic) 1883 s=s+my_listdesc(ic)%nbnivest 1884 end do 1885 vpf=0 1886 !calcul de la vraissemblance de plein frere 1887 do kd=1,na 1888 vpf=vpf+dot_product(matmul(V(kd,:),VCI),V(kd,:)) 1889 end do 1890 vmere=vmere+pbr*dexp(-0.5d0*vpf) 1891 1892 ! on increment 1893 ok=.true. 1894 do iq=1,nqtl 1895 if (ok) then 1896 if ((ig(iq) < ngenom(current_chr(iq),jm+1))) then 1897 ig(iq)=ig(iq)+1 1898 ok=.false. 1899 end if 1900 end if 1901 end do 1902 end do ! iig 1903 1904 1905 if (vmere == 0) then 1906 f=INIFINY_REAL_VALUE 1907 else 1908 !finallement, la vraissemblance de la famille ip/jm: 1909 f=f-dlog(vmere)+dble(kd2-kd1+1)*0.5*log(determ) 1910 end if 1911 ! end do 1912 ! print *,ip,jm,f 1913 ! stop 1914 1915 end subroutine likelihood_ncar_hn_family
likelihood_ncar_hn_family_LU
[ Top ] [ m_qtlmap_incidence_multi ] [ Subroutines ]
NAME
likelihood_ncar_hn_family_LU
DESCRIPTION
NOTE
multivariate analysis
(indicecar-1)*np +ip : SIG ip,indice_caracetere
ncar*np + (ic-1)*ncar+jc : RHO (ic,jc) residual correlation between two traits....
ncar*np + ncar*ncar + nbnivest
SOURCE
2976 subroutine likelihood_ncar_hn_family_LU(ip,jm,n,x,f,iuser,user) 2977 integer , intent(in) :: ip,jm,n 2978 real (kind=dp) ,dimension(n), intent(in) :: x 2979 real (kind=dp) , intent(inout) :: f 2980 integer , dimension(1), intent(inout) :: iuser 2981 real (kind=dp) ,dimension(1), intent(inout) :: user 2982 2983 !local 2984 real(kind=dp) ,dimension(ncar,ncar) :: VCI 2985 real(kind=dp) :: determ 2986 ! real(kind=dp),dimension(ncar,ncar) :: VCITemp 2987 real(kind=dp),dimension(dataset%nkd_max_by_fam,ncar) :: V 2988 2989 2990 ! 2991 integer :: jjm,ig(my_listdesc(1)%nqtl),ifem,indf,indm,iq,ngg,iig,z,irank,ic,jc 2992 real(kind=dp) :: effm,pbr,vmere,vpf 2993 integer :: kd1,kd2,jj,nqtl,nbnivest,s,kd,na 2994 logical :: ok,valide 2995 2996 f=0.d0 2997 nqtl=my_listdesc(1)%nqtl 2998 2999 call get_inv_residual_covariance_matrix_LU(ip,n,x,vci,determ,ok) 3000 3001 ! do ip=1,np 3002 ! VCITemp = VCI(ip,:,:) 3003 ! do jm=nmp(ip)+1,nmp(ip+1) 3004 3005 jjm=jm-nmp(ip) 3006 3007 if (dataset%lSires(ip)%full_sib(jjm)%lastKd>0) then 3008 kd1=dataset%lSires(ip)%full_sib(jjm)%firstKd 3009 kd2=dataset%lSires(ip)%full_sib(jjm)%lastKd 3010 na = kd2-kd1+1 3011 effm=dble(na) 3012 else 3013 return 3014 end if 3015 3016 vmere=0.d0 3017 if ( estime_multi(jm) ) ifem=iam_multi(repfem(jm)) 3018 !ngg : le nombre de genotype possible sur tous les qtls... 3019 ngg=1 3020 do iq=1,nqtl 3021 ig(iq)=ngenom(current_chr(iq),jm)+1 3022 ngg=ngg*(ngenom(current_chr(iq),jm+1)-ngenom(current_chr(iq),jm)) 3023 end do 3024 !pour toutes les combinaisons possibles des genotypes 3025 do iig=1,ngg 3026 pbr=1 3027 !on modifie la matrice d incidence pour les n qtl 3028 do iq=1,nqtl 3029 do ic=1,ncar 3030 indm=my_listdesc(ic)%ntniv_qtlsires(iq)+ip-1 3031 if ( estime_multi(jm) ) then 3032 !Si la mere est estimable, on place dans la matrice les 3033 !pdds dam et pdds sires (si ces effets sont estimables) 3034 indf=my_listdesc(ic)%ntniv_qtldams(iq)+ifem-1 3035 if ( my_listdesc(ic)%vecsol(indf)) then 3036 indf=my_listdesc(ic)%corr_niv_nivb(indf) ! vrai correspondance dans le tableau si estimable 3037 My_XINCREDUITMUL(kd1:kd2,indf,ic)=dataset%lSires(ip)%full_sib(jjm)%pmt(iq,(ig(iq)-ngenom(current_chr(iq),jm)),:) 3038 end if 3039 if ( my_listdesc(ic)%vecsol(indm) ) then 3040 indm=my_listdesc(ic)%corr_niv_nivb(indm) ! vrai correspondance dans le tableau si estimable 3041 My_XINCREDUITMUL(kd1:kd2,indm,ic)=dataset%lSires(ip)%full_sib(jjm)%ppt(iq,(ig(iq)-ngenom(current_chr(iq),jm)),:) 3042 end if 3043 else 3044 !la mere n est pas estimable, on place seulement les pdds males 3045 if ( my_listdesc(ic)%vecsol(indm) ) then 3046 indm=my_listdesc(ic)%corr_niv_nivb(indm) ! vrai correspondance dans le tableau si estimable 3047 My_XINCREDUITMUL(kd1:kd2,indm,ic)=dataset%lSires(ip)%ppt(iq,kd1:kd2) 3048 end if 3049 end if 3050 end do ! ic 3051 !print *,ig(iq)-ngenom(current_chr(iq),jm),ig(iq) 3052 pbr=pbr*probg(current_chr(iq),ig(iq)) 3053 end do ! iq 3054 ! Dans la famille du pere ip, on calcul le Y (performance auquel on soustrait l'ensemble des effets) 3055 ! pour chaque caractere 3056 s=np*ncar+ncar*(ncar-1)/2 3057 do ic=1,ncar 3058 nbnivest=my_listdesc(ic)%nbnivest 3059 call matmul_incidence(kd1,kd2,my_listdesc(ic),nd,ntnivmaxtotal,& 3060 My_XINCREDUITMUL(:,:,ic),X(s+1:s+nbnivest),dataset%nkd_max_by_fam,1,V(:,ic)) 3061 V(:na,ic) = my_listdesc(ic)%dataset%Y(ic,kd1:kd2)-V(:na,ic) 3062 s=s+my_listdesc(ic)%nbnivest 3063 end do 3064 vpf=0 3065 !calcul de la vraissemblance de plein frere 3066 do kd=1,na 3067 vpf=vpf+dot_product(matmul(V(kd,:),VCI),V(kd,:)) 3068 end do 3069 vmere=vmere+pbr*dexp(-0.5d0*vpf) 3070 3071 ! on increment 3072 ok=.true. 3073 do iq=1,nqtl 3074 if (ok) then 3075 if ((ig(iq) < ngenom(current_chr(iq),jm+1))) then 3076 ig(iq)=ig(iq)+1 3077 ok=.false. 3078 end if 3079 end if 3080 end do 3081 end do ! iig 3082 3083 3084 if (vmere == 0) then 3085 f=INIFINY_REAL_VALUE 3086 else 3087 !finallement, la vraissemblance de la famille ip/jm: 3088 f=f-dlog(vmere)+dble(kd2-kd1+1)*0.5*log(determ) 3089 end if 3090 ! end do 3091 ! end do 3092 3093 end subroutine likelihood_ncar_hn_family_LU
likelihood_ncar_hn_family_withcd
[ Top ] [ m_qtlmap_incidence_multi ] [ Subroutines ]
NAME
likelihood_ncar_hn_family_withcd
DESCRIPTION
NOTE
multivariate analysis
(indicecar-1)*np +ip : SIG ip,indice_caracetere
ncar*np + (ic-1)*ncar+jc : RHO (ic,jc) residual correlation between two traits....
ncar*np + ncar*ncar + nbnivest
SOURCE
2140 subroutine likelihood_ncar_hn_family_withcd(ip,jm,n,x,f,iuser,user) 2141 integer , intent(in) :: ip,jm,n 2142 real (kind=dp) ,dimension(n), intent(in) :: x 2143 real (kind=dp) , intent(inout) :: f 2144 integer , dimension(1), intent(inout) :: iuser 2145 real (kind=dp) ,dimension(1), intent(inout) :: user 2146 2147 2148 real(kind=dp),dimension(ncar,ncar,dataset%nkd_max_by_fam) :: VCI 2149 real(kind=dp),dimension(dataset%nkd_max_by_fam) :: determ 2150 real(kind=dp),dimension(dataset%nkd_max_by_fam,ncar) :: V 2151 ! 2152 integer :: jjm,ig(my_listdesc(1)%nqtl),ifem,indf,indm,iq,ngg,iig,z,irank,ic,jc 2153 real(kind=dp) :: effm,pbr,vmere,vpf 2154 integer :: kd1,kd2,jj,nqtl,nbnivest,s,kd,kkk,na 2155 logical :: ok,valide 2156 2157 determ=0.d0 2158 kd=1 2159 do kkk=ndm(jm)+1,ndm(jm+1) 2160 if ( count(presentc(:,kkk)) == ncar ) then 2161 call get_inv_residual_covariance_matrix_cd(ip,kkk,n,x,vci(:,:,kd),determ(kd),ok) 2162 if ( .not. ok ) then 2163 f=INIFINY_REAL_VALUE 2164 return 2165 end if 2166 kd = kd + 1 2167 end if 2168 end do 2169 2170 f=0.d0 2171 nqtl=my_listdesc(1)%nqtl 2172 2173 jjm=jm-nmp(ip) 2174 2175 if (dataset%lSires(ip)%full_sib(jjm)%lastKd>0) then 2176 kd1=dataset%lSires(ip)%full_sib(jjm)%firstKd 2177 kd2=dataset%lSires(ip)%full_sib(jjm)%lastKd 2178 na = kd2-kd1+1 2179 effm=dble(na) 2180 else 2181 return 2182 end if 2183 2184 vmere=0.d0 2185 if ( estime_multi(jm) ) ifem=iam_multi(repfem(jm)) 2186 !ngg : le nombre de genotype possible sur tous les qtls... 2187 ngg=1 2188 do iq=1,nqtl 2189 ig(iq)=ngenom(current_chr(iq),jm)+1 2190 ngg=ngg*(ngenom(current_chr(iq),jm+1)-ngenom(current_chr(iq),jm)) 2191 end do 2192 !pour toutes les combinaisons possibles des genotypes 2193 do iig=1,ngg 2194 pbr=1 2195 !on modifie la matrice d incidence pour les n qtl 2196 do iq=1,nqtl 2197 do ic=1,ncar 2198 indm=my_listdesc(ic)%ntniv_qtlsires(iq)+ip-1 2199 if ( estime_multi(jm) ) then 2200 !Si la mere est estimable, on place dans la matrice les 2201 !pdds dam et pdds sires (si ces effets sont estimables) 2202 indf=my_listdesc(ic)%ntniv_qtldams(iq)+ifem-1 2203 if ( my_listdesc(ic)%vecsol(indf)) then 2204 indf=my_listdesc(ic)%corr_niv_nivb(indf) ! vrai correspondance dans le tableau si estimable 2205 My_XINCREDUITMUL(kd1:kd2,indf,ic)=dataset%lSires(ip)%full_sib(jjm)%pmt(iq,(ig(iq)-ngenom(current_chr(iq),jm)),:) 2206 end if 2207 if ( my_listdesc(ic)%vecsol(indm) ) then 2208 indm=my_listdesc(ic)%corr_niv_nivb(indm) ! vrai correspondance dans le tableau si estimable 2209 My_XINCREDUITMUL(kd1:kd2,indm,ic)=dataset%lSires(ip)%full_sib(jjm)%ppt(iq,(ig(iq)-ngenom(current_chr(iq),jm)),:) 2210 end if 2211 else 2212 !la mere n est pas estimable, on place seulement les pdds males 2213 if ( my_listdesc(ic)%vecsol(indm) ) then 2214 indm=my_listdesc(ic)%corr_niv_nivb(indm) ! vrai correspondance dans le tableau si estimable 2215 My_XINCREDUITMUL(kd1:kd2,indm,ic)=dataset%lSires(ip)%ppt(iq,kd1:kd2) 2216 end if 2217 end if 2218 end do ! ic 2219 !print *,ig(iq)-ngenom(current_chr(iq),jm),ig(iq) 2220 pbr=pbr*probg(current_chr(iq),ig(iq)) 2221 end do ! iq 2222 ! Dans la famille du pere ip, on calcul le Y (performance auquel on soustrait l'ensemble des effets) 2223 ! pour chaque caractere 2224 s=np*ncar+ncar*(ncar-1)/2 2225 do ic=1,ncar 2226 nbnivest=my_listdesc(ic)%nbnivest 2227 call matmul_incidence(kd1,kd2,my_listdesc(ic),nd,ntnivmaxtotal,& 2228 My_XINCREDUITMUL(:,:,ic),X(s+1:s+nbnivest),dataset%nkd_max_by_fam,1,V(:,ic)) 2229 V(:na,ic) = my_listdesc(ic)%dataset%Y(ic,kd1:kd2)-V(:na,ic) 2230 s=s+my_listdesc(ic)%nbnivest 2231 end do 2232 vpf=0 2233 !calcul de la vraissemblance de plein frere 2234 do kd=1,na 2235 vpf=vpf+dot_product(matmul(V(kd,:),vci(:,:,kd)),V(kd,:)) 2236 end do 2237 vmere=vmere+pbr*dexp(-0.5d0*vpf) 2238 ! on increment 2239 ok=.true. 2240 do iq=1,nqtl 2241 if (ok) then 2242 if ((ig(iq) < ngenom(current_chr(iq),jm+1))) then 2243 ig(iq)=ig(iq)+1 2244 ok=.false. 2245 end if 2246 end if 2247 end do 2248 end do ! iig 2249 2250 2251 if (vmere == 0) then 2252 f=INIFINY_REAL_VALUE 2253 else 2254 !finallement, la vraissemblance de la famille ip/jm: 2255 f=f-dlog(vmere)+0.5*sum(log(determ(:na))) 2256 end if 2257 ! print *,f 2258 2259 end subroutine likelihood_ncar_hn_family_withcd
likelihood_ncar_hn_LU
[ Top ] [ m_qtlmap_incidence_multi ] [ Subroutines ]
NAME
likelihood_ncar_hn_LU
DESCRIPTION
NOTE
multivariate analysis
(indicecar-1)*np +ip : SIG ip,indice_caracetere
ncar*np + (ic-1)*ncar+jc : RHO (ic,jc) residual correlation between two traits....
ncar*np + ncar*ncar + nbnivest
SOURCE
3107 subroutine likelihood_ncar_hn_LU(n,x,f,iuser,user) 3108 integer , intent(in) :: n 3109 real (kind=dp) ,dimension(n), intent(in) :: x 3110 real (kind=dp) , intent(inout) :: f 3111 integer , dimension(1), intent(inout) :: iuser 3112 real (kind=dp) ,dimension(1), intent(inout) :: user 3113 3114 !local 3115 real(kind=dp) ,dimension(ncar,ncar) :: VCI 3116 real(kind=dp) :: determ 3117 real(kind=dp),dimension(dataset%nkd_max_by_fam,ncar) :: V 3118 3119 3120 ! 3121 integer :: jjm,ig(my_listdesc(1)%nqtl),ifem,indf,indm,iq,ngg,iig,z,irank,ic,jc 3122 real(kind=dp) :: effm,pbr,vmere,vpf 3123 integer :: kd1,kd2,jj,nqtl,nbnivest,s,kd,na,ip,jm 3124 logical :: ok,valide 3125 3126 f=0.d0 3127 nqtl=my_listdesc(1)%nqtl 3128 3129 do ip=1,np 3130 call get_inv_residual_covariance_matrix_LU(ip,n,x,vci,determ,ok) 3131 do jm=nmp(ip)+1,nmp(ip+1) 3132 3133 jjm=jm-nmp(ip) 3134 3135 if (dataset%lSires(ip)%full_sib(jjm)%lastKd>0) then 3136 kd1=dataset%lSires(ip)%full_sib(jjm)%firstKd 3137 kd2=dataset%lSires(ip)%full_sib(jjm)%lastKd 3138 na = kd2-kd1+1 3139 effm=dble(na) 3140 else 3141 cycle 3142 end if 3143 3144 vmere=0.d0 3145 if ( estime_multi(jm) ) ifem=iam_multi(repfem(jm)) 3146 !ngg : le nombre de genotype possible sur tous les qtls... 3147 ngg=1 3148 do iq=1,nqtl 3149 ig(iq)=ngenom(current_chr(iq),jm)+1 3150 ngg=ngg*(ngenom(current_chr(iq),jm+1)-ngenom(current_chr(iq),jm)) 3151 end do 3152 !pour toutes les combinaisons possibles des genotypes 3153 do iig=1,ngg 3154 pbr=1 3155 !on modifie la matrice d incidence pour les n qtl 3156 do iq=1,nqtl 3157 do ic=1,ncar 3158 indm=my_listdesc(ic)%ntniv_qtlsires(iq)+ip-1 3159 if ( estime_multi(jm) ) then 3160 !Si la mere est estimable, on place dans la matrice les 3161 !pdds dam et pdds sires (si ces effets sont estimables) 3162 indf=my_listdesc(ic)%ntniv_qtldams(iq)+ifem-1 3163 if ( my_listdesc(ic)%vecsol(indf)) then 3164 indf=my_listdesc(ic)%corr_niv_nivb(indf) ! vrai correspondance dans le tableau si estimable 3165 My_XINCREDUITMUL(kd1:kd2,indf,ic)=dataset%lSires(ip)%full_sib(jjm)%pmt(iq,(ig(iq)-ngenom(current_chr(iq),jm)),:) 3166 end if 3167 if ( my_listdesc(ic)%vecsol(indm) ) then 3168 indm=my_listdesc(ic)%corr_niv_nivb(indm) ! vrai correspondance dans le tableau si estimable 3169 My_XINCREDUITMUL(kd1:kd2,indm,ic)=dataset%lSires(ip)%full_sib(jjm)%ppt(iq,(ig(iq)-ngenom(current_chr(iq),jm)),:) 3170 end if 3171 else 3172 !la mere n est pas estimable, on place seulement les pdds males 3173 if ( my_listdesc(ic)%vecsol(indm) ) then 3174 indm=my_listdesc(ic)%corr_niv_nivb(indm) ! vrai correspondance dans le tableau si estimable 3175 My_XINCREDUITMUL(kd1:kd2,indm,ic)=dataset%lSires(ip)%ppt(iq,kd1:kd2) 3176 end if 3177 end if 3178 end do ! ic 3179 !print *,ig(iq)-ngenom(current_chr(iq),jm),ig(iq) 3180 pbr=pbr*probg(current_chr(iq),ig(iq)) 3181 end do ! iq 3182 ! Dans la famille du pere ip, on calcul le Y (performance auquel on soustrait l'ensemble des effets) 3183 ! pour chaque caractere 3184 s=np*ncar+ncar*(ncar-1)/2 3185 do ic=1,ncar 3186 nbnivest=my_listdesc(ic)%nbnivest 3187 call matmul_incidence(kd1,kd2,my_listdesc(ic),nd,ntnivmaxtotal,& 3188 My_XINCREDUITMUL(:,:,ic),X(s+1:s+nbnivest),dataset%nkd_max_by_fam,1,V(:,ic)) 3189 V(:na,ic) = my_listdesc(ic)%dataset%Y(ic,kd1:kd2)-V(:na,ic) 3190 s=s+my_listdesc(ic)%nbnivest 3191 end do 3192 vpf=0 3193 !calcul de la vraissemblance de plein frere 3194 do kd=1,na 3195 vpf=vpf+dot_product(matmul(V(kd,:),VCI),V(kd,:)) 3196 end do 3197 vmere=vmere+pbr*dexp(-0.5d0*vpf) 3198 3199 ! on increment 3200 ok=.true. 3201 do iq=1,nqtl 3202 if (ok) then 3203 if ((ig(iq) < ngenom(current_chr(iq),jm+1))) then 3204 ig(iq)=ig(iq)+1 3205 ok=.false. 3206 end if 3207 end if 3208 end do 3209 end do ! iig 3210 3211 3212 if (vmere == 0) then 3213 f=INIFINY_REAL_VALUE 3214 else 3215 !finallement, la vraissemblance de la famille ip/jm: 3216 f=f-dlog(vmere)+dble(kd2-kd1+1)*0.5*log(determ) 3217 end if 3218 end do 3219 end do 3220 3221 end subroutine likelihood_ncar_hn_LU
likelihood_ncar_hn_LU_family_withcd
[ Top ] [ m_qtlmap_incidence_multi ] [ Subroutines ]
NAME
likelihood_ncar_hn_LU_family_withcd
DESCRIPTION
NOTE
multivariate analysis
(indicecar-1)*np +ip : SIG ip,indice_caracetere
ncar*np + (ic-1)*ncar+jc : RHO (ic,jc) residual correlation between two traits....
ncar*np + ncar*ncar + nbnivest
SOURCE
3235 subroutine likelihood_ncar_hn_LU_family_withcd(ip,jm,n,x,f,iuser,user) 3236 integer , intent(in) :: ip,jm,n 3237 real (kind=dp) ,dimension(n), intent(in) :: x 3238 real (kind=dp) , intent(inout) :: f 3239 integer , dimension(1), intent(inout) :: iuser 3240 real (kind=dp) ,dimension(1), intent(inout) :: user 3241 3242 3243 real(kind=dp),dimension(ncar,ncar,dataset%nkd_max_by_fam) :: VCI 3244 real(kind=dp),dimension(dataset%nkd_max_by_fam) :: determ 3245 real(kind=dp),dimension(dataset%nkd_max_by_fam,ncar) :: V 3246 ! 3247 integer :: jjm,ig(my_listdesc(1)%nqtl),ifem,indf,indm,iq,ngg,iig,z,irank,ic,jc 3248 real(kind=dp) :: effm,pbr,vmere,vpf 3249 integer :: kd1,kd2,jj,nqtl,nbnivest,s,kd,kkk,na 3250 logical :: ok,valide 3251 3252 determ=0.d0 3253 kd=1 3254 do kkk=ndm(jm)+1,ndm(jm+1) 3255 if ( count(presentc(:,kkk)) == ncar ) then 3256 call get_inv_residual_covariance_matrix_LU_cd(ip,kkk,n,x,vci(:,:,kd),determ(kd),ok) 3257 ! print *,kd,determ(kd) 3258 if ( .not. ok ) then 3259 f=INIFINY_REAL_VALUE 3260 return 3261 end if 3262 kd = kd + 1 3263 end if 3264 end do 3265 3266 3267 f=0.d0 3268 nqtl=my_listdesc(1)%nqtl 3269 3270 jjm=jm-nmp(ip) 3271 3272 if (dataset%lSires(ip)%full_sib(jjm)%lastKd>0) then 3273 kd1=dataset%lSires(ip)%full_sib(jjm)%firstKd 3274 kd2=dataset%lSires(ip)%full_sib(jjm)%lastKd 3275 na = kd2-kd1+1 3276 effm=dble(na) 3277 else 3278 return 3279 end if 3280 3281 vmere=0.d0 3282 if ( estime_multi(jm) ) ifem=iam_multi(repfem(jm)) 3283 !ngg : le nombre de genotype possible sur tous les qtls... 3284 ngg=1 3285 do iq=1,nqtl 3286 ig(iq)=ngenom(current_chr(iq),jm)+1 3287 ngg=ngg*(ngenom(current_chr(iq),jm+1)-ngenom(current_chr(iq),jm)) 3288 end do 3289 !pour toutes les combinaisons possibles des genotypes 3290 do iig=1,ngg 3291 pbr=1 3292 !on modifie la matrice d incidence pour les n qtl 3293 do iq=1,nqtl 3294 do ic=1,ncar 3295 indm=my_listdesc(ic)%ntniv_qtlsires(iq)+ip-1 3296 if ( estime_multi(jm) ) then 3297 !Si la mere est estimable, on place dans la matrice les 3298 !pdds dam et pdds sires (si ces effets sont estimables) 3299 indf=my_listdesc(ic)%ntniv_qtldams(iq)+ifem-1 3300 if ( my_listdesc(ic)%vecsol(indf)) then 3301 indf=my_listdesc(ic)%corr_niv_nivb(indf) ! vrai correspondance dans le tableau si estimable 3302 My_XINCREDUITMUL(kd1:kd2,indf,ic)=dataset%lSires(ip)%full_sib(jjm)%pmt(iq,(ig(iq)-ngenom(current_chr(iq),jm)),:) 3303 end if 3304 if ( my_listdesc(ic)%vecsol(indm) ) then 3305 indm=my_listdesc(ic)%corr_niv_nivb(indm) ! vrai correspondance dans le tableau si estimable 3306 My_XINCREDUITMUL(kd1:kd2,indm,ic)=dataset%lSires(ip)%full_sib(jjm)%ppt(iq,(ig(iq)-ngenom(current_chr(iq),jm)),:) 3307 end if 3308 else 3309 !la mere n est pas estimable, on place seulement les pdds males 3310 if ( my_listdesc(ic)%vecsol(indm) ) then 3311 indm=my_listdesc(ic)%corr_niv_nivb(indm) ! vrai correspondance dans le tableau si estimable 3312 My_XINCREDUITMUL(kd1:kd2,indm,ic)=dataset%lSires(ip)%ppt(iq,kd1:kd2) 3313 end if 3314 end if 3315 end do ! ic 3316 !print *,ig(iq)-ngenom(current_chr(iq),jm),ig(iq) 3317 pbr=pbr*probg(current_chr(iq),ig(iq)) 3318 end do ! iq 3319 ! Dans la famille du pere ip, on calcul le Y (performance auquel on soustrait l'ensemble des effets) 3320 ! pour chaque caractere 3321 s=np*ncar+ncar*(ncar-1)/2 3322 do ic=1,ncar 3323 nbnivest=my_listdesc(ic)%nbnivest 3324 call matmul_incidence(kd1,kd2,my_listdesc(ic),nd,ntnivmaxtotal,& 3325 My_XINCREDUITMUL(:,:,ic),X(s+1:s+nbnivest),dataset%nkd_max_by_fam,1,V(:,ic)) 3326 V(:na,ic) = my_listdesc(ic)%dataset%Y(ic,kd1:kd2)-V(:na,ic) 3327 s=s+my_listdesc(ic)%nbnivest 3328 end do 3329 vpf=0 3330 !calcul de la vraissemblance de plein frere 3331 do kd=1,na 3332 vpf=vpf+dot_product(matmul(V(kd,:),vci(:,:,kd)),V(kd,:)) 3333 end do 3334 vmere=vmere+pbr*dexp(-0.5d0*vpf) 3335 ! on increment 3336 ok=.true. 3337 do iq=1,nqtl 3338 if (ok) then 3339 if ((ig(iq) < ngenom(current_chr(iq),jm+1))) then 3340 ig(iq)=ig(iq)+1 3341 ok=.false. 3342 end if 3343 end if 3344 end do 3345 end do ! iig 3346 3347 3348 if (vmere == 0) then 3349 f=INIFINY_REAL_VALUE 3350 else 3351 !finallement, la vraissemblance de la famille ip/jm: 3352 f=f-dlog(vmere)+0.5*sum(log(determ(:na))) 3353 end if 3354 ! print *,f 3355 3356 end subroutine likelihood_ncar_hn_LU_family_withcd
model_optim_multi_family
[ Top ] [ m_qtlmap_incidence_multi ] [ Subroutines ]
NAME
model_optim_multi_family
DESCRIPTION
NOTE
SOURCE
1481 subroutine model_optim_multi_family(xinc,listDesc,workstruct,sigsquareEstime,rhoiestim,Bestim,& 1482 performPrecision,FUNCT_PART,FUNCT_PART_CENSURE,tConf,tempForConfusion) 1483 real (kind=dp) , dimension(ncar,nd,ntnivmaxtotal), intent(in) :: xinc 1484 type(INCIDENCE_GEN_STRUCT) ,intent(inout) :: workstruct 1485 real (kind=dp) , dimension(ncar,np) ,intent(out) :: sigsquareEstime 1486 real (kind=dp) , dimension(ncar,ntnivmaxtotal) ,intent(out) :: Bestim 1487 real (kind=dp) , dimension(ncar,ncar) ,intent(out) :: rhoiestim 1488 type(INCIDENCE_TYPE) , dimension(ncar) , intent(inout),target :: listDesc 1489 logical ,intent(in) :: performPrecision,tConf 1490 real (kind=dp),dimension(ntnivmaxtotal,ntnivmaxtotal,ncar) ,intent(out) :: tempForConfusion 1491 external :: FUNCT_PART,FUNCT_PART_CENSURE 1492 1493 integer :: j,i,ip,ifail,ibound,npar 1494 ! real(kind=dp) ,dimension(np+ntnivmax) :: par 1495 real (kind=dp) , dimension(ntnivmaxtotal,ntnivmaxtotal) :: XX 1496 real (kind=dp) , dimension(ntnivmaxtotal,ntnivmaxtotal) :: triang 1497 integer :: iuser(1),ix,kd1,kd2,jm,jjm,ic,s,jc,iic 1498 real (kind=dp) ,dimension(1) :: user 1499 logical ,dimension(ncar,ntnivmaxtotal) :: lastvecsol 1500 real(kind=dp) :: vci_t(ncar,ncar,np),determ_t(np),fp(np),fm(nm),r,tempLin(ncar*ncar) 1501 1502 !Filter for the mimimization - dim (np,nm,nd) 1503 logical , dimension(np,nm,ncar*np+ncar*(ncar-1)/2+ncar*ntnivmaxtotal) :: filter_inc 1504 !Filter for calcul on the vci matrix : dim np, 1505 logical , dimension(np,ncar*np+ncar*(ncar-1)/2+ncar*ntnivmaxtotal) ,target :: filter_vci 1506 logical , dimension(:,:) , pointer :: p_filter_vci 1507 1508 1509 iuser=0 1510 user=0 1511 my_listDesc=>listDesc 1512 dataset => listDesc(1)%dataset 1513 allocate (my_xincreduitmul(nd,ntnivmaxtotal,ncar)) 1514 1515 !filter car initialisation 1516 filter_vci=.false. 1517 do ip=1,np 1518 do ic=1,ncar 1519 filter_vci(ip,(ic-1)*np+ip)=.true. 1520 end do 1521 end do 1522 1523 filter_vci(:,ncar*np+1:ncar*np+ncar*(ncar-1)/2)=.true. 1524 1525 1526 my_xincreduitmul=0.d0 1527 do ic=1,ncar 1528 lastvecsol(ic,:my_listDesc(ic)%ntniv)=my_listDesc(ic)%vecsol(:my_listDesc(ic)%ntniv) 1529 ! create X'.X matrix from incidence matrix 1530 call model_XT_X(xinc(ic,:,:),my_listDesc(ic),XX) 1531 ! Check all parameters to remove from the estimation 1532 call estim_cholesky(XX,my_listDesc(ic),ntnivmaxtotal,triang) 1533 ! compute the precision of each parameter 1534 if (performPrecision) call get_precision(XX,tempForConfusion(:,:,ic),my_listDesc(ic)) 1535 call set_corrxinc(xinc(ic,:,:),my_listDesc(ic),my_xincreduitmul(:,:,ic)) 1536 !optimisation : on sauvegarde les index des elements non nul de la matrice reduite 1537 call fill_nonull_elements(my_listDesc(ic),size(my_xincreduitmul,1),size(my_xincreduitmul,2),my_xincreduitmul(:,:,ic)) 1538 !call debug_write_incidence(xinc(ic,:,:),my_listDesc(ic)) 1539 end do 1540 1541 ! Optimisation de la vraisemblance a la position dx 1542 ifail=1 1543 ibound=0 1544 npar=ncar*np+ncar*(ncar-1)/2 1545 1546 j=ncar*np+ncar*(ncar-1)/2 1547 do ic=1,ncar 1548 npar=npar+my_listDesc(ic)%nbnivest 1549 1550 if (count(lastvecsol(ic,:))>0) then ! autre que initialisation 1551 do i=1,my_listDesc(ic)%ntniv 1552 if(my_listDesc(ic)%vecsol(i))then 1553 j=j+1 1554 if(.not. lastvecsol(ic,i)) my_listDesc(ic)%par(j)=0.d0 1555 end if 1556 end do 1557 else ! initialisation des correlations phenotypiques 1558 s=ncar*np 1559 do iic=2,ncar 1560 do jc=1,iic-1 1561 s=s+1 1562 my_listDesc(ic)%par(s)=log( (1+RhoP(iic,jc)) / (1-(RhoP(iic,jc)))) 1563 ! print *, my_listDesc(ic)%par(s) 1564 end do 1565 end do 1566 end if 1567 end do 1568 ! allocate (filter_inc(np,npar)) 1569 s=np*ncar+ncar*(ncar-1)/2 1570 filter_inc=.true. 1571 filter_inc(:,:,1:np*ncar)=.false. 1572 filter_inc(:,:,np*ncar+1:s)=.true. !rho 1573 do ic=1,ncar 1574 do ip=1,np 1575 filter_inc(ip,nmp(ip)+1:nmp(ip+1),(ic-1)*np+ip)=.true. 1576 jjm=0 1577 do jm=nmp(ip)+1,nmp(ip+1) 1578 jjm=jjm+1 1579 if (dataset%lSires(ip)%full_sib(jjm)%lastKd>0) then 1580 kd1=dataset%lSires(ip)%full_sib(jjm)%firstkd 1581 kd2=dataset%lSires(ip)%full_sib(jjm)%lastkd 1582 filter_inc(ip,jm,s+1:s+my_listDesc(ic)%nbnivest)=any(my_xincreduitmul(kd1:kd2,:my_listDesc(ic)%nbnivest,ic)/=0.d0,dim=1) 1583 else 1584 ! print *,'pas de desc....' 1585 filter_inc(ip,jm,s+1:s+my_listDesc(ic)%nbnivest)=.false. 1586 end if 1587 end do 1588 end do 1589 s=s+my_listDesc(ic)%nbnivest 1590 end do 1591 1592 if ( datasetUser%na>0 ) then 1593 !prise en compte des donnees censures 1594 call minimizing_funct_family(npar,ibound,FUNCT_PART_CENSURE,filter_inc,fm,fp,& 1595 my_listDesc(1)%borni,my_listDesc(1)%borns,my_listDesc(1)%par,listDesc(1)%fmax,iuser,user,ifail) 1596 else 1597 p_filter_vci=>filter_vci 1598 !pas de donnee censure=> optimisation pour le calcul de l'inverse de matrice de covariance 1599 call minimizing_funct_family_multi(npar,ibound,FUNCT_PART,filter_inc,p_filter_vci,vci_t,determ_t,fm,fp,& 1600 my_listDesc(1)%borni,my_listDesc(1)%borns,my_listDesc(1)%par,listDesc(1)%fmax,iuser,user,ifail) 1601 end if 1602 1603 Bestim=0.d0 1604 s=np*ncar+ncar*(ncar-1)/2 1605 !getting standart deviation 1606 do ic=1,ncar 1607 do ip=1,np 1608 sigsquareEstime(ic,ip)=my_listDesc(1)%par((ic-1)*np+ip)*my_listDesc(1)%par((ic-1)*np+ip) 1609 end do 1610 !The solution 1611 Bestim(ic,:my_listDesc(ic)%nbnivest)=my_listDesc(1)%par(s+1:s+my_listDesc(ic)%nbnivest) 1612 s=s+my_listDesc(ic)%nbnivest 1613 end do 1614 1615 rhoiestim=0.d0 1616 s=ncar*np 1617 do ic=2,ncar 1618 do jc=1,ic-1 1619 s=s+1 1620 ! *** RESCALE OFI *** 1621 !rhoiestim(jc,ic)=my_listDesc(1)%par(s) 1622 r = my_listDesc(1)%par(s) 1623 rhoiestim(jc,ic)=(dexp(r)-1.d0)/(dexp(r)+1.d0) 1624 rhoiestim(ic,jc)=rhoiestim(jc,ic) 1625 end do 1626 end do 1627 1628 i=0 1629 do ic=2,ncar 1630 do jc=1,ic-1 1631 i=i+1 1632 tempLin(i)=rhoiestim(ic,jc) 1633 end do 1634 end do 1635 1636 i=0 1637 do jc=1,ncar-1 1638 do ic=jc+1,ncar 1639 i = i + 1 1640 rhoiestim(ic,jc) = tempLin(i) 1641 rhoiestim(jc,ic) = rhoiestim(ic,jc) 1642 end do 1643 end do 1644 1645 deallocate (my_xincreduitmul) 1646 1647 end subroutine model_optim_multi_family
model_optim_multi_h0
[ Top ] [ m_qtlmap_incidence_multi ] [ Subroutines ]
NAME
model_optim_multi_h0
DESCRIPTION
NOTE
SOURCE
1390 subroutine model_optim_multi_h0(xinc,listDesc,workstruct,sigsquareEstime,rhoiestim,Bestim,performPrecision) 1391 real (kind=dp) , dimension(ncar,nd,ntnivmaxtotal), intent(in) :: xinc 1392 type(INCIDENCE_GEN_STRUCT) ,intent(inout) :: workstruct 1393 real (kind=dp) , dimension(ncar,ntnivmaxtotal) ,intent(out) :: Bestim 1394 real (kind=dp) , dimension(ncar,ncar) ,intent(out) :: rhoiestim 1395 real (kind=dp) , dimension(ncar,np) ,intent(out) :: sigsquareEstime 1396 type(INCIDENCE_TYPE) , dimension(ncar) , intent(inout) :: listDesc 1397 logical ,intent(in) :: performPrecision 1398 1399 real(kind=dp) , dimension(ntnivmaxtotal,ntnivmaxtotal,ncar) :: xxx 1400 1401 1402 call model_optim_multi_family(xinc,listDesc,workstruct,sigsquareEstime,rhoiestim,Bestim,& 1403 performPrecision,likelihood_ncar_h0_family,likelihood_ncar_h0_family_withcd,.true.,xxx) 1404 1405 !stop 1406 !workstruct%sigsquareEstime=osig*osig 1407 1408 end subroutine model_optim_multi_h0
model_optim_multi_h0_LU
[ Top ] [ m_qtlmap_incidence_multi ] [ Subroutines ]
NAME
model_optim_multi_h0_LU
DESCRIPTION
NOTE
SOURCE
2276 subroutine model_optim_multi_h0_LU(xinc,listDesc,workstruct,sigsquareEstime,rhoiestim,Bestim,performPrecision) 2277 real (kind=dp) , dimension(ncar,nd,ntnivmaxtotal), intent(in) :: xinc 2278 type(INCIDENCE_GEN_STRUCT) ,intent(inout) :: workstruct 2279 real (kind=dp) , dimension(ncar,ntnivmaxtotal) ,intent(out) :: Bestim 2280 real (kind=dp) , dimension(ncar,ncar) ,intent(out) :: rhoiestim 2281 real (kind=dp) , dimension(ncar,np) ,intent(out) :: sigsquareEstime 2282 type(INCIDENCE_TYPE) , dimension(ncar) , intent(inout) :: listDesc 2283 logical ,intent(in) :: performPrecision 2284 2285 real(kind=dp) , dimension(ntnivmaxtotal,ntnivmaxtotal,ncar) :: xxx 2286 2287 2288 call model_optim_multi_LU(xinc,listDesc,workstruct,sigsquareEstime,rhoiestim,Bestim,& 2289 performPrecision,likelihood_ncar_h0_family_LU,likelihood_ncar_h0_LU_family_withcd,.true.,xxx) 2290 2291 ! call model_optim_multi_LU(xinc,listDesc,workstruct,sigsquareEstime,rhoiestim,Bestim,& 2292 ! performPrecision,likelihood_ncar_h0_LU,likelihood_ncar_h0_LU,.true.,xxx) 2293 2294 2295 end subroutine model_optim_multi_h0_LU
model_optim_multi_hn
[ Top ] [ m_qtlmap_incidence_multi ] [ Subroutines ]
NAME
model_optim_multi_hn
DESCRIPTION
NOTE
SOURCE
1419 subroutine model_optim_multi_hn(xinc,listDesc,curPos,workstruct,sigsquareEstime,rhoiestim,Bestim,& 1420 performPrecision,tConf,tempForConfusion,invlrt) 1421 real (kind=dp) , dimension(ncar,nd,ntnivmaxtotal), intent(in) :: xinc 1422 type(INCIDENCE_GEN_STRUCT) ,intent(inout) :: workstruct 1423 type(POSITION_LRT_INCIDENCE) ,intent(inout) :: curPos 1424 real (kind=dp) , dimension(ncar,ntnivmaxtotal) ,intent(out) :: Bestim 1425 real (kind=dp) , dimension(ncar,ncar) ,intent(out) :: rhoiestim 1426 real (kind=dp) , dimension(ncar,np) ,intent(out) :: sigsquareEstime 1427 type(INCIDENCE_TYPE) , dimension(ncar) , intent(inout) :: listDesc 1428 logical ,intent(in) :: performPrecision,tConf 1429 real (kind=dp),dimension(ntnivmaxtotal,ntnivmaxtotal,ncar) ,intent(out) :: tempForConfusion 1430 logical ,intent(in) :: invlrt 1431 1432 integer :: i,hypothesis 1433 character(len=LEN_L) :: dx 1434 1435 allocate (current_chr(workstruct%nqtl)) 1436 current_chr=curPos%listChr 1437 1438 call model_optim_multi_family(xinc,listDesc,workstruct,sigsquareEstime,rhoiestim,Bestim,& 1439 performPrecision,likelihood_ncar_hn_family,likelihood_ncar_hn_family_withcd,tConf,tempForConfusion) 1440 ! print *,"F:",workstruct%fmax 1441 1442 1443 if ((listDesc(1)%fmax == INIFINY_REAL_VALUE) ) then 1444 dx="" 1445 do i=1,workstruct%nqtl-1 1446 dx=trim(dx)//trim(str(absi(curPos%listChr(i),curPos%listN(i))))//"," 1447 end do 1448 1449 dx=trim(dx)//str(absi(curPos%listChr(workstruct%nqtl),curPos%listN(workstruct%nqtl))) 1450 call log_mess("dx ["//trim(dx)// "]. Can not optimize likelihood....The start point is reinitializing",WARNING_DEF) 1451 call init_startpoint(workStruct,listDesc(1)) 1452 curPos%lrt=0.d0 1453 listDesc(1)%fmax=0.d0 1454 return 1455 end if 1456 1457 !compute LRT 1458 if (invLrt) then 1459 do hypothesis=1,workstruct%nqtl 1460 curPos%lrt(hypothesis)=-2.d0*(workstruct%fnqtl(hypothesis)-listDesc(1)%fmax) 1461 end do 1462 else 1463 do hypothesis=1,workstruct%nqtl 1464 curPos%lrt(hypothesis)=-2.d0*(listDesc(1)%fmax-workstruct%fnqtl(hypothesis)) 1465 end do 1466 end if 1467 1468 deallocate (current_chr) 1469 1470 end subroutine model_optim_multi_hn
model_optim_multi_hn_LU
[ Top ] [ m_qtlmap_incidence_multi ] [ Subroutines ]
NAME
model_optim_multi_hn_LU
DESCRIPTION
NOTE
SOURCE
2306 subroutine model_optim_multi_hn_LU(xinc,listDesc,curPos,workstruct,sigsquareEstime,rhoiestim,Bestim,& 2307 performPrecision,tConf,tempForConfusion,invlrt) 2308 real (kind=dp) , dimension(ncar,nd,ntnivmaxtotal), intent(in) :: xinc 2309 type(INCIDENCE_GEN_STRUCT) ,intent(inout) :: workstruct 2310 type(POSITION_LRT_INCIDENCE) ,intent(inout) :: curPos 2311 real (kind=dp) , dimension(ncar,ntnivmaxtotal) ,intent(out) :: Bestim 2312 real (kind=dp) , dimension(ncar,ncar) ,intent(out) :: rhoiestim 2313 real (kind=dp) , dimension(ncar,np) ,intent(out) :: sigsquareEstime 2314 type(INCIDENCE_TYPE) , dimension(ncar) , intent(inout) :: listDesc 2315 logical ,intent(in) :: performPrecision,tConf 2316 real (kind=dp),dimension(ntnivmaxtotal,ntnivmaxtotal,ncar) ,intent(out) :: tempForConfusion 2317 logical ,intent(in) :: invlrt 2318 2319 integer :: i,hypothesis 2320 character(len=LEN_L) :: dx 2321 2322 allocate (current_chr(workstruct%nqtl)) 2323 current_chr=curPos%listChr 2324 2325 call model_optim_multi_LU(xinc,listDesc,workstruct,sigsquareEstime,rhoiestim,Bestim,& 2326 performPrecision,likelihood_ncar_hn_family_LU,likelihood_ncar_hn_LU_family_withcd,tConf,tempForConfusion) 2327 2328 ! call model_optim_multi_LU(xinc,listDesc,workstruct,sigsquareEstime,rhoiestim,Bestim,& 2329 ! performPrecision,likelihood_ncar_hn_LU,likelihood_ncar_hn_LU,tConf,tempForConfusion) 2330 ! print *,"F:",workstruct%fmax 2331 2332 2333 if ((listDesc(1)%fmax == INIFINY_REAL_VALUE) ) then 2334 dx="" 2335 do i=1,workstruct%nqtl-1 2336 dx=trim(dx)//trim(str(absi(curPos%listChr(i),curPos%listN(i))))//"," 2337 end do 2338 2339 dx=trim(dx)//str(absi(curPos%listChr(workstruct%nqtl),curPos%listN(workstruct%nqtl))) 2340 call log_mess("dx ["//trim(dx)// "]. Can not optimize likelihood....The start point is reinitializing",WARNING_DEF) 2341 call init_startpoint(workStruct,listDesc(1)) 2342 curPos%lrt=0.d0 2343 listDesc(1)%fmax=0.d0 2344 return 2345 end if 2346 2347 !compute LRT 2348 if (invLrt) then 2349 do hypothesis=1,workstruct%nqtl 2350 curPos%lrt(hypothesis)=-2.d0*(workstruct%fnqtl(hypothesis)-listDesc(1)%fmax) 2351 end do 2352 else 2353 do hypothesis=1,workstruct%nqtl 2354 curPos%lrt(hypothesis)=-2.d0*(listDesc(1)%fmax-workstruct%fnqtl(hypothesis)) 2355 end do 2356 end if 2357 2358 deallocate (current_chr) 2359 2360 end subroutine model_optim_multi_hn_LU
model_optim_multi_LU
[ Top ] [ m_qtlmap_incidence_multi ] [ Subroutines ]
NAME
model_optim_multi_LU
DESCRIPTION
NOTE
SOURCE
2371 subroutine model_optim_multi_LU(xinc,listDesc,workstruct,sigsquareEstime,rhoiestim,Bestim,& 2372 performPrecision,FUNCT_PART,FUNCT_PART_CENSURE,tConf,tempForConfusion) 2373 real (kind=dp) , dimension(ncar,nd,ntnivmaxtotal), intent(in) :: xinc 2374 type(INCIDENCE_GEN_STRUCT) ,intent(inout) :: workstruct 2375 real (kind=dp) , dimension(ncar,np) ,intent(out) :: sigsquareEstime 2376 real (kind=dp) , dimension(ncar,ntnivmaxtotal) ,intent(out) :: Bestim 2377 real (kind=dp) , dimension(ncar,ncar) ,intent(out) :: rhoiestim 2378 type(INCIDENCE_TYPE) , dimension(ncar) , intent(inout),target :: listDesc 2379 logical ,intent(in) :: performPrecision,tConf 2380 real (kind=dp),dimension(ntnivmaxtotal,ntnivmaxtotal,ncar) ,intent(out) :: tempForConfusion 2381 external :: FUNCT_PART,FUNCT_PART_CENSURE 2382 2383 integer :: j,i,ip,ifail,ibound,npar 2384 ! real(kind=dp) ,dimension(np+ntnivmax) :: par 2385 real (kind=dp) , dimension(ntnivmaxtotal,ntnivmaxtotal) :: XX 2386 real (kind=dp) , dimension(ntnivmaxtotal,ntnivmaxtotal) :: triang 2387 integer :: iuser(1),ix,kd1,kd2,jm,k,jjm,ic,s,jc 2388 real (kind=dp) ,dimension(1) :: user 2389 logical ,dimension(ncar,ntnivmaxtotal) :: lastvecsol 2390 real(kind=dp) :: vci_t(ncar,ncar,np),determ_t(np),fp(np),fm(nm),r 2391 real(kind=dp),dimension(ncar,ncar) :: mat_L,V 2392 real(kind=dp),dimension(np,ncar) :: mat_H 2393 real(kind=dp) :: sigtemp 2394 2395 !Filter for the mimimization - dim (np,nm,nd) 2396 logical , dimension(np,nm,ncar*np+ncar*(ncar-1)/2+ncar*ntnivmaxtotal) :: filter_inc 2397 2398 iuser=0 2399 user=0 2400 my_listDesc=>listDesc 2401 dataset => listDesc(1)%dataset 2402 allocate (my_xincreduitmul(nd,ntnivmaxtotal,ncar)) 2403 2404 my_xincreduitmul=0.d0 2405 do ic=1,ncar 2406 lastvecsol(ic,:my_listDesc(ic)%ntniv)=my_listDesc(ic)%vecsol(:my_listDesc(ic)%ntniv) 2407 ! create X'.X matrix from incidence matrix 2408 call model_XT_X(xinc(ic,:,:),my_listDesc(ic),XX) 2409 ! Check all parameters to remove from the estimation 2410 call estim_cholesky(XX,my_listDesc(ic),ntnivmaxtotal,triang) 2411 ! compute the precision of each parameter 2412 if (performPrecision) call get_precision(XX,tempForConfusion(:,:,ic),my_listDesc(ic)) 2413 call set_corrxinc(xinc(ic,:,:),my_listDesc(ic),my_xincreduitmul(:,:,ic)) 2414 !optimisation : on sauvegarde les index des elements non nul de la matrice reduite 2415 call fill_nonull_elements(my_listDesc(ic),size(my_xincreduitmul,1),size(my_xincreduitmul,2),my_xincreduitmul(:,:,ic)) 2416 !call debug_write_incidence(xinc(ic,:,:),my_listDesc(ic)) 2417 end do 2418 2419 ! Optimisation de la vraisemblance a la position dx 2420 ifail=1 2421 ibound=0 2422 npar=ncar*np+ncar*(ncar-1)/2 2423 2424 j=ncar*np+ncar*(ncar-1)/2 2425 do ic=1,ncar 2426 npar=npar+my_listDesc(ic)%nbnivest 2427 2428 if (count(lastvecsol(ic,:))>0) then ! autre que initialisation 2429 do i=1,my_listDesc(ic)%ntniv 2430 if(my_listDesc(ic)%vecsol(i))then 2431 j=j+1 2432 if(.not. lastvecsol(ic,i)) my_listDesc(ic)%par(j)=0.d0 2433 end if 2434 end do 2435 end if 2436 end do 2437 2438 filter_inc=.true. 2439 2440 ! allocate (filter_inc(np,npar)) 2441 s=np*ncar+ncar*(ncar-1)/2 2442 ! filter_inc=.true. 2443 ! filter_inc(:,:,1:np*ncar)=.false. 2444 ! filter_inc(:,:,np*ncar+1:s)=.true. !rho 2445 filter_inc(:,:,1:(np-1)*ncar)=.false. 2446 do ic=1,ncar 2447 do ip=1,np 2448 if ( ip>1) then 2449 filter_inc(ip,:,(ip-2)*ncar+1:(ip-1)*ncar)=.true. !rho 2450 endif 2451 filter_inc(ip,nmp(ip)+1:nmp(ip+1),(ic-1)*np+ip)=.true. 2452 jjm=0 2453 do jm=nmp(ip)+1,nmp(ip+1) 2454 jjm=jjm+1 2455 if (dataset%lSires(ip)%full_sib(jjm)%lastKd>0) then 2456 kd1=dataset%lSires(ip)%full_sib(jjm)%firstkd 2457 kd2=dataset%lSires(ip)%full_sib(jjm)%lastkd 2458 filter_inc(ip,jm,s+1:s+my_listDesc(ic)%nbnivest)=any(my_xincreduitmul(kd1:kd2,:my_listDesc(ic)%nbnivest,ic)/=0.d0,dim=1) 2459 else 2460 ! print *,'pas de desc....' 2461 filter_inc(ip,jm,s+1:s+my_listDesc(ic)%nbnivest)=.false. 2462 end if 2463 end do 2464 end do 2465 s=s+my_listDesc(ic)%nbnivest 2466 end do 2467 2468 ! print *,"PAR*** H" 2469 ! print *,my_listDesc(1)%par(1:(np-1)*ncar) 2470 ! print *,"PAR*** L" 2471 ! print *,my_listDesc(1)%par((np-1)*ncar+1:(np-1)*ncar+ncar*(ncar+1)/2) 2472 ! print *,"PAR*** B" 2473 ! print *,my_listDesc(1)%par((np-1)*ncar+ncar*(ncar+1)/2+1:(np-1)*ncar+ncar*(ncar+1)/2+my_listDesc(1)%ntniv) 2474 2475 if ( datasetUser%na>0 ) then 2476 !prise en compte des donnees censures 2477 call minimizing_funct_family(npar,ibound,FUNCT_PART_CENSURE,filter_inc,fm,fp,& 2478 my_listDesc(1)%borni,my_listDesc(1)%borns,my_listDesc(1)%par,listDesc(1)%fmax,iuser,user,ifail) 2479 ! call minimizing_funct(npar,ibound,FUNCT_PART_CENSURE,& 2480 ! my_listDesc(1)%borni,my_listDesc(1)%borns,my_listDesc(1)%par,listDesc(1)%fmax,iuser,user,ifail) 2481 else 2482 !pas de donnee censure=> optimisation pour le calcul de l'inverse de matrice de covariance 2483 call minimizing_funct_family(npar,ibound,FUNCT_PART,filter_inc,fm,fp,& 2484 my_listDesc(1)%borni,my_listDesc(1)%borns,my_listDesc(1)%par,listDesc(1)%fmax,iuser,user,ifail) 2485 ! call minimizing_funct(npar,ibound,FUNCT_PART,& 2486 ! my_listDesc(1)%borni,my_listDesc(1)%borns,my_listDesc(1)%par,listDesc(1)%fmax,iuser,user,ifail) 2487 end if 2488 2489 Bestim=0.d0 2490 2491 2492 ! on fixe pour ip=1 H 2493 mat_H(1,:)=1.d0 2494 2495 do ip=2,np 2496 do jc=1,ncar 2497 mat_H(ip,jc)=my_listDesc(1)%par((ip-2)*ncar+jc) 2498 end do 2499 end do 2500 2501 mat_L=0.d0 2502 k=0 2503 do ic=1,ncar 2504 do jc=1,ic 2505 k=k+1 2506 mat_L(ic,jc) = my_listDesc(1)%par(ncar*(np-1)+k) 2507 end do 2508 end do 2509 2510 V = matmul(mat_L,transpose(mat_L)) 2511 2512 s=np*ncar+ncar*(ncar-1)/2 2513 !getting standart deviation 2514 do ic=1,ncar 2515 do ip=1,np 2516 sigsquareEstime(ic,ip)=V(ic,ic)*mat_H(ip,ic)*mat_H(ip,ic) 2517 end do 2518 !The solution 2519 Bestim(ic,:my_listDesc(ic)%nbnivest)=my_listDesc(1)%par(s+1:s+my_listDesc(ic)%nbnivest) 2520 s=s+my_listDesc(ic)%nbnivest 2521 end do 2522 2523 rhoiestim=0.d0 2524 2525 do ic=1,ncar 2526 do jc=1,ic 2527 rhoiestim(ic,jc)=V(ic,jc) / ( dsqrt( V(ic,ic)*V(jc,jc)) ) 2528 rhoiestim(jc,ic)=rhoiestim(ic,jc) 2529 end do 2530 end do 2531 2532 deallocate (my_xincreduitmul) 2533 2534 end subroutine model_optim_multi_LU
opti_0qtl_multi
[ Top ] [ m_qtlmap_incidence_multi ] [ Subroutines ]
NAME
opti_0qtl_multi
DESCRIPTION
NOTE
SOURCE
1111 subroutine opti_0qtl_multi( incsol , & ! 1112 workstruct, & 1113 maxqtl,FUNCT_MODEL,resol_lu) ! maximum qtl to find 1114 1115 integer ,intent(in) :: maxqtl 1116 type(INCIDENCE_GEN_STRUCT) ,intent(inout) :: workstruct 1117 type(TYPE_INCIDENCE_SOLUTION) ,intent(out) :: incsol 1118 logical ,intent(in) :: resol_lu 1119 external :: FUNCT_MODEL ! function with a model specific 1120 1121 integer :: nkd,ip,i,listnteff(1) 1122 1123 real (kind=dp) ,dimension(:,:,:),pointer :: xinc 1124 real (kind=dp) , dimension(:,:) ,pointer :: Bestim,rhoiestim 1125 type(INCIDENCE_TYPE) ,dimension(ncar) :: listDesc 1126 type(POSITION_LRT_INCIDENCE) :: curPos 1127 real (kind=dp) :: sigsquareEstime(ncar,np),f 1128 integer :: k,ic,jc 1129 1130 !initialisation of incidence matrix 1131 call init_incidence_multi(0,xinc,listDesc,workstruct%type_incidence,resol_lu) 1132 dataset=>listDesc(1)%dataset 1133 allocate (Bestim(ncar,ntnivmaxtotal)) 1134 allocate (rhoiestim(ncar,ncar)) 1135 Bestim=0.d0 1136 rhoiestim=0.d0 1137 sigsquareEstime=0.d0 1138 workstruct%nqtl=0 1139 Bestim=0.d0 1140 call add_general_mean_multi(xinc,listDesc) 1141 workstruct%eff_general_mean=1 1142 !add polygenic effect 1143 call add_polygenic_effect_multi(xinc,listDesc) 1144 !add fixed effect and covariate 1145 call add_effcov_multi(xinc,listDesc,0,0) 1146 1147 ! call model_optim_multi_h0(xinc,listDesc,workstruct,sigsquareEstime,rhoiestim,Bestim,.true.) 1148 call FUNCT_MODEL(xinc,listDesc,workstruct,sigsquareEstime,rhoiestim,Bestim,.true.) 1149 call set_solution_multi(0,workstruct,sigsquareEstime,rhoiestim,Bestim,listDesc,incsol,1,listnteff) 1150 1151 if (size(workstruct%sigsquare,1)>=1) then 1152 do ip=1,np 1153 do ic=1,ncar 1154 workstruct%sigsquare(1,ip,ic)=sigsquareEstime(ic,ip) 1155 end do 1156 end do 1157 end if 1158 1159 if (size(workstruct%fnqtl,1)>=1) then 1160 workstruct%fnqtl(1)=listDesc(1)%fmax 1161 end if 1162 1163 if (associated(workstruct%rhoi)) then 1164 workstruct%rhoi(1,:,:) = rhoiestim 1165 end if 1166 1167 !call model analysis 1168 !call FUNCT_MODEL(xinc,incidenceDesc,workstruct,Bestim,.true.) 1169 ! call set_solution(ic,0,workstruct,Bestim,incidenceDesc,incsol,1,listnteff) 1170 1171 call end_incidence_multi(listDesc) 1172 deallocate (xinc) 1173 deallocate (Bestim) 1174 deallocate (rhoiestim) 1175 end subroutine opti_0qtl_multi
opti_nqtl_multi
[ Top ] [ m_qtlmap_incidence_multi ] [ Subroutines ]
NAME
opti_nqtl_multi
DESCRIPTION
NOTE
SOURCE
1241 subroutine opti_nqtl_multi( nqtl, & ! under hypothesis NQTL 1242 workStruct, & ! variance estimated under NQTL-1 hypothesis 1243 incsol , & ! incidence solution (estimation of each effect) 1244 lrtsol , & ! maximum lrt finding at a position 1245 FUNCT_MODEL, & ! function with a model specific 1246 maxqtl,resol_lu) 1247 1248 integer ,intent(in) :: nqtl,maxqtl 1249 type(INCIDENCE_GEN_STRUCT) ,intent(inout) :: workstruct 1250 type(TYPE_INCIDENCE_SOLUTION) ,intent(out) :: incsol 1251 type(TYPE_LRT_SOLUTION) ,intent(out) :: lrtsol 1252 logical ,intent(in) :: resol_lu 1253 external :: FUNCT_MODEL 1254 1255 1256 integer :: i 1257 1258 real (kind=dp) ,dimension(:,:,:),pointer :: xxx,xinc ! incidence matrix and temp for testing nuisance effect 1259 type(INCIDENCE_TYPE) ,dimension(ncar) :: listDesc ! description incidence matrix 1260 type(POSITION_LRT_INCIDENCE) :: curPos ! the position with additional info for the model 1261 integer :: n,j,chr,k,ic,jc,ip 1262 1263 ! solution of the estimation 1264 real (kind=dp) , dimension(:,:) ,pointer :: Bestim,rhoiestim 1265 real(kind=dp) :: sigsquareEstime(ncar,np) 1266 1267 if ( nqtl <= 0 ) then 1268 call stop_application("Error dev: can not call opti_nqtl_linear with nqtl:"//trim(str(nqtl))) 1269 end if 1270 1271 !Position Allocation 1272 call init_position(nqtl,curPos) 1273 1274 !initialisation of incidence matrix 1275 call init_incidence_multi(nqtl,xinc,listDesc,workstruct%type_incidence,resol_lu) 1276 1277 dataset=>listDesc(1)%dataset 1278 1279 allocate (Bestim(ncar,ntnivmaxtotal)) 1280 allocate (rhoiestim(ncar,ncar)) 1281 allocate (xxx(ncar,ntnivmaxtotal,ntnivmaxtotal)) 1282 1283 if ( .not. resol_lu ) call init_startpoint(workstruct,listDesc(1)) 1284 1285 Bestim=0.d0 1286 rhoiestim=0.d0 1287 sigsquareEstime=0.d0 1288 1289 do i=1,nqtl 1290 n=0 1291 j=i 1292 do chr=1,nchr 1293 n=n+get_ilong(chr) 1294 if (i<=n) exit 1295 j=1 1296 end do 1297 1298 if (chr>nchr) call stop_application("Not enough sampling point to detect QTLs") 1299 1300 curPos%listChr(i)=1 ! attention si pas assez d echantillonage on doit passer au chromosome suivant.... 1301 curPos%listN(i)=i 1302 end do 1303 1304 ! call build_incidence_matrix(workstruct,curPos,xinc,incidenceDesc,0,0,0) 1305 Bestim=0.d0 1306 call add_general_mean_multi(xinc,listDesc) 1307 workstruct%eff_general_mean=1 1308 1309 workstruct%listnteff(1)=2 1310 do i=1,workstruct%nqtl 1311 !add qtl effect/interaction at position n to estim 1312 call add_qtleffect_multi(workstruct%listnteff(i),i,curPos%listChr(i),curPos%listN(i),& 1313 xinc,listDesc,0,workstruct%linear) 1314 if (i<workstruct%nqtl) workstruct%listnteff(i+1)=listDesc(1)%nteff+1 1315 end do 1316 !add polygenic effect 1317 call add_polygenic_effect_multi(xinc,listDesc) 1318 !add fixed effect and covariate 1319 call add_effcov_multi(xinc,listDesc,0,0) 1320 1321 !initalizing lrtsolution 1322 call new(nqtl,lrtsol) 1323 1324 lrtsol%lrtmax=-INIFINY_REAL_VALUE 1325 1326 ! call gen_loop_opti_nqtl(1,nqtl,curPos,workstruct,xinc,& 1327 ! incidenceDesc,Bestim,lrtsol,.true.,FUNCT_MODEL) 1328 !call gen_loop_opti_nqtl_multi(1,nqtl,curPos,workstruct,& 1329 ! xinc,listDesc,sigsquareEstime,rhoiestim,Bestim,lrtsol,.false.,FUNCT_MODEL) 1330 1331 1332 call gen_opti_nqtl_multi(nqtl,curPos,workstruct,xinc,listDesc,sigsquareEstime,rhoiestim,Bestim,lrtsol,FUNCT_MODEL) 1333 1334 ! Compute precision at the maximum LRT in position posx 1335 do i=1,nqtl 1336 call change_qtleffect_multi(workstruct%listnteff(i),i,lrtsol%chrmax(i-1),lrtsol%nxmax(i-1),& 1337 xinc,listDesc,0) 1338 curPos%listChr(i)=lrtsol%chrmax(i-1) 1339 curPos%listN(i)=lrtsol%nxmax(i-1) 1340 end do 1341 1342 ! call debug_write_incidence(xinc,incidenceDesc) 1343 call FUNCT_MODEL(xinc,listDesc,curPos,workstruct,sigsquareEstime,rhoiestim,Bestim,& 1344 .true.,workstruct%performConfusion,xxx,.false.) 1345 1346 ! print *,"****************************************************************" 1347 if ( workstruct%performConfusion) then 1348 call confusion_multi_type1(listDesc,xxx,workstruct) 1349 end if 1350 if ( workstruct%performTestNuis ) then 1351 ! call test_nuisances(ic,incidenceDesc%nbnivest,lrtsol,curPos,workstruct,FUNCT_MODEL) 1352 end if 1353 call set_solution_multi(nqtl,workstruct,sigsquareEstime,rhoiestim,Bestim,listDesc,incsol,nqtl,workstruct%listnteff) 1354 1355 if (size(workstruct%sigsquare,1)>nqtl) then 1356 do ip=1,np 1357 do ic=1,ncar 1358 workstruct%sigsquare(nqtl+1,ip,ic)=sigsquareEstime(ic,ip) 1359 end do 1360 end do 1361 end if 1362 1363 if (size(workstruct%fnqtl,1)>nqtl) then 1364 workstruct%fnqtl(nqtl+1)=listDesc(1)%fmax 1365 end if 1366 1367 if (size(workstruct%rhoi,1)>nqtl) then 1368 workstruct%rhoi(nqtl+1,:,:) = rhoiestim 1369 end if 1370 1371 call end_incidence_multi(listDesc) 1372 1373 deallocate (xxx) 1374 deallocate (xinc) 1375 deallocate (Bestim) 1376 deallocate (rhoiestim) 1377 call end_position(curPos) 1378 end subroutine opti_nqtl_multi
pdd_at_position_multi
[ Top ] [ m_qtlmap_incidence_multi ] [ Subroutines ]
NAME
pdd_at_position_multi
DESCRIPTION
NOTE
SOURCE
863 subroutine pdd_at_position_multi(dataset,iq,chr,n) 864 type(DATASET_TYPE) , intent(inout) :: dataset 865 integer ,intent(in) :: iq,chr,n 866 867 integer :: ip,jm,kd,ig,kkd,igg,ifem,kkkd,kds 868 real (kind=dp) :: ppt,pmt 869 870 do ip=1,np 871 dataset%lSires(ip)%ppt(iq,:)=0.d0 872 ifem=0 873 kkkd=0 874 do jm=nmp(ip)+1,nmp(ip+1) 875 ifem=ifem+1 876 if (dataset%lSires(ip)%full_sib(ifem)%firstkd>=0) then 877 dataset%lSires(ip)%full_sib(ifem)%ppt(iq,:,:)=0.d0 878 dataset%lSires(ip)%full_sib(ifem)%pmt(iq,:,:)=0.d0 879 igg=0 880 kds=kkkd 881 ! print *,"ip,jm,firstkd:",ip,jm,lSires(ip)%half_sib%firstKd 882 do ig=ngenom(chr,jm)+1,ngenom(chr,jm+1) 883 igg=igg+1 884 kkd=0 885 kkkd=kds 886 do kd=ngend(chr,ig)+1,ngend(chr,ig+1) 887 if(count(presentc(:,ndesc(chr,kd)))==ncar) then 888 kkd=kkd+1 889 kkkd=kkkd+1 890 if( estime_multi(jm) )then 891 dataset%lSires(ip)%full_sib(ifem)%ppt(iq,igg,kkd)=& 892 -pdd(chr,kd,1,n)-pdd(chr,kd,2,n)+pdd(chr,kd,3,n)+pdd(chr,kd,4,n) 893 dataset%lSires(ip)%full_sib(ifem)%pmt(iq,igg,kkd)=& 894 -pdd(chr,kd,1,n)+pdd(chr,kd,2,n)-pdd(chr,kd,3,n)+pdd(chr,kd,4,n) 895 else 896 dataset%lSires(ip)%ppt(iq,dataset%lSires(ip)%half_sib%firstKd+kkkd-1)=& 897 -pdd(chr,kd,1,n)+pdd(chr,kd,3,n) 898 end if 899 end if 900 end do !! kd 901 end do !! ig 902 end if 903 end do ! jm 904 end do ! ip 905 end subroutine pdd_at_position_multi