m_qtlmap_incidence
[ Top ] [ ANALYSE ] [ Modules ]
NAME
m_qtlmap_incidence
DESCRIPTION
NOTES
SEE ALSO
DATASET_TYPE
[ Top ] [ m_qtlmap_incidence ] [ Types ]
NAME
DATASET_TYPE
DESCRIPTION
description of the user dataset
SOURCE
183 type DATASET_TYPE 184 !Index trait 185 integer :: ncar = 0 186 !Struct of Families (missing data are removes) with pdd for full and half sib families 187 type(SireFamily_type) , dimension(:) , pointer :: lSires => NULL() 188 !Number of animals referenced in the contingence matrix (missing data removed) 189 integer :: nkd = 0 190 !Maximum number recorded in a family full sib 191 integer :: nkd_max_by_fam = 0 192 !Performance without missing data : size : nkd 193 real(kind=dp) , dimension(:,:), pointer :: Y => NULL() 194 integer , dimension(:,:), pointer :: YDISCRETORD => NULL() 195 !Censored data without missing data : size nkd 196 real(kind=dp) , dimension(:,:), pointer :: CD => NULL() 197 !Covariance of traits (H0) dim : ncar,ncar 198 real(kind=dp) , dimension(:,:), pointer :: cov => NULL() 199 !Residual correlation (H0) dim : ncar,ncar 200 real(kind=dp) , dimension(:,:), pointer :: rhoi => NULL() 201 !Standart deviation (H0), dim:ncar 202 real(kind=dp) , dimension(:) , pointer :: sig => NULL() 203 !Mean (H0), dim:ncar 204 real(kind=dp) , dimension(:) , pointer :: xmu => NULL() 205 206 end type DATASET_TYPE
family_type
[ Top ] [ m_qtlmap_incidence ] [ Types ]
NAME
family_type
DESCRIPTION
Basic type with the proba of transmission of a family
SOURCE
149 type family_type 150 integer :: firstKd = -1 ! correspond a l index du premier descendant (equivalent ligne mat incidence) 151 integer :: lastKd = -2 ! correspond a l index du dernier descedant (equivalent ligne mat incidence) 152 real(kind=dp) , dimension(:,:,:),pointer :: ppt => null()!proba de trsanmition recu par le pere 153 real(kind=dp) , dimension(:,:,:),pointer :: pmt => null()! proba de transmision recu par la mere si estim 154 ! integer ,dimension(:) ,pointer :: ntniv_dam !indice dans la matrice d incidence correspondant a l effet pere et le qtl 155 end type family_type
INCIDENCE_GEN_STRUCT
[ Top ] [ m_qtlmap_incidence ] [ Types ]
NAME
INCIDENCE_GEN_STRUCT
DESCRIPTION
structure pour developper des methodes generiques
ce type est un agrega des structures (non commune au mode lineaire et optim) necessaires pour la solution de la vraissemblance a un point n
SOURCE
281 type INCIDENCE_GEN_STRUCT 282 !!!!! initializing by the caller !!!!! 283 !************************************** 284 integer :: type_model = -1 ! linear homo,hetero,optim,... 285 integer :: type_incidence = -1 ! classic, traits as co, cox,ld 286 logical :: performConfusion = .false. 287 logical :: performTestNuis = .false. 288 logical :: linear = .false. 289 logical :: hdam = .false. 290 ! logical :: nopoly = .false. 291 integer :: nqtl 292 293 integer :: eff_general_mean = -1 294 !id nteff inside xinc 295 integer ,dimension(:), pointer :: listnteff => null() 296 !id nteff inside xinc 297 integer ,dimension(:), pointer :: listnteffhap => null() 298 !linear information 299 !------------------ 300 !0-->maxiqtl,np 301 real (kind=dp) ,dimension(:,:,:) ,pointer :: sigsquare => null() 302 !maxqtl,ncar,ncar 303 real (kind=dp) ,dimension(:,:,:) ,pointer :: rhoi => null() 304 305 !F0,F1,etc... dim:nqtl 306 real (kind=dp) ,dimension(:) ,pointer :: fnqtl => null() 307 !ADDITIONAL OUTPUTS INFORMATIONS 308 !------------------------------- 309 ! alert confusion (fill if performConfusion is true ) 310 type(CORR_ALERT_TYPE) ,dimension(:),pointer :: alertQtl => null() 311 ! size of alertQtl 312 integer :: nalert = 0 313 ! Maximum correlation finding between qtl and other effect 314 real (kind=dp) :: corrmax 315 316 ! list of nuisance effect and test (fill if performTestNuis is true) 317 type(TEST_NUISANCES_TYPE) ,dimension(:),pointer :: listtestnuis => null() 318 ! size of listtestnuis 319 integer :: ntest = 0 320 ! index to each trait to init as a covariate 321 integer :: startTraitAsCov,endTraitAsCov 322 end type INCIDENCE_GEN_STRUCT
INCIDENCE_TYPE
[ Top ] [ m_qtlmap_incidence ] [ Types ]
NAME
INCIDENCE_TYPE
DESCRIPTION
description of an effect enumered in the INCIDENCE_TYPE
SOURCE
216 type INCIDENCE_TYPE 217 !Index trait 218 integer :: ic = 0 219 !The associated dataset 220 type(DATASET_TYPE) ,pointer :: dataset => NULL() 221 !number of qtl to estim 222 integer :: nqtl = 0 223 ! Number of level 224 integer :: ntniv = 0 225 ! Number of effect 226 integer :: nteff = 0 227 ! tab of effect to print in the eqtl format 228 logical ,dimension(20) :: eqtl_print = .true. 229 integer :: ntnivmax=0 230 !Maximum of effect defined in the model : General mean, polygenic sire,polygenic dam, qtl sire,.... 231 integer :: nteffmax=0 232 ! Number of level which can be estimate 233 integer :: nbnivest = 0 234 ! Number of qtl*interaction level by sire or dam 235 integer :: ntlev = 0 236 ! Number of level for fixed effect 237 integer :: nbniv = 0 238 !description of each effect 239 type(DESC_EFFECT_TYPE),dimension(:),pointer :: desc => NULL() 240 !estimability of each level ntniv 241 logical , dimension(:) , pointer :: vecsol => NULL() 242 !precision of each level estim 243 real (kind=dp) ,dimension(:),pointer :: precis => NULL() 244 !dim:nqtl : referenced qtl sire ntniv in the contingence matrix 245 integer ,dimension(:),pointer :: ntniv_qtlsires => NULL() 246 !dim:nqtl : referenced qtl dams ntniv in the contingence matrix 247 integer ,dimension(:),pointer :: ntniv_qtldams => NULL() 248 !dim:ntnivmax correspondance si estimable de la colonne dans la matrice xinc reduite 249 integer ,dimension(:),pointer :: corr_niv_nivb => NULL() 250 !optim struc 251 real (kind=dp) , dimension(:) ,pointer :: borni => null() 252 real (kind=dp) , dimension(:) ,pointer :: borns => null()! les bornes de l optimisation 253 real (kind=dp) , dimension(:) ,pointer :: par => null() 254 real(kind=dp) :: fmax = 0.d0 255 256 !optimisation qui fait office de corniv dans le modlin originale... 257 ! on stocke pour chque kd, les indices des elements des elements non nul.... 258 integer , dimension(:,:) ,pointer :: nonull => null() 259 integer , dimension(:) ,pointer :: nbnonull => null()! nombre d element non nul 260 261 end type INCIDENCE_TYPE
POSITION_LRT_INCIDENCE
[ Top ] [ m_qtlmap_incidence ] [ Types ]
NAME
POSITION_LRT_INCIDENCE
DESCRIPTION
structure pour developper des methodes generiques
ce type est un agrega des structures (non commune au mode lineaire et optim) necessaires pour la solution de la vraissemblance a un point n
SOURCE
333 type POSITION_LRT_INCIDENCE 334 !Information bind to the current position 335 !---------------------------------------- 336 !nqtl 337 integer ,dimension(:) ,pointer :: listChr => null() 338 integer ,dimension(:) ,pointer :: listN => null() 339 !output info fill by the model 340 real(kind=dp) ,dimension(:) ,pointer :: lrt => null() 341 342 end type POSITION_LRT_INCIDENCE
SireFamily_type
[ Top ] [ m_qtlmap_incidence ] [ Types ]
NAME
SireFamily_type
DESCRIPTION
Complex type with the proba of transmission of a sire family
SOURCE
165 type SireFamily_type 166 type(family_type) :: half_sib ! famille complete de demi-frere (ligne ref incidence) 167 type(family_type) , dimension(:),pointer :: full_sib => null()! famille de plein frere (ligne ref incidence) 168 real(kind=dp) , dimension(:),pointer :: sig0 => null()!variance du caractere dans la famille de pere 169 real(kind=dp) , dimension(:),pointer :: xmu0 => null()!moyenne du caractere dans la famille de pere 170 real(kind=dp) , dimension(:,:),pointer :: ppt => null()!proba de transmition recu par le pere 171 !integer ,dimension(:) ,pointer :: ntniv_sire !indice dans la matrice d incidence correspondant a l effet pere et le qtl 172 end type SireFamily_type
INCIDENCE_INTERACTION
[ Top ] [ m_qtlmap_incidence ] [ Constants ]
NAME
INCIDENCE_INTERACTION
DESCRIPTION
Constant id for construction type of the contingence matrix: <DEV>
NOTES
build_incidence_matrix
INCIDENCE_TRAIT_COV
[ Top ] [ m_qtlmap_incidence ] [ Constants ]
NAME
INCIDENCE_TRAIT_COV
DESCRIPTION
Constant id for construction type of the contingence matrix: * General mean * QTL effect * Polygenic effect * Others continue trait as a covariate * Covariate * Fixed Effect
NOTES
build_incidence_matrix
INCIDENCE_TYPE_CLASSIC
[ Top ] [ m_qtlmap_incidence ] [ Constants ]
NAME
INCIDENCE_TYPE_CLASSIC
DESCRIPTION
Constant id for construction type of the contingence matrix: * General mean * Polygenic effect * QTL effect * Covariate * Fixed Effect
NOTES
build_incidence_matrix
INCIDENCE_TYPE_LD
[ Top ] [ m_qtlmap_incidence ] [ Constants ]
NAME
INCIDENCE_TYPE_LD
DESCRIPTION
Constant id for construction type of the contingence matrix: * General mean * Haplotype effect * Polygenic effect * Others continue trait as a covariate * Covariate * Fixed Effect
NOTES
build_incidence_matrix
INCIDENCE_TYPE_LDLA
[ Top ] [ m_qtlmap_incidence ] [ Constants ]
NAME
INCIDENCE_TYPE_LDLA
DESCRIPTION
Constant id for construction type of the contingence matrix: * General mean * Haplotype effect * QTL effect * Polygenic effect * Others continue trait as a covariate * Covariate * Fixed Effect
NOTES
build_incidence_matrix
LINEAR_HETEROSCEDASTIC
[ Top ] [ m_qtlmap_incidence ] [ Constants ]
NAME
LINEAR_HETEROSCEDASTIC
DESCRIPTION
Constant id for linear heteroscedastic analysis
NOTES
LINEAR_HOMOSCEDASTIC
[ Top ] [ m_qtlmap_incidence ] [ Constants ]
NAME
LINEAR_HOMOSCEDASTIC
DESCRIPTION
Constant id for linear homoscedastic analysis
NOTES
MAX_KD_PRINT
[ Top ] [ m_qtlmap_incidence ] [ Constants ]
NAME
MAX_KD_PRINT
DESCRIPTION
Constant : maximum of progeny to print for the debug mode
NOTES
debug_write_incidence
OPTIM_HETEROSCEDASTIC
[ Top ] [ m_qtlmap_incidence ] [ Constants ]
NAME
OPTIM_HETEROSCEDASTIC
DESCRIPTION
Constant id for optimized heteroscedastic analysis
NOTES
likelyprob
[ Top ] [ m_qtlmap_incidence ] [ Variables ]
NAME
likelyprob
DESCRIPTION
give the index the likely genotype dam DIMENSION * nchr * np * nm
NOTES
init_analyse_linear
add_effcov
[ Top ] [ m_qtlmap_incidence ] [ Subroutines ]
NAME
add_effcov
DESCRIPTION
Add the nuisances effects inside a contingence matrix with a description * An effect is adding for each covariate and fixed effect
INPUTS
meff : fixed effect to not added in the contingence matrix (if meff <= 0, all fixed effects are adding) mcov : covariate effect to not added in the contingence matrix (if mcov <= 0, all covariates effects are adding)
OUTPUTS
xinc : The contingence matrix incidenceDesc : The description of contingence matrix
NOTES
SOURCE
1927 subroutine add_effcov(xinc,incidenceDesc,meff,mcov) 1928 type(INCIDENCE_TYPE) , intent(inout) :: incidenceDesc 1929 real (kind=dp) , dimension(nd,incidenceDesc%ntnivmax) , intent(inout) :: xinc 1930 integer , intent(in) :: meff ! indice on effet fixed list to remove (<=0 otherwise) 1931 integer , intent(in) :: mcov ! indice on covariate list to remove (<=0 otherwise) 1932 integer :: nkd,ntniv,ip,jm,kd,ntnifix 1933 integer :: nbef,nbco,nteff,jef,jcov,ieff 1934 integer :: nbniv,nivd,i,j,nbef2,nbco2 1935 1936 call log_mess("incidence : ** Add eff and cov ** ",DEBUG_DEF) 1937 1938 nbef=listModelTrait(incidenceDesc%ic)%nbfe 1939 nbco=listModelTrait(incidenceDesc%ic)%nbco 1940 1941 nbef2=nbef 1942 if ( meff /= 0 ) nbef2=nbef-1 1943 nbco2=nbco 1944 if ( mcov /=0 ) nbco2=nbco-1 1945 1946 1947 if (nbef2+nbco2 <=0)return 1948 1949 nbniv=0 1950 ieff=incidenceDesc%nteff 1951 1952 if ( nbef2 > 0 ) then 1953 ieff=ieff+1 1954 incidenceDesc%eqtl_print(ieff)=.false. 1955 do i=1,nbef 1956 if ( i == meff ) cycle 1957 nbniv=nlev(listModelTrait(incidenceDesc%ic)%indexFixedEffect(i))+nbniv 1958 end do 1959 incidenceDesc%desc(ieff)%isVar = .false. 1960 incidenceDesc%desc(ieff)%name="Fixed effects" 1961 incidenceDesc%desc(ieff)%start=incidenceDesc%ntniv+1 1962 incidenceDesc%desc(ieff)%end=incidenceDesc%ntniv+nbniv 1963 incidenceDesc%desc(ieff)%haveSubDesc=.true. 1964 1965 !init of xinc 1966 xinc(:,incidenceDesc%desc(ieff)%start:incidenceDesc%desc(ieff)%end)=0.d0 1967 1968 allocate ( incidenceDesc%desc(ieff)%listSubDesc(nbniv) ) 1969 1970 j=0 1971 do jef=1,nbef 1972 if ( jef == meff ) cycle 1973 do i=1,nlev(listModelTrait(incidenceDesc%ic)%indexFixedEffect(jef)) 1974 j=j+1 1975 incidenceDesc%desc(ieff)%listSubDesc(j)%name=trim(namefix(listModelTrait(incidenceDesc%ic)%indexFixedEffect(jef)))& 1976 //' level '//trim(str(i)) 1977 incidenceDesc%desc(ieff)%listSubDesc(j)%start=incidenceDesc%ntniv+j 1978 incidenceDesc%desc(ieff)%listSubDesc(j)%end=incidenceDesc%ntniv+j 1979 end do 1980 end do 1981 incidenceDesc%nteff = incidenceDesc%nteff+1 1982 end if 1983 1984 if ( nbco2 > 0 ) then 1985 ieff=ieff+1 1986 incidenceDesc%eqtl_print(ieff)=.false. 1987 incidenceDesc%desc(ieff)%isVar = .false. 1988 incidenceDesc%desc(ieff)%name="Covariates" 1989 incidenceDesc%desc(ieff)%start=incidenceDesc%ntniv+nbniv+1 1990 incidenceDesc%desc(ieff)%end=incidenceDesc%ntniv+nbniv+nbco2 1991 incidenceDesc%desc(ieff)%haveSubDesc=.true. 1992 1993 !init of xinc 1994 xinc(:,incidenceDesc%desc(ieff)%start:incidenceDesc%desc(ieff)%end)=0.d0 1995 1996 allocate ( incidenceDesc%desc(ieff)%listSubDesc(nbco2) ) 1997 1998 j=0 1999 do jcov=1,nbco 2000 if ( jcov == mcov ) cycle 2001 j=j+1 2002 incidenceDesc%desc(ieff)%listSubDesc(j)%name=trim(namecov(listModelTrait(incidenceDesc%ic)%indexCovariate(jcov))) 2003 incidenceDesc%desc(ieff)%listSubDesc(j)%start=incidenceDesc%ntniv+nbniv+j 2004 incidenceDesc%desc(ieff)%listSubDesc(j)%end=incidenceDesc%ntniv+nbniv+j 2005 end do 2006 incidenceDesc%nteff = incidenceDesc%nteff+1 2007 end if 2008 2009 2010 nkd = 0 2011 nbniv=0 2012 do ip=1,np 2013 do jm=nmp(ip)+1,nmp(ip+1) 2014 do kd=ndm(jm)+1,ndm(jm+1) 2015 if ( .not. presentc(incidenceDesc%ic,kd) ) cycle 2016 ntniv=incidenceDesc%ntniv 2017 nkd=nkd+1 2018 ! 2019 ! autres effets fixes 2020 ! 2021 nbniv=0 2022 do jef=1,nbef 2023 if ( jef == meff ) cycle ! this effet fixed is not adding 2024 ! incidenceDesc%nivdir(kd,nteff+jef)=ntniv+nbniv+nivx(kd,modele(ic,3+jef)) 2025 nivd=ntniv+nbniv+nivx(kd,listModelTrait(incidenceDesc%ic)%indexFixedEffect(jef)) 2026 xinc(nkd,nivd)=1 2027 nbniv=nbniv+nlev(listModelTrait(incidenceDesc%ic)%indexFixedEffect(jef)) 2028 end do 2029 2030 ntnifix=ntniv+nbniv 2031 !covariate 2032 j=0 2033 do jcov=1,nbco 2034 if ( jcov == mcov ) cycle ! this covariate is not adding 2035 j=j+1 2036 ! incidenceDesc%covdir(kd,jcov)=covar(kd,modele(ic,nbt+jcov)) 2037 xinc(nkd,ntnifix+j)=covar(kd,listModelTrait(incidenceDesc%ic)%indexCovariate(jcov)) 2038 end do 2039 !ntniv=ntnifix+nbco 2040 !nbt=nbt+nbco 2041 end do 2042 end do 2043 end do 2044 2045 incidenceDesc%ntniv=incidenceDesc%ntniv+nbniv+nbco2 2046 end subroutine add_effcov
add_general_mean
[ Top ] [ m_qtlmap_incidence ] [ Subroutines ]
NAME
add_general_mean
DESCRIPTION
Add a general mean effect inside a contingence matrix with a description
OUTPUTS
xinc : The contingence matrix incidenceDesc : The description of contingence matrix
NOTES
SOURCE
1250 subroutine add_general_mean(xinc,incidenceDesc) 1251 type(INCIDENCE_TYPE) , intent(inout) :: incidenceDesc 1252 real (kind=dp) , dimension(nd,incidenceDesc%ntnivmax) , intent(inout) :: xinc 1253 integer :: ntniv 1254 call log_mess("incidence : ** Add General mean ** ",DEBUG_DEF) 1255 incidenceDesc%ntniv = incidenceDesc%ntniv+1 1256 1257 ntniv = incidenceDesc%ntniv 1258 incidenceDesc%nteff = incidenceDesc%nteff+1 1259 incidenceDesc%desc(incidenceDesc%nteff)%name="General Mean" 1260 incidenceDesc%desc(incidenceDesc%nteff)%start=incidenceDesc%nteff 1261 incidenceDesc%desc(incidenceDesc%nteff)%end=incidenceDesc%nteff 1262 incidenceDesc%desc(incidenceDesc%nteff)%isVar = .false. 1263 xinc(:,ntniv) = 1 1264 !incidenceDesc%nivdir(:,incidenceDesc%nteff) = 1 1265 end subroutine add_general_mean
add_haplotype_effect
[ Top ] [ m_qtlmap_incidence ] [ Subroutines ]
NAME
add_haplotype_effect
DESCRIPTION
Add a haplotype effect inside a contingence matrix with a description An effect is adding by haplotype finded in the population at the specific position The number of haplotype is growing depending hdam optional variable.
INPUTS
nteff : index of the effect
numhap : index of the haplotype effect
chr : chromosome localisation
n : the position
hdam :
start : using in a inialize context
OUTPUTS
xinc : The contingence matrix incidenceDesc : The description of contingence matrix
NOTES
SOURCE
1495 subroutine add_haplotype_effect(nteff,numhap,chr,n,xinc,incidenceDesc,hdam,start) 1496 integer ,intent(in) :: nteff,numhap,chr,n 1497 type(INCIDENCE_TYPE) ,intent(inout) :: incidenceDesc 1498 real (kind=dp) , dimension(nd,incidenceDesc%ntnivmax) ,intent(inout) :: xinc 1499 logical ,intent(in) :: hdam,start 1500 1501 !local 1502 integer :: ip,jm,kd,ntniv,nkd,j_gam,i_haplo,j,i_gam,geno,nteffstart 1503 real (kind=dp) ,dimension(:,:),allocatable :: pp_ldla,pm_ldla 1504 real (kind=dp) ,dimension(nd,NB_HAPLO_PRIOR) :: pb_haplo 1505 real (kind=dp) :: su 1506 1507 allocate (pp_ldla(nd,2)) 1508 allocate (pm_ldla(nd,2)) 1509 1510 ntniv = incidenceDesc%desc(nteff-1)%end 1511 nkd=0 1512 1513 call set_haplo_for_ldla(chr,absi(chr,n),n,.true.,hdam) 1514 1515 if (start) then 1516 nteffstart = incidenceDesc%nteff 1517 incidenceDesc%nteff = incidenceDesc%nteff+1 1518 incidenceDesc%ntniv = ntniv + 2**longhap 1519 else 1520 nteffstart = nteff - 1 1521 end if 1522 1523 incidenceDesc%desc(nteffstart+1)%isVar = .true. 1524 incidenceDesc%desc(nteffstart+1)%name="Haplotypes effects["//trim(str(numhap))//"]" 1525 incidenceDesc%desc(nteffstart+1)%start=ntniv+1 1526 incidenceDesc%desc(nteffstart+1)%end=ntniv+2**longhap ! nb_haplo_reduit 1527 1528 incidenceDesc%desc(nteffstart+1)%haveSubDesc=.true. 1529 if ( associated(incidenceDesc%desc(nteffstart+1)%listSubDesc) ) & 1530 deallocate (incidenceDesc%desc(nteffstart+1)%listSubDesc) 1531 1532 allocate (incidenceDesc%desc(nteffstart+1)%listSubDesc(nb_haplo_reduit)) 1533 1534 do j=1,nb_haplo_reduit 1535 incidenceDesc%desc(nteffstart+1)%listSubDesc(j)%name=trim(name_haplo_reduit(j))//"(race="//trim(race_haplo_reduit(j)) //& 1536 ", freq="//trim(str(pb_haplo_reduit(j)))//")" 1537 incidenceDesc%desc(nteffstart+1)%listSubDesc(j)%start=incidenceDesc%desc(nteffstart+1)%start+j-1 1538 incidenceDesc%desc(nteffstart+1)%listSubDesc(j)%end=incidenceDesc%desc(nteffstart+1)%start+j-1 1539 end do 1540 1541 xinc(:,ntniv+1:ntniv+2**longhap) = 0.d0 1542 1543 do ip=1,np 1544 do jm=nmp(ip)+1,nmp(ip+1) 1545 !if ( linear ) then 1546 call getprobld_lin(chr,n,incidenceDesc%ic,ip,jm,ndm(jm)+1,ndm(jm+1),pp_ldla,pm_ldla) 1547 ! else 1548 ! call getglobalprob(chr,n,incidenceDesc%ic,ip,jm,ndm(jm)+1,ndm(jm+1),pp_ldla,pm_ldla) 1549 !end if 1550 do kd=ndm(jm)+1,ndm(jm+1) 1551 if ( presentc(incidenceDesc%ic,kd) ) then 1552 nkd=nkd+1 1553 1554 do i_haplo=1,nb_haplo_reduit 1555 pb_haplo(kd,i_haplo)=0.d0 1556 end do 1557 1558 ! haplotype du pere 1559 !if(hsire)then 1560 do j=1,2 1561 su = dble(sum(pb_haplo_reduit(num_haplo_pere(ip,j,1:nb_gam_pere(ip,j))))) 1562 if ( su == 0 ) cycle 1563 do i_gam=1,nb_gam_pere(ip,j) 1564 pb_haplo(kd,num_haplo_pere(ip,j,i_gam))= pb_haplo(kd,num_haplo_pere(ip,j,i_gam)) & 1565 +pp_ldla(kd,j) * ( pb_haplo_reduit(num_haplo_pere(ip,j,i_gam)) / su ) 1566 end do !i_gam 1567 end do !j 1568 1569 1570 1571 !end if ! opt_model 1572 ! haplotype de la mere si phase connue 1573 if(hdam)then 1574 if(estime(incidenceDesc%ic,jm)) then 1575 geno = likelyprob(chr,ip,jm) 1576 do j=1,2 1577 su = dble(sum(pb_haplo_reduit(num_haplo_mere(geno,j,1:nb_gam_mere(geno,j))))) 1578 if ( su == 0 ) cycle 1579 do i_gam=1,nb_gam_mere(geno,j) 1580 pb_haplo(kd,num_haplo_mere(geno,j,i_gam)) = pb_haplo(kd,num_haplo_mere(geno,j,i_gam)) & 1581 +pm_ldla(kd,j) * ( pb_haplo_reduit(num_haplo_mere(geno,j,i_gam)) / su ) 1582 end do !i_gam 1583 end do !j 1584 ! haplotype de la mere si phase inconnue 1585 else 1586 ! do j_gam=1,nb_gam(kd) 1587 ! su = dble(sum(pb_haplo_reduit(num_haplo_desc(kd,j_gam,1:nb_gam_desc(kd,j_gam))))) 1588 ! if ( su == 0 ) cycle 1589 ! do i_gam=1,nb_gam_desc(kd,j_gam) 1590 ! pb_haplo(kd,num_haplo_desc(kd,j_gam,i_gam)) = pb_haplo(kd,num_haplo_desc(kd,j_gam,i_gam)) & 1591 ! +prob_gam(kd,j_gam) * (pb_haplo_reduit(num_haplo_desc(kd,j_gam,i_gam)) /su) 1592 ! end do !i_gam 1593 ! end do !j_gam 1594 1595 ! write(111,*)'n,kd,nb_gam_desc(kd)',n,kd,nb_gam_desc(kd) 1596 ! do i_gam=1,nb_gam_desc(kd) 1597 ! pb_haplo(kd,num_haplo_desc(kd,i_gam))=pb_haplo(kd,num_haplo_desc(kd,i_gam)) + & 1598 ! pb_haplo_desc(kd,i_gam) 1599 do i_gam=1,nb_haplo_reduit 1600 pb_haplo(kd,i_gam)=pb_haplo(kd,i_gam) + pb_haplo_desc(kd,i_gam) 1601 end do !i_gam 1602 1603 1604 end if !estime 1605 end if !opt_model (hdam) 1606 ! 1607 ! ligne d'incidence 1608 ! 1609 do i_haplo=1,nb_haplo_reduit 1610 xinc(nkd,ntniv+i_haplo)=pb_haplo(kd,i_haplo) 1611 end do 1612 end if !presentc 1613 end do 1614 end do 1615 end do 1616 1617 !ntniv=ntniv+nb_haplo_reduit 1618 ! nteff=nteff+nb_haplo_reduit 1619 1620 deallocate (pp_ldla) 1621 deallocate (pm_ldla) 1622 !ajout ofi 1623 !call free_internal_struct_ldla 1624 1625 end subroutine add_haplotype_effect
add_other_continue_traits_as_cov
[ Top ] [ m_qtlmap_incidence ] [ Subroutines ]
NAME
add_other_continue_traits_as_cov
DESCRIPTION
Add traits (that are not analysing in the current session) as covariate effects inside a contingence matrix with a description * An effect is adding for each covariate and fixed effect
INPUTS
mcov : covariate effect to not added in the contingence matrix (if mcov <= 0, all covariates effects are adding)
OUTPUTS
xinc : The contingence matrix incidenceDesc : The description of contingence matrix
NOTES
SOURCE
2064 subroutine add_other_continue_traits_as_cov(xinc,incidenceDesc,mcov,startTraitAsCov,endTraitAsCov) 2065 type(INCIDENCE_TYPE) , intent(inout) :: incidenceDesc 2066 real (kind=dp) , dimension(nd,incidenceDesc%ntnivmax) , intent(inout) :: xinc 2067 integer , intent(in) :: startTraitAsCov,endTraitAsCov,mcov ! indice on covariate list to remove (<=0 otherwise) 2068 2069 integer :: nbntniv,i,j,jj,nkd,ieff,ip,jm,kd,jcov 2070 2071 ieff=incidenceDesc%nteff 2072 ieff=ieff+1 2073 2074 nbntniv=0 2075 do i=startTraitAsCov,endTraitAsCov 2076 if ( i == incidenceDesc%ic )cycle 2077 if ( natureY(i) /= 'r') cycle 2078 nbntniv=nbntniv+1 2079 end do 2080 if ( nbntniv <= 0 ) return 2081 incidenceDesc%eqtl_print(ieff)=.false. 2082 incidenceDesc%desc(ieff)%isVar = .false. 2083 incidenceDesc%desc(ieff)%name="Traits as covariates" 2084 incidenceDesc%desc(ieff)%start=incidenceDesc%ntniv+1 2085 incidenceDesc%desc(ieff)%end=incidenceDesc%ntniv+nbntniv 2086 incidenceDesc%desc(ieff)%haveSubDesc=.true. 2087 2088 !init of xinc 2089 xinc(:,incidenceDesc%desc(ieff)%start:incidenceDesc%desc(ieff)%end)=0.d0 2090 2091 allocate ( incidenceDesc%desc(ieff)%listSubDesc(nbntniv) ) 2092 2093 jj=0;j=0 2094 do jcov=startTraitAsCov,endTraitAsCov 2095 if ( jcov == incidenceDesc%ic )cycle ! the current trait... 2096 if ( natureY(jcov) /= 'r') cycle ! can not be treat... 2097 jj=jj+1 2098 if ( jj == mcov ) cycle 2099 j=j+1 2100 incidenceDesc%desc(ieff)%listSubDesc(j)%name=trim(carac(jcov)) 2101 incidenceDesc%desc(ieff)%listSubDesc(j)%start=incidenceDesc%ntniv+j 2102 incidenceDesc%desc(ieff)%listSubDesc(j)%end=incidenceDesc%ntniv+j 2103 2104 nkd=0 2105 do ip=1,np 2106 do jm=nmp(ip)+1,nmp(ip+1) 2107 do kd=ndm(jm)+1,ndm(jm+1) 2108 if ( .not. presentc(incidenceDesc%ic,kd) ) cycle 2109 nkd=nkd+1 2110 if ( presentc(jcov,kd) ) then 2111 xinc(nkd,incidenceDesc%desc(ieff)%listSubDesc(j)%start)= y(jcov,kd) 2112 else 2113 xinc(nkd,incidenceDesc%desc(ieff)%listSubDesc(j)%start)= 0.d0 2114 end if 2115 end do !kd 2116 end do !jm 2117 end do !ip 2118 end do ! jcov 2119 incidenceDesc%nteff = incidenceDesc%nteff+1 2120 incidenceDesc%ntniv=incidenceDesc%ntniv+nbntniv 2121 2122 end subroutine add_other_continue_traits_as_cov 2123 2124 !remove the last effect describe by the INCIDENCE_TYPE 2125 subroutine remove_last_effect(incidenceDesc,xinc) 2126 type(INCIDENCE_TYPE) , intent(inout) :: incidenceDesc 2127 real (kind=dp) , dimension(nd,incidenceDesc%ntnivmax) , intent(inout) :: xinc 2128 2129 xinc(:,incidenceDesc%desc(incidenceDesc%nteff)%start:incidenceDesc%desc(incidenceDesc%nteff)%end)=0.d0 2130 incidenceDesc%ntniv=incidenceDesc%desc(incidenceDesc%nteff)%start-1 2131 ! deallocation of sub effect 2132 if (associated(incidenceDesc%desc(incidenceDesc%nteff)%listSubDesc)) then 2133 deallocate ( incidenceDesc%desc(incidenceDesc%nteff)%listSubDesc ) 2134 end if 2135 2136 incidenceDesc%nteff=incidenceDesc%nteff-1 2137 end subroutine remove_last_effect
add_polygenic_effect
[ Top ] [ m_qtlmap_incidence ] [ Subroutines ]
NAME
add_polygenic_effect
DESCRIPTION
Add a polygenic effect inside a contingence matrix with a description * An effect is adding for each Sire family * An effect is adding for each full sib family (depends ndmin variable)
OUTPUTS
xinc : The contingence matrix incidenceDesc : The description of contingence matrix
NOTES
SOURCE
1282 subroutine add_polygenic_effect(xinc,incidenceDesc) 1283 type(INCIDENCE_TYPE) , intent(inout) :: incidenceDesc 1284 real (kind=dp) , dimension(nd,incidenceDesc%ntnivmax) , intent(inout) :: xinc 1285 integer :: ip,jm,kd,nt,j,fem,nteff,nkd 1286 call log_mess("incidence : ** Add Polygenic effect to estimate **",DEBUG_DEF) 1287 nt=incidenceDesc%ntniv 1288 nteff = incidenceDesc%nteff 1289 1290 incidenceDesc%desc(incidenceDesc%nteff+1)%name="Sire polygenic effects" 1291 incidenceDesc%desc(incidenceDesc%nteff+1)%start=nt+1 1292 incidenceDesc%desc(incidenceDesc%nteff+1)%end=nt+np 1293 incidenceDesc%desc(incidenceDesc%nteff+1)%haveSubDesc=.true. 1294 incidenceDesc%desc(incidenceDesc%nteff+1)%isVar = .false. 1295 allocate (incidenceDesc%desc(incidenceDesc%nteff+1)%listSubDesc(np)) 1296 do ip=1,np 1297 incidenceDesc%desc(incidenceDesc%nteff+1)%listSubDesc(ip)%name="Sire "//trim(pere(ip)) 1298 incidenceDesc%desc(incidenceDesc%nteff+1)%listSubDesc(ip)%start=nt+ip 1299 incidenceDesc%desc(incidenceDesc%nteff+1)%listSubDesc(ip)%end=nt+ip 1300 end do 1301 incidenceDesc%nteff = incidenceDesc%nteff+1 1302 incidenceDesc%ntniv = incidenceDesc%ntniv+np 1303 1304 if ( count(estime(incidenceDesc%ic,:))>0) then 1305 incidenceDesc%desc(incidenceDesc%nteff+1)%name="Dam polygenic effects" 1306 incidenceDesc%desc(incidenceDesc%nteff+1)%start=incidenceDesc%desc(incidenceDesc%nteff)%end+1 1307 incidenceDesc%desc(incidenceDesc%nteff+1)%end=nt+np+count(estime(incidenceDesc%ic,:))!nbfemLoc(ic) 1308 incidenceDesc%desc(incidenceDesc%nteff+1)%haveSubDesc=.true. 1309 incidenceDesc%desc(incidenceDesc%nteff+1)%isVar = .false. 1310 allocate (incidenceDesc%desc(incidenceDesc%nteff+1)%listSubDesc(count(estime(incidenceDesc%ic,:)))) 1311 fem=0 1312 do ip=1,np 1313 do jm=nmp(ip)+1,nmp(ip+1) 1314 if (estime(incidenceDesc%ic,jm) ) then 1315 fem=fem+1 1316 incidenceDesc%desc(incidenceDesc%nteff+1)%listSubDesc(fem)%name="Dam "//trim(mere(jm))//" [Sire "//trim(pere(ip))//"]" 1317 incidenceDesc%desc(incidenceDesc%nteff+1)%listSubDesc(fem)%start=nt+np+fem 1318 incidenceDesc%desc(incidenceDesc%nteff+1)%listSubDesc(fem)%end=nt+np+fem 1319 end if 1320 end do 1321 end do 1322 incidenceDesc%nteff = incidenceDesc%nteff+1 1323 incidenceDesc%ntniv = incidenceDesc%ntniv+count(estime(incidenceDesc%ic,:))!nbfemLoc(ic) 1324 else 1325 1326 end if 1327 j=0 1328 1329 !initialisation 1330 xinc(:,nt+1:nt+np+count(estime(incidenceDesc%ic,:nm))) = 0 1331 1332 nkd=0 1333 do ip=1,np 1334 xinc(:,nt+ip) = 0 1335 do jm=nmp(ip)+1,nmp(ip+1) 1336 do kd=ndm(jm)+1,ndm(jm+1) 1337 if ( presentc(incidenceDesc%ic,kd) ) then 1338 nkd =nkd +1 1339 xinc(nkd,nt+ip) = 1 1340 ! incidenceDesc%nivdir(kd,nteff+1)= nt+ip 1341 if (estime(incidenceDesc%ic,jm)) then 1342 fem=count(estime(incidenceDesc%ic,:jm)) ! on compte le nombre estimable dans la famille de pere ip 1343 xinc(nkd,nt+np+fem) = 1 1344 ! incidenceDesc%nivdir(kd,nteff+2)= nt+np+fem 1345 end if 1346 end if 1347 end do 1348 end do 1349 end do 1350 end subroutine add_polygenic_effect
add_qtleffect
[ Top ] [ m_qtlmap_incidence ] [ Subroutines ]
NAME
add_qtleffect
DESCRIPTION
Add a qtl effect inside a contingence matrix with a description * An effect is adding for each Sire family * An effect is adding for each dam with enough progenies (ndmin)
INPUTS
nteff : index of the effect
numqtl : index of the qtl effect
chr : chromosome localisation
n : the position
mint : index of the interaction to remove (otherwise mint <= 0)
linear : use the likely probability genotype dam or not
OUTPUTS
xinc : The contingence matrix incidenceDesc : The description of contingence matrix
NOTES
SOURCE
1375 subroutine add_qtleffect(nteff,numqtl,chr,n,xinc,incidenceDesc,mint,linear) 1376 integer ,intent(in) :: nteff,numqtl,chr,n 1377 type(INCIDENCE_TYPE) ,intent(inout) :: incidenceDesc 1378 real (kind=dp) , dimension(nd,incidenceDesc%ntnivmax) ,intent(inout) :: xinc 1379 integer ,intent(in) :: mint ! index of interaction to remove from incidence 1380 logical ,intent(in) :: linear 1381 1382 integer :: ntniv,nteffstart,ip,fem,jm,l,s,u,nbint,jef 1383 logical ,dimension(nm) :: femInit 1384 1385 call log_mess("incidence : ** Add Qtl effect to estimate **",DEBUG_DEF) 1386 1387 nteffstart = incidenceDesc%nteff 1388 incidenceDesc%nteff = incidenceDesc%nteff+1 1389 if ( count(estime(incidenceDesc%ic,:)) > 0 ) incidenceDesc%nteff = incidenceDesc%nteff+1 1390 1391 nbint=listModelTrait(incidenceDesc%ic)%nbint(numqtl) !modele(incidenceDesc%ic,3) 1392 incidenceDesc%ntlev=1 1393 1394 do jef=1,nbint 1395 if ( mint == jef) cycle ! this interaction is not adding 1396 incidenceDesc%ntlev=incidenceDesc%ntlev*& 1397 nlev(listModelTrait(incidenceDesc%ic)%indexFixedEffectWithInteraction(numqtl,jef)) 1398 end do 1399 1400 incidenceDesc%desc(nteffstart+1)%isVar = .true. 1401 incidenceDesc%desc(nteffstart+1)%name="Sire QTL effects ["//trim(str(numqtl))//"]" 1402 incidenceDesc%desc(nteffstart+1)%start=incidenceDesc%ntniv+1 1403 incidenceDesc%ntniv_qtlsires(numqtl)=incidenceDesc%desc(nteffstart+1)%start 1404 incidenceDesc%desc(nteffstart+1)%end=incidenceDesc%ntniv+incidenceDesc%ntlev*np 1405 1406 incidenceDesc%desc(nteffstart+1)%haveSubDesc=.true. 1407 allocate (incidenceDesc%desc(nteffstart+1)%listSubDesc(np*incidenceDesc%ntlev)) 1408 u=0 1409 do ip=1,np 1410 l=0 1411 do s=incidenceDesc%ntniv+(ip-1)*incidenceDesc%ntlev+1,incidenceDesc%ntniv+ip*incidenceDesc%ntlev 1412 l=l+1 1413 u=u+1 1414 incidenceDesc%desc(nteffstart+1)%listSubDesc(u)%name="Sire "//trim(pere(ip))//" Level "//trim(str(l)) 1415 incidenceDesc%desc(nteffstart+1)%listSubDesc(u)%start=s 1416 incidenceDesc%desc(nteffstart+1)%listSubDesc(u)%end=s 1417 end do 1418 end do 1419 1420 if ( count(estime(incidenceDesc%ic,:)) > 0 ) then 1421 incidenceDesc%desc(nteffstart+2)%isVar = .true. 1422 incidenceDesc%desc(nteffstart+2)%name="Dam QTL effects ["//trim(str(numqtl))//"]" 1423 incidenceDesc%desc(nteffstart+2)%start=incidenceDesc%desc(nteffstart+1)%end+1 1424 incidenceDesc%ntniv_qtldams(numqtl)=incidenceDesc%desc(nteffstart+2)%start 1425 femInit=.false. 1426 l=0 1427 do jm=1,nm 1428 if (estime(incidenceDesc%ic,jm) ) then 1429 fem=iam(incidenceDesc%ic,repfem(jm)) 1430 if (femInit(fem)) then 1431 cycle 1432 end if 1433 femInit(fem)=.true. 1434 l=l+1 1435 end if 1436 end do 1437 1438 allocate (incidenceDesc%desc(nteffstart+2)%listSubDesc(l*incidenceDesc%ntlev)) 1439 incidenceDesc%desc(nteffstart+2)%end=incidenceDesc%ntniv+incidenceDesc%ntlev*(np+l) 1440 incidenceDesc%desc(nteffstart+2)%haveSubDesc=.true. 1441 !fem=0 1442 femInit=.false. 1443 u=0 1444 do jm=1,nm 1445 if (estime(incidenceDesc%ic,jm) ) then 1446 fem=iam(incidenceDesc%ic,repfem(jm)) 1447 if (femInit(fem)) then 1448 cycle 1449 end if 1450 femInit(fem)=.true. 1451 l=0 1452 do s=incidenceDesc%desc(nteffstart+2)%start+(fem-1)*incidenceDesc%ntlev,& 1453 incidenceDesc%desc(nteffstart+2)%start-1+fem*incidenceDesc%ntlev 1454 l=l+1 1455 u=u+1 1456 !fem=fem+1 1457 incidenceDesc%desc(nteffstart+2)%listSubDesc(u)%name="Dam "//trim(mere(jm))//" Level "//trim(str(l)) 1458 incidenceDesc%desc(nteffstart+2)%listSubDesc(u)%start=s 1459 incidenceDesc%desc(nteffstart+2)%listSubDesc(u)%end=s 1460 end do 1461 end if 1462 end do 1463 1464 incidenceDesc%ntniv = incidenceDesc%desc(nteffstart+2)%end 1465 else 1466 incidenceDesc%ntniv = incidenceDesc%desc(nteffstart+1)%end 1467 end if 1468 1469 call change_qtleffect(nteff,numqtl,chr,n,xinc,incidenceDesc,mint,linear) 1470 1471 end subroutine add_qtleffect
add_qtlinteraction
[ Top ] [ m_qtlmap_incidence ] [ Subroutines ]
NAME
add_qtlinteraction
DESCRIPTION
<DEV> modele interaction
NOTES
SOURCE
1786 subroutine add_qtlinteraction(nteff,chr1,chr2,n1,n2,xinc,incidenceDesc) 1787 integer ,intent(in) :: nteff,chr1,chr2,n1,n2 1788 type(INCIDENCE_TYPE) ,intent(inout) :: incidenceDesc 1789 real (kind=dp) , dimension(nd,incidenceDesc%ntnivmax) ,intent(inout) :: xinc 1790 1791 integer :: ntniv,nteffstart,ip,fem,jm,start,nkd 1792 logical ,dimension(nfem) :: femInit 1793 1794 call log_mess("incidence : ** Add Qtl Interaction effect to estimate **",DEBUG_DEF) 1795 1796 nteffstart = incidenceDesc%nteff 1797 incidenceDesc%nteff = incidenceDesc%nteff+1 1798 1799 incidenceDesc%desc(nteffstart+1)%name="QTL Interaction effects" 1800 incidenceDesc%desc(nteffstart+1)%start=incidenceDesc%ntniv+1 1801 incidenceDesc%desc(nteffstart+1)%end=incidenceDesc%ntniv 1802 1803 incidenceDesc%desc(nteffstart+1)%haveSubDesc=.true. 1804 1805 allocate ( incidenceDesc%desc(nteffstart+1)%listSubDesc(np+count(estime(incidenceDesc%ic,:))) ) 1806 1807 start=incidenceDesc%desc(nteffstart+1)%start 1808 1809 do ip=1,np 1810 ! on estime forcement les effets des haplotype pere QQ Qq qQ qq 1811 incidenceDesc%desc(nteffstart+1)%listSubDesc(ip)%name="Sire "//trim(pere(ip)) 1812 incidenceDesc%desc(nteffstart+1)%listSubDesc(ip)%start=start 1813 incidenceDesc%desc(nteffstart+1)%listSubDesc(ip)%end=start+4-1 1814 start=start+4 1815 incidenceDesc%desc(nteffstart+1)%end=incidenceDesc%desc(nteffstart+1)%end+4 1816 end do 1817 1818 fem=0 1819 do ip=1,np 1820 do jm=nmp(ip)+1,nmp(ip+1) 1821 ! on estime les effets combines pere et mere QQWW QQWw QQww QQwW..... 1822 if ( estime(incidenceDesc%ic,jm)) then 1823 fem=fem+1 1824 incidenceDesc%desc(nteffstart+1)%listSubDesc(np+fem)%name="Sire "//trim(pere(ip))//'- Dam '//trim(mere(jm)) 1825 incidenceDesc%desc(nteffstart+1)%listSubDesc(np+fem)%start=start 1826 incidenceDesc%desc(nteffstart+1)%listSubDesc(np+fem)%end=start+16-1 1827 start=start+16 1828 end if 1829 end do 1830 end do 1831 incidenceDesc%desc(nteffstart+1)%end=incidenceDesc%desc(nteffstart+1)%start+start 1832 1833 incidenceDesc%ntniv = start-1 1834 1835 call change_interaction_effect(nteff,chr1,chr2,n1,n2,xinc,incidenceDesc) 1836 1837 end subroutine add_qtlinteraction
build_incidence_matrix
[ Top ] [ m_qtlmap_incidence ] [ Subroutines ]
NAME
build_incidence_matrix
DESCRIPTION
Build the contingence matrix
INPUTS
hypothesis : the hypothesis
workstruct :
curPos : the initial position
meff : index of the fixed effect to no add
mcov : index of the covariate to no add
mint : index of the fixed effect-qtl interaction to not add
numqtl : number of qtl to test confusion with the mint interaction
OUTPUTS
incidenceDesc : description of the incidence
xinc : the contingence matrix to build
NOTES
SOURCE
3298 subroutine build_incidence_matrix(workstruct,curPos,xinc,incidenceDesc,meff,mcov,mint,numqtl) 3299 integer ,intent(in) :: meff,mcov,mint,numqtl 3300 type(POSITION_LRT_INCIDENCE) ,intent(in) :: curPos 3301 type(INCIDENCE_GEN_STRUCT) ,intent(inout) :: workstruct 3302 type(INCIDENCE_TYPE) ,intent(inout) :: incidenceDesc 3303 real (kind=dp) , dimension(nd,incidenceDesc%ntnivmax) ,intent(inout) :: Xinc 3304 3305 !local 3306 integer :: i,ic2,mint2 3307 3308 SELECT CASE(workstruct%type_incidence) 3309 CASE (INCIDENCE_TYPE_CLASSIC,INCIDENCE_INTERACTION) 3310 !add general mean to estim 3311 if (natureY(incidenceDesc%ic)/='i') then 3312 call add_general_mean(xinc,incidenceDesc) 3313 workstruct%eff_general_mean=1 3314 if (workstruct%nqtl>0) workstruct%listnteff(1)=2 3315 else 3316 workstruct%eff_general_mean=-1 3317 if (workstruct%nqtl>0) workstruct%listnteff(1)=1 3318 end if 3319 do i=1,workstruct%nqtl 3320 mint2=0 3321 !add qtl effect/interaction at position n to estim 3322 if ( numqtl == i ) mint2=mint 3323 call add_qtleffect(workstruct%listnteff(i),i,curPos%listChr(i),curPos%listN(i),& 3324 xinc,incidenceDesc,mint2,workstruct%linear) 3325 if (i<workstruct%nqtl) workstruct%listnteff(i+1)=incidenceDesc%nteff+1 3326 end do 3327 !add polygenic effect 3328 ! if(.not.workstruct%nopoly) call add_polygenic_effect(xinc,incidenceDesc) 3329 call add_polygenic_effect(xinc,incidenceDesc) 3330 !add fixed effect and covariate 3331 call add_effcov(xinc,incidenceDesc,meff,mcov) 3332 3333 CASE(INCIDENCE_TRAIT_COV) 3334 if (natureY(incidenceDesc%ic)/='r') then 3335 call stop_application(" ** INCIDENCE_TRAIT_COV : manage only continue trait ** ") 3336 end if 3337 if ( workstruct%startTraitAsCov == 0 .or. workstruct%endTraitAsCov ==0 ) then 3338 call stop_application(" ** INCIDENCE_TRAIT_COV : init workstruct%startTraitAsCov and workstruct%endTraitAsCov ** ") 3339 end if 3340 3341 call add_general_mean(xinc,incidenceDesc) 3342 workstruct%eff_general_mean=1 3343 if (workstruct%nqtl>0) workstruct%listnteff(1)=2 3344 do i=1,workstruct%nqtl 3345 mint2=0 3346 !add qtl effect/interaction at position n to estim 3347 if ( numqtl == i ) mint2=mint 3348 call add_qtleffect(workstruct%listnteff(i),i,curPos%listChr(i),curPos%listN(i),& 3349 xinc,incidenceDesc,mint2,workstruct%linear) 3350 if (i<workstruct%nqtl) workstruct%listnteff(i+1)=incidenceDesc%nteff+1 3351 end do 3352 !add polygenic effect 3353 ! if(.not.workstruct%nopoly) call add_polygenic_effect(xinc,incidenceDesc) 3354 call add_polygenic_effect(xinc,incidenceDesc) 3355 !add others traits as covariate 3356 call add_other_continue_traits_as_cov(xinc,incidenceDesc,mcov,& 3357 workstruct%startTraitAsCov,workstruct%endTraitAsCov) 3358 !add fixed effect and covariate 3359 call add_effcov(xinc,incidenceDesc,meff,mcov) 3360 3361 ! call debug_write_incidence(xinc,incidenceDesc) 3362 3363 CASE (INCIDENCE_TYPE_LDLA) 3364 !add general mean to estim 3365 call add_general_mean(xinc,incidenceDesc) 3366 workstruct%eff_general_mean=1 3367 !haplotype effect 3368 do i=1,workstruct%nqtl 3369 workstruct%listnteffhap(i)=i+1 3370 call add_haplotype_effect(workstruct%listnteffhap(i),i,curPos%listChr(i),curPos%listN(i),& 3371 xinc,incidenceDesc,workstruct%hdam,.true.) 3372 end do 3373 !qtl effect 3374 do i=1,workstruct%nqtl 3375 mint2=0 3376 workstruct%listnteff(i)=incidenceDesc%nteff+1 3377 !add qtl effect/interaction at position n to estim 3378 if ( numqtl == i ) mint2=mint 3379 call add_qtleffect(workstruct%listnteff(i),i,curPos%listChr(i),curPos%listN(i),& 3380 xinc,incidenceDesc,mint2,workstruct%linear) 3381 end do 3382 !add polygenic effect 3383 3384 ! if(.not.workstruct%nopoly) call add_polygenic_effect(xinc,incidenceDesc) 3385 call add_polygenic_effect(xinc,incidenceDesc) 3386 3387 !add fixed effect and covariate 3388 call add_effcov(xinc,incidenceDesc,meff,mcov) 3389 CASE (INCIDENCE_TYPE_LD) 3390 !add general mean to estim 3391 if (natureY(incidenceDesc%ic)/='i') then 3392 call add_general_mean(xinc,incidenceDesc) 3393 workstruct%eff_general_mean=1 3394 do i=1,workstruct%nqtl 3395 workstruct%listnteffhap(i)=i+1 3396 call add_haplotype_effect(workstruct%listnteffhap(i),i,curPos%listChr(i),curPos%listN(i),& 3397 xinc,incidenceDesc,workstruct%hdam,.true.) 3398 end do 3399 !add polygenic effect 3400 3401 ! if(.not.workstruct%nopoly) call add_polygenic_effect(xinc,incidenceDesc) 3402 call add_polygenic_effect(xinc,incidenceDesc) 3403 3404 !add fixed effect and covariate 3405 call add_effcov(xinc,incidenceDesc,meff,mcov) 3406 end if 3407 3408 CASE DEFAULT 3409 call stop_application("build_incidence_matrix:Unknown data type description for model ["//& 3410 trim(str(workstruct%type_incidence))//"]") 3411 END SELECT 3412 end subroutine build_incidence_matrix
change_interaction_effect
[ Top ] [ m_qtlmap_incidence ] [ Subroutines ]
NAME
change_interaction_effect
DESCRIPTION
<DEV> modele interaction
NOTES
SOURCE
1847 subroutine change_interaction_effect(nteff,chr1,chr2,n1,n2,xinc,incidenceDesc) 1848 integer ,intent(in) :: nteff,chr1,chr2,n1,n2 1849 type(INCIDENCE_TYPE) ,intent(inout) :: incidenceDesc 1850 real (kind=dp) , dimension(nd,incidenceDesc%ntnivmax) ,intent(inout) :: xinc 1851 1852 !local 1853 integer :: ip,jm,kd,ntniv,kkd,ig1,ig2,kd2,fem,nkd 1854 1855 1856 ntniv = incidenceDesc%desc(nteff)%start -1 1857 ! print *,"NTNIV:",ntniv,' NTEFF:',nteff 1858 nkd=0 1859 1860 do ip=1,np 1861 do jm=nmp(ip)+1,nmp(ip+1) 1862 ig1=likelyprob(chr1,ip,jm) 1863 do kd=ngend(chr1,ig1)+1,ngend(chr1,ig1+1) 1864 kkd = ndesc(chr1,kd) 1865 if(presentc(incidenceDesc%ic,kkd)) then 1866 !il faut rechercher le kd associe au chr2.... 1867 do kd2=ngend(chr1,ig1)+1,ngend(chr1,ig1+1) 1868 if ( kkd == ndesc(chr2,kd2)) exit 1869 end do 1870 nkd=nkd+1 1871 !naif pour l instant.... 1872 ! effet haplotype pere 1873 ! QW, Qw qW qw 1874 xinc(nkd,ntniv+(ip-1)*4+1) = ( pdd(chr1,kd,1,n1)+ pdd(chr1,kd,2,n1)) * (pdd(chr2,kd2,1,n2)+pdd(chr2,kd2,2,n2)) 1875 xinc(nkd,ntniv+(ip-1)*4+2) = ( pdd(chr1,kd,1,n1)+ pdd(chr1,kd,2,n1)) * (pdd(chr2,kd2,3,n2)+pdd(chr2,kd2,4,n2)) 1876 xinc(nkd,ntniv+(ip-1)*4+3) = ( pdd(chr1,kd,3,n1)+ pdd(chr1,kd,4,n1)) * (pdd(chr2,kd2,1,n2)+pdd(chr2,kd2,2,n2)) 1877 xinc(nkd,ntniv+(ip-1)*4+4) = ( pdd(chr1,kd,3,n1)+ pdd(chr1,kd,4,n1)) * (pdd(chr2,kd2,3,n2)+pdd(chr2,kd2,4,n2)) 1878 1879 if ( estime (incidenceDesc%ic,jm) ) then 1880 ! effet haplotype combine QQWW QqWW qQWW qqWW QQwW QqwW qQwW qqwW QQWw QqWw qQWw qqWw QQww Qqww qQww qqww 1881 fem=count(estime(incidenceDesc%ic,:jm)) 1882 xinc(nkd,ntniv+np*4+(fem-1)*16+1) = pdd(chr1,kd,1,n1) * pdd(chr2,kd2,1,n2) ! QQWW 1883 xinc(nkd,ntniv+np*4+(fem-1)*16+2) = pdd(chr1,kd,3,n1) * pdd(chr2,kd2,1,n2) ! QqWW 1884 xinc(nkd,ntniv+np*4+(fem-1)*16+3) = pdd(chr1,kd,2,n1) * pdd(chr2,kd2,1,n2) ! qQWW 1885 xinc(nkd,ntniv+np*4+(fem-1)*16+4) = pdd(chr1,kd,4,n1) * pdd(chr2,kd2,1,n2) ! qqWW 1886 1887 xinc(nkd,ntniv+np*4+(fem-1)*16+5) = pdd(chr1,kd,1,n1) * pdd(chr2,kd2,2,n2) ! QQwW 1888 xinc(nkd,ntniv+np*4+(fem-1)*16+6) = pdd(chr1,kd,3,n1) * pdd(chr2,kd2,2,n2) ! QqwW 1889 xinc(nkd,ntniv+np*4+(fem-1)*16+7) = pdd(chr1,kd,2,n1) * pdd(chr2,kd2,2,n2) ! qQwW 1890 xinc(nkd,ntniv+np*4+(fem-1)*16+8) = pdd(chr1,kd,4,n1) * pdd(chr2,kd2,2,n2) ! qqwW 1891 1892 xinc(nkd,ntniv+np*4+(fem-1)*16+9) = pdd(chr1,kd,1,n1) * pdd(chr2,kd2,3,n2) ! QQWw 1893 xinc(nkd,ntniv+np*4+(fem-1)*16+10) = pdd(chr1,kd,3,n1) * pdd(chr2,kd2,3,n2) ! QqWw 1894 xinc(nkd,ntniv+np*4+(fem-1)*16+11) = pdd(chr1,kd,2,n1) * pdd(chr2,kd2,3,n2) ! qQWw 1895 xinc(nkd,ntniv+np*4+(fem-1)*16+12) = pdd(chr1,kd,4,n1) * pdd(chr2,kd2,3,n2) ! qqWw 1896 1897 xinc(nkd,ntniv+np*4+(fem-1)*16+13) = pdd(chr1,kd,1,n1) * pdd(chr2,kd2,4,n2) ! QQww 1898 xinc(nkd,ntniv+np*4+(fem-1)*16+14) = pdd(chr1,kd,3,n1) * pdd(chr2,kd2,4,n2) ! Qqww 1899 xinc(nkd,ntniv+np*4+(fem-1)*16+15) = pdd(chr1,kd,2,n1) * pdd(chr2,kd2,4,n2) ! qQww 1900 xinc(nkd,ntniv+np*4+(fem-1)*16+16) = pdd(chr1,kd,4,n1) * pdd(chr2,kd2,4,n2) ! qqww 1901 end if 1902 end if 1903 end do ! kd 1904 end do ! jm 1905 end do ! ip 1906 1907 1908 end subroutine change_interaction_effect
change_qtleffect
[ Top ] [ m_qtlmap_incidence ] [ Subroutines ]
NAME
change_qtleffect
DESCRIPTION
Change the contingence matrix (probabilities of transmission information) during the chromosome scan .
INPUTS
nteff : index of the effect
iq : index of the qtl effect
chr : chromosome localisation
n : the position
mint : index of the interaction to remove (otherwise mint <= 0)
linear : use the likely probability genotype dam or not
OUTPUTS
xinc : The contingence matrix incidenceDesc : The description of contingence matrix
NOTES
SOURCE
1647 subroutine change_qtleffect(nteff,iq,chr,n,xinc,incidenceDesc,mint,linear) 1648 integer ,intent(in) :: nteff,iq,chr,n 1649 type(INCIDENCE_TYPE) ,intent(inout) :: incidenceDesc 1650 real (kind=dp) , dimension(nd,incidenceDesc%ntnivmax) ,intent(inout) :: xinc 1651 integer ,intent(in) :: mint ! index of interaction to remove from incidence 1652 logical ,intent(in) :: linear 1653 1654 !local 1655 integer :: nbint,jef,nivint,ntlev,ip,jm,kd,ntt,km,ntniv,nkd 1656 real (kind=dp) ,dimension(:),allocatable:: pp,pm 1657 1658 allocate (pp(nd)) 1659 allocate (pm(nd)) 1660 1661 nbint=listModelTrait(incidenceDesc%ic)%nbint(iq) 1662 ntniv = incidenceDesc%desc(nteff)%start -1 1663 nkd=0 1664 1665 do ip=1,np 1666 do jm=nmp(ip)+1,nmp(ip+1) 1667 if ( linear ) then 1668 call getprob_lin(chr,n,incidenceDesc%ic,ip,jm,ndm(jm)+1,ndm(jm+1),pp,pm) 1669 else 1670 call getglobalprob(chr,n,incidenceDesc%ic,ip,jm,ndm(jm)+1,ndm(jm+1),pp,pm) 1671 end if 1672 do kd=ndm(jm)+1,ndm(jm+1) 1673 if ( presentc(incidenceDesc%ic,kd) ) then 1674 nkd=nkd+1 1675 1676 ! si des effets fixes sont en interaction avec l'allele paternel au qtl 1677 ! il faut estimer les effets qtl pour chacun des niveaux concernes 1678 ! 1679 ! on commence donc par creer un effet composites regroupant l'ensemble de 1680 ! ces effets 1681 ! 1682 nivint=1 1683 ntlev=1 1684 !nt=ntniv 1685 do jef=1,nbint 1686 if ( mint == jef) cycle ! this interaction is not adding 1687 nivint=nivint+ntlev*(nivx(kd,listModelTrait(incidenceDesc%ic)%indexFixedEffectWithInteraction(iq,jef))-1) 1688 ntlev=ntlev*nlev(listModelTrait(incidenceDesc%ic)%indexFixedEffectWithInteraction(iq,jef)) 1689 end do 1690 ntt=ntniv+nivint+ntlev*(ip-1) 1691 xinc(nkd,ntt)=pp(kd) 1692 1693 ! incidenceDesc%nivdir(kd,nteff+1)=ntt 1694 ! print *,'nt:',nivdir(kd,nteff+1),' vp:',pp(kd),' nteff:',nteff,' kd:',kd 1695 1696 if(estime(incidenceDesc%ic,jm)) then 1697 km=iam(incidenceDesc%ic,repfem(jm)) 1698 ntt= ntniv+ntlev*np+nivint+ntlev*(km-1) 1699 xinc(nkd,ntt)=pm(kd) 1700 ! incidenceDesc%nivdir(kd,nteff+2)=ntt 1701 ! print *,'nt:',nivdir(kd,nteff+2),' vp:',pm(kd) 1702 end if 1703 end if 1704 end do 1705 end do 1706 end do 1707 1708 deallocate (pp) 1709 deallocate (pm) 1710 ! a modifier 1711 ! --> on devrait integerer l intialisation des pdds dans la boucle precedente 1712 ! attention les interaction qtl ne sont pas encore pris en compte 1713 if (.not. linear) then 1714 call pdd_at_position(incidenceDesc%ic,incidenceDesc%dataset,iq,chr,n) 1715 end if 1716 end subroutine change_qtleffect
confusion_between_qtl_and_effect
[ Top ] [ m_qtlmap_incidence ] [ Subroutines ]
NAME
confusion_between_qtl_and_effect
DESCRIPTION
Detect a confusion betwwen a qtl effect and the other effect described by the model.
INPUTS
incidenceDesc : the description of the incidence matrix.
xcor : the correlation matrix
index_effect : the index qtl effect to test with the other effect
OUTPUTS
sizealert : number of alert detected
alertCorrQtl : the correlation
corrmax : the maximum correlation finded succesively
NOTES
xxx come from the get_precision function.
SOURCE
2420 subroutine confusion_between_qtl_and_effect(incidenceDesc,xcor,index_effect,& 2421 sizealert,alertCorrQtl,nalert,corrmax) 2422 type(INCIDENCE_TYPE) ,intent(in) :: incidenceDesc 2423 real (kind=dp) ,dimension(incidenceDesc%ntnivmax,incidenceDesc%ntnivmax) ,intent(in) :: xcor 2424 integer ,intent(in) :: index_effect 2425 integer ,intent(in) :: sizealert 2426 type(CORR_ALERT_TYPE) ,dimension(sizealert) ,intent(out) :: alertCorrQtl 2427 integer ,intent(out) :: nalert 2428 real (kind=dp) ,intent(inout) :: corrmax 2429 2430 integer :: numqtl,ip,jm,iqtl,jniv,ntlev,id,s,e,i,fem 2431 logical ,dimension(:),allocatable :: femInit 2432 2433 nalert=0 2434 allocate (femInit(maxval(iam))) 2435 ! 2436 ! on caclule les correlations entre les colonnes de X'X des effets qtl pere et de 2437 ! tous les autres effets 2438 ! 2439 do numqtl=1,incidenceDesc%nqtl 2440 do ip=1,np 2441 ntlev=0 2442 !print *,"ip:",ip 2443 s=incidenceDesc%ntniv_qtlsires(numqtl)+(ip-1)*incidenceDesc%ntlev 2444 e=s+incidenceDesc%ntlev-1 2445 do iqtl=s,e 2446 ! print *,"iqtl:",iqtl 2447 ntlev=ntlev+1 2448 2449 if (incidenceDesc%vecsol(iqtl))then 2450 ! 2451 ! test d'une confusion avec l'effect concerné 2452 ! 2453 2454 do jniv=incidenceDesc%desc(index_effect)%start,incidenceDesc%desc(index_effect)%end 2455 if (abs(xcor(iqtl,jniv))>=corrmax) corrmax=xcor(iqtl,jniv) 2456 ! si la correlation est superieure a corMax, on imprime une alerte 2457 2458 if(dabs(xcor(iqtl,jniv)).ge.THRES_CONFUSION)then 2459 nalert=nalert+1 2460 alertCorrQtl(nalert)%qtl=numqtl 2461 alertCorrQtl(nalert)%ip=ip 2462 alertCorrQtl(nalert)%corr=xcor(iqtl,jniv) 2463 alertCorrQtl(nalert)%ntlev=ntlev 2464 id=jniv-incidenceDesc%desc(index_effect)%start+1 2465 alertCorrQtl(nalert)%name_effect=incidenceDesc%desc(index_effect)%name 2466 alertCorrQtl(nalert)%name_level=incidenceDesc%desc(index_effect)%listSubDesc(id)%name 2467 end if 2468 end do ! jniv 2469 end if 2470 end do ! iqtl 2471 ! 2472 ! fin pour le niveau du qtl du pere 2473 ! 2474 end do ! ip 2475 2476 femInit=.false. 2477 do jm=1,nm 2478 if (estime(incidenceDesc%ic,jm) ) then 2479 fem=iam(incidenceDesc%ic,repfem(jm)) 2480 if (femInit(fem)) then 2481 cycle 2482 end if 2483 femInit(fem)=.true. 2484 ! print *,'jm:',jm 2485 ntlev=0 2486 s=incidenceDesc%ntniv_qtldams(numqtl)+(fem-1)*incidenceDesc%ntlev 2487 e=s+incidenceDesc%ntlev-1 2488 do iqtl=s,e 2489 ! print *,"iqtl:",iqtl 2490 ntlev=ntlev+1 2491 2492 if (incidenceDesc%vecsol(iqtl))then 2493 ! 2494 ! test d'une confusion avec l'effect concerné 2495 ! 2496 do jniv=incidenceDesc%desc(index_effect)%start,incidenceDesc%desc(index_effect)%end 2497 ! print *,'jniv:',jniv 2498 ! print *,abs(xcor(iqtl,jniv)) 2499 if (abs(xcor(iqtl,jniv))>=corrmax) corrmax=xcor(iqtl,jniv) 2500 ! si la correlation est superieure a corMax, on imprime une alerte 2501 ! 2502 if(dabs(xcor(iqtl,jniv)).ge.THRES_CONFUSION)then 2503 nalert=nalert+1 2504 alertCorrQtl(nalert)%qtl=numqtl 2505 alertCorrQtl(nalert)%jm=jm 2506 alertCorrQtl(nalert)%corr=xcor(iqtl,jniv) 2507 alertCorrQtl(nalert)%ntlev=ntlev 2508 id=jniv-incidenceDesc%desc(index_effect)%start+1 2509 alertCorrQtl(nalert)%name_effect=incidenceDesc%desc(index_effect)%name 2510 alertCorrQtl(nalert)%name_level=incidenceDesc%desc(index_effect)%listSubDesc(id)%name 2511 end if 2512 end do ! jniv 2513 end if 2514 end do ! iqtl 2515 2516 ! 2517 ! fin pour le niveau du qtl du pere 2518 ! 2519 end if 2520 end do ! jm 2521 end do ! numqtl 2522 deallocate (femInit) 2523 2524 end subroutine confusion_between_qtl_and_effect
confusion_type1
[ Top ] [ m_qtlmap_incidence ] [ Subroutines ]
NAME
confusion_type1
DESCRIPTION
INPUTS
incidenceDesc : the description of the incidence matrix.
xxx :
OUTPUTS
workstruct :
NOTES
xxx come from the get_precision function.
SOURCE
2541 subroutine confusion_type1(incidenceDesc,xxx,workstruct) 2542 type(INCIDENCE_TYPE) ,intent(in) :: incidenceDesc 2543 real (kind=dp) ,dimension(incidenceDesc%ntnivmax,incidenceDesc%ntnivmax) ,intent(in) :: xxx 2544 type(INCIDENCE_GEN_STRUCT) ,intent(inout) :: workstruct 2545 2546 real (kind=dp) ,dimension(incidenceDesc%ntnivmax,incidenceDesc%ntnivmax) :: xcor 2547 type(CORR_ALERT_TYPE) ,dimension(:),allocatable :: alerts 2548 type(CORR_ALERT_TYPE) ,dimension(:),allocatable :: as 2549 integer :: ialert,sizealert,eff,i,j,starteff 2550 2551 call get_correlations(incidenceDesc,xxx,xcor) 2552 workstruct%corrmax=0.d0 2553 2554 sizealert=( np+count(estime(incidenceDesc%ic,:)) )*incidenceDesc%ntlev*incidenceDesc%nbnivest 2555 2556 allocate (alerts(sizealert)) 2557 allocate (as(sizealert*incidenceDesc%nqtl)) 2558 2559 workstruct%nalert=0 2560 starteff = workstruct%listnteff(workstruct%nqtl)+1 2561 if (count(estime(incidenceDesc%ic,:))>0) starteff = starteff+1 2562 do eff=starteff,incidenceDesc%nteff 2563 call log_mess("computing confusion beetween qtls and:"//incidenceDesc%desc(eff)%name,VERBOSE_DEF) 2564 call confusion_between_qtl_and_effect(incidenceDesc,xcor,& 2565 eff,sizealert,alerts,ialert,workstruct%corrmax) 2566 if (ialert>0) then 2567 ! print *,"ialert:",alerts(1)%corr,alerts(2)%corr,alerts(3)%corr 2568 as(workstruct%nalert+1:workstruct%nalert+ialert)=alerts(:ialert) 2569 workstruct%nalert=workstruct%nalert+ialert 2570 end if 2571 end do 2572 2573 allocate (workstruct%alertQtl(workstruct%nalert)) 2574 workstruct%alertQtl=as(:workstruct%nalert) 2575 2576 deallocate (alerts) 2577 deallocate (as) 2578 2579 end subroutine confusion_type1
copy_dataset
[ Top ] [ m_qtlmap_incidence ] [ Subroutines ]
NAME
copy_dataset
DESCRIPTION
Allocate new arrays and copy a variable of type DATASET_TYPE.
INPUTS
dataset : the original variable
OUTPUTS
out_dataset : the copy of dataset
NOTES
SOURCE
894 subroutine copy_dataset(dataset,out_dataset) 895 type(DATASET_TYPE) , intent(in) :: dataset 896 type(DATASET_TYPE) , intent(out) :: out_dataset 897 integer :: ip,jm,ifem,s1,s2,s3 898 899 out_dataset%nkd = dataset%nkd 900 out_dataset%nkd_max_by_fam = dataset%nkd_max_by_fam 901 902 allocate (out_dataset%lSires(np)) 903 904 do ip=1,np 905 out_dataset%lSires(ip)%half_sib%firstKd = dataset%lSires(ip)%half_sib%firstKd 906 out_dataset%lSires(ip)%half_sib%lastKd = dataset%lSires(ip)%half_sib%lastKd 907 908 if ( associated(dataset%lSires(ip)%ppt) ) then 909 s1=size(dataset%lSires(ip)%ppt,1) 910 s2=size(dataset%lSires(ip)%ppt,2) 911 allocate (out_dataset%lSires(ip)%ppt(s1,s2)) 912 out_dataset%lSires(ip)%ppt = dataset%lSires(ip)%ppt 913 end if 914 915 if ( associated(dataset%lSires(ip)%full_sib) ) then 916 917 nfem = size(dataset%lSires(ip)%full_sib) 918 allocate ( out_dataset%lSires(ip)%full_sib(nfem) ) 919 do ifem=1,nfem 920 out_dataset%lSires(ip)%full_sib(ifem)%firstKd = dataset%lSires(ip)%full_sib(ifem)%firstKd 921 out_dataset%lSires(ip)%full_sib(ifem)%lastKd = dataset%lSires(ip)%full_sib(ifem)%lastKd 922 if ( associated(dataset%lSires(ip)%full_sib(ifem)%ppt) ) then 923 s1=size(dataset%lSires(ip)%full_sib(ifem)%ppt,1) 924 s2=size(dataset%lSires(ip)%full_sib(ifem)%ppt,2) 925 s3=size(dataset%lSires(ip)%full_sib(ifem)%ppt,3) 926 allocate ( out_dataset%lSires(ip)%full_sib(ifem)%ppt(s1,s2,s3) ) 927 out_dataset%lSires(ip)%full_sib(ifem)%ppt = dataset%lSires(ip)%full_sib(ifem)%ppt 928 end if 929 if ( associated(dataset%lSires(ip)%full_sib(ifem)%pmt) ) then 930 s1=size(dataset%lSires(ip)%full_sib(ifem)%pmt,1) 931 s2=size(dataset%lSires(ip)%full_sib(ifem)%pmt,2) 932 s3=size(dataset%lSires(ip)%full_sib(ifem)%pmt,3) 933 allocate ( out_dataset%lSires(ip)%full_sib(ifem)%pmt(s1,s2,s3) ) 934 out_dataset%lSires(ip)%full_sib(ifem)%pmt = dataset%lSires(ip)%full_sib(ifem)%pmt 935 end if 936 end do 937 938 end if 939 940 if ( associated(dataset%lSires(ip)%sig0) ) then 941 allocate (out_dataset%lSires(ip)%sig0(size(dataset%lSires(ip)%sig0))) 942 out_dataset%lSires(ip)%sig0 = dataset%lSires(ip)%sig0 943 end if 944 945 if ( associated(dataset%lSires(ip)%xmu0) ) then 946 allocate (out_dataset%lSires(ip)%xmu0(size(dataset%lSires(ip)%xmu0))) 947 out_dataset%lSires(ip)%xmu0 = dataset%lSires(ip)%xmu0 948 end if 949 950 end do 951 952 953 if (associated (dataset%Y)) then 954 s1=size(dataset%Y,1) 955 s2=size(dataset%Y,2) 956 allocate(out_dataset%Y(s1,s2)) 957 allocate(out_dataset%CD(s1,s2)) 958 out_dataset%Y = dataset%Y 959 out_dataset%CD = dataset%CD 960 else 961 s1=size(dataset%YDISCRETORD,1) 962 s2=size(dataset%YDISCRETORD,2) 963 allocate(out_dataset%YDISCRETORD(s1,s2)) 964 out_dataset%YDISCRETORD = dataset%YDISCRETORD 965 end if 966 967 if (associated (dataset%rhoi) ) then 968 s1=size(dataset%rhoi,1) 969 s2=size(dataset%rhoi,2) 970 allocate(out_dataset%rhoi(s1,s2)) 971 out_dataset%rhoi = dataset%rhoi 972 end if 973 974 if (associated (dataset%cov) ) then 975 s1=size(dataset%cov,1) 976 s2=size(dataset%cov,2) 977 allocate(out_dataset%cov(s1,s2)) 978 out_dataset%cov = dataset%cov 979 end if 980 981 if (associated (dataset%sig) ) then 982 s1=size(dataset%sig,1) 983 allocate(out_dataset%sig(s1)) 984 out_dataset%sig = dataset%sig 985 end if 986 987 if (associated (dataset%xmu) ) then 988 s1=size(dataset%xmu,1) 989 allocate(out_dataset%xmu(s1)) 990 out_dataset%xmu = dataset%xmu 991 end if 992 993 end subroutine copy_dataset
copy_incidence_desc
[ Top ] [ m_qtlmap_incidence ] [ Subroutines ]
NAME
copy_incidence_desc
DESCRIPTION
Allocate partial arrays and copy a variable of type DATASET_TYPE. Static informations are referenced from the original variable in_inc.
INPUTS
in_inc : the original variable
OUTPUTS
out_inc : the copy of incidence description
NOTES
SOURCE
1010 subroutine copy_incidence_desc (in_inc,out_inc) 1011 type(INCIDENCE_TYPE) , intent(in) :: in_inc 1012 type(INCIDENCE_TYPE) , intent(out) :: out_inc 1013 integer :: ip,ifem,jm,s1,s2,n,i 1014 1015 out_inc%ic = in_inc%ic 1016 out_inc%dataset => in_inc%dataset ! attention on ne recopie pas le jeu de donnees qui est considerer comme statique et non modifiable 1017 out_inc%nqtl = in_inc%nqtl 1018 out_inc%ntniv = in_inc%ntniv 1019 out_inc%nteff = in_inc%nteff 1020 out_inc%eqtl_print = in_inc%eqtl_print 1021 out_inc%nbnivest = in_inc%nbnivest 1022 out_inc%ntlev = in_inc%ntlev 1023 out_inc%nbniv = in_inc%nbniv 1024 1025 if ( associated (in_inc%desc) ) then 1026 allocate (out_inc%desc(size(in_inc%desc))) 1027 do i=1,size(in_inc%desc) 1028 call copy(in_inc%desc(i),out_inc%desc(i)) 1029 end do 1030 else 1031 out_inc%desc => null() 1032 end if 1033 1034 allocate(out_inc%vecsol(size(in_inc%vecsol))) 1035 allocate(out_inc%precis(size(in_inc%precis))) 1036 if ( out_inc%nqtl > 0) then 1037 allocate (out_inc%ntniv_qtlsires(out_inc%nqtl)) 1038 allocate (out_inc%ntniv_qtldams(out_inc%nqtl)) 1039 out_inc%ntniv_qtlsires=0 1040 out_inc%ntniv_qtldams=0 1041 end if 1042 allocate (out_inc%corr_niv_nivb(size(in_inc%corr_niv_nivb))) 1043 1044 out_inc%vecsol = in_inc%vecsol 1045 out_inc%precis = in_inc%precis 1046 out_inc%ntniv_qtlsires = in_inc%ntniv_qtlsires 1047 out_inc%ntniv_qtldams = in_inc%ntniv_qtldams 1048 out_inc%corr_niv_nivb = in_inc%corr_niv_nivb 1049 1050 allocate (out_inc%dataset) 1051 !call init_dataset_1car(out_inc%ic,out_inc%nqtl,out_inc%dataset) 1052 allocate(out_inc%dataset%lSires(np)) 1053 do ip=1,np 1054 allocate(out_inc%dataset%lSires(ip)%full_sib(nmp(ip+1)-nmp(ip))) 1055 ifem=0 1056 do jm=nmp(ip)+1,nmp(ip+1) 1057 ifem=ifem+1 1058 out_inc%dataset%lSires(ip)%full_sib(ifem)%firstKd=in_inc%dataset%lSires(ip)%full_sib(ifem)%firstKd 1059 out_inc%dataset%lSires(ip)%full_sib(ifem)%lastKd=in_inc%dataset%lSires(ip)%full_sib(ifem)%lastKd 1060 1061 if ( associated(in_inc%dataset%lSires(ip)%full_sib(ifem)%ppt) ) then 1062 n=size(in_inc%dataset%lSires(ip)%full_sib(ifem)%ppt,1) 1063 s1=size(in_inc%dataset%lSires(ip)%full_sib(ifem)%ppt,2) 1064 s2=size(in_inc%dataset%lSires(ip)%full_sib(ifem)%ppt,3) 1065 allocate(out_inc%dataset%lSires(ip)%full_sib(ifem)%ppt(n,s1,s2)) 1066 allocate(out_inc%dataset%lSires(ip)%full_sib(ifem)%pmt(n,s1,s2)) 1067 out_inc%dataset%lSires(ip)%full_sib(ifem)%ppt=in_inc%dataset%lSires(ip)%full_sib(ifem)%ppt 1068 out_inc%dataset%lSires(ip)%full_sib(ifem)%pmt=in_inc%dataset%lSires(ip)%full_sib(ifem)%pmt 1069 end if 1070 1071 end do 1072 1073 out_inc%dataset%lSires(ip)%half_sib%firstKd=in_inc%dataset%lSires(ip)%half_sib%firstKd 1074 out_inc%dataset%lSires(ip)%half_sib%lastKd=in_inc%dataset%lSires(ip)%half_sib%lastKd 1075 if ( associated(in_inc%dataset%lSires(ip)%ppt) ) then 1076 n=size(in_inc%dataset%lSires(ip)%ppt,1) 1077 s1=in_inc%dataset%lSires(ip)%half_sib%firstKd 1078 s2=in_inc%dataset%lSires(ip)%half_sib%lastKd 1079 allocate(out_inc%dataset%lSires(ip)%ppt(n,s1:s2)) 1080 out_inc%dataset%lSires(ip)%ppt=in_inc%dataset%lSires(ip)%ppt 1081 end if 1082 1083 out_inc%dataset%lSires(ip)%sig0=>in_inc%dataset%lSires(ip)%sig0 1084 out_inc%dataset%lSires(ip)%xmu0=>in_inc%dataset%lSires(ip)%xmu0 1085 1086 end do 1087 1088 out_inc%dataset%nkd = in_inc%dataset%nkd 1089 out_inc%dataset%nkd_max_by_fam = in_inc%dataset%nkd_max_by_fam 1090 out_inc%ntnivmax =in_inc%ntnivmax 1091 out_inc%nteffmax =in_inc%nteffmax 1092 1093 out_inc%dataset%YDISCRETORD => in_inc%dataset%YDISCRETORD 1094 out_inc%dataset%Y => in_inc%dataset%Y 1095 out_inc%dataset%CD => in_inc%dataset%CD 1096 out_inc%dataset%cov => in_inc%dataset%cov 1097 out_inc%dataset%rhoi => in_inc%dataset%rhoi 1098 1099 ! --- Fin copy dataset 1100 1101 if (associated(in_inc%nonull) ) then 1102 allocate (out_inc%nonull(size(in_inc%nonull,1),size(in_inc%nonull,2))) 1103 allocate (out_inc%nbnonull(size(in_inc%nbnonull))) 1104 out_inc%nonull = in_inc%nonull 1105 out_inc%nbnonull = in_inc%nbnonull 1106 end if 1107 1108 !les bornes sont sense ne pas change... 1109 out_inc%borni=>in_inc%borni 1110 out_inc%borns=>in_inc%borns 1111 1112 allocate (out_inc%par(size(in_inc%par))) 1113 out_inc%par = in_inc%par 1114 1115 end subroutine copy_incidence_desc
debug_write_incidence
[ Top ] [ m_qtlmap_incidence ] [ Subroutines ]
NAME
debug_write_incidence
DESCRIPTION
Print the contingence matrix and additionnal information. The user have the possibility to print a part of the contingence matrix
INPUTS
incidenceDesc : the description of contingence matrix xinc : the contingence matrix
NOTES
end_analyse_linear
[ Top ] [ m_qtlmap_incidence ] [ Subroutines ]
NAME
end_analyse_linear
DESCRIPTION
close the linear mode
NOTES
SOURCE
454 subroutine end_analyse_linear 455 deallocate (likelyprob) 456 end subroutine end_analyse_linear
end_dataset
[ Top ] [ m_qtlmap_incidence ] [ Subroutines ]
NAME
end_dataset
DESCRIPTION
- release the variable dataset of type DATASET_TYPE
OUTPUTS
dataset : the user datatset (DATASET_TYPE)
NOTES
SOURCE
845 subroutine end_dataset(dataset) 846 type(DATASET_TYPE) , intent(inout) :: dataset 847 integer :: ip,jm,ifem 848 849 do ip=1,np 850 if ( associated(dataset%lSires(ip)%ppt) ) deallocate (dataset%lSires(ip)%ppt) 851 if ( associated(dataset%lSires(ip)%sig0) ) deallocate (dataset%lSires(ip)%sig0) 852 if ( associated(dataset%lSires(ip)%xmu0) ) deallocate (dataset%lSires(ip)%xmu0) 853 ifem=0 854 do jm=nmp(ip)+1,nmp(ip+1) 855 ifem=ifem+1 856 if (dataset%lSires(ip)%full_sib(ifem)%firstkd<0) cycle 857 if ( associated(dataset%lSires(ip)%full_sib(ifem)%ppt) ) & 858 deallocate(dataset%lSires(ip)%full_sib(ifem)%ppt) 859 if ( associated(dataset%lSires(ip)%full_sib(ifem)%pmt) ) & 860 deallocate(dataset%lSires(ip)%full_sib(ifem)%pmt) 861 end do 862 deallocate(dataset%lSires(ip)%full_sib) 863 end do 864 deallocate (dataset%lSires) 865 866 if (associated (dataset%Y)) then 867 deallocate(dataset%Y) 868 deallocate(dataset%CD) 869 else 870 deallocate(dataset%YDISCRETORD) 871 end if 872 873 if (associated (dataset%rhoi) ) deallocate(dataset%rhoi) 874 if (associated (dataset%cov) ) deallocate(dataset%cov) 875 if (associated (dataset%sig) ) deallocate(dataset%sig) 876 if (associated (dataset%xmu) ) deallocate(dataset%xmu) 877 878 879 end subroutine end_dataset
end_incidence
[ Top ] [ m_qtlmap_incidence ] [ Subroutines ]
NAME
end_incidence
DESCRIPTION
release the incidence structure (INCIDENCE_TYPE)
INPUTS
incidenceDesc : the description of the contingence matrix
NOTES
SOURCE
563 subroutine end_incidence(incidenceDesc) 564 type(INCIDENCE_TYPE) , intent(inout) :: incidenceDesc 565 integer :: i,ip,jm,ifem 566 567 call log_mess("end_incidence for incidence...",DEBUG_DEF) 568 !release structure recursive of description of each effect 569 do i=1,incidenceDesc%nteff 570 if ( incidenceDesc%desc(i)%haveSubDesc) then 571 deallocate (incidenceDesc%desc(i)%listSubDesc) 572 end if 573 end do 574 !release all alocated array 575 deallocate(incidenceDesc%desc) 576 deallocate(incidenceDesc%vecsol) 577 deallocate(incidenceDesc%precis) 578 deallocate (incidenceDesc%corr_niv_nivb) 579 580 if ( incidenceDesc%nqtl > 0) then 581 deallocate (incidenceDesc%ntniv_qtlsires) 582 deallocate (incidenceDesc%ntniv_qtldams) 583 end if 584 585 if (associated(incidenceDesc%dataset)) then 586 call end_dataset(incidenceDesc%dataset) 587 deallocate (incidenceDesc%dataset) 588 end if 589 590 if (associated(incidenceDesc%par) ) then 591 deallocate (incidenceDesc%borni) 592 deallocate (incidenceDesc%borns) 593 deallocate (incidenceDesc%par) 594 end if 595 596 if ( associated (incidenceDesc%nonull) ) deallocate(incidenceDesc%nonull) 597 if ( associated (incidenceDesc%nbnonull) ) deallocate(incidenceDesc%nbnonull) 598 599 end subroutine end_incidence
end_position
[ Top ] [ m_qtlmap_incidence ] [ Subroutines ]
NAME
end_position
DESCRIPTION
release position of type POSITION_LRT_INCIDENCE
OUTPUTS
position : the variable poisiton
NOTES
SOURCE
1228 subroutine end_position (position) 1229 type(POSITION_LRT_INCIDENCE) , intent(inout) :: position 1230 1231 deallocate (position%lrt) 1232 deallocate (position%listChr) 1233 deallocate (position%listN) 1234 1235 end subroutine end_position
estim_cholesky
[ Top ] [ m_qtlmap_incidence ] [ Subroutines ]
NAME
estim_cholesky
DESCRIPTION
Estimability with a decomposition of cholesky method
X'.X = L . L*
where L is a lower triangular matrix with strictly positive diagonal entries, and L* denotes the conjugate transpose of L
output : incidenceDesc%vecsol : array of boolean : true if estim , false otherwise
incidenceDesc%nbnivest : number of effect to estimate according the decompostion
INPUTS
XX : X'.X
sizeTriangComp : size of the vector result
OUTPUTS
incidenceDesc%vecsol : a logical array with the estimability of each level incidence matrix incidenceDesc%nbnivest : the number of estimable level triangComp : X' V-1 X
NOTES
triangComp is using in the resolution for the linear context
SOURCE
2834 subroutine estim_cholesky(XX,incidenceDesc,sizeTriangComp,triangComp) 2835 type(INCIDENCE_TYPE) , intent(inout) :: incidenceDesc 2836 real (kind=dp) , dimension(incidenceDesc%ntnivmax,incidenceDesc%ntnivmax), intent(in) :: XX 2837 integer ,intent(in) :: sizeTriangComp 2838 real (kind=dp) , dimension(sizeTriangComp,sizeTriangComp),intent(out) :: triangComp 2839 integer :: i,j,k,ief,jef 2840 real (kind=dp) , dimension(incidenceDesc%ntniv,incidenceDesc%ntniv) :: triang 2841 2842 ! application de la d�composition de choleski pour d�terminer les contraintes 2843 ! 2844 ! Dans le vecteur vecsol les effets � estimer sont indiqu�s � TRUE et ceux 2845 ! qui ne sont pas estimables � FALSE 2846 ! 2847 ! La matrice triang est la matrice triangulaire sup�rieure M telle que 2848 ! M'M = X'X 2849 ! 2850 if ( sizeTriangComp < incidenceDesc%ntniv ) then 2851 call stop_application("Devel error in m_incidence:estim_cholesky") 2852 end if 2853 2854 incidenceDesc%vecsol=.true. 2855 triang=0.d0 2856 triangComp=0.d0 2857 do j=1,incidenceDesc%ntniv 2858 triang(j,j)=xx(j,j) 2859 do k=1,j-1 2860 triang(j,j)=triang(j,j)- triang(j,k)*triang(j,k) 2861 end do 2862 if(triang(j,j).gt. SEUIL_CHO) then 2863 incidenceDesc%nbnivest=incidenceDesc%nbnivest+1 2864 triang(j,j)=dsqrt(triang(j,j)) 2865 do i=j+1,incidenceDesc%ntniv 2866 triang(i,j)=xx(i,j) 2867 do k=1,j-1 2868 triang(i,j)=triang(i,j)-triang(i,k)*triang(j,k) 2869 end do 2870 triang(i,j)=triang(i,j)/triang(j,j) 2871 triang(j,i)=triang(i,j) 2872 end do 2873 else 2874 incidenceDesc%vecsol(j)=.false. 2875 end if 2876 end do 2877 2878 i=0 2879 incidenceDesc%nbnivest=0 2880 do ief=1,incidenceDesc%ntniv 2881 if(incidenceDesc%vecsol(ief))then 2882 incidenceDesc%nbnivest=incidenceDesc%nbnivest+1 2883 i=i+1 2884 j=0 2885 do jef=1,incidenceDesc%ntniv 2886 if (incidenceDesc%vecsol(jef)) then 2887 j=j+1 2888 triangComp(i,j) = triang(ief,jef) 2889 end if 2890 end do 2891 ! corniv(ief)=nbnivest 2892 end if 2893 end do 2894 2895 end subroutine estim_cholesky
fill_nonull_elements
[ Top ] [ m_qtlmap_incidence ] [ Subroutines ]
NAME
fill_nonull_elements
DESCRIPTION
Compacte a sparse matrice.
INPUTS
xincred : the reduction contingence matrix
OUTPUTS
incidenceDesc%nonull : the compacted array with elements not null. incidenceDesc%nbnonull : index in the xincred of no null element by line
NOTES
SOURCE
2625 subroutine fill_nonull_elements(incidenceDesc,s1,s2,xincred) 2626 integer ,intent(in) :: s1,s2 2627 type(INCIDENCE_TYPE) ,intent(inout) :: incidenceDesc 2628 real (kind=dp) ,dimension(s1,s2) ,intent(in) :: xincred 2629 integer :: i,j,k 2630 2631 do i=1,incidenceDesc%dataset%nkd 2632 k=0 2633 do j=1,incidenceDesc%nbnivest 2634 if ( xincred(i,j) /= 0.d0 ) then 2635 k=k+1 2636 incidenceDesc%nonull(i,k)=j 2637 end if 2638 end do 2639 incidenceDesc%nbnonull(i)=k 2640 end do 2641 2642 end subroutine fill_nonull_elements
gen_opti_nqtl
[ Top ] [ m_qtlmap_incidence ] [ Subroutines ]
NAME
gen_opti_nqtl
DESCRIPTION
scan the genome under the hypothesis nqtl
INPUTS
nqtl : the hypothesis
workstruct : residual variance and maximum Likelihood under H(nqtl-1)
FUNCT_MODEL : the function for the resolution (model_lin_h0,model_lin_hn,model_optim_h0,model_optim_hn)
xinc : the contingence matrix
OUTPUTS
curPosMax : maximum position finded
workstruct : residual variance and maximum Likelihood under H(nqtl)
incidenceDescMax : the description of the contingence matrix at the maximum position
sigsquareEstimeMax : the residual variance
BestimMax : the solution
lrtsol : LRT : curves and maximum under H(nqtl)
NOTES
SOURCE
4434 subroutine gen_opti_nqtl( nqtl, & 4435 curPosMax, & 4436 workstruct,& 4437 xincOrig, & 4438 incidenceDescMax, & 4439 sigsquareEstimeMax, & 4440 BestimMax, & 4441 lrtsol, & 4442 FUNCT_MODEL) 4443 integer , intent(in) :: nqtl ! under hypothesis nqtl 4444 type(POSITION_LRT_INCIDENCE) ,intent(inout) :: curPosMax 4445 type(INCIDENCE_GEN_STRUCT) ,intent(inout) :: workstruct 4446 type(INCIDENCE_TYPE) , intent(inout) :: incidenceDescMax ! description of the contingence matrix 4447 real (kind=dp) , dimension(nd,incidenceDescMax%ntnivmax) , target, intent(inout) :: xincOrig ! contingence matrix 4448 real (kind=dp) , dimension(np) , intent(inout) :: sigsquareEstimeMax ! estimation of parameter on the maximum finded 4449 real (kind=dp) , dimension(incidenceDescMax%ntnivmax) , intent(inout) :: BestimMax ! estimation of parameter on the maximum finded 4450 type(TYPE_LRT_SOLUTION) ,intent(inout) :: lrtsol 4451 external :: FUNCT_MODEL 4452 4453 4454 real (kind=dp) , dimension(incidenceDescMax%ntnivmax) :: Bestim 4455 real (kind=dp) , dimension(:,:) , pointer :: xincSub 4456 real (kind=dp) , dimension(np) :: sigsquareEstime 4457 4458 integer :: iip,ifem,ii,chr,ip,iqtl,i,nlinReal,status 4459 logical hypotheseOne, hypotheseTwo,ok 4460 real (kind=dp) ,dimension(:,:),pointer :: xxx 4461 4462 integer :: nGL,nGLTotal,nGLFact,InGL,nlinMax,nlin,nmax(nqtl),chmax(nqtl),nlinValide 4463 integer , dimension(:,:) ,allocatable :: lchr,lnpos 4464 real(kind=dp) , dimension(:,:) ,allocatable :: lrt_save,bestim_save,sigsq_save 4465 type(POSITION_LRT_INCIDENCE) :: curPos 4466 type(INCIDENCE_TYPE) :: incidenceDesc 4467 logical :: loopFalse 4468 4469 hypotheseOne = ( lrtsol%nqtl == 1 ) 4470 hypotheseTwo = ( lrtsol%nqtl == 2 ) 4471 4472 nGL = 0 4473 4474 !Nombre de point sur le groupe de liaison 4475 !on cherche le nombre de point dependant de l hyp nqtl du numbre de point sur le groupe de liason 4476 4477 do chr=1,nchr 4478 nGL=nGL + get_npo(chr) 4479 end do 4480 4481 nGLTotal=1 4482 do iqtl=1,nqtl 4483 nGLTotal=nGLTotal*nGL 4484 end do 4485 4486 curPosMax%listN(nqtl)=0 4487 curPosMax%listN(1:nqtl-1)=1 4488 curPosMax%listChr=1 4489 4490 allocate (lchr(nGLTotal,nqtl)) 4491 allocate (lnpos(nGLTotal,nqtl)) 4492 4493 nlinValide=0 4494 do nlin=1,nGLTotal 4495 i=nqtl 4496 ok=.false. 4497 do while ( (.not. ok) .and. (i>=1)) 4498 curPosMax%listN(i) = curPosMax%listN(i)+1 4499 if ( curPosMax%listN(i) > get_npo(curPosMax%listChr(i)) ) then 4500 if ( nchr > curPosMax%listChr(i) ) then 4501 curPosMax%listN(i)=1 4502 curPosMax%listChr(i)=curPosMax%listChr(i)+1 4503 ok=.true. 4504 else 4505 curPosMax%listN(i)=1 4506 curPosMax%listChr(i)=1 4507 i=i-1 4508 end if 4509 else 4510 ok = .true. 4511 end if 4512 end do 4513 ok=.true. 4514 do iqtl=2,nqtl 4515 ok = ok .and. ( curPosMax%listN(iqtl-1)<=curPosMax%listN(iqtl) ) 4516 end do 4517 4518 if ( ok ) then 4519 nlinValide=nlinValide+1 4520 lchr(nlinValide,:)=curPosMax%listChr 4521 lnpos(nlinValide,:)=curPosMax%listN 4522 end if 4523 end do 4524 nGLTotal=nlinValide 4525 4526 allocate (bestim_save(nGLTotal,incidenceDescMax%ntnivmax)) 4527 allocate (sigsq_save(nGLTotal,np)) 4528 allocate (lrt_save(nGLTotal,nqtl)) 4529 lrt_save=0.d0 4530 4531 !call debug_write_incidence(xinc,incidenceDescMax) 4532 4533 !$OMP PARALLEL DEFAULT(SHARED) & 4534 !$OMP PRIVATE(xincSub,curPos,i,iqtl,ok,sigsquareEstime,Bestim,xxx,ip,iip,ifem,incidenceDesc,loopFalse) 4535 4536 allocate (xincSub(nd,incidenceDescMax%ntnivmax),stat=status) 4537 xincSub=xincOrig 4538 if ( status /= 0 ) call stop_application("Can not allocate incidence matrix ** not enough memory ** ") 4539 call init_position (nqtl,curPos) 4540 allocate (xxx(incidenceDescMax%ntnivmax,incidenceDescMax%ntnivmax),stat=status) 4541 if ( status /= 0 ) call stop_application("Can not allocate xxx array ** not enough memory ** ") 4542 !initialisation of contingence matrix 4543 call copy_incidence_desc (incidenceDescMax,incidenceDesc) 4544 4545 curPos%lrt(nqtl)=0 4546 curPos%listN(nqtl)=0 4547 curPos%listN(1:nqtl-1)=1 4548 curPos%listChr=1 4549 4550 !$OMP DO 4551 do nlin=1,nGLTotal 4552 call log_mess( trim(str((float(nlin)/float(nGLTotal))*100.d0))//"%", DEBUG_DEF ) 4553 curPos%listN=lnpos(nlin,:) 4554 curPos%listChr=lchr(nlin,:) 4555 4556 !traitement specifique pour le LD 4557 if ( workstruct%type_incidence == INCIDENCE_TYPE_LD .or. workstruct%type_incidence == INCIDENCE_TYPE_LDLA) then 4558 i=1 4559 loopFalse = .false. 4560 do while (i<=nqtl .and. .not. loopFalse ) 4561 if (absi(curPos%listChr(i),curPos%listN(i))<posi(curPos%listChr(i),1+longhap/2) ) & 4562 loopFalse = .true. 4563 if (absi(curPos%listChr(i),curPos%listN(i))>posi(curPos%listChr(i),nmk(curPos%listChr(i))-longhap/2)) & 4564 loopFalse = .true. 4565 i = i + 1 4566 end do 4567 4568 if ( loopFalse ) cycle 4569 do i=1,nqtl 4570 call add_haplotype_effect(workstruct%listnteffhap(i),i,curPos%listChr(i),curPos%listN(i),& 4571 xincSub,incidenceDesc,workstruct%hdam,.false.) 4572 end do 4573 end if 4574 4575 if ( workstruct%type_incidence /= INCIDENCE_TYPE_LD ) then 4576 do iqtl=1,nqtl 4577 call change_qtleffect(workstruct%listnteff(iqtl),iqtl,curPos%listChr(iqtl),curPos%listN(iqtl),& 4578 xincSub,incidenceDesc,0,workstruct%linear) 4579 end do 4580 end if 4581 4582 ! call debug_write_incidence(xinc,incidenceDesc) 4583 call FUNCT_MODEL(xincSub,incidenceDesc,curPos,workstruct,sigsquareEstime,Bestim,.false.,.false.,xxx,.false.) 4584 4585 !save value 4586 bestim_save(nlin,:incidenceDesc%nbnivest)=Bestim(:incidenceDesc%nbnivest) 4587 sigsq_save(nlin,:np)=sigsquareEstime(:np) 4588 lrt_save(nlin,:)=curPos%lrt 4589 ! print *,curPos%listN,":",curPos%lrt(nqtl),'-->nlin:',nlin!,OMP_GET_THREAD_NUM() 4590 4591 if ( hypotheseOne ) then 4592 lrtsol%lrt1(curPos%listChr(1),curPos%listN(1))=curPos%lrt(1) 4593 iip=1 4594 do ip=1,np 4595 if ( incidenceDesc%vecsol(1+ip)) then 4596 iip=iip+1 4597 ! ATTENTION VALIDE DANS LE CAS SANS INTERACTION AVEC QTL 4598 lrtsol%pater_eff(curPos%listChr(1),ip,curPos%listN(1))=Bestim(iip) 4599 end if 4600 end do 4601 ifem=0 4602 do ii=1,nm 4603 if ( estime(incidenceDesc%ic,ii)) then 4604 ifem=ifem+1 4605 if ( incidenceDesc%vecsol(np+1+ifem) ) then 4606 iip=iip+1 4607 ! ATTENTION VALIDE DANS LE CAS SANS INTERACTION AVEC QTL 4608 lrtsol%mater_eff(curPos%listChr(1),ifem,curPos%listN(1))=Bestim(iip) 4609 end if 4610 end if 4611 end do 4612 end if 4613 if ( hypotheseTwo ) then 4614 lrtsol%lrt1_2(curPos%listChr(1),curPos%listChr(2),curPos%listN(1),curPos%listN(2))=curPos%lrt(2) 4615 lrtsol%lrt0_2(curPos%listChr(1),curPos%listChr(2),curPos%listN(1),curPos%listN(2))=curPos%lrt(1) 4616 end if 4617 ! end if !ok 4618 end do ! nlin 4619 !$OMP END DO 4620 call end_position (curPos) 4621 deallocate (xxx) 4622 call release_copy_incidence_desc(incidenceDesc) 4623 deallocate (xincSub) 4624 !$OMP END PARALLEL 4625 4626 nlinMax=1 4627 lrtsol%lrtmax(0)=lrt_save(nlinMax,nqtl) 4628 nmax=1 4629 chmax=1 4630 4631 call init_position (nqtl,curPos) 4632 4633 curPos%listN(nqtl)=0 4634 curPos%listN(1:nqtl-1)=1 4635 curPos%listChr=1 4636 4637 do nlin=1,nGLTotal 4638 curPos%listN=lnpos(nlin,:) 4639 curPos%listChr=lchr(nlin,:) 4640 4641 if(lrtsol%lrtmax(nqtl-1) < lrt_save(nlin,nqtl)) then 4642 nmax=curPos%listN 4643 chmax=curPos%listChr 4644 nlinMax=nlin 4645 lrtsol%lrtmax=lrt_save(nlin,:) 4646 end if 4647 end do 4648 4649 do iqtl=1,nqtl 4650 lrtsol%nxmax(iqtl-1)=nmax(iqtl) 4651 lrtsol%chrmax(iqtl-1)=chmax(iqtl) 4652 BestimMax = bestim_save(nlinMax,:) 4653 sigsquareEstimeMax=sigsq_save(nlinMax,:) 4654 curPosMax%listN=curPos%listN 4655 curPosMax%listChr=curPos%listChr 4656 end do 4657 4658 call end_position (curPos) 4659 4660 deallocate (bestim_save) 4661 deallocate (sigsq_save) 4662 deallocate (lrt_save) 4663 deallocate (lchr) 4664 deallocate (lnpos) 4665 4666 end subroutine gen_opti_nqtl
gen_opti_nqtl_cuda
[ Top ] [ m_qtlmap_incidence ] [ Subroutines ]
NAME
gen_opti_nqtl
DESCRIPTION
scan the genome under the hypothesis nqtl
INPUTS
nqtl : the hypothesis
workstruct : residual variance and maximum Likelihood under H(nqtl-1)
FUNCT_MODEL : the function for the resolution (model_lin_h0,model_lin_hn,model_optim_h0,model_optim_hn)
xinc : the contingence matrix
OUTPUTS
curPosMax : maximum position finded
workstruct : residual variance and maximum Likelihood under H(nqtl)
incidenceDescMax : the description of the contingence matrix at the maximum position
sigsquareEstimeMax : the residual variance
BestimMax : the solution
lrtsol : LRT : curves and maximum under H(nqtl)
NOTES
SOURCE
3941 subroutine gen_opti_nqtl_cuda( nsim, & 3942 nqtl, & 3943 curPosMax, & 3944 workstruct,& 3945 xinc, & 3946 incidenceDescMax, & 3947 sigsquareEstimeMax, & 3948 BestimMax, & 3949 lrtsol,& 3950 ySIMUL,& 3951 sizeF, & 3952 AllsigQtlMoinUn,& 3953 listLrtSolSim) 3954 integer , intent(in) :: nsim,nqtl ! under hypothesis nqtl 3955 type(POSITION_LRT_INCIDENCE) ,intent(inout) :: curPosMax 3956 type(INCIDENCE_GEN_STRUCT) ,target ,intent(inout) :: workstruct 3957 type(INCIDENCE_TYPE) , intent(inout) ,target :: incidenceDescMax ! description of the contingence matrix 3958 real (kind=dp) , dimension(nd,incidenceDescMax%ntnivmax) ,target, intent(inout) :: xinc ! contingence matrix 3959 real (kind=dp) , dimension(np) , intent(inout) :: sigsquareEstimeMax ! estimation of parameter on the maximum finded 3960 real (kind=dp) , dimension(incidenceDescMax%ntnivmax) , intent(inout) :: BestimMax ! estimation of parameter on the maximum finded 3961 type(TYPE_LRT_SOLUTION) ,intent(inout) :: lrtsol 3962 type(TYPE_LRT_SOLUTION) , intent(inout) ,dimension(nsim,ncar,nqtl) :: listLrtSolSim 3963 real (kind=dp), dimension (nsim+1,nd) , intent(in) :: ySIMUL 3964 integer ,dimension(np) , intent(in) :: sizeF 3965 real (kind=dp), dimension (nsim+1,np) , intent(inout) :: AllsigQtlMoinUn 3966 3967 3968 real (kind=dp) , dimension(incidenceDescMax%ntnivmax) :: Bestim 3969 real (kind=dp) , dimension(np) :: sigsquareEstime 3970 3971 integer :: iip,ifem,ii,chr,ip,iqtl,i,nlinReal,nsim2,ieff,stat_ok,j 3972 logical hypotheseOne, hypotheseTwo,ok 3973 real (kind=dp) ,dimension(:,:),pointer :: xxx 3974 3975 integer :: nGL,nGLTotal,nGLFact,InGL,nlinMax,nlin,nmax(nqtl),chmax(nqtl),nlinValide 3976 type(WORK_CUDA) :: w_cuda 3977 real(kind=dp) , dimension(:,:) ,allocatable :: lrt_save,bestim_save,sigsq_save 3978 type(POSITION_LRT_INCIDENCE) :: curPos 3979 integer , dimension(:),allocatable ,target :: corrLevel 3980 3981 integer :: isim,c,do,u,nteffS,nteffE,nkd,jm,kd,ikd 3982 integer :: n,kd1,kd2,checkIndice,mode 3983 integer ,dimension(:) ,pointer :: allmaxpos 3984 real (kind=dp) :: seuil_cholesky 3985 3986 #ifdef CUDA_SP 3987 real ,dimension(:,:,:) ,pointer :: allLrt 3988 real ,dimension(:,:) ,pointer :: allmaxlrt 3989 real ,dimension(:,:,:) ,pointer :: allbestim,allosig 3990 3991 real ,dimension (:,:) ,pointer :: AllsigQtlMoinUnR 3992 real ,dimension (:,:) ,pointer :: ySIMULR 3993 real ,dimension (:) ,pointer :: CDR 3994 #else 3995 real (kind=dp) ,dimension(:,:,:) ,pointer :: allLrt 3996 real (kind=dp) ,dimension(:,:) ,pointer :: allmaxlrt 3997 real (kind=dp) ,dimension(:,:,:) ,pointer :: allbestim,allosig 3998 #endif 3999 4000 4001 !******************************************************************* 4002 4003 4004 hypotheseOne = ( lrtsol%nqtl == 1 ) 4005 hypotheseTwo = ( lrtsol%nqtl == 2 ) 4006 4007 nGL = 0 4008 4009 !Nombre de point sur le groupe de liaison 4010 !on cherche le nombre de point dependant de l hyp nqtl du numbre de point sur le groupe de liason 4011 4012 do chr=1,nchr 4013 nGL=nGL + get_npo(chr) 4014 end do 4015 4016 nGLTotal=1 4017 do iqtl=1,nqtl 4018 nGLTotal=nGLTotal*nGL 4019 end do 4020 4021 curPosMax%listN(nqtl)=0 4022 curPosMax%listN(1:nqtl-1)=1 4023 curPosMax%listChr=1 4024 4025 ! structure pour transporter l info du code C cuda a Fortran 4026 allocate (w_cuda%lchr(nGLTotal,nqtl)) 4027 allocate (w_cuda%lnpos(nGLTotal,nqtl)) 4028 w_cuda%workstruct=>workstruct 4029 w_cuda%xinc=>xinc 4030 w_cuda%incidenceDesc=>incidenceDescMax 4031 4032 4033 nlinValide=0 4034 do nlin=1,nGLTotal 4035 i=nqtl 4036 ok=.false. 4037 do while ( (.not. ok) .and. (i>=1)) 4038 curPosMax%listN(i) = curPosMax%listN(i)+1 4039 if ( curPosMax%listN(i) > get_npo(curPosMax%listChr(i)) ) then 4040 if ( nchr > curPosMax%listChr(i) ) then 4041 curPosMax%listN(i)=1 4042 curPosMax%listChr(i)=curPosMax%listChr(i)+1 4043 ok=.true. 4044 else 4045 curPosMax%listN(i)=1 4046 curPosMax%listChr(i)=1 4047 i=i-1 4048 end if 4049 else 4050 ok = .true. 4051 end if 4052 end do 4053 ok=.true. 4054 4055 do iqtl=2,nqtl 4056 ok = ok .and. ( curPosMax%listN(iqtl-1)<=curPosMax%listN(iqtl) ) 4057 end do 4058 4059 if ( workstruct%type_incidence == INCIDENCE_TYPE_LD .or. workstruct%type_incidence == INCIDENCE_TYPE_LDLA ) then 4060 ok = ok .and. (& 4061 .not. (absi(curPosMax%listChr(1),curPosMax%listN(1))<& 4062 posi(curPosMax%listChr(1),1+longhap/2)) ) & 4063 .and. & 4064 .not. (absi(curPosMax%listChr(nqtl),curPosMax%listN(nqtl))>& 4065 posi(curPosMax%listChr(nqtl),nmk(curPosMax%listChr(nqtl))-longhap/2)) 4066 end if 4067 4068 if ( ok ) then 4069 nlinValide=nlinValide+1 4070 w_cuda%lchr(nlinValide,:)=curPosMax%listChr 4071 w_cuda%lnpos(nlinValide,:)=curPosMax%listN 4072 end if 4073 end do 4074 4075 nGLTotal=nlinValide 4076 4077 4078 nsim2 = nsim + 1 4079 4080 ! allocate (allvecsol(nGLTotal,incidenceDesc%ntniv)) 4081 allocate (allbestim(nGLTotal,nsim2,incidenceDescMax%ntniv),stat=stat_ok) 4082 if ( stat_ok /= 0 ) then 4083 call stop_application("** 1/6 -> Not enough host memory ** number of position to test:"//trim(str(nGLTotal))) 4084 endif 4085 allocate (allosig(nGLTotal,nsim2,np),stat=stat_ok) 4086 if ( stat_ok /= 0 ) then 4087 call stop_application("** 2/6 -> Not enough host memory ** number of position to test:"//trim(str(nGLTotal))) 4088 endif 4089 allocate (allLrt(nGLTotal,nsim2,nqtl),stat=stat_ok) 4090 if ( stat_ok /= 0 ) then 4091 call stop_application("** 3/6 -> Not enough host memory ** number of position to test:"//trim(str(nGLTotal))) 4092 endif 4093 4094 allocate (allmaxlrt(nsim2,nqtl),stat=stat_ok) 4095 if ( stat_ok /= 0 ) then 4096 call stop_application("** 4/6 -> Not enough host memory ** number of position to test:"//trim(str(nGLTotal))) 4097 endif 4098 4099 allocate (allmaxpos(nsim2),stat=stat_ok) 4100 if ( stat_ok /= 0 ) then 4101 call stop_application("** 5/6 -> Not enough host memory ** number of position to test:"//trim(str(nGLTotal))) 4102 endif 4103 allocate (corrLevel(incidenceDescMax%ntniv),stat=stat_ok) 4104 if ( stat_ok /= 0 ) then 4105 call stop_application("** 6/6 -> Not enough host memory ** number of position to test:"//trim(str(nGLTotal))) 4106 endif 4107 4108 w_cuda%corrLevel=>corrLevel 4109 4110 u=1 4111 w_cuda%nLevelFix = 0 4112 !les effets fixe a la position sont en premier 4113 do ieff=1,incidenceDescMax%nteff 4114 if ( incidenceDescMax%desc(ieff)%isVar ) cycle 4115 ! print *,"fix:",ieff,incidenceDescMax%desc(ieff)%start,incidenceDescMax%desc(ieff)%end 4116 do i=incidenceDescMax%desc(ieff)%start,incidenceDescMax%desc(ieff)%end 4117 corrLevel(u)=i-1 4118 u=u+1 4119 ! print *,i,"->",corrLevel(i) 4120 w_cuda%nLevelFix = w_cuda%nLevelFix + 1 4121 end do 4122 end do 4123 4124 w_cuda%nLevelVar = 0 4125 !puis les effets variables 4126 do ieff=1,incidenceDescMax%nteff 4127 if ( .not. incidenceDescMax%desc(ieff)%isVar ) cycle 4128 ! print *,"VAR:",ieff,incidenceDescMax%desc(ieff)%start,incidenceDescMax%desc(ieff)%end 4129 do i=incidenceDescMax%desc(ieff)%start,incidenceDescMax%desc(ieff)%end 4130 corrLevel(u)=i-1 4131 u=u+1 4132 w_cuda%nLevelVar = w_cuda%nLevelVar + 1 4133 ! print *,i,"->",corrLevel(i) 4134 end do 4135 end do 4136 4137 call log_mess("call CUDA analysis.....",INFO_DEF) 4138 4139 mode = 0 4140 if (workstruct%type_model==LINEAR_HETEROSCEDASTIC) mode=1 4141 4142 #ifdef CUDA_SP 4143 4144 allocate (AllsigQtlMoinUnR(nsim2,np)) 4145 AllsigQtlMoinUnR = AllsigQtlMoinUn 4146 allocate (ySIMULR(nsim2,incidenceDescMax%dataset%nkd)) 4147 allocate (CDR(incidenceDescMax%dataset%nkd)) 4148 4149 CDR=incidenceDescMax%dataset%CD(1,:) 4150 ySIMULR=ySIMUL(:,:incidenceDescMax%dataset%nkd) 4151 seuil_cholesky=0.3 4152 4153 call log_mess("** simple precision ** ",INFO_DEF) 4154 call cuda_homoscedatic_resolv_genome(mode, & 4155 nqtl, & 4156 AllsigQtlMoinUnR, & 4157 w_cuda%nLevelFix, & 4158 w_cuda%nLevelVar, & 4159 corrLevel, & 4160 w_cuda, & 4161 get_xinc_fixed_effect, & 4162 get_partial_xinc, & 4163 ySIMULR, & 4164 CDR, & 4165 nsim2, & 4166 nd, & 4167 incidenceDescMax%dataset%nkd, & 4168 np, & 4169 sizeF, & 4170 incidenceDescMax%ntnivmax, & 4171 nGLTotal, & 4172 seuil_cholesky, & 4173 allbestim, & 4174 allosig, & 4175 allLrt, & 4176 allmaxlrt, & 4177 allmaxpos ) 4178 4179 AllsigQtlMoinUn = AllsigQtlMoinUnR 4180 deallocate (AllsigQtlMoinUnR) 4181 deallocate (ySIMULR) 4182 deallocate (CDR) 4183 4184 #else 4185 seuil_cholesky=SEUIL_CHO 4186 call log_mess("** double precision ** ",INFO_DEF) 4187 call cuda_homoscedatic_resolv_genome(mode, & 4188 nqtl, & 4189 AllsigQtlMoinUn, & 4190 w_cuda%nLevelFix, & 4191 w_cuda%nLevelVar, & 4192 corrLevel, & 4193 w_cuda, & 4194 get_xinc_fixed_effect, & 4195 get_partial_xinc, & 4196 ySIMUL, & 4197 incidenceDescMax%dataset%CD, & 4198 nsim2, & 4199 nd, & 4200 incidenceDescMax%dataset%nkd, & 4201 np, & 4202 sizeF, & 4203 incidenceDescMax%ntnivmax, & 4204 nGLTotal, & 4205 seuil_cholesky, & 4206 allbestim, & 4207 allosig, & 4208 allLrt, & 4209 allmaxlrt, & 4210 allmaxpos ) 4211 4212 4213 #endif 4214 4215 ! print *,'Fortran ALLOSIG ' 4216 ! print *,allosig(1,1,:) 4217 ! 4218 ! print *,'Fortran ALLOBEST ' 4219 ! print *,allbestim(1,1,:) 4220 ! 4221 ! print *,'Fortran MAXLRT ' 4222 ! print *,allmaxlrt(:10) 4223 ! print *,'Fortran MAXPOS ' 4224 ! print *,allmaxpos(:10) 4225 4226 isim=1 4227 4228 if ( workstruct%nqtl == 1 ) then 4229 do nlin=1,nGLTotal 4230 lrtsol%lrt1(w_cuda%lchr(nlin,1),w_cuda%lnpos(nlin,1))=allLrt(nlin,isim,1) 4231 end do 4232 end if 4233 if ( workstruct%nqtl == 2 ) then 4234 lrtsol%lrt0_2 = 0.d0 ! pas implemente.... 4235 do nlin=1,nGLTotal 4236 lrtsol%lrt1_2(w_cuda%lchr(nlin,1),w_cuda%lchr(nlin,2),w_cuda%lnpos(nlin,1),w_cuda%lnpos(nlin,2))=allLrt(nlin,isim,2) 4237 end do 4238 end if 4239 4240 if ( allmaxpos(isim) > 0 ) then 4241 sigsquareEstimeMax=allosig(allmaxpos(isim),isim,:) 4242 incidenceDescMax%vecsol=.false. 4243 4244 j=1 4245 do i=1,incidenceDescMax%ntniv 4246 if ( allbestim(allmaxpos(isim),isim,i) /= 0.d0 ) then 4247 BestimMax(j) = allbestim(allmaxpos(isim),isim,i) 4248 incidenceDescMax%vecsol(i)=.true. 4249 j=j+1 4250 end if 4251 end do 4252 4253 AllsigQtlMoinUn(isim,:) = allosig(allmaxpos(isim),isim,:) 4254 4255 do iqtl=1,nqtl 4256 lrtsol%nxmax(iqtl-1)=w_cuda%lnpos(allmaxpos(isim),iqtl) 4257 lrtsol%chrmax(iqtl-1)=w_cuda%lchr(allmaxpos(isim),iqtl) 4258 lrtsol%lrtmax(iqtl-1)=allmaxlrt(isim,iqtl) 4259 curPosMax%listN=w_cuda%lnpos(allmaxpos(isim),iqtl) 4260 curPosMax%listChr=w_cuda%lchr(allmaxpos(isim),iqtl) 4261 end do 4262 ! print *,lrtsol%lrtmax 4263 end if 4264 4265 do isim=1,nsim 4266 if ( allmaxpos(isim+1) <= 0 .or. allmaxpos(isim+1) > nGLTotal) cycle 4267 AllsigQtlMoinUn(isim+1,:) = allosig(allmaxpos(isim+1),isim+1,:) 4268 allocate (listLrtSolSim(isim,incidenceDescMax%ic,nqtl)%lrtmax(0:nqtl-1)) 4269 allocate (listLrtSolSim(isim,incidenceDescMax%ic,nqtl)%nxmax(0:nqtl-1)) 4270 allocate (listLrtSolSim(isim,incidenceDescMax%ic,nqtl)%chrmax(0:nqtl-1)) 4271 4272 ! listLrtSolSim(isim,incidenceDescMax%ic,nqtl)%lrtmax=allmaxlrt(isim+1,nqtl) 4273 4274 do iqtl=1,nqtl 4275 listLrtSolSim(isim,incidenceDescMax%ic,nqtl)%nxmax(iqtl-1)=w_cuda%lnpos(allmaxpos(isim+1),iqtl) 4276 listLrtSolSim(isim,incidenceDescMax%ic,nqtl)%chrmax(iqtl-1)=w_cuda%lchr(allmaxpos(isim+1),iqtl) 4277 listLrtSolSim(isim,incidenceDescMax%ic,nqtl)%lrtmax(iqtl-1)=allmaxlrt(isim,iqtl) 4278 end do 4279 end do 4280 4281 deallocate (corrLevel) 4282 ! deallocate (xinc2) 4283 ! deallocate (allvecsol) 4284 deallocate (allbestim) 4285 deallocate (allosig) 4286 deallocate (w_cuda%lchr) 4287 deallocate (w_cuda%lnpos) 4288 deallocate (allLrt) 4289 deallocate (allmaxlrt) 4290 deallocate (allmaxpos) 4291 4292 end subroutine gen_opti_nqtl_cuda
get_correlations
[ Top ] [ m_qtlmap_incidence ] [ Subroutines ]
NAME
get_correlations
DESCRIPTION
Compute the correlation between level of the incidence matrix.
INPUTS
incidenceDesc : the description of the incidence matrix.
xxx : buffer array provide by get_precision function
OUTPUTS
xcor : the correlation
NOTES
xxx come from the get_precision function.
SOURCE
2374 subroutine get_correlations(incidenceDesc,xxx,xcor) 2375 type(INCIDENCE_TYPE) ,intent(in) :: incidenceDesc 2376 real (kind=dp) ,dimension(incidenceDesc%ntnivmax,incidenceDesc%ntnivmax) ,intent(in) :: xxx 2377 real (kind=dp) ,dimension(incidenceDesc%ntnivmax,incidenceDesc%ntnivmax) ,intent(out) :: xcor 2378 integer ::ief,jef,i,j 2379 2380 xcor=0.d0 2381 ief=0 2382 do i=1,incidenceDesc%ntniv 2383 if(incidenceDesc%vecsol(i))then 2384 ief=ief+1 2385 jef=0 2386 do j=1,incidenceDesc%ntniv 2387 if(incidenceDesc%vecsol(j))then 2388 jef=jef+1 2389 if (xxx(ief,ief) == 0 .or. xxx(jef,jef) == 0) then 2390 xcor(i,j)=0.d0 2391 else 2392 xcor(i,j)=xxx(ief,jef)/dsqrt(xxx(ief,ief)*xxx(jef,jef)) 2393 end if 2394 end if 2395 end do 2396 end if 2397 end do 2398 end subroutine get_correlations
get_incidence_reduit
[ Top ] [ m_qtlmap_incidence ] [ Subroutines ]
NAME
get_incidence_reduit
DESCRIPTION
INPUTS
incidenceDesc : the description of the incidence matrix.
xinc : the contingence matrix
OUTPUTS
xincred : the incontingence matrixithout level which are not estimable (the reduction contingence matrix)
NOTES
SOURCE
2595 subroutine get_incidence_reduit(incidenceDesc,xinc,xincred) 2596 type(INCIDENCE_TYPE) ,intent(in) :: incidenceDesc 2597 real (kind=dp) ,dimension(nd,incidenceDesc%ntnivmax) ,intent(in) :: xinc 2598 real (kind=dp) ,dimension(incidenceDesc%dataset%nkd,incidenceDesc%nbnivest) ,intent(inout) :: xincred 2599 2600 integer :: i,j 2601 j=0 2602 do i=1,incidenceDesc%ntniv 2603 if ( incidenceDesc%vecsol(i) ) then 2604 j=j+1 2605 xincred(:,j)=xinc(:incidenceDesc%dataset%nkd,i) 2606 end if 2607 end do 2608 2609 end subroutine get_incidence_reduit
get_precision
[ Top ] [ m_qtlmap_incidence ] [ Subroutines ]
NAME
get_precision
DESCRIPTION
Compute the precision for each level of the incidence matrix.
INPUTS
xx : X'.X ,X contingence matrix xxx : buffer array jm : index of dam nd1 : first index of progenies of family (ip,jm) nd2 : last index of progenies of family (ip,jm) chr : chromosome localisation n : the position
OUTPUTS
incidenceDesc%precis : Precision of each level of the incidence matrix.
NOTES
ginv1
SOURCE
2315 subroutine get_precision(xx,xxx,incidenceDesc) 2316 type(INCIDENCE_TYPE) ,intent(inout) :: incidenceDesc 2317 real (kind=dp) ,dimension(incidenceDesc%ntnivmax,incidenceDesc%ntnivmax) , intent(in) :: xx 2318 real (kind=dp) ,dimension(incidenceDesc%ntnivmax,incidenceDesc%ntnivmax) , intent(out) :: xxx 2319 2320 integer :: iniv,kniv,jniv,lniv,irank 2321 real (kind=dp) :: ups 2322 2323 ! 2324 ! Reduction de la matrice d'incidence 2325 ! 2326 xxx(:incidenceDesc%ntniv,:incidenceDesc%ntniv)=0.d0 2327 2328 kniv=0 2329 do iniv =1,incidenceDesc%ntniv 2330 if(incidenceDesc%vecsol(iniv)) then 2331 kniv=kniv+1 2332 lniv=0 2333 do jniv=1,incidenceDesc%ntniv 2334 if(incidenceDesc%vecsol(jniv)) then 2335 lniv=lniv+1 2336 xxx(kniv,lniv)=xx(iniv,jniv) 2337 end if 2338 end do 2339 end if 2340 end do 2341 ! 2342 ! inversion de la matrice d'incidence reduite pour calcul des precisions 2343 ! 2344 2345 ups=1.do-15 2346 call ginv1(xxx,lniv,size(xxx,1),ups,irank) 2347 2348 jniv=0 2349 do iniv=1,incidenceDesc%ntniv 2350 incidenceDesc%precis(iniv)=9999.d0 2351 if(incidenceDesc%vecsol(iniv)) then 2352 jniv=jniv+1 2353 incidenceDesc%precis(iniv)=xxx(jniv,jniv) 2354 end if 2355 end do 2356 2357 end subroutine get_precision
getglobalprob
[ Top ] [ m_qtlmap_incidence ] [ Subroutines ]
NAME
getglobalprob
DESCRIPTION
Get the probabilities to receive a qtleffect in a general context
INPUTS
ic : index of the trait
ip : index of sire
jm : index of dam
nd1 : first index of progenies of family (ip,jm)
nd2 : last index of progenies of family (ip,jm)
chr : chromosome localisation
n : the position
OUTPUTS
probp :
probm :
NOTES
SOURCE
2263 subroutine getglobalprob(chr,n,ic,ip,jm,nd1,nd2,probp,probm) 2264 integer ,intent(in) :: chr,n,ic,ip,jm,nd1,nd2 2265 real (kind=dp) ,dimension(nd) ,intent(out) :: probp,probm 2266 2267 integer :: kd,ig,kkd 2268 real (kind=dp) :: ppt,pmt 2269 probp=0.d0 2270 probm=0.d0 2271 do kkd=nd1,nd2 2272 do ig=ngenom(chr,jm)+1,ngenom(chr,jm+1) 2273 do kd=ngend(chr,ig)+1,ngend(chr,ig+1) 2274 if ( kkd == ndesc(chr,kd)) then 2275 if(presentc(ic,ndesc(chr,kkd))) then 2276 if( estime(ic,jm) )then 2277 ppt=-pdd(chr,kd,1,n)-pdd(chr,kd,2,n)+pdd(chr,kd,3,n)+pdd(chr,kd,4,n) 2278 probp(kkd)=probp(kkd)+probg(chr,ig)*ppt 2279 pmt=-pdd(chr,kd,1,n)+pdd(chr,kd,2,n)-pdd(chr,kd,3,n)+pdd(chr,kd,4,n) 2280 probm(kkd)=probm(kkd)+probg(chr,ig)*pmt 2281 else 2282 ppt=-pdd(chr,kd,1,n)+pdd(chr,kd,3,n) 2283 probp(kkd)=ppt 2284 end if 2285 exit 2286 end if 2287 end if 2288 end do !! kd 2289 end do !! ig 2290 end do !! kkd 2291 2292 end subroutine getglobalprob
getprob_lin
[ Top ] [ m_qtlmap_incidence ] [ Subroutines ]
NAME
getprob_lin
DESCRIPTION
Get the probabilities to receive qtl in a linear context depending the probabilities of transmission
INPUTS
ic : index of the trait
ip : index of sire
jm : index of dam
nd1 : first index of progenies of family (ip,jm)
nd2 : last index of progenies of family (ip,jm)
chr : chromosome localisation
n : the position
OUTPUTS
probp : The probabilities to receive the qtl from the sire
probm : The probabilities to receive the qtl from the dam
NOTES
SOURCE
2159 subroutine getprob_lin(chr,n,ic,ip,jm,nd1,nd2,probp,probm) 2160 integer ,intent(in) :: chr,n,ic,ip,jm,nd1,nd2 2161 real (kind=dp) ,dimension(nd) ,intent(out) :: probp,probm 2162 2163 integer :: kd,ig,kkd 2164 real (kind=dp) :: p 2165 2166 probp=0.d0 2167 probm=0.d0 2168 ig=likelyprob(chr,ip,jm) 2169 do kkd=nd1,nd2 2170 do kd=ngend(chr,ig)+1,ngend(chr,ig+1) 2171 if ( kkd == ndesc(chr,kd)) then 2172 if(presentc(ic,ndesc(chr,kkd))) then 2173 if( estime(ic,jm) )then 2174 probp(kkd)=-pdd(chr,kd,1,n)-pdd(chr,kd,2,n)+pdd(chr,kd,3,n)+pdd(chr,kd,4,n) 2175 probm(kkd)=-pdd(chr,kd,1,n)+pdd(chr,kd,2,n)-pdd(chr,kd,3,n)+pdd(chr,kd,4,n) 2176 exit 2177 else 2178 probp(kkd)=-pdd(chr,kd,1,n)+pdd(chr,kd,3,n) 2179 exit 2180 end if 2181 end if 2182 end if 2183 end do 2184 end do 2185 2186 end subroutine getprob_lin
getprobld_lin
[ Top ] [ m_qtlmap_incidence ] [ Subroutines ]
NAME
getprobld_lin
DESCRIPTION
Get the probabilities to receive an haplotype
INPUTS
ic : index of the trait
ip : index of sire
jm : index of dam
nd1 : first index of progenies of family (ip,jm)
nd2 : last index of progenies of family (ip,jm)
chr : chromosome localisation
n : the position
OUTPUTS
pp_ldla : pm_ldla :
NOTES
SOURCE
2209 !attention pp_ldla,pm_ldla de taille nd car lineaire (1 seul genotype pour toute les meres) 2210 subroutine getprobld_lin(chr,n,ic,ip,jm,nd1,nd2,pp_ldla,pm_ldla) 2211 integer ,intent(in) :: chr,n,ic,ip,jm,nd1,nd2 2212 real (kind=dp) ,dimension(nd,2) ,intent(out) :: pp_ldla,pm_ldla 2213 2214 integer :: kd,ig,kkd 2215 real (kind=dp) :: p 2216 2217 pp_ldla=0.d0 2218 pm_ldla=0.d0 2219 ig=likelyprob(chr,ip,jm) 2220 do kkd=nd1,nd2 2221 do kd=ngend(chr,ig)+1,ngend(chr,ig+1) 2222 if ( kkd == ndesc(chr,kd)) then 2223 if(presentc(ic,kkd)) then 2224 if( estime(ic,jm) )then 2225 pp_ldla(kkd,1)=pp_ldla(kkd,1)+probg(chr,ig)*(pdd(chr,kd,1,n)+pdd(chr,kd,2,n)) 2226 pp_ldla(kkd,2)=pp_ldla(kkd,2)+probg(chr,ig)*(pdd(chr,kd,3,n)+pdd(chr,kd,4,n)) 2227 pm_ldla(kkd,1)=probg(chr,ig)*(pdd(chr,kd,1,n)+pdd(chr,kd,3,n)) 2228 pm_ldla(kkd,2)=probg(chr,ig)*(pdd(chr,kd,2,n)+pdd(chr,kd,4,n)) 2229 exit 2230 else 2231 pp_ldla(kkd,1)=pdd(chr,kd,1,n) 2232 pp_ldla(kkd,2)=pdd(chr,kd,3,n) 2233 exit 2234 end if 2235 end if 2236 end if 2237 end do 2238 end do 2239 2240 end subroutine getprobld_lin
init_analyse_linear
[ Top ] [ m_qtlmap_incidence ] [ Subroutines ]
NAME
init_analyse_linear
DESCRIPTION
initialize the linear mode
* likelyprob is building from probg array
NOTES
SOURCE
413 subroutine init_analyse_linear() 414 integer :: chr,ip,jm,ngeno1,ngeno2,igprob,ig,ic 415 416 allocate (likelyprob(nchr,np,nm)) 417 418 likelyprob = 0 419 do chr=1,nchr 420 do ip=1,np 421 do jm=nmp(ip)+1,nmp(ip+1) 422 ngeno1=ngenom(chr,jm)+1 423 ngeno2=ngenom(chr,jm+1) 424 igprob=ngeno1 425 426 do ig=ngeno1+1,ngeno2 427 if (probg(chr,igprob)<probg(chr,ig))then 428 igprob = ig 429 end if 430 end do 431 likelyprob(chr,ip,jm)=igprob 432 if ( probg(chr,igprob) <1.d0 ) then 433 call log_mess("*************************************",WARNING_DEF) 434 call log_mess("CAREFULLY : Linear model use the most likely genotype",WARNING_DEF) 435 call log_mess("Chromosome: "//trim(ch(chr)),WARNING_DEF) 436 call log_mess("Sire Family ["//trim(pere(ip))//"]",WARNING_DEF) 437 call log_mess("Genotype prob for dam ["//trim(mere(jm))//"] :"//str(probg(chr,igprob)),WARNING_DEF) 438 call log_mess("*************************************",WARNING_DEF) 439 end if 440 end do !jm 441 end do ! ip 442 end do ! chr 443 end subroutine init_analyse_linear
init_dataset_1car
[ Top ] [ m_qtlmap_incidence ] [ Subroutines ]
NAME
init_dataset_1car
DESCRIPTION
* build information family without unkwnown trait (trait and data censured/ discrete trait)
(algorithm dot not use presentc and estime arrays)
build an subtype array : lSire of type SireFamily_type
* variance and mean for each family (sig0, xmu0)
* allocation of probabilities of transmission arrays (by sire family and full family)
INPUTS
ic : index of the trait
nqtl : number of qtl (hypothesis)
OUTPUTS
dataset : the user datatset (DATASET_TYPE)
NOTES
SOURCE
698 subroutine init_dataset_1car(ic,nqtl,dataset) 699 type(DATASET_TYPE) , intent(inout) :: dataset 700 integer , intent(in) :: ic,nqtl 701 logical ,dimension(nd) :: pres 702 integer :: ip,jm,kd,ifem,kdd,s1,s2,chr,iq,kkd,effp 703 real(kind=dp) :: somyp 704 logical :: isContinue 705 706 707 call log_mess("build_family for incidence...",DEBUG_DEF) 708 709 allocate(dataset%lSires(np)) 710 kdd=0 711 712 do ip=1,np 713 allocate(dataset%lSires(ip)%full_sib(nmp(ip+1)-nmp(ip))) 714 ifem=0 715 do jm=nmp(ip)+1,nmp(ip+1) 716 ifem=ifem+1 717 kd=ndm(jm)+1 718 isContinue = presentc(ic,kd) 719 do while ( kd<=ndm(jm+1) .and. .not. isContinue ) 720 kd=kd+1 721 if ( kd <= nd ) then 722 isContinue = presentc(ic,kd) 723 end if 724 end do 725 726 if (kd>ndm(jm+1)) cycle 727 kdd=kdd+1 728 ! on a trouve le premier kd valide du pere ip dans la matrice d incidence 729 dataset%lSires(ip)%full_sib(ifem)%firstKd=kdd 730 ! print *,trim(pere(ip)),"-",trim(mere(jm)),ndm(jm)+1,ndm(jm+1),count(presentc(ic,kd+1:ndm(jm+1))) 731 ! print *,kd+1,ndm(jm+1) 732 733 kdd=kdd + count(presentc(ic,kd+1:ndm(jm+1))) 734 dataset%lSires(ip)%full_sib(ifem)%lastKd=kdd 735 call log_mess("Full Sib ["//trim(pere(ip))//"-"//trim(mere(jm))//"] :"//& 736 trim(str(dataset%lSires(ip)%full_sib(ifem)%firstKd))//"->"& 737 //trim(str(dataset%lSires(ip)%full_sib(ifem)%lastKd)),VERBOSE_DEF) 738 s1=0 739 do chr=1,nchr 740 s1=max(s1,ngenom(chr,jm+1)-ngenom(chr,jm)) 741 end do 742 s2=dataset%lSires(ip)%full_sib(ifem)%lastKd-dataset%lSires(ip)%full_sib(ifem)%firstKd+1 743 if ( nqtl >= 1) then 744 allocate(dataset%lSires(ip)%full_sib(ifem)%ppt(nqtl,s1,s2)) 745 allocate(dataset%lSires(ip)%full_sib(ifem)%pmt(nqtl,s1,s2)) 746 end if 747 end do 748 ifem=0 749 do jm=nmp(ip)+1,nmp(ip+1) 750 ifem=ifem+1 751 if (dataset%lSires(ip)%full_sib(ifem)%firstKd<0) cycle 752 dataset%lSires(ip)%half_sib%firstKd= dataset%lSires(ip)%full_sib(ifem)%firstKd 753 exit 754 end do 755 ifem=nmp(ip+1)-nmp(ip)+1 756 do jm=nmp(ip+1),nmp(ip)+1,-1 757 ifem=ifem-1 758 if (dataset%lSires(ip)%full_sib(ifem)%lastKd<0) cycle 759 dataset%lSires(ip)%half_sib%lastKd=dataset%lSires(ip)%full_sib(ifem)%lastKd 760 exit 761 end do 762 763 call log_mess("Half Sib ["//trim(pere(ip))//"] :"//& 764 trim(str(dataset%lSires(ip)%half_sib%firstKd))//"->"//trim(str(dataset%lSires(ip)%half_sib%lastKd)),VERBOSE_DEF) 765 !s1=lSires(ip)%half_sib%lastKd-lSires(ip)%half_sib%firstKd+1 766 if ( nqtl >= 1) then 767 allocate(dataset%lSires(ip)%ppt(nqtl,dataset%lSires(ip)%half_sib%firstKd:dataset%lSires(ip)%half_sib%lastKd)) 768 end if 769 end do ! ip 770 771 isContinue = natureY(ic) /= 'i' 772 !pour ne pas a tester la presence des caracteres on garde ce qui est seulement present 773 if ( isContinue ) then 774 allocate (dataset%Y(1,count(presentc(ic,:)))) 775 allocate (dataset%CD(1,count(presentc(ic,:)))) 776 else 777 allocate (dataset%YDISCRETORD(1,count(presentc(ic,:)))) 778 end if 779 780 dataset%nkd=0 781 dataset%nkd_max_by_fam=0 782 783 do ip=1,np 784 do jm=nmp(ip)+1,nmp(ip+1) 785 dataset%nkd_max_by_fam=max(dataset%nkd_max_by_fam,count(presentc(ic,ndm(jm)+1:ndm(jm+1)))) 786 do kd=ndm(jm)+1,ndm(jm+1) 787 if (presentc(ic,kd)) then 788 dataset%nkd=dataset%nkd+1 789 if ( isContinue ) then 790 dataset%Y(1,dataset%nkd)=y(ic,kd) 791 dataset%CD(1,dataset%nkd)=cd(ic,kd) 792 else 793 dataset%YDISCRETORD(1,dataset%nkd)=ydiscretord(ic,kd) 794 end if 795 end if 796 end do !kd 797 end do ! jm 798 end do ! ip 799 800 do ip=1,np 801 allocate(dataset%lSires(ip)%sig0(1)) 802 allocate(dataset%lSires(ip)%xmu0(1)) 803 effp=0 804 somyp=0 805 do jm=nmp(ip)+1,nmp(ip+1) 806 do kd=ndm(jm)+1,ndm(jm+1) 807 if (presentc(ic,kd)) then 808 effp=effp+1 809 somyp=somyp+y(ic,kd)*cd(ic,kd) 810 end if 811 end do 812 end do 813 814 if (effp == 0.d0) then 815 call stop_application('Father ['//trim(pere(ip))//'] has got no child with trait value ['//trim(carac(ic))//"]") 816 end if 817 !moyenne du caractere dans la famille du pere ip 818 dataset%lSires(ip)%xmu0(1)=somyp/dble(effp) 819 820 !ecart type du caractere dans la famille du pere ip 821 dataset%lSires(ip)%sig0(1)=0.d0 822 do jm=nmp(ip)+1,nmp(ip+1) 823 do kd=ndm(jm)+1,ndm(jm+1) 824 if (presentc(ic,kd)) then 825 dataset%lSires(ip)%sig0(1)=dataset%lSires(ip)%sig0(1)+(y(ic,kd)-& 826 dataset%lSires(ip)%xmu0(1))*(y(ic,kd)-dataset%lSires(ip)%xmu0(1)) 827 end if 828 end do !kd 829 end do !jm 830 dataset%lSires(ip)%sig0(1)=dsqrt(dataset%lSires(ip)%sig0(1)/dble(effp-1)) 831 end do !ip 832 end subroutine init_dataset_1car
init_incidence
[ Top ] [ m_qtlmap_incidence ] [ Subroutines ]
NAME
init_incidence
DESCRIPTION
initilise the incidence structure (INCIDENCE_TYPE)
INPUTS
ic : index of the trait
nqtl : number of hypothesis to test
type_incidence : type of contingence matrix construction (INCIDENCE_TYPE_CLASSIC, INCIDENCE_TYPE_LD, INCIDENCE_TYPE_LDLA, INCIDENCE_TRAIT_COV)
dataset : use dataset (type DATASET_TYPE)
OUTPUTS
incidenceDesc : the description of the contingence matrix
NOTES
see copy_dataset,init_dataset_1car if dataset are not created init_dataset_1car is called otherwise copy_dataset
SOURCE
476 subroutine init_incidence(ic,nqtl,incidenceDesc,type_incidence,dataset) 477 integer , intent(in) :: ic,type_incidence 478 integer , intent(in) :: nqtl 479 type(INCIDENCE_TYPE) , intent(inout) :: incidenceDesc 480 type(DATASET_TYPE) , intent(in),optional :: dataset 481 integer :: npar 482 incidenceDesc%nqtl=nqtl 483 call set_ntnivmax(ic,nqtl,incidenceDesc%ntnivmax,incidenceDesc%nteffmax,incidenceDesc%ntlev,incidenceDesc%nbniv) 484 485 if ( type_incidence==INCIDENCE_INTERACTION ) then 486 !! pour l instant que 2 qtl... 487 incidenceDesc%ntnivmax=incidenceDesc%ntnivmax+np*4+count(estime(ic,:))*16 488 incidenceDesc%nteffmax=incidenceDesc%nteffmax+1 489 end if 490 491 if ( type_incidence==INCIDENCE_TRAIT_COV ) then 492 incidenceDesc%ntnivmax=incidenceDesc%ntnivmax+ncar 493 incidenceDesc%nteffmax=incidenceDesc%nteffmax+1 494 end if 495 496 incidenceDesc%ic = ic 497 incidenceDesc%ntniv = 0 498 incidenceDesc%nbnivest = 0 499 incidenceDesc%nteff = 0 500 501 call log_mess("********************* INIT *****************************",DEBUG_DEF) 502 call log_mess("NTNIVMAX:"//str(incidenceDesc%ntnivmax),DEBUG_DEF) 503 call log_mess("NTEFFMAX:"//str(incidenceDesc%nteffmax),DEBUG_DEF) 504 call log_mess("NTLEVQTL / Reproductor:"//str(incidenceDesc%ntlev),DEBUG_DEF) 505 call log_mess("NBLEV FIXED EFFECT:"//str(incidenceDesc%nbniv),DEBUG_DEF) 506 call log_mess("*********************************************************",DEBUG_DEF) 507 508 allocate(incidenceDesc%desc(incidenceDesc%nteffmax)) 509 allocate(incidenceDesc%vecsol(incidenceDesc%ntnivmax)) 510 incidenceDesc%vecsol=.false. 511 ! allocate(incidenceDesc%corniv(ntnivmax)) 512 allocate(incidenceDesc%precis(incidenceDesc%ntnivmax)) 513 incidenceDesc%precis=0.d0 514 if ( nqtl > 0) then 515 allocate (incidenceDesc%ntniv_qtlsires(nqtl)) 516 allocate (incidenceDesc%ntniv_qtldams(nqtl)) 517 incidenceDesc%ntniv_qtlsires=0 518 incidenceDesc%ntniv_qtldams=0 519 end if 520 allocate (incidenceDesc%corr_niv_nivb(incidenceDesc%ntnivmax)) 521 522 allocate (incidenceDesc%dataset) 523 if ( present(dataset) ) then 524 call copy_dataset(dataset,incidenceDesc%dataset) 525 else 526 call init_dataset_1car(ic,nqtl,incidenceDesc%dataset) 527 end if 528 529 npar = np+incidenceDesc%ntnivmax 530 531 if ( nmod(ic) > 1 ) npar = npar + nmod(ic)-1 532 533 allocate (incidenceDesc%borni(npar)) 534 allocate (incidenceDesc%borns(npar)) 535 allocate (incidenceDesc%par(npar)) 536 537 incidenceDesc%borni(1:np)=SIG_MIN 538 incidenceDesc%borns(1:np)=SIG_MAX 539 incidenceDesc%borni(np+1:)=DEFAULT_PARAM_MIN 540 incidenceDesc%borns(np+1:)=DEFAULT_PARAM_MAX 541 542 if (nmod(ic)>1) then 543 incidenceDesc%borni(np+incidenceDesc%ntnivmax+1:)=SIG_MIN 544 incidenceDesc%borns(np+incidenceDesc%ntnivmax+1:)=DEFAULT_PARAM_MAX 545 end if 546 547 incidenceDesc%par=0.d0 548 549 allocate (incidenceDesc%nonull(incidenceDesc%dataset%nkd,incidenceDesc%ntnivmax)) 550 allocate (incidenceDesc%nbnonull(incidenceDesc%dataset%nkd)) 551 end subroutine init_incidence
init_position
[ Top ] [ m_qtlmap_incidence ] [ Subroutines ]
NAME
init_position
DESCRIPTION
Initialize a position for a multi qtl hypothesis test (POSITION_LRT_INCIDENCE), e.g : * allocation of a LRT array value (comparaison with Hi-1, Hi-2,...) * allocation of the chromosome list (genome position of all qtl) * allocation of a discrete position list (index of absi array for each qtl position)
INPUTS
nqtl : the number of qtl to test
OUTPUTS
position : the variable poisiton
NOTES
SOURCE
1185 subroutine init_position (nqtl,position) 1186 integer , intent(in) :: nqtl 1187 type(POSITION_LRT_INCIDENCE) , intent(inout) :: position 1188 1189 integer :: i,j,n,chr 1190 1191 allocate (position%lrt(nqtl)) 1192 allocate (position%listChr(nqtl)) 1193 allocate (position%listN(nqtl)) 1194 1195 position%lrt=0.d0; 1196 1197 do i=1,nqtl 1198 n=0 1199 j=i 1200 do chr=1,nchr 1201 n=n+get_ilong(chr) 1202 if (i<=n) exit 1203 j=1 1204 end do 1205 1206 if (chr>nchr) call stop_application("Not enough sampling point to detect QTLs") 1207 1208 position%listChr(i)=1 ! attention si pas assez d echantillonage on doit passer au chromosome suivant.... 1209 position%listN(i)=i 1210 do while (absi(1,position%listN(i))<posi(1,1+longhap/2)) 1211 position%listN(i)=position%listN(i)+1 1212 end do 1213 end do 1214 1215 end subroutine init_position
init_workstruct
[ Top ] [ m_qtlmap_incidence ] [ Subroutines ]
NAME
init_workstruct
DESCRIPTION
initialize the workstruct structure (INCIDENCE_GEN_STRUCT) The following information is given for each position on the chromosome : * variance family * residual correlation matrix (multi trait analysis case) * value of the likelihood
INPUTS
type_incidence : type of contingence matrix construction (INCIDENCE_TYPE_CLASSIC, INCIDENCE_TYPE_LD, INCIDENCE_TYPE_LDLA, INCIDENCE_TRAIT_COV) type_resolution : type of resolution (LINEAR_HOMOSCEDASTIC, LINEAR_HETEROSCEDASTIC, OPTIM_HETEROSCEDASTIC ) nqtl : number of qtl (hypothesis) pc : need to compute the confusion step pt : need to compute the precision step multi : boolean for multi trait analysis
OUTPUTS
workstruct : buffer structure to save information about result among the chromosome
NOTES
SOURCE
622 subroutine init_workstruct(workstruct,type_incidence,type_resolution,nqtl,pc,pt,multi) 623 type(INCIDENCE_GEN_STRUCT) ,intent(inout) :: workstruct 624 integer ,intent(in) :: type_incidence,type_resolution,nqtl 625 logical ,intent(in) :: pc,pt,multi 626 627 workstruct%type_model=type_resolution 628 workstruct%type_incidence=type_incidence 629 630 workstruct%linear=.false. 631 if (workstruct%type_model==LINEAR_HOMOSCEDASTIC .or. workstruct%type_model==LINEAR_HETEROSCEDASTIC) then 632 workstruct%linear=.true. 633 end if 634 workstruct%nqtl=nqtl 635 workstruct%performConfusion=pc 636 workstruct%performTestNuis=pt 637 638 if (multi) then 639 allocate (workstruct%sigsquare(nqtl,np,ncar)) 640 allocate (workstruct%rhoi(nqtl,ncar,ncar)) 641 else 642 allocate (workstruct%sigsquare(nqtl,np,1)) 643 end if 644 645 allocate (workstruct%listnteff(nqtl)) 646 if (workstruct%type_incidence==INCIDENCE_TYPE_LD .or. workstruct%type_incidence==INCIDENCE_TYPE_LDLA) then 647 allocate (workstruct%listnteffhap(nqtl)) 648 end if 649 allocate (workstruct%fnqtl(nqtl)) 650 651 end subroutine init_workstruct
matmul_incidence
[ Top ] [ m_qtlmap_incidence ] [ Subroutines ]
NAME
matmul_incidence
DESCRIPTION
multiplication matrix routine using the compacted sparse matrix structure (incidenceDesc%nonull,incidenceDesc%nbnonull)
INPUTS
kd1 : first line index
kd2 : last line index
incidenceDesc%nonull : the compacted array with elements not null.
incidenceDesc%nbnonull : index in the xincred of no null element by line
sizeInXinc1 : size of the first dimension of xinc
sizeInXinc2 : size of the second dimension of xinc
in_X : the vector to multiplied
SizeOutVec : size of the result vector
startRes : first index in the result vector to set result
OUTPUTS
out_vec : the result vector
NOTES
SOURCE
2666 subroutine matmul_incidence(kd1,kd2,incidenceDesc,sizeInXinc1,sizeInXinc2,in_xinc,in_X,SizeOutVec,startRes,out_vec) 2667 integer ,intent(in) :: kd1,kd2,SizeOutVec,startRes,sizeInXinc1,sizeInXinc2 2668 type(INCIDENCE_TYPE) ,intent(in) :: incidenceDesc 2669 real (kind=dp) ,dimension(sizeInXinc1,sizeInXinc2) ,intent(in) :: in_xinc 2670 real (kind=dp) ,dimension(incidenceDesc%nbnivest) ,intent(in) :: in_X 2671 real (kind=dp) ,dimension(SizeOutVec) ,intent(out) :: out_vec 2672 2673 integer :: i,j,k 2674 out_vec(:) = 0.d0 2675 k=startRes 2676 do i=kd1,kd2 2677 do j=1,incidenceDesc%nbnonull(i) 2678 out_vec(k) = out_vec(k) + in_X(incidenceDesc%nonull(i,j))*in_xinc(i,incidenceDesc%nonull(i,j)) 2679 end do 2680 k=k+1 2681 end do 2682 2683 end subroutine matmul_incidence
model_XT_V_X
[ Top ] [ m_qtlmap_incidence ] [ Subroutines ]
NAME
model_XT_V_X
DESCRIPTION
Build X' V-1 X : X : the contingence matrix V-1 : diagonal matrix / VInv(k,k) = var(i)^2 / CD(i,j,k) , i sire,j dam ,k progeny
INPUTS
incidenceDesc : the description of the contingence matrix xinc : the contingence matrix VInv : diagonal matrix / VInv(k,k) = var(i)^2 / CD(i,j,k) , i sire,j dam ,k progeny
OUTPUTS
XVX : X' V-1 X
NOTES
Context heteroscedastic
SOURCE
2786 subroutine model_XT_V_X(xinc,incidenceDesc,VInv,XVX) 2787 type(INCIDENCE_TYPE) ,intent(in) :: incidenceDesc 2788 real (kind=dp) , dimension(nd,incidenceDesc%ntnivmax), intent(in) :: xinc 2789 real (kind=dp) , dimension(incidenceDesc%dataset%nkd,incidenceDesc%dataset%nkd) , intent(in):: VInv 2790 real (kind=dp) , dimension(incidenceDesc%ntnivmax,incidenceDesc%ntnivmax), intent(out) :: XVX 2791 2792 2793 real (kind=dp) , dimension(incidenceDesc%ntniv,incidenceDesc%dataset%nkd) :: XVXTemp 2794 real (kind=dp) , dimension(incidenceDesc%dataset%nkd,incidenceDesc%ntniv) :: temp 2795 2796 integer :: i,j,k 2797 ! 2798 ! construction de X' V-1 X 2799 ! 2800 XVX=0.d0 2801 !stop 2802 temp = xinc(:incidenceDesc%dataset%nkd,:incidenceDesc%ntniv) 2803 2804 XVXTemp=matmul(transpose(temp),VInv) 2805 XVX(:incidenceDesc%ntniv,:incidenceDesc%ntniv)=matmul(XVXTemp,temp) 2806 2807 end subroutine model_XT_V_X
model_XT_X
[ Top ] [ m_qtlmap_incidence ] [ Subroutines ]
NAME
model_XT_X
DESCRIPTION
X contingence matrix Build X' X : Model linear homoscedastic without cd or Non linear model (modlin) Dimension of X : ND line NTNIV Column Mu sire1 ...sireN ..MatEff1 Mateff Animal1 1 .... 1 AnimalN 1
INPUTS
xinc : the contingence matrix
OUTPUTS
XX : X' X
NOTES
Context homoscedastic
SOURCE
2709 subroutine model_XT_X(xinc,incidenceDesc,XX) 2710 type(INCIDENCE_TYPE) ,intent(in) :: incidenceDesc 2711 real (kind=dp) , dimension(:,:), intent(in) :: xinc 2712 real (kind=dp) , dimension(incidenceDesc%ntnivmax,incidenceDesc%ntnivmax), intent(out) :: XX 2713 2714 real (kind=dp) , dimension(incidenceDesc%dataset%nkd,incidenceDesc%ntniv) :: temp 2715 integer :: i 2716 ! print *,'-->',ntniv 2717 temp = xinc(:incidenceDesc%dataset%nkd,:incidenceDesc%ntniv) 2718 ! 2719 ! construction de X'X 2720 ! 2721 XX=0.d0 2722 XX(:incidenceDesc%ntniv,:incidenceDesc%ntniv)=matmul(transpose(temp),temp) 2723 ! print *,'' 2724 ! print *,' --XT_X-- ' 2725 ! do i=1,incidenceDesc%ntniv 2726 ! print *,XX(i,:incidenceDesc%ntniv) 2727 ! end do 2728 2729 end subroutine model_XT_X
opti_0qtl
[ Top ] [ m_qtlmap_incidence ] [ Subroutines ]
NAME
opti_0qtl
DESCRIPTION
compute the LRT under H0 1 trait
INPUTS
ic : index of the trait
FUNCT_MODEL : the function for the resolution (model_lin_h0,model_lin_hn,model_optim_h0,model_optim_hn)
OUTPUTS
incsol : the solution of the estimation
workstruct : residual variance and maximum Likelihood under H0
NOTES
SOURCE
3430 subroutine opti_0qtl( ic, & ! the trait 3431 incsol , & ! 3432 workstruct, & 3433 FUNCT_MODEL) ! function with a model specific 3434 integer ,intent(in) :: ic 3435 type(INCIDENCE_GEN_STRUCT) ,intent(inout) :: workstruct 3436 type(TYPE_INCIDENCE_SOLUTION) ,intent(out) :: incsol 3437 external :: FUNCT_MODEL ! function with a model specific 3438 3439 integer :: nkd,ip,i,listnteff(1) 3440 3441 real (kind=dp) ,dimension(:,:),pointer :: xinc,xxx 3442 real (kind=dp) , dimension(:) ,pointer :: Bestim 3443 real (kind=dp) , dimension(np) :: sigsquareEstime 3444 type(INCIDENCE_TYPE) :: incidenceDesc 3445 type(POSITION_LRT_INCIDENCE) :: curPos 3446 3447 !initialisation of contingence matrix 3448 call init_incidence(ic,0,incidenceDesc,workstruct%type_incidence) 3449 workstruct%nqtl=0 3450 allocate (Bestim(incidenceDesc%ntnivmax)) 3451 allocate (xinc(nd,incidenceDesc%ntnivmax)) 3452 xinc=0.d0 3453 Bestim=0.d0 3454 call build_incidence_matrix(workstruct,curPos,xinc,incidenceDesc,0,0,0,0) 3455 !call debug_write_incidence(xinc,incidenceDesc) 3456 !call model analysis 3457 call FUNCT_MODEL(xinc,incidenceDesc,workstruct,sigsquareEstime,Bestim,.true.) 3458 call set_solution(0,workstruct,sigsquareEstime,Bestim,incidenceDesc,incsol,1,listnteff) 3459 workstruct%sigsquare(1,:,1)=sigsquareEstime 3460 workstruct%fnqtl(1)=incidenceDesc%fmax 3461 3462 call end_incidence(incidenceDesc) 3463 deallocate (xinc) 3464 deallocate (Bestim) 3465 3466 end subroutine opti_0qtl
opti_0qtl_cuda
[ Top ] [ m_qtlmap_incidence ] [ Subroutines ]
NAME
opti_0qtl
DESCRIPTION
compute the LRT under H0 1 trait
INPUTS
ic : index of the trait
FUNCT_MODEL : the function for the resolution (model_lin_h0,model_lin_hn,model_optim_h0,model_optim_hn)
OUTPUTS
incsol : the solution of the estimation
workstruct : residual variance and maximum Likelihood under H0
NOTES
SOURCE
3514 subroutine opti_0qtl_cuda( nsim, & 3515 ic, & ! the trait 3516 incsol , & ! 3517 workstruct, & ! 3518 ySIMUL, & ! 3519 sizeF, & ! 3520 oAllsig, & 3521 FUNCT_MODEL) ! function with a model specific 3522 integer ,intent(in) :: nsim,ic 3523 type(INCIDENCE_GEN_STRUCT) ,target ,intent(inout) :: workstruct 3524 type(TYPE_INCIDENCE_SOLUTION) ,intent(out) :: incsol 3525 external :: FUNCT_MODEL ! function with a model specific 3526 real (kind=dp), dimension (nsim+1,nd) , intent(in) :: ySIMUL 3527 integer ,dimension(np) , intent(out) :: sizeF 3528 real (kind=dp), dimension (nsim+1,np) , intent(out) :: oAllsig 3529 3530 integer :: nkd,ip,i,listnteff(1),kd1,kd2,nLevelFix,nLevelVar,j 3531 3532 real (kind=dp) ,dimension(:,:),pointer :: xinc,xxx 3533 real (kind=dp) , dimension(:) ,pointer :: Bestim 3534 real (kind=dp) , dimension(np) :: sigsquareEstime 3535 type(INCIDENCE_TYPE) ,target :: incidenceDesc 3536 type(POSITION_LRT_INCIDENCE) :: curPos 3537 integer , dimension(:),allocatable ,target :: corrLevel 3538 integer :: nqtl,u,nsim2,nGLTotal,ikd,jm,kd 3539 3540 integer ,dimension(:),pointer :: allmaxpos 3541 real (kind=dp) :: seuil_cholesky 3542 3543 #ifdef CUDA_SP 3544 real , dimension(:,:,:) ,pointer :: allLrt 3545 real , dimension(:,:) ,pointer :: allmaxlrt 3546 real , dimension(:,:) ,pointer :: isigsq 3547 real , dimension(:,:,:) ,pointer :: allbestim 3548 3549 real , dimension (:,:) ,pointer :: ySIMULR 3550 real , dimension (:) ,pointer :: CDR 3551 real , dimension (:,:) ,pointer :: oAllsigR 3552 #else 3553 real (kind=dp) ,dimension(:,:,:),pointer :: allLrt 3554 real (kind=dp) ,dimension(:,:) ,pointer :: allmaxlrt 3555 real (kind=dp) , dimension(:,:) ,pointer :: isigsq 3556 real (kind=dp) , dimension(:,:,:) ,pointer :: allbestim 3557 #endif 3558 3559 3560 type(WORK_CUDA) :: w_cuda 3561 integer :: mode 3562 3563 3564 !initialisation of contingence matrix 3565 call init_incidence(ic,0,incidenceDesc,workstruct%type_incidence) 3566 workstruct%nqtl=0 3567 allocate (Bestim(incidenceDesc%ntnivmax)) 3568 allocate (xinc(nd,incidenceDesc%ntnivmax)) 3569 3570 xinc=0.d0 3571 Bestim=0.d0 3572 call build_incidence_matrix(workstruct,curPos,xinc,incidenceDesc,0,0,0,0) 3573 3574 w_cuda%xinc=>xinc 3575 w_cuda%incidenceDesc=>incidenceDesc 3576 w_cuda%workstruct=>workstruct 3577 3578 !call debug_write_incidence(xinc,incidenceDesc) 3579 !call model analysis 3580 !call FUNCT_MODEL(xinc,incidenceDesc,workstruct,sigsquareEstime,Bestim,.true.) 3581 3582 3583 do ip=1,np 3584 kd1=incidenceDesc%dataset%lSires(ip)%half_sib%firstKd 3585 kd2=incidenceDesc%dataset%lSires(ip)%half_sib%lastkd 3586 if ( kd2>kd1) then 3587 sizeF(ip) = kd2-kd1+1 3588 else 3589 sizeF(ip) = 0 3590 end if 3591 end do 3592 3593 nLevelFix=incidenceDesc%ntniv 3594 nLevelVar=0 3595 3596 w_cuda%nLevelFix=nLevelFix 3597 w_cuda%nLevelVar=nLevelVar 3598 3599 nqtl=0 3600 nGLTotal=1 3601 nsim2 = nsim+1 3602 3603 allocate (corrLevel(incidenceDesc%ntniv)) 3604 ! allocate (allvecsol(nGLTotal,incidenceDesc%ntniv)) 3605 allocate (allbestim(nGLTotal,nsim2,incidenceDesc%ntniv)) 3606 3607 w_cuda%corrLevel=>corrLevel 3608 3609 ! ** Normalement pas utilise sous H0 ** 3610 ! ************************************* 3611 allocate (isigsq(nsim2,np)) 3612 ! INITIALISER SEULEMENT POUR LE CAS HETEROSCEDASTIC 3613 !--------------------------------------------------- 3614 ! do ip=1,np 3615 ! isigsq(:,ip)=incidenceDesc%dataset%lSires(ip)%sig0(1) 3616 ! end do 3617 3618 ! FIN HETERO 3619 !--------------------------------------------------- 3620 3621 allocate (allLrt(nGLTotal,1,1)) 3622 allocate (allmaxlrt(1,1)) 3623 allocate (allmaxpos(1)) 3624 3625 do u=1,incidenceDesc%ntniv 3626 corrLevel(u)=u-1 3627 end do 3628 3629 mode = 0 3630 if (workstruct%type_model==LINEAR_HETEROSCEDASTIC) mode=1 3631 3632 #ifdef CUDA_SP 3633 allocate (ySIMULR(nsim2,incidenceDesc%dataset%nkd)) 3634 allocate (CDR(incidenceDesc%dataset%nkd)) 3635 CDR(:)=incidenceDesc%dataset%CD(1,:) 3636 ySIMULR=ySIMUL(:,:incidenceDesc%dataset%nkd) 3637 allocate (oAllsigR(nsim2,np)) 3638 seuil_cholesky=0.3 3639 3640 call cuda_homoscedatic_resolv_genome(mode, & 3641 nqtl, & 3642 isigsq, & 3643 nLevelFix, & 3644 nLevelVar, & 3645 corrLevel, & 3646 w_cuda, & 3647 get_xinc_fixed_effect, & 3648 get_partial_xinc0, & 3649 ySIMULR, & 3650 CDR, & 3651 nsim2, & 3652 nd, & 3653 incidenceDesc%dataset%nkd, & 3654 np, & 3655 sizeF, & 3656 incidenceDesc%ntnivmax, & 3657 nGLTotal, & 3658 seuil_cholesky, & 3659 allbestim, & 3660 oAllsigR, & 3661 allLrt, & 3662 allmaxlrt, & 3663 allmaxpos ) 3664 deallocate (ySIMULR) 3665 deallocate (CDR) 3666 oAllsig = oAllsigR 3667 deallocate (oAllsigR) 3668 #else 3669 seuil_cholesky=SEUIL_CHO 3670 call cuda_homoscedatic_resolv_genome(mode, & 3671 nqtl, & 3672 isigsq, & 3673 nLevelFix, & 3674 nLevelVar, & 3675 corrLevel, & 3676 w_cuda, & 3677 get_xinc_fixed_effect, & 3678 get_partial_xinc0, & 3679 ySIMUL, & 3680 incidenceDesc%dataset%CD, & 3681 nsim2, & 3682 nd, & 3683 incidenceDesc%dataset%nkd, & 3684 np, & 3685 sizeF, & 3686 incidenceDesc%ntnivmax, & 3687 nGLTotal, & 3688 seuil_cholesky, & 3689 allbestim, & 3690 oAllsig, & 3691 allLrt, & 3692 allmaxlrt, & 3693 allmaxpos ) 3694 3695 #endif 3696 3697 sigsquareEstime = oAllsig(1,:) 3698 3699 incidenceDesc%vecsol = .false. 3700 j=1 3701 do i=1,incidenceDesc%ntniv 3702 if ( allbestim(1,1,i) /= 0.d0 ) then 3703 Bestim(j) = allbestim(1,1,i) 3704 incidenceDesc%vecsol(i) = .true. 3705 j=j+1 3706 end if 3707 end do 3708 3709 !call FUNCT_MODEL(xinc,incidenceDesc,workstruct,sigsquareEstime,Bestim,.true.) 3710 call set_solution(0,workstruct,sigsquareEstime,Bestim,incidenceDesc,incsol,1,listnteff) 3711 workstruct%sigsquare(1,:,1)=sigsquareEstime 3712 workstruct%fnqtl(1)=incidenceDesc%fmax 3713 !stop 3714 3715 call end_incidence(incidenceDesc) 3716 3717 deallocate (xinc) 3718 deallocate (Bestim) 3719 deallocate (corrLevel) 3720 deallocate (isigsq) 3721 ! deallocate (allvecsol) 3722 deallocate (allbestim) 3723 deallocate (allLrt) 3724 deallocate (allmaxlrt) 3725 deallocate (allmaxpos) 3726 3727 end subroutine opti_0qtl_cuda
opti_nqtl
[ Top ] [ m_qtlmap_incidence ] [ Subroutines ]
NAME
opti_nqtl
DESCRIPTION
compute the LRT under H(nqtl) 1 trait nqtl>=1
INPUTS
nqtl : the hypothesis
ic : index of the trait
workstruct : residual variance and maximum Likelihood under H(nqtl-1)
FUNCT_MODEL : the function for the resolution (model_lin_h0,model_lin_hn,model_optim_h0,model_optim_hn)
OUTPUTS
workstruct : residual variance and maximum Likelihood under H(nqtl)
incsol : the solution of the estimation at the maximum finded
lrtsol : LRT : curves and maximum under H(nqtl)
NOTES
use the iterative scan gen_opti_nqtl function
SOURCE
4313 subroutine opti_nqtl( nqtl, & ! under hypothesis NQTL 4314 ic, & ! the trait 4315 workStruct, & ! variance / likelihood estimated under NQTL-1 hypothesis 4316 incsol , & ! incidence solution (estimation of each effect) 4317 lrtsol , & ! maximum lrt finding at a position 4318 FUNCT_MODEL ) ! function with a model specific 4319 4320 integer ,intent(in) :: nqtl,ic 4321 type(INCIDENCE_GEN_STRUCT) ,intent(inout) :: workstruct 4322 type(TYPE_INCIDENCE_SOLUTION) ,intent(out) :: incsol 4323 type(TYPE_LRT_SOLUTION) ,intent(out) :: lrtsol 4324 external :: FUNCT_MODEL 4325 4326 4327 integer :: i 4328 4329 real (kind=dp) ,dimension(:,:),pointer :: xinc ,xxx ! incidence matrix and temp for testing nuisance effect 4330 real (kind=dp) , dimension(:) ,pointer :: Bestim ! solution of the estimation 4331 real (kind=dp) , dimension(np) :: sigsquareEstime 4332 type(INCIDENCE_TYPE) :: incidenceDesc ! description incidence matrix 4333 type(POSITION_LRT_INCIDENCE) :: curPos ! the position with additional info for the model 4334 integer :: n,j,chr 4335 4336 if ( nqtl <= 0 ) then 4337 call stop_application("Error dev: can not call opti_nqtl_linear with nqtl:"//trim(str(nqtl))) 4338 end if 4339 4340 !Position Allocation 4341 call init_position (nqtl,curPos) 4342 4343 !initialisation of incidence matrix 4344 call init_incidence(ic,nqtl,incidenceDesc,workstruct%type_incidence) 4345 4346 allocate (xinc(nd,incidenceDesc%ntnivmax)) 4347 xinc=0.d0 4348 allocate (Bestim(incidenceDesc%ntnivmax)) 4349 allocate (xxx(incidenceDesc%ntnivmax,incidenceDesc%ntnivmax)) 4350 4351 Bestim=0.d0 4352 4353 call build_incidence_matrix(workstruct,curPos,xinc,incidenceDesc,0,0,0,0) 4354 4355 call FUNCT_MODEL(xinc,incidenceDesc,curPos,workstruct,sigsquareEstime,Bestim,.false.,.false.,xxx,.false.) 4356 4357 !initalizing lrtsolution 4358 call new(nqtl,lrtsol) 4359 4360 lrtsol%lrtmax=curPos%lrt 4361 lrtsol%chrmax=curPos%listChr 4362 lrtsol%nxmax=curPos%listN 4363 4364 call gen_opti_nqtl(nqtl,curPos,workstruct,xinc,incidenceDesc,sigsquareEstime,Bestim,lrtsol,FUNCT_MODEL) 4365 4366 if ( workstruct%type_incidence == INCIDENCE_TYPE_LD .or. workstruct%type_incidence == INCIDENCE_TYPE_LDLA ) then 4367 do i=1,nqtl 4368 call add_haplotype_effect(workstruct%listnteffhap(i),i,lrtsol%chrmax(i-1),lrtsol%nxmax(i-1),& 4369 xinc,incidenceDesc,workstruct%hdam,.false.) 4370 end do 4371 end if 4372 4373 if ( workstruct%type_incidence /= INCIDENCE_TYPE_LD ) then 4374 ! Compute precision at the maximum LRT in position posx 4375 do i=1,nqtl 4376 call change_qtleffect(workstruct%listnteff(i),i,lrtsol%chrmax(i-1),lrtsol%nxmax(i-1),& 4377 xinc,incidenceDesc,0,workstruct%linear) 4378 end do 4379 end if 4380 4381 do i=1,nqtl 4382 curPos%listChr(i)=lrtsol%chrmax(i-1) 4383 curPos%listN(i)=lrtsol%nxmax(i-1) 4384 end do 4385 ! call debug_write_incidence(xinc,incidenceDesc) 4386 call FUNCT_MODEL(xinc,incidenceDesc,curPos,workstruct,sigsquareEstime,Bestim,.true.,workstruct%performConfusion,xxx,.false.) 4387 if ( workstruct%performConfusion) then 4388 call confusion_type1(incidenceDesc,xxx,workstruct) 4389 end if 4390 4391 if ( workstruct%performTestNuis ) then 4392 call test_nuisances(ic,incidenceDesc,lrtsol,curPos,workstruct,sigsquareEstime,FUNCT_MODEL) 4393 end if 4394 4395 call set_solution(nqtl,workstruct,sigsquareEstime,Bestim,incidenceDesc,incsol,nqtl,workstruct%listnteff) 4396 4397 if (size(workstruct%fnqtl,1)>nqtl) then 4398 workstruct%fnqtl(nqtl+1)=incidenceDesc%fmax 4399 end if 4400 4401 if (size(workstruct%sigsquare,1)>nqtl) then 4402 workstruct%sigsquare(nqtl+1,:,1)=sigsquareEstime 4403 end if 4404 4405 call end_incidence(incidenceDesc) 4406 deallocate (xinc) 4407 deallocate (xxx) 4408 deallocate (Bestim) 4409 call end_position(curPos) 4410 4411 end subroutine opti_nqtl
opti_nqtl_cuda
[ Top ] [ m_qtlmap_incidence ] [ Subroutines ]
NAME
opti_nqtl
DESCRIPTION
compute the LRT under H(nqtl) 1 trait nqtl>=1
INPUTS
nqtl : the hypothesis
ic : index of the trait
workstruct : residual variance and maximum Likelihood under H(nqtl-1)
FUNCT_MODEL : the function for the resolution (model_lin_h0,model_lin_hn,model_optim_h0,model_optim_hn)
OUTPUTS
workstruct : residual variance and maximum Likelihood under H(nqtl)
incsol : the solution of the estimation at the maximum finded
lrtsol : LRT : curves and maximum under H(nqtl)
NOTES
use the iterative scan gen_opti_nqtl function
SOURCE
3747 subroutine opti_nqtl_cuda( nsim, & 3748 nqtl, & ! under hypothesis NQTL 3749 ic, & ! the trait 3750 workStruct, & ! variance / likelihood estimated under NQTL-1 hypothesis 3751 incsol , & ! incidence solution (estimation of each effect) 3752 lrtsol , & ! maximum lrt finding at a position 3753 ySIMUL , & 3754 listLrtSolSim, & 3755 sizeF, & 3756 AllsigQtlMoinUn, & 3757 FUNCT_MODEL) 3758 3759 integer ,intent(in) :: nsim,nqtl,ic 3760 type(INCIDENCE_GEN_STRUCT) ,intent(inout) :: workstruct 3761 type(TYPE_INCIDENCE_SOLUTION) ,intent(out) :: incsol 3762 type(TYPE_LRT_SOLUTION) ,intent(out) :: lrtsol 3763 external :: FUNCT_MODEL 3764 type(TYPE_LRT_SOLUTION) , intent(inout) ,dimension(nsim,ncar,nqtl) :: listLrtSolSim 3765 real (kind=dp), dimension (nsim+1,nd) , intent(in) :: ySIMUL 3766 integer ,dimension(np) , intent(in) :: sizeF 3767 real (kind=dp), dimension (nsim+1,np) , intent(inout) :: AllsigQtlMoinUn 3768 3769 3770 integer :: i 3771 3772 real (kind=dp) ,dimension(:,:),pointer :: xinc ,xxx ! incidence matrix and temp for testing nuisance effect 3773 real (kind=dp) , dimension(:) ,pointer :: Bestim ! solution of the estimation 3774 real (kind=dp) , dimension(np) :: sigsquareEstime 3775 type(INCIDENCE_TYPE) :: incidenceDesc ! description incidence matrix 3776 type(POSITION_LRT_INCIDENCE) :: curPos ! the position with additional info for the model 3777 integer :: n,j,chr 3778 3779 if ( nqtl <= 0 ) then 3780 call stop_application("Error dev: can not call opti_nqtl_linear with nqtl:"//trim(str(nqtl))) 3781 end if 3782 3783 !Position Allocation 3784 call init_position (nqtl,curPos) 3785 3786 !initialisation of incidence matrix 3787 call init_incidence(ic,nqtl,incidenceDesc,workstruct%type_incidence) 3788 3789 allocate (xinc(nd,incidenceDesc%ntnivmax)) 3790 xinc=0.d0 3791 allocate (Bestim(incidenceDesc%ntnivmax)) 3792 allocate (xxx(incidenceDesc%ntnivmax,incidenceDesc%ntnivmax)) 3793 3794 Bestim=0.d0 3795 3796 call build_incidence_matrix(workstruct,curPos,xinc,incidenceDesc,0,0,0,0) 3797 3798 !initalizing lrtsolution 3799 call new(nqtl,lrtsol) 3800 3801 lrtsol%lrtmax=curPos%lrt 3802 lrtsol%chrmax=curPos%listChr 3803 lrtsol%nxmax=curPos%listN 3804 3805 call gen_opti_nqtl_cuda(nsim,nqtl,curPos,workstruct,xinc,incidenceDesc,sigsquareEstime,& 3806 Bestim,lrtsol,ySIMUL,sizeF,AllsigQtlMoinUn,listLrtSolSim) 3807 3808 if ( workstruct%type_incidence == INCIDENCE_TYPE_LD .or. workstruct%type_incidence == INCIDENCE_TYPE_LDLA ) then 3809 do i=1,nqtl 3810 call add_haplotype_effect(workstruct%listnteffhap(i),i,lrtsol%chrmax(i-1),lrtsol%nxmax(i-1),& 3811 xinc,incidenceDesc,workstruct%hdam,.false.) 3812 end do 3813 end if 3814 3815 if ( workstruct%type_incidence /= INCIDENCE_TYPE_LD ) then 3816 ! Compute precision at the maximum LRT in position posx 3817 do i=1,nqtl 3818 call change_qtleffect(workstruct%listnteff(i),i,lrtsol%chrmax(i-1),lrtsol%nxmax(i-1),& 3819 xinc,incidenceDesc,0,workstruct%linear) 3820 end do 3821 end if 3822 3823 do i=1,nqtl 3824 curPos%listChr(i)=lrtsol%chrmax(i-1) 3825 curPos%listN(i)=lrtsol%nxmax(i-1) 3826 end do 3827 ! call debug_write_incidence(xinc,incidenceDesc) 3828 ! call FUNCT_MODEL(xinc,incidenceDesc,curPos,workstruct,sigsquareEstime,Bestim,.true.,workstruct%performConfusion,xxx,.false.) 3829 ! if ( workstruct%performConfusion) then 3830 ! call confusion_type1(incidenceDesc,xxx,workstruct) 3831 ! end if 3832 ! 3833 ! if ( workstruct%performTestNuis ) then 3834 ! call test_nuisances(ic,incidenceDesc,lrtsol,curPos,workstruct,sigsquareEstime,FUNCT_MODEL) 3835 ! end if 3836 call set_solution(nqtl,workstruct,sigsquareEstime,Bestim,incidenceDesc,incsol,nqtl,workstruct%listnteff) 3837 3838 if (size(workstruct%fnqtl,1)>nqtl) then 3839 workstruct%fnqtl(nqtl+1)=incidenceDesc%fmax 3840 end if 3841 3842 if (size(workstruct%sigsquare,1)>nqtl) then 3843 workstruct%sigsquare(nqtl+1,:,1)=sigsquareEstime 3844 end if 3845 3846 call end_incidence(incidenceDesc) 3847 deallocate (xinc) 3848 deallocate (xxx) 3849 deallocate (Bestim) 3850 call end_position(curPos) 3851 3852 end subroutine opti_nqtl_cuda
pdd_at_position
[ Top ] [ m_qtlmap_incidence ] [ Subroutines ]
NAME
pdd_at_position
DESCRIPTION
Compute the probalities to receive the qtl from the sire
INPUTS
ic : index of the trait
iq : index of the qtl
chr : chromosome localisation
n : the position
OUTPUTS
dataset : The user dataset
NOTES
SOURCE
1735 subroutine pdd_at_position(ic,dataset,iq,chr,n) 1736 type (DATASET_TYPE) , intent(inout) :: dataset 1737 integer ,intent(in) :: ic,iq,chr,n 1738 1739 integer :: ip,jm,kd,ig,kkd,igg,ifem,kkkd,kds 1740 real (kind=dp) :: ppt,pmt 1741 1742 do ip=1,np 1743 dataset%lSires(ip)%ppt(iq,:)=0.d0 1744 ifem=0 1745 kkkd=0 1746 do jm=nmp(ip)+1,nmp(ip+1) 1747 ifem=ifem+1 1748 if (dataset%lSires(ip)%full_sib(ifem)%firstkd<0) cycle 1749 dataset%lSires(ip)%full_sib(ifem)%ppt(iq,:,:)=0.d0 1750 dataset%lSires(ip)%full_sib(ifem)%pmt(iq,:,:)=0.d0 1751 igg=0 1752 kds=kkkd 1753 ! print *,"ip,jm,firstkd:",ip,jm,lSires(ip)%half_sib%firstKd 1754 do ig=ngenom(chr,jm)+1,ngenom(chr,jm+1) 1755 igg=igg+1 1756 kkd=0 1757 kkkd=kds 1758 do kd=ngend(chr,ig)+1,ngend(chr,ig+1) 1759 if(presentc(ic,ndesc(chr,kd))) then 1760 kkd=kkd+1 1761 kkkd=kkkd+1 1762 if( estime(ic,jm) )then 1763 dataset%lSires(ip)%full_sib(ifem)%ppt(iq,igg,kkd)=& 1764 -pdd(chr,kd,1,n)-pdd(chr,kd,2,n)+pdd(chr,kd,3,n)+pdd(chr,kd,4,n) 1765 dataset%lSires(ip)%full_sib(ifem)%pmt(iq,igg,kkd)=& 1766 -pdd(chr,kd,1,n)+pdd(chr,kd,2,n)-pdd(chr,kd,3,n)+pdd(chr,kd,4,n) 1767 else 1768 dataset%lSires(ip)%ppt(iq,dataset%lSires(ip)%half_sib%firstKd+kkkd-1)=& 1769 -pdd(chr,kd,1,n)+pdd(chr,kd,3,n) 1770 end if 1771 end if 1772 end do !! kd 1773 end do !! ig 1774 end do ! jm 1775 end do ! ip 1776 end subroutine pdd_at_position
release_copy_incidence_desc
[ Top ] [ m_qtlmap_incidence ] [ Subroutines ]
NAME
release_copy_incidence_desc
DESCRIPTION
release only arrays created during a copy (see copy_incidence_desc)
OUTPUTS
copy : the copy of the incidence description
NOTES
SOURCE
1129 subroutine release_copy_incidence_desc (copy) 1130 type(INCIDENCE_TYPE) , intent(inout) :: copy 1131 integer :: i,ip,jm 1132 1133 if ( associated (copy%desc) ) then 1134 do i=1,size(copy%desc) 1135 call release(copy%desc(i)) 1136 end do 1137 deallocate (copy%desc) 1138 end if 1139 1140 deallocate (copy%vecsol) 1141 deallocate (copy%precis) 1142 1143 if ( copy%nqtl > 0 ) then 1144 deallocate (copy%ntniv_qtlsires) 1145 deallocate (copy%ntniv_qtldams) 1146 end if 1147 1148 deallocate (copy%corr_niv_nivb) 1149 ! call end_dataset(copy%dataset) 1150 do ip=1,np 1151 do jm=1,size(copy%dataset%lSires(ip)%full_sib) 1152 if ( associated(copy%dataset%lSires(ip)%full_sib(jm)%ppt) ) then 1153 deallocate ( copy%dataset%lSires(ip)%full_sib(jm)%ppt ) 1154 deallocate ( copy%dataset%lSires(ip)%full_sib(jm)%pmt ) 1155 end if 1156 end do 1157 deallocate ( copy%dataset%lSires(ip)%full_sib) 1158 if ( associated(copy%dataset%lSires(ip)%ppt) ) deallocate(copy%dataset%lSires(ip)%ppt ) 1159 end do 1160 deallocate ( copy%dataset%lSires ) 1161 deallocate ( copy%dataset ) 1162 1163 deallocate (copy%par) 1164 if (associated(copy%nonull) ) deallocate (copy%nonull) 1165 if (associated(copy%nbnonull) ) deallocate (copy%nbnonull) 1166 1167 end subroutine release_copy_incidence_desc
release_ws
[ Top ] [ m_qtlmap_incidence ] [ Subroutines ]
NAME
release_ws
DESCRIPTION
release the workstruct structure (INCIDENCE_GEN_STRUCT)
INPUTS
workstruct : buffer structure to save information about result among the chromosome
NOTES
SOURCE
663 subroutine release_ws(workstruct) 664 type(INCIDENCE_GEN_STRUCT) ,intent(inout) :: workstruct 665 666 deallocate(workstruct%sigsquare) 667 deallocate (workstruct%listnteff) 668 deallocate (workstruct%fnqtl) 669 if (associated(workstruct%rhoi)) deallocate(workstruct%rhoi) 670 if (associated(workstruct%alertQtl)) deallocate(workstruct%alertQtl) 671 if (associated(workstruct%listtestnuis)) deallocate(workstruct%listtestnuis) 672 if (associated(workstruct%listnteffhap)) deallocate(workstruct%listnteffhap) 673 674 end subroutine release_ws
set_corrxinc
[ Top ] [ m_qtlmap_incidence ] [ Subroutines ]
NAME
set_corrxinc
DESCRIPTION
remove parameter which are not estimable from the original contingence matrix Sample ( 1 1 0 ) ( T ) ( 1 1 ) X'X = ( 1 1 0 ) Vecsol ( T ) ==> ( 1 1 ) ( 1 0 1 ) ( F ) ( 1 0 )
INPUTS
incidenceDesc : the description of contingence matrix xinc : the contingence matrix
OUTPUTS
XXOUT : xinc with column parameter not estimable removed
NOTES
SOURCE
3040 subroutine set_corrxinc(Xinc,incidenceDesc,XXOUT) 3041 type(INCIDENCE_TYPE) , intent(inout) :: incidenceDesc 3042 real (kind=dp) , dimension(:,:), intent(in) :: Xinc 3043 real (kind=dp) , dimension(nd,incidenceDesc%ntnivmax) , intent(out) :: XXOUT 3044 integer :: i,inb 3045 3046 inb=0 3047 incidenceDesc%corr_niv_nivb=-1 3048 do i=1,incidenceDesc%ntniv 3049 if(incidenceDesc%vecsol(i)) then 3050 inb=inb+1 3051 incidenceDesc%corr_niv_nivb(i)=inb 3052 XXOUT(:,inb)=Xinc(:,i) 3053 end if 3054 end do 3055 3056 end subroutine set_corrxinc
set_MR
[ Top ] [ m_qtlmap_incidence ] [ Subroutines ]
NAME
set_MR
DESCRIPTION
reduction Matrix :
Mr = Xr' Vestim -1 Xr
INPUTS
XX : X'.X
incidenceDesc : the description of contingence matrix
sizeTriangComp : size of the vector result
OUTPUTS
Mr : the number of estimable level
NOTES
This function are not using....Mr was using for the resolution of Bestim. We are
using the LU resolution (with L lower triangular matrix from cholesky decomposition)
SOURCE
2918 subroutine set_MR(XX,incidenceDesc,Mr) 2919 type(INCIDENCE_TYPE) , intent(in) :: incidenceDesc 2920 real (kind=dp) , dimension(incidenceDesc%ntnivmax,incidenceDesc%ntnivmax), intent(in) :: XX 2921 real (kind=dp) , dimension(incidenceDesc%ntnivmax,incidenceDesc%ntnivmax) , intent(out) :: Mr 2922 2923 integer :: ki,kj,i,j,ix 2924 ! 2925 ! compactage de triang 2926 ! 2927 Mr=0.d0 2928 ki=0 2929 2930 do i=1,incidenceDesc%ntniv 2931 if(incidenceDesc%vecsol(i)) then 2932 ki=ki+1 2933 kj=0 2934 do j=1,incidenceDesc%ntniv 2935 if(incidenceDesc%vecsol(j))then 2936 kj=kj+1 2937 Mr(ki,kj)=XX(i,j) 2938 end if 2939 end do 2940 end if 2941 end do 2942 2943 end subroutine set_MR
set_RHS
[ Top ] [ m_qtlmap_incidence ] [ Subroutines ]
NAME
set_RHS
DESCRIPTION
Build reduction of matrix : RHS = Xr' Xr
INPUTS
incidenceDesc : the description of contingence matrix xincreduit : the reduction of the contingence matrix
OUTPUTS
RHS : Xr' Xr
NOTES
SOURCE
2959 subroutine set_RHS(xincreduit,incidenceDesc,RHS) 2960 type(INCIDENCE_TYPE) , intent(in) :: incidenceDesc 2961 real (kind=dp) , dimension(nd,incidenceDesc%ntnivmax), intent(in) :: xincreduit 2962 real (kind=dp) , dimension(incidenceDesc%ntnivmax) , intent(out) :: RHS 2963 2964 real (kind=dp) , dimension(incidenceDesc%dataset%nkd,incidenceDesc%nbnivest) :: temp 2965 2966 RHS=0.d0 2967 temp = xincreduit(:incidenceDesc%dataset%nkd,:incidenceDesc%nbnivest) 2968 RHS(:incidenceDesc%nbnivest)=matmul(transpose(temp),incidenceDesc%dataset%Y(1,:)) 2969 2970 end subroutine set_RHS
set_RHS_V
[ Top ] [ m_qtlmap_incidence ] [ Subroutines ]
NAME
set_RHS
DESCRIPTION
Let the diagonal matrix V -1 = var(i)^2 / CD(i,j,k) , i sire,j dam ,k progeny. Build reduction of matrix : RHS = Xr' V -1 Xr
INPUTS
incidenceDesc : the description of contingence matrix xincreduit : the reduction of the contingence matrix VInv : var(i)^2 / CD(i,j,k) , i sire,j dam ,k progeny.
OUTPUTS
RHS : Xr' V -1 Xr
NOTES
SOURCE
2990 subroutine set_RHS_V(xincreduit,VInv,incidenceDesc,RHS) 2991 type(INCIDENCE_TYPE) , intent(in) :: incidenceDesc 2992 real (kind=dp) , dimension(nd,incidenceDesc%ntnivmax), intent(in) :: xincreduit 2993 real (kind=dp) , dimension(incidenceDesc%dataset%nkd,incidenceDesc%dataset%nkd) , intent(in):: VInv 2994 real (kind=dp) , dimension(incidenceDesc%ntnivmax) , intent(out) :: RHS 2995 2996 real (kind=dp) , dimension(incidenceDesc%dataset%nkd,incidenceDesc%nbnivest) :: temp 2997 real (kind=dp) , dimension(incidenceDesc%dataset%nkd)::II 2998 real (kind=dp) , dimension(incidenceDesc%nbnivest,incidenceDesc%dataset%nkd) :: RHSTemp 2999 integer :: i,kd 3000 3001 RHS=0.d0 3002 temp = xincreduit(:incidenceDesc%dataset%nkd,:incidenceDesc%nbnivest) 3003 3004 RHSTemp=matmul(transpose(temp),VInv) 3005 RHS(:incidenceDesc%nbnivest)=matmul(RHSTemp,incidenceDesc%dataset%Y(1,:)) 3006 !! II=1 3007 3008 ! RHS(:incidenceDesc%nbnivest)=matmul(transpose(temp),incidenceDesc%dataset%Y(1,:)) 3009 !! do i=1,incidenceDesc%nbnivest 3010 ! do kd=1,incidenceDesc%dataset%nkd 3011 !! RHS(i)= incidenceDesc%dataset%Y(1,incidenceDesc%dataset%nkd) 3012 !!! end do 3013 !! end do 3014 3015 end subroutine set_RHS_V
set_solution
[ Top ] [ m_qtlmap_incidence ] [ Subroutines ]
NAME
set_solution
DESCRIPTION
Fill incsol the solution from the estimimation give in parameters
INPUTS
hypothesis : the hypothesis
workstruct :
sigsquareEstime : residual variance
Bestim : solution
incidenceDesc : description of the incidence
lennteffQtl : length of listnteffQtl
listnteffQtl : index of each qtl effect inside the contingence matrix
OUTPUTS
incsol : description of the solution
NOTES
SOURCE
3178 subroutine set_solution(hypothesis,workstruct,sigsquareEstime,Bestim,incidenceDesc,incsol,lennteffQtl,listnteffQtl) 3179 integer , intent(in) :: hypothesis 3180 type(INCIDENCE_GEN_STRUCT) ,intent(in) :: workstruct 3181 real (kind=dp) , dimension(np) , intent(in) :: sigsquareEstime 3182 type(INCIDENCE_TYPE) , intent(in) :: incidenceDesc 3183 real (kind=dp) , dimension(incidenceDesc%ntnivmax) , intent(in) :: Bestim ! estimation of parameter on the maximum finded 3184 type(TYPE_INCIDENCE_SOLUTION) ,intent(inout) :: incsol 3185 integer , intent(in) ::lennteffQtl 3186 integer , dimension(lennteffQtl) , intent(in) ::listnteffQtl 3187 integer :: i,j,k,ip,maxNbPar,g,nb,INDEX,ntniv,nt 3188 3189 incsol%hypothesis=hypothesis 3190 call log_mess( "" , VERBOSE_DEF) 3191 call log_mess( "-------------------------------", VERBOSE_DEF) 3192 call log_mess( "Estimation of parameters under "//trim(str(incsol%hypothesis)), VERBOSE_DEF) 3193 call log_mess( "-------------------------------", VERBOSE_DEF) 3194 call log_mess( "rank of incidence:"//trim(str(incidenceDesc%ntniv)), VERBOSE_DEF) 3195 call log_mess( "rank max to resolve :"//trim(str(incidenceDesc%ntnivmax)), VERBOSE_DEF) 3196 3197 allocate (incsol%sig(1,np)) 3198 allocate (incsol%groupeName(incidenceDesc%nteff)) 3199 allocate (incsol%eqtl_print(incidenceDesc%nteff)) 3200 allocate (incsol%nbParameterGroup(incidenceDesc%nteff)) 3201 if (hypothesis > 0) then 3202 allocate (incsol%qtl_groupeName(1,hypothesis)) 3203 3204 if ( lennteffQtl < hypothesis ) then 3205 call stop_application("Devel error : bad init of array ** listnteff **") 3206 end if 3207 3208 do i=1,hypothesis 3209 incsol%qtl_groupeName(1,i)=listnteffQtl(i) 3210 end do 3211 3212 end if 3213 3214 maxNbPar=1 3215 do i=1,incidenceDesc%nteff 3216 if ( incidenceDesc%desc(i)%haveSubDesc) then 3217 nb=0 3218 do j=1,size(incidenceDesc%desc(i)%listSubDesc) 3219 nb = nb + incidenceDesc%desc(i)%listSubDesc(j)%end - incidenceDesc%desc(i)%listSubDesc(j)%start + 1 3220 end do 3221 maxNbPar = max(maxNbPar,nb) 3222 end if 3223 end do 3224 3225 allocate (incsol%parameterName(incidenceDesc%nteff,maxNbPar)) 3226 allocate (incsol%paramaterValue(incidenceDesc%nteff,maxNbPar)) 3227 allocate (incsol%parameterVecsol(incidenceDesc%nteff,maxNbPar)) 3228 allocate (incsol%parameterPrecis(incidenceDesc%nteff,maxNbPar)) 3229 3230 do ip=1,np 3231 incsol%sig(1,ip)=sqrt(sigsquareEstime(ip))*sigt(incidenceDesc%ic) 3232 end do 3233 3234 nt=0 3235 do i=1,incidenceDesc%nteff 3236 incsol%eqtl_print(i) = incidenceDesc%eqtl_print(i) 3237 incsol%groupeName(i) = incidenceDesc%desc(i)%name 3238 if ( incidenceDesc%desc(i)%haveSubDesc) then 3239 incsol%nbParameterGroup(i)=0 3240 do j=1,size(incidenceDesc%desc(i)%listSubDesc) 3241 INDEX=0 3242 do k=incidenceDesc%desc(i)%listSubDesc(j)%start,incidenceDesc%desc(i)%listSubDesc(j)%end 3243 INDEX=INDEX+1 3244 incsol%nbParameterGroup(i) = incsol%nbParameterGroup(i) +1 3245 g=incsol%nbParameterGroup(i) 3246 incsol%parameterName(i,g)=incidenceDesc%desc(i)%listSubDesc(j)%name//" "//adjustl(STR(INDEX)) 3247 ntniv=incidenceDesc%desc(i)%listSubDesc(j)%start+INDEX -1 3248 incsol%parameterVecsol(i,g)=incidenceDesc%vecsol(ntniv) 3249 if ( incsol%parameterVecsol(i,g) ) then 3250 nt = nt + 1 3251 incsol%paramaterValue(i,g)= Bestim(nt)*sigt(incidenceDesc%ic) 3252 incsol%parameterPrecis(i,g) = incidenceDesc%precis(ntniv) 3253 end if 3254 end do 3255 end do 3256 else 3257 incsol%nbParameterGroup(i)=1 3258 incsol%parameterName(i,1)=incidenceDesc%desc(i)%name 3259 incsol%parameterVecsol(i,1)=incidenceDesc%vecsol(incidenceDesc%desc(i)%start) 3260 if ( incsol%parameterVecsol(i,1) ) then 3261 nt = nt + 1 3262 incsol%paramaterValue(i,1) = Bestim(nt)*sigt(incidenceDesc%ic) 3263 incsol%parameterPrecis(i,1) = incidenceDesc%precis(incidenceDesc%desc(i)%start) 3264 end if 3265 end if 3266 end do 3267 3268 !specific treatment for general mean 3269 if ( workstruct%eff_general_mean > 0 ) then 3270 if ( incsol%parameterVecsol(workstruct%eff_general_mean,1) ) then 3271 incsol%paramaterValue(workstruct%eff_general_mean,1) = & 3272 incsol%paramaterValue(workstruct%eff_general_mean,1)+xmut(incidenceDesc%ic) 3273 end if 3274 end if 3275 end subroutine set_solution
set_VInv
[ Top ] [ m_qtlmap_incidence ] [ Subroutines ]
NAME
set_VInv
DESCRIPTION
Build V and VInv diagonal matrix with :
V(k,k) = CD(i,j,k) / var(i)^2
VInv(k,k) = var(i)^2 / CD(i,j,k)
i sire,j dam ,k progeny
INPUTS
incidenceDesc : description of the incidence matrix
OUTPUTS
V : CD(i,j,k) / var(i)^2 , i sire,j dam ,k progeny.
VInv : var(i)^2 / CD(i,j,k) , i sire,j dam ,k progeny.
NOTES
Context heteroscedastic
SOURCE
2749 subroutine set_VInv( incidenceDesc,sig, VInv ) 2750 type(INCIDENCE_TYPE) , intent(in) :: incidenceDesc 2751 real (kind=dp) , dimension(np), intent(in) :: sig 2752 real (kind=dp) , dimension(incidenceDesc%dataset%nkd,incidenceDesc%dataset%nkd), intent(out) :: VInv 2753 integer :: l,ip,jm,kd 2754 real (kind=dp) :: sigsquare 2755 VInv=0.d0 2756 2757 do ip=1,np 2758 sigsquare = sig(ip)*sig(ip) 2759 do kd=incidenceDesc%dataset%lSires(ip)%half_sib%firstKd,incidenceDesc%dataset%lSires(ip)%half_sib%lastKd 2760 VInv(kd,kd)=incidenceDesc%dataset%cd(1,kd)/sigsquare 2761 end do 2762 2763 end do 2764 2765 end subroutine set_VInv
test_nuisances
[ Top ] [ m_qtlmap_incidence ] [ Subroutines ]
NAME
test_nuisances
DESCRIPTION
Test des differents effets de nuisance du modele par une LRT compare a une chi2
INPUTS
ic : index of the trait
incidenceDescOrigine : description of the contingence matrix at the maximum
curPosMax : maximum position
sigsquareEstime : the residual variance at the maximum
FUNCT_MODEL : the function for the resolution (model_lin_h0,model_lin_hn,model_optim_h0,model_optim_hn)
OUTPUTS
workstruct : residual variance and maximum Likelihood under H(nqtl)
NOTES
SOURCE
4685 subroutine test_nuisances(ic,incidenceDescOrigine,lrtsol,curPosMax,workstruct,sigsquareEstime,FUNCT_MODEL) 4686 integer ,intent(in) :: ic!,nbnivestmax 4687 type(INCIDENCE_TYPE) ,intent(in) :: incidenceDescOrigine 4688 type(TYPE_LRT_SOLUTION) ,intent(in) :: lrtsol 4689 type(POSITION_LRT_INCIDENCE) ,intent(in) :: curPosMax 4690 type(INCIDENCE_GEN_STRUCT) ,intent(inout) :: workstruct 4691 real (kind=dp) , dimension(np) ,intent(in) :: sigsquareEstime 4692 external :: FUNCT_MODEL 4693 4694 real (kind=dp) ,dimension(:,:),pointer :: xinc,xxx 4695 real (kind=dp) , dimension(:) ,pointer :: Bestim,SIGest 4696 type(INCIDENCE_TYPE) :: incidenceDesc 4697 real (kind=dp) :: prob,sigEstSave(np),fmaxSave 4698 integer :: itest,i,j,nbreduit,nbef,nbco,nkd,ip,ifail,nqtl,nbnivestmax,nbeffnuis,numqtl 4699 4700 nbnivestmax = incidenceDescOrigine%nbnivest 4701 nbef=listModelTrait(ic)%nbfe 4702 nbco=listModelTrait(ic)%nbco 4703 4704 ! nbint=listModelTrait(incidenceDesc%ic)%nbint(numqtl) 4705 nqtl=workstruct%nqtl 4706 nbeffnuis = nbef+nbco 4707 do numqtl=1,nqtl 4708 nbeffnuis = nbeffnuis + listModelTrait(ic)%nbint(numqtl) 4709 end do 4710 4711 if (nbeffnuis==0) then 4712 workstruct%ntest = 0 4713 return 4714 end if 4715 4716 sigEstSave=workstruct%sigsquare(nqtl,:,1) 4717 fmaxSave=workstruct%fnqtl(nqtl) 4718 4719 workstruct%sigsquare(workstruct%nqtl,:,1)=sigsquareEstime 4720 !workstruct%fnqtl(nqtl)=incidenceDesc%fmax 4721 4722 workstruct%ntest=nbeffnuis 4723 allocate (workstruct%listtestnuis(workstruct%ntest)) 4724 itest=0 4725 4726 !initialisation of contingence matrix 4727 call init_incidence(ic,nqtl,incidenceDesc,workstruct%type_incidence,incidenceDescOrigine%dataset) 4728 4729 allocate (xinc(nd,incidenceDescOrigine%ntnivmax)) 4730 xinc=0.d0 4731 allocate (Bestim(incidenceDescOrigine%ntnivmax)) 4732 allocate (SIGest(np)) 4733 allocate (xxx(incidenceDescOrigine%ntnivmax,incidenceDescOrigine%ntnivmax)) 4734 4735 !on laisse seulement les variances comme point de depart 4736 incidenceDesc%par(np+1:)=0.d0 4737 4738 Bestim=0.d0 4739 ! workStruct% 4740 4741 do i=1,nbef+nbco 4742 if ( i <= nbef ) then 4743 call build_incidence_matrix(workstruct,curPosMax,xinc,incidenceDesc,i,0,0,0) 4744 end if 4745 if ( i > nbef .and. i <= nbef+nbco ) then 4746 call build_incidence_matrix(workstruct,curPosMax,xinc,incidenceDesc,0,i-nbef,0,0) 4747 end if 4748 4749 !call debug_write_incidence(xinc,incidenceDesc) 4750 !call model analysis 4751 4752 call FUNCT_MODEL(xinc,incidenceDesc,curPosMax,workstruct,SIGest,Bestim,.false.,.false.,xxx,.true.) 4753 4754 ifail=0 4755 if ( i > nbef .and. i <= nbco ) then 4756 nbreduit=1 4757 else 4758 nbreduit = nbnivestmax - incidenceDesc%nbnivest 4759 end if 4760 4761 prob=MATH_QTLMAP_G01ECF('U',curPosMax%lrt(workstruct%nqtl),dble(nbreduit),ifail) 4762 4763 itest=itest+1 4764 workStruct%listtestnuis(itest)%directeffect=.true. 4765 if ( i <= nbef ) then 4766 workStruct%listtestnuis(itest)%name=trim(namefix(listModelTrait(ic)%indexFixedEffect(i))) 4767 else 4768 workStruct%listtestnuis(itest)%name=trim(namefix(listModelTrait(ic)%indexCovariate(i-nbef))) 4769 end if 4770 workStruct%listtestnuis(itest)%df=nbreduit 4771 workStruct%listtestnuis(itest)%lrt=curPosMax%lrt(workstruct%nqtl) 4772 workStruct%listtestnuis(itest)%pvalue=prob 4773 4774 do while(incidenceDesc%nteff>0) 4775 call remove_last_effect(incidenceDesc,xinc) 4776 end do 4777 end do 4778 4779 do numqtl=1,nqtl 4780 do i=1,listModelTrait(incidenceDesc%ic)%nbint(numqtl) 4781 call build_incidence_matrix(workstruct,curPosMax,xinc,incidenceDesc,0,0,i,numqtl) 4782 !call debug_write_incidence(xinc,incidenceDesc) 4783 !call model analysis 4784 4785 call FUNCT_MODEL(xinc,incidenceDesc,curPosMax,workstruct,SIGest,Bestim,.false.,.false.,xxx,.true.) 4786 4787 ifail=0 4788 nbreduit = nbnivestmax - incidenceDesc%nbnivest 4789 4790 prob=MATH_QTLMAP_G01ECF('U',curPosMax%lrt(workstruct%nqtl),dble(nbreduit),ifail) 4791 4792 itest=itest+1 4793 workStruct%listtestnuis(itest)%directeffect=.false. 4794 workStruct%listtestnuis(itest)%name=trim(namefix(listModelTrait(ic)%indexFixedEffectWithInteraction(numqtl,i)))//& 4795 " QTL["//trim(str(numqtl))//"]" 4796 workStruct%listtestnuis(itest)%df=nbreduit 4797 workStruct%listtestnuis(itest)%lrt=curPosMax%lrt(workstruct%nqtl) 4798 workStruct%listtestnuis(itest)%pvalue=prob 4799 4800 do while(incidenceDesc%nteff>0) 4801 call remove_last_effect(incidenceDesc,xinc) 4802 end do 4803 end do 4804 end do 4805 4806 workstruct%sigsquare(nqtl,:,1)=sigEstSave 4807 workstruct%fnqtl(nqtl)=fmaxSave 4808 4809 call end_incidence(incidenceDesc) 4810 4811 deallocate (SIGest) 4812 deallocate (xinc) 4813 deallocate (Bestim) 4814 deallocate (xxx) 4815 4816 end subroutine test_nuisances