m_qtlmap_incidence_optim
[ Top ] [ ANALYSE ] [ Modules ]
NAME
m_qtlmap_incidence_optim
DESCRIPTION
NOTES
SEE ALSO
grand
[ Top ] [ m_qtlmap_incidence_optim ] [ Constants ]
NAME
grand
DESCRIPTION
Constant used by the discrete data analysis optimization
NOTES
SOURCE
40 real (kind=dp) ,parameter :: grand = 1.d6
current_chr
[ Top ] [ m_qtlmap_incidence_optim ] [ Variables ]
NAME
current_chr
DESCRIPTION
The currents chromosomes analysed (one for each qtl added in the contigence)
DIMENSIONS
nd, nbnivest
NOTES
used by the likelihood function
dloggrand
[ Top ] [ m_qtlmap_incidence_optim ] [ Variables ]
NAME
dloggrand
DESCRIPTION
Constant used by the discrete data analysis optimization (log of grand)
NOTES
my_desc
[ Top ] [ m_qtlmap_incidence_optim ] [ Variables ]
NAME
my_desc
DESCRIPTION
pointer of the incidence description to acces into the likelihood function (objective function to optimized)
NOTES
pointer of INCIDENCE_TYPE
my_xincreduit
[ Top ] [ m_qtlmap_incidence_optim ] [ Variables ]
NAME
my_xincreduit
DESCRIPTION
the reduction of the originale incidence matrix (without level which are not estimable according the cholesky decomposition)
DIMENSIONS
nd, nbnivest
NOTES
used by the likelihood function
likelihood_discret_h0_family
[ Top ] [ m_qtlmap_incidence_optim ] [ Subroutines ]
NAME
likelihood_discret_h0_family
DESCRIPTION
Calculus of the likelihood under H0 (without qtl effect) -- discrete data
INPUTS
ip : the index of sire jm : the index of dam n : number of parameter x : the values of parameters iuser : integer array (defined by the user) user : real array (defined by the user)
OUTPUTS
f : the value of the likelihood
NOTES
SOURCE
275 subroutine likelihood_discret_h0_family(ip,jm,n,x,f,iuser,user) 276 277 implicit none 278 integer , intent(in) :: ip,jm,n 279 ! 280 ! Tableaux dimensionnes selon n, le nombre de parametres a estimer 281 ! 282 real (kind=dp) ,dimension(n), intent(in) :: x 283 real (kind=dp) ,intent(inout) :: f 284 integer ,dimension(1), intent(in) :: iuser 285 real (kind=dp) ,dimension(1), intent(in) :: user 286 287 ! Tableaux dimensionnes selon nm, le nombre de meres 288 real (kind=dp) ,dimension(nmod(my_desc%ic)) :: threshold 289 290 !modif mkw 291 ! declaration de l comme entier et on a plus besoin de la variable vpf 292 293 integer :: i, ifail 294 integer :: neffet,ief,ifem,kd,kd1,kd2 295 real (kind=dp) :: vpf,tig 296 real(kind=dp),dimension(my_desc%dataset%nkd) :: V 297 !mkw 298 ! 299 neffet=np+my_desc%nbnivest 300 301 !****************************************************************************** 302 ! 303 ifail=1 304 305 threshold(1)=x(neffet+1) 306 do i=2,nmod(my_desc%ic)-1 307 threshold(i)=threshold(i-1)+x(neffet+i) 308 end do 309 f=0 310 ifem=jm-nmp(ip) 311 !V = matmul(my_xincreduit(:my_desc%dataset%nkd,:My_DESC%nbnivest),X(np+1:neffet)) 312 kd1=my_desc%dataset%lSires(ip)%full_sib(ifem)%firstkd 313 kd2=my_desc%dataset%lSires(ip)%full_sib(ifem)%lastkd 314 315 if (kd2<0) then 316 return 317 end if 318 319 call matmul_incidence(kd1,kd2,my_desc,nd,my_desc%ntnivmax,my_xincreduit,X(np+1:neffet),my_desc%dataset%nkd,kd1,V) 320 321 tig=grand 322 if (x(ip)> 0)tig=1.d0/x(ip) 323 324 if (my_desc%dataset%lSires(ip)%full_sib(ifem)%lastKd<0) then 325 return 326 end if 327 328 do kd=kd1,kd2 329 if (my_desc%dataset%ydiscretord(1,kd)==1) then 330 vpf=MATH_QTLMAP_G01EAF('L',((threshold(1)-v(kd))*tig),ifail) 331 else 332 if (my_desc%dataset%ydiscretord(1,kd)==nmod(my_desc%ic)) then 333 vpf=1-MATH_QTLMAP_G01EAF('L',(threshold(nmod(my_desc%ic)-1)-v(kd))*tig,ifail) 334 else 335 vpf=MATH_QTLMAP_G01EAF('L',(threshold(my_desc%dataset%ydiscretord(1,kd))-v(kd))*tig,ifail)-& 336 MATH_QTLMAP_G01EAF('L',(threshold(my_desc%dataset%ydiscretord(1,kd)-1)-v(kd))*tig,ifail) 337 end if 338 end if 339 340 if (vpf <= 0) then 341 f=INIFINY_REAL_VALUE 342 else 343 f=f-dlog(vpf) 344 end if 345 end do ! kd 346 347 end subroutine likelihood_discret_h0_family
likelihood_discret_hn_family
[ Top ] [ m_qtlmap_incidence_optim ] [ Subroutines ]
NAME
likelihood_discret_hn_family
DESCRIPTION
Calculus of the likelihood under H(nqtl) (with nqtl qtl effect) -- discrete data
INPUTS
ip : the index of sire jm : the index of dam n : number of parameter x : the values of parameters iuser : integer array (defined by the user) user : real array (defined by the user)
OUTPUTS
f : the value of the likelihood
NOTES
SOURCE
367 subroutine likelihood_discret_hn_family(ip,jm,n,x,f,iuser,user) 368 369 implicit none 370 371 integer , intent(in) :: n,ip,jm 372 ! 373 ! Tableaux dimensionnes selon n, le nombre de parametres a estimer 374 ! 375 real (kind=dp) ,dimension(n), intent(in) :: x 376 real (kind=dp) , intent(inout) :: f 377 integer , dimension(1), intent(in) :: iuser 378 real (kind=dp) ,dimension(1), intent(in) :: user 379 real (kind=dp) ,dimension(nmod(my_desc%ic)) :: threshold 380 381 ! Tableaux dimensionnes selon nm, le nombre de meres 382 real (kind=dp) ,dimension(nm) :: effm 383 384 385 ! 386 ! Divers 387 integer :: neffet,kd,kkd,ilev,lambda 388 integer :: ntnivmax,neffmax, m, i, j, m1, temp, k, ii, ifail 389 390 391 real (kind=dp) :: sig,var,vmere,vpf,wpf,tig 392 real (kind=dp) :: pbr 393 394 integer :: jjm,ig(my_desc%nqtl),ifem,indf,indm,iq,ngg,iig 395 real(kind=dp),dimension(my_desc%dataset%nkd) :: V 396 integer :: kd1,kd2,jj 397 logical :: ok 398 399 !****************************************************************************** 400 ! 401 neffet=np+my_desc%nbnivest 402 ifail=1 403 404 threshold(1)=x(neffet+1) 405 do i=2,nmod(my_desc%ic)-1 406 threshold(i)=threshold(i-1)+x(neffet+i) 407 end do 408 f=0 409 ok=.false. 410 411 tig=grand 412 if (x(ip)> 0)tig=1.d0/x(ip) 413 jjm=jm-nmp(ip) 414 415 if (my_desc%dataset%lSires(ip)%full_sib(jjm)%lastKd>0) then 416 vmere=0.d0 417 if (estime(my_desc%ic,jm)) ifem=iam(my_desc%ic,repfem(jm)) 418 kd1=my_desc%dataset%lSires(ip)%full_sib(jjm)%firstKd 419 kd2=my_desc%dataset%lSires(ip)%full_sib(jjm)%lastKd 420 ! igg=0 421 !ngg : le nombre de genotype possible sur tous les qtls... 422 ngg=1 423 do iq=1,my_desc%nqtl 424 ig(iq)=ngenom(current_chr(iq),jm)+1 425 ngg=ngg*(ngenom(current_chr(iq),jm+1)-ngenom(current_chr(iq),jm)) 426 end do 427 do iig=1,ngg 428 ! igg=0 429 pbr=1 430 !on modifie la matrice d incidence pour les n qtl 431 do iq=1,my_desc%nqtl 432 indm=my_desc%ntniv_qtlsires(iq)+ip-1 433 if ( estime(my_desc%ic,jm) ) then 434 !Si la mere est estimable, on place dans la matrice les 435 !pdds dam et pdds sires (si ces effets sont estimables) 436 indf=my_desc%ntniv_qtldams(iq)+ifem-1 437 if ( My_desc%vecsol(indf)) then 438 indf=my_desc%corr_niv_nivb(indf) ! vrai correspondance dans le tableau si estimable 439 my_xincreduit(kd1:kd2,indf)=my_desc%dataset%lSires(ip)%full_sib(jjm)%pmt(iq,ig(iq)-ngenom(current_chr(iq),jm),:) 440 end if 441 if ( My_desc%vecsol(indm) ) then 442 indm=my_desc%corr_niv_nivb(indm) ! vrai correspondance dans le tableau si estimable 443 my_xincreduit(kd1:kd2,indm)=my_desc%dataset%lSires(ip)%full_sib(jjm)%ppt(iq,ig(iq)-ngenom(current_chr(iq),jm),:) 444 end if 445 else 446 !la mere n est pas estimable, on place seulement les pdds males 447 if ( My_desc%vecsol(indm) ) then 448 indm=my_desc%corr_niv_nivb(indm) ! vrai correspondance dans le tableau si estimable 449 my_xincreduit(kd1:kd2,indm)=my_desc%dataset%lSires(ip)%ppt(iq,kd1:kd2) 450 end if 451 end if 452 pbr=pbr*probg(current_chr(iq),ig(iq)) 453 end do ! iq 454 455 !V(kd1:kd2) = matmul(my_xincreduit(kd1:kd2,:My_DESC%nbnivest),X(np+1:neffet)) 456 call matmul_incidence(kd1,kd2,my_desc,nd,my_desc%ntnivmax,my_xincreduit,X(np+1:neffet),my_desc%dataset%nkd,kd1,V) 457 458 vpf=grand 459 do kkd=kd1,kd2 460 if (my_desc%dataset%ydiscretord(1,kkd)==1) then 461 wpf=MATH_QTLMAP_G01EAF('L',((threshold(1)-v(kkd))*tig),ifail) 462 else 463 if (my_desc%dataset%ydiscretord(1,kkd)==nmod(my_desc%ic)) then 464 wpf=1.d0-MATH_QTLMAP_G01EAF('L',(threshold(nmod(my_desc%ic)-1)-v(kkd))*tig,ifail) 465 else 466 wpf=MATH_QTLMAP_G01EAF('L',(threshold(my_desc%dataset%ydiscretord(1,kkd))-v(kkd))*tig,ifail)-& 467 MATH_QTLMAP_G01EAF('L',(threshold(my_desc%dataset%ydiscretord(1,kkd)-1)-v(kkd))*tig,ifail) 468 end if 469 end if 470 vpf=vpf*wpf 471 end do 472 473 vmere=vmere+pbr*vpf 474 ! on increment 475 do iq=1,my_desc%nqtl 476 ! print *,"test:",ig(iq),ngenom(current_chr(iq),jm),iig,ngg 477 if ((ig(iq) < ngenom(current_chr(iq),jm+1))) then 478 ig(iq)=ig(iq)+1 479 exit 480 end if 481 end do 482 end do !ig 483 484 if (vmere <= 0) then 485 f=INIFINY_REAL_VALUE 486 return 487 endif 488 489 f=dloggrand-dlog(vmere) 490 end if 491 492 return 493 end subroutine likelihood_discret_hn_family
likelihood_h0_family
[ Top ] [ m_qtlmap_incidence_optim ] [ Subroutines ]
NAME
likelihood_h0_family
DESCRIPTION
Calculus of the likelihood under H0 (without qtl effect) -- continue data
INPUTS
ip : the index of sire jm : the index of dam n : number of parameter x : the values of parameters iuser : integer array (defined by the user) user : real array (defined by the user)
OUTPUTS
f : the value of the likelihood
NOTES
SOURCE
109 subroutine likelihood_h0_family(ip,jm,n,x,f,iuser,user) 110 integer , intent(in) :: ip,jm,n 111 real (kind=dp) ,dimension(n), intent(in) :: x 112 real (kind=dp) , intent(inout) :: f 113 integer , dimension(1), intent(in) :: iuser 114 real (kind=dp) ,dimension(1), intent(in) :: user 115 116 real (kind=dp) :: var,sig,vmere,vpf 117 integer :: ifem,kd1,kd2 118 real(kind=dp),dimension(my_desc%dataset%nkd) :: V 119 120 f = 0.d0 121 V = 0.d0 122 ifem=jm-nmp(ip) 123 124 kd2=my_desc%dataset%lSires(ip)%full_sib(ifem)%lastkd 125 126 if (kd2<0) then 127 return 128 end if 129 kd1=my_desc%dataset%lSires(ip)%full_sib(ifem)%firstkd 130 131 ! V = matmul(my_xincreduit(:my_desc%dataset%nkd,:My_DESC%nbnivest),X(np+1:)) 132 call matmul_incidence(kd1,kd2,my_desc,nd,my_desc%ntnivmax,my_xincreduit,X(np+1:),my_desc%dataset%nkd,kd1,V) 133 V(kd1:kd2) = my_desc%dataset%Y(1,kd1:kd2)-V(kd1:kd2) 134 V(kd1:kd2) = V(kd1:kd2)*V(kd1:kd2)*my_desc%dataset%CD(1,kd1:kd2) 135 sig=x(ip) 136 var=sig*sig 137 vpf = sum(V(kd1:kd2)) 138 vmere=dexp(-0.5d0*vpf/var) 139 if (vmere == 0) then 140 f=INIFINY_REAL_VALUE 141 return 142 end if 143 144 f=-dlog(vmere)+(kd2-kd1+1)*dlog(sig) 145 146 end subroutine likelihood_h0_family
likelihood_hn_family
[ Top ] [ m_qtlmap_incidence_optim ] [ Subroutines ]
NAME
likelihood_hn_family
DESCRIPTION
Calculus of the likelihood under H(nqtl) (with nqtl qtl effect) -- continue data
INPUTS
ip : the index of sire jm : the index of dam n : number of parameter x : the values of parameters iuser : integer array (defined by the user) user : real array (defined by the user)
OUTPUTS
f : the value of the likelihood
NOTES
SOURCE
166 subroutine likelihood_hn_family(ip,jm,n,x,f,iuser,user) 167 ! use OMP_LIB 168 169 integer , intent(in) :: ip,jm,n 170 real (kind=dp) ,dimension(n), intent(in) :: x 171 real (kind=dp) , intent(inout) :: f 172 integer , dimension(1), intent(inout) :: iuser 173 real (kind=dp) ,dimension(1), intent(inout) :: user 174 175 real (kind=dp) :: var,sig,vmere,vpf,pbr,effm 176 integer :: jjm,ig(my_desc%nqtl),ifem,indf,indm,iq,ngg,iig,z 177 real(kind=dp),dimension(my_desc%dataset%nkd) :: V 178 integer :: kd1,kd2,jj 179 logical :: ok 180 181 f=0.d0 182 sig=x(ip) 183 var=sig*sig 184 ! print *,'ip:',ip,' jm:',jm,' nmp(ip):',nmp(ip),' nmp(ip+1)',nmp(ip+1) 185 jjm=jm-nmp(ip) 186 ! print *,'jjm:',jjm 187 if (my_desc%dataset%lSires(ip)%full_sib(jjm)%lastKd>0) then 188 kd1=my_desc%dataset%lSires(ip)%full_sib(jjm)%firstKd 189 kd2=my_desc%dataset%lSires(ip)%full_sib(jjm)%lastKd 190 effm=dble(kd2-kd1+1) 191 vmere=0.d0 192 193 if (estime(my_desc%ic,jm)) ifem=iam(my_desc%ic,repfem(jm)) 194 !ngg : le nombre de genotype possible sur tous les qtls... 195 ngg=1 196 do iq=1,my_desc%nqtl 197 ig(iq)=ngenom(current_chr(iq),jm)+1 198 ngg=ngg*(ngenom(current_chr(iq),jm+1)-ngenom(current_chr(iq),jm)) 199 end do 200 !pour toutes les combinaisons possibles des genotypes 201 do iig=1,ngg 202 pbr=1 203 !on modifie la matrice d incidence pour les n qtl 204 do iq=1,my_desc%nqtl 205 indm=my_desc%ntniv_qtlsires(iq)+ip-1 206 if ( estime(my_desc%ic,jm) ) then 207 !Si la mere est estimable, on place dans la matrice les 208 !pdds dam et pdds sires (si ces effets sont estimables) 209 indf=my_desc%ntniv_qtldams(iq)+ifem-1 210 if ( My_desc%vecsol(indf)) then 211 indf=my_desc%corr_niv_nivb(indf) ! vrai correspondance dans le tableau si estimable 212 my_xincreduit(kd1:kd2,indf)=my_desc%dataset%lSires(ip)%full_sib(jjm)%pmt(iq,(ig(iq)-ngenom(current_chr(iq),jm)),:) 213 end if 214 if ( My_desc%vecsol(indm) ) then 215 indm=my_desc%corr_niv_nivb(indm) ! vrai correspondance dans le tableau si estimable 216 my_xincreduit(kd1:kd2,indm)=my_desc%dataset%lSires(ip)%full_sib(jjm)%ppt(iq,(ig(iq)-ngenom(current_chr(iq),jm)),:) 217 end if 218 else 219 !la mere n est pas estimable, on place seulement les pdds males 220 if ( My_desc%vecsol(indm) ) then 221 indm=my_desc%corr_niv_nivb(indm) ! vrai correspondance dans le tableau si estimable 222 my_xincreduit(kd1:kd2,indm)=my_desc%dataset%lSires(ip)%ppt(iq,kd1:kd2) 223 end if 224 end if 225 !print *,ig(iq)-ngenom(current_chr(iq),jm),ig(iq) 226 pbr=pbr*probg(current_chr(iq),ig(iq)) 227 end do ! iq 228 229 !V(kd1:kd2) = matmul(my_xincreduit(kd1:kd2,:My_DESC%nbnivest),X(np+1:)) 230 call matmul_incidence(kd1,kd2,my_desc,nd,my_desc%ntnivmax,my_xincreduit,X(np+1:),my_desc%dataset%nkd,kd1,V) 231 V(kd1:kd2) = my_desc%dataset%Y(1,kd1:kd2) - V(kd1:kd2) 232 V(kd1:kd2) = V(kd1:kd2)*V(kd1:kd2)*my_desc%dataset%CD(1,kd1:kd2) 233 vpf = sum(V(kd1:kd2)) 234 235 vmere=vmere+pbr*dexp(-0.5d0*vpf/var) 236 ! on increment 237 ok=.true. 238 do iq=1,my_desc%nqtl 239 if (ok) then 240 if ((ig(iq) < ngenom(current_chr(iq),jm+1))) then 241 ig(iq)=ig(iq)+1 242 ok=.false. 243 end if 244 end if 245 end do 246 end do ! iig 247 248 if (vmere == 0) then 249 f=INIFINY_REAL_VALUE 250 else 251 f=-dlog(vmere)+effm*dlog(sig) 252 end if 253 end if 254 255 end subroutine likelihood_hn_family
model_optim_family
[ Top ] [ m_qtlmap_incidence_optim ] [ Subroutines ]
NAME
model_optim_family
DESCRIPTION
Generic model for an optimization of a likelihood (discrete and continue data)
process :
* Cholesky decomposition (to have the estimability of each level of the contingence)
* build a matrix without no-null element to improve the multiplication matrix during the likelihood calculus
* compute threhold in discrete data case
* build a filter to improve the optimization (see m_qtlmap_optimization)
INPUTS
xinc : the contigence matrix incidenceDesc : description of the contingence and incidence matrix performPrecision : indicates if the precision have to be computed FUNCT_PART : the likelihood function to optimize tConf : indicates if the caller want to keep a buffer array in order to compute the confusion between parameter
OUTPUTS
f : the maximum value of the likelihood sigsquareEstime : residual variance Bestim : solution values of each level of the contingence matrix tempForConfusion : a buffer array for the confusion calculus
NOTES
SOURCE
524 subroutine model_optim_family(xinc,incidenceDesc,f,osig,sigsquareEstime,Bestim,& 525 performPrecision,FUNCT_PART,tConf,tempForConfusion) 526 type(INCIDENCE_TYPE),target , intent(inout) :: incidenceDesc 527 real (kind=dp) , dimension(nd,incidenceDesc%ntnivmax), intent(in) :: xinc 528 real (kind=dp) ,intent(in) :: f 529 ! type(INCIDENCE_GEN_STRUCT) ,intent(in) :: workstruct 530 real (kind=dp) , dimension(np) ,intent(out) :: osig 531 real (kind=dp) , dimension(np) ,intent(out) :: sigsquareEstime 532 real (kind=dp) , dimension(incidenceDesc%ntnivmax) ,intent(out) :: Bestim 533 logical ,intent(in) :: performPrecision,tConf 534 real (kind=dp),dimension(incidenceDesc%ntnivmax,incidenceDesc%ntnivmax) ,intent(out) :: tempForConfusion 535 external :: FUNCT_PART 536 537 integer :: j,i,ip,ifail,ibound,npar 538 ! real(kind=dp) ,dimension(np+ntnivmax) :: par 539 real (kind=dp) , dimension(incidenceDesc%ntnivmax,incidenceDesc%ntnivmax) :: XX,xxx 540 real (kind=dp) , dimension(incidenceDesc%ntniv,incidenceDesc%ntniv) :: triang 541 integer :: iuser(1),ix,kd1,kd2,jm,jjm 542 real (kind=dp) :: user(1) 543 logical ,dimension(incidenceDesc%ntnivmax) :: lastvecsol 544 logical , dimension(:,:,:),pointer :: filter_inc 545 real (kind=dp) ,dimension(np) :: fp 546 real (kind=dp) ,dimension(nm) :: fm 547 548 my_desc=>incidenceDesc 549 lastvecsol=my_desc%vecsol 550 ! create X'.X matrix from incidence matrix 551 call model_XT_X(xinc,my_desc,XX) 552 ! Check all parameters to remove from the estimation 553 call estim_cholesky(XX,my_desc,my_desc%ntniv,triang) 554 ! compute the precision of each parameter 555 if (performPrecision) call get_precision(XX,xxx,my_desc) 556 if (tConf) tempForConfusion=xxx 557 allocate (my_xincreduit(nd,incidenceDesc%ntnivmax)) 558 my_xincreduit=0.d0 559 call set_corrxinc(xinc,my_desc,my_xincreduit) 560 !optimisation : on sauvegarde les index des elements non nul de la matrice reduite 561 call fill_nonull_elements(my_desc,size(my_xincreduit,1),size(my_xincreduit,2),my_xincreduit) 562 ! Optimisation de la vraisemblance a la position dx 563 ifail=1 564 ibound=0 565 npar=np+my_desc%nbnivest 566 567 if (count(lastvecsol)>0) then ! autre que initialisation 568 j=np 569 do i=1,my_desc%ntniv 570 if(my_desc%vecsol(i))then 571 j=j+1 572 if(.not. lastvecsol(i)) incidenceDesc%par(j)=0.d0 573 end if 574 end do 575 else 576 ! print *,'initialisation......' 577 incidenceDesc%par=0.d0 578 do ip=1,np 579 incidenceDesc%par(ip)=incidenceDesc%dataset%lSires(ip)%sig0(1) 580 end do 581 end if 582 583 if ( nmod(incidenceDesc%ic) > 1) then 584 npar=npar + nmod(incidenceDesc%ic)-1 585 incidenceDesc%par(np+my_desc%nbnivest+1)=seuil(incidenceDesc%ic,1) 586 do ix=2,nmod(incidenceDesc%ic)-1 587 incidenceDesc%par(np+my_desc%nbnivest+ix)=seuil(incidenceDesc%ic,ix)-seuil(incidenceDesc%ic,ix-1) 588 end do 589 end if 590 !call debug_write_incidence(xinc,incidenceDesc) 591 !Preparation de l optimisation.... 592 !call init_optim(incidenceDesc%ic) 593 594 allocate (filter_inc(np,nm,npar)) 595 filter_inc=.true. 596 filter_inc(:,:,1:np+1)=.false. 597 do ip=1,np 598 filter_inc(ip,nmp(ip)+1:nmp(ip+1),ip)=.true. 599 jjm=0 600 do jm=nmp(ip)+1,nmp(ip+1) 601 jjm=jjm+1 602 if (my_desc%dataset%lSires(ip)%full_sib(jjm)%lastKd>0) then 603 kd1=my_desc%dataset%lSires(ip)%full_sib(jjm)%firstkd 604 kd2=my_desc%dataset%lSires(ip)%full_sib(jjm)%lastkd 605 filter_inc(ip,jm,np+1:np+my_desc%nbnivest)=any(my_xincreduit(kd1:kd2,:my_desc%nbnivest)/=0.d0,dim=1) 606 else 607 ! print *,'pas de desc....' 608 filter_inc(ip,jm,np+1:np+my_desc%nbnivest)=.false. 609 end if 610 end do 611 end do 612 ! filter_inc=.true. 613 614 !stop 615 616 call minimizing_funct_family(npar,ibound,FUNCT_PART,filter_inc,fm,fp,incidenceDesc%borni,incidenceDesc%borns,& 617 incidenceDesc%par,f,iuser,user,ifail) 618 619 deallocate (filter_inc) 620 621 !getting standart deviation 622 do ip=1,np 623 osig(ip)=incidenceDesc%par(ip) 624 sigsquareEstime(ip)=osig(ip)*osig(ip) 625 end do 626 !The solution 627 Bestim(:my_desc%nbnivest)=incidenceDesc%par(np+1:np+my_desc%nbnivest) 628 629 deallocate (my_xincreduit) 630 ! 631 end subroutine model_optim_family
model_optim_h0
[ Top ] [ m_qtlmap_incidence_optim ] [ Subroutines ]
NAME
model_optim_h0
DESCRIPTION
INPUTS
xinc : the contigence matrix incidenceDesc : description of the contingence and incidence matrix performPrecision : indicates if the precision of each parameter have to be computed
INPUTS/OUTPUTS
workstruct : result information about current point (not available in this function) and the specific method to used (linear or not, homo,heteroscedastic)
OUTPUTS
sigsquareEstime : residual variance Bestim : solution values of each level of the contingence matrix
NOTES
SOURCE
654 subroutine model_optim_h0(xinc,incidenceDesc,workstruct,sigsquareEstime,Bestim,performPrecision) 655 type(INCIDENCE_TYPE) , intent(inout) :: incidenceDesc 656 real (kind=dp) , dimension(nd,incidenceDesc%ntnivmax), intent(in) :: xinc 657 type(INCIDENCE_GEN_STRUCT) , intent(in) :: workstruct 658 real (kind=dp) , dimension(np) , intent(out) :: sigsquareEstime 659 real (kind=dp) , dimension(incidenceDesc%ntnivmax) , intent(out) :: Bestim 660 logical , intent(in) :: performPrecision 661 662 real (kind=dp) ,dimension(np) :: osig 663 real(kind=dp) , dimension(incidenceDesc%ntnivmax,incidenceDesc%ntnivmax) :: xxx 664 665 integer :: ip 666 667 if ( natureY(incidenceDesc%ic) == 'r') then 668 call model_optim_family(xinc,incidenceDesc,incidenceDesc%fmax,& 669 osig,sigsquareEstime,Bestim,performPrecision,likelihood_h0_family,.false.,xxx) 670 671 else if ( natureY(incidenceDesc%ic) == 'i' ) then 672 dloggrand= dlog(grand) 673 call model_optim_family(xinc,incidenceDesc,& 674 incidenceDesc%fmax,osig,sigsquareEstime,Bestim,performPrecision,likelihood_discret_h0_family,.false.,xxx) 675 end if 676 677 end subroutine model_optim_h0
model_optim_hn
[ Top ] [ m_qtlmap_incidence_optim ] [ Subroutines ]
NAME
model_optim_hn
DESCRIPTION
INPUTS
xinc : the contigence matrix incidenceDesc : description of the contingence and incidence matrix curPos : information about the current position performPrecision : indicates if the precision of each parameter have to be computed tConf : indicates if the caller want to keep a buffer array in order to compute the confusion between parameter invlrt : true => 2* (F0max - Fnqtlmax), false => 2* ( Fnqtlmax - F0max )
INPUTS/OUTPUTS
workstruct : result information about current point and the specific method to used (linear or not, homo,heteroscedastic)
OUTPUTS
sigsquareEstime : residual variance Bestim : solution values of each level of the contingence matrix tempForConfusion : a buffer array for the confusion calculus
NOTES
SOURCE
703 subroutine model_optim_hn(xinc,incidenceDesc,curPos,workstruct,sigsquareEstime,Bestim,& 704 performPrecision,tConf,tempForConfusion,invlrt) 705 type(INCIDENCE_TYPE) , intent(inout) :: incidenceDesc 706 real (kind=dp) , dimension(nd,incidenceDesc%ntnivmax), intent(in) :: xinc 707 type(POSITION_LRT_INCIDENCE) ,intent(inout) :: curPos 708 type(INCIDENCE_GEN_STRUCT) ,intent(inout) :: workstruct 709 real (kind=dp) , dimension(incidenceDesc%ntnivmax) ,intent(out) :: Bestim 710 real (kind=dp) , dimension(np) ,intent(out) :: sigsquareEstime 711 logical ,intent(in) :: performPrecision,tConf 712 real (kind=dp),dimension(incidenceDesc%ntnivmax,incidenceDesc%ntnivmax) ,intent(out) :: tempForConfusion 713 logical ,intent(in) :: invlrt 714 715 real (kind=dp) :: osig(np),t1 716 integer :: hypothesis,i 717 character(len=LEN_L) :: dx 718 719 allocate (current_chr(workstruct%nqtl)) 720 current_chr=curPos%listChr 721 722 incidenceDesc%par(:np)=sqrt(workstruct%sigsquare(workstruct%nqtl,:,1)) 723 724 if ( natureY(incidenceDesc%ic) == 'r') then 725 call model_optim_family(xinc,incidenceDesc,& 726 incidenceDesc%fmax,osig,sigsquareEstime,Bestim,performPrecision,likelihood_hn_family,tConf,tempForConfusion) 727 else if ( natureY(incidenceDesc%ic) == 'i') then 728 dloggrand= dlog(grand) 729 call model_optim_family(xinc,incidenceDesc,& 730 incidenceDesc%fmax,osig,sigsquareEstime,Bestim,performPrecision,likelihood_discret_hn_family,tConf,tempForConfusion) 731 732 end if 733 734 if ((incidenceDesc%fmax == INIFINY_REAL_VALUE) ) then 735 dx="" 736 do i=1,workstruct%nqtl-1 737 dx=trim(dx)//trim(str(absi(curPos%listChr(i),curPos%listN(i))))//"," 738 end do 739 740 dx=trim(dx)//str(absi(curPos%listChr(workstruct%nqtl),curPos%listN(workstruct%nqtl))) 741 call log_mess("dx ["//trim(dx)// "]. Can not optimize likelihood....The start point is reinitializing",WARNING_DEF) 742 incidenceDesc%par=0.d0 743 curPos%lrt=0.d0 744 incidenceDesc%fmax=0.d0 745 return 746 end if 747 748 !compute LRT 749 if (invLrt) then 750 do hypothesis=1,workstruct%nqtl 751 curPos%lrt(hypothesis)=-2.d0*(workstruct%fnqtl(hypothesis)-incidenceDesc%fmax) 752 end do 753 else 754 do hypothesis=1,workstruct%nqtl 755 curPos%lrt(hypothesis)=-2.d0*(incidenceDesc%fmax-workstruct%fnqtl(hypothesis)) 756 end do 757 end if 758 759 deallocate (current_chr) 760 761 end subroutine model_optim_hn