animal model - gauss seidel, iteration on data

김시동 박사의 프로그램

 

animal model

 

gauss siedel, iteration on data 방식으로 해 구하기

 

각 effect별로 동시에 해를 구함.animal effect를 다룰 때에는 animal을 자손으로 조상 순으로 처리.

 

자료

 

1 12 30
2 12 32
3 11 38
4 12 33
5 11 40
6 12 28
7 12 34
8 11 37
9 12 31
10 11 39
1st - animal effect2nd - fixed effect3rd - observation

 

위 자료를 tayler.d 로 저장

 

혈통

 

1 0 0
2 0 0
3 0 0
4 0 0
5 0 0
6 3 2
7 3 4
8 5 0
9 5 0
10 8 7
1st - animal2nd - sire3rd - dam

 

위 혈통을 tayler.p 로 저장

 

프로그램

 

!
! Solve Animal Model
! GS iteration on data
! 2009. 4. 18
! programmed by sidong kim
!
program gs_iod
integer :: no_col, it, if, itr, iostat
real, dimension(100) :: data
real (kind=4), dimension(1000) :: sol, diag, rhs
character(128) :: dfn, pfn, solfn
integer, dimension(10) :: level_s, level_e
integer :: idx, idx2, rpos, noped
integer, dimension(100,4) :: ped ! for pedigree
real :: alpha
real, dimension(0:2) :: dd
integer :: offset dd(0) = 1.0
dd(1) = 4.0/3.0
dd(2) = 2.0 write(*,fmt='(A)', advance="no") 'Data file name :'
read *, dfn
write(*,fmt='(A)', advance="no") 'Solution file name :'
read *, solfn
write(*,fmt='(A)', advance="no") 'Pedigree file name :'
read *, pfn
write(*,fmt='(A)', advance="no") 'Number of iteration :'
read *, itr
write(*,fmt='(A)', advance="no") 'Number of column :'
read *, no_col
write(*,fmt='(A)', advance="no") 'Position of random effect :'
read *, rpos
write(*,fmt='(A)', advance="no") 'Alpha = '
read *, alpha open(unit=10, file=dfn, status="old")
open(unit=11, file=solfn, status="replace")
open(unit=12, file=pfn, status="old") ! read pedigree file and store it to ped array
ped = 0 ! initialize
ii = 1
do
read(unit=12, fmt=*, iostat=iostat) ped(ii,:3) ! read animal, sire and dam
if (iostat /= 0) exit
if(ped(ii,2) == 0 .and. ped(ii,3) == 0) then ! both parents are unknown
ped(ii,4) = 0
else if(ped(ii,2) == 0 .or. ped(ii,3) == 0) then ! one parent is known
ped(ii,4) = 1
else ! both parents are known
ped(ii,4) = 2
end if
ii = ii + 1
end do
close(12)
noped = ii - 1 ! number of animals in pedigree ! read data file
! check level
level_s = 1
level_e = 1
do
read(unit=10, fmt=*,iostat=iostat) data(:no_col) ! read one data line, animal, fixed effect, observation
if(iostat /=0) exit
! check level of effects
do if = 1, no_col - 1
if(data(if) > level_e(if)) level_e(if) = data(if)
end do
end do
! check start and end level of each effect
do if = 1, no_col - 2
level_s(if+1) = level_e(if)+1
end do print *, level_s
print *, level_e

 

! solve
sol = 0
do it = 1, itr ! iteration
do if = 1, no_col - 1 ! get solution for each effect
rewind 10 ! rewind data file
rhs = 0
diag = 0
do ! as fixed effect
read(unit=10, fmt=*, iostat=iostat) data(:no_col)
if(iostat /= 0) exit
idx = data(if)
rhs(idx) = rhs(idx) + data(no_col)
do ii = 1, no_col-1
if(if == ii) then ! diagonal
diag(idx) = diag(idx) + 1.0
else ! off diagonal
idx2= data(ii)
rhs(idx) = rhs(idx) - sol(idx2)
end if
end do
end do ! get solution for effect, if
if(if == rpos) then ! it's animal effect
do ii = noped, 1, -1 ! process yongest to oldest
offset = ped(ii,1)
! offset = level_s(if) + ii - 1 ! i-th animal's diagonal
diag(offset) = diag(offset) + dd(ped(ii,4))*alpha
! i-th animal's offdiagonal
if(ped(ii,2) /=0) & ! if sire exists
rhs(offset) = rhs(offset) + dd(ped(ii,4))/2.0*alpha*sol(ped(ii,2))
if(ped(ii,3) /=0) & ! if dam exists
rhs(offset) = rhs(offset) + dd(ped(ii,4))/2.0*alpha*sol(ped(ii,3)) ! get new solution
sol(offset) = rhs(offset) / diag(offset) ! adjust the rest pedigree
if(ped(ii,2) /=0) then
! s,s
diag(ped(ii,2)) = diag(ped(ii,2)) + dd(ped(ii,4))/4.0*alpha
! s,i
rhs(ped(ii,2)) = rhs(ped(ii,2)) + dd(ped(ii,4))/2.0*alpha*sol(ii)
! s,d
if(ped(ii,3) /=0) &
rhs(ped(ii,2)) = rhs(ped(ii,2)) - dd(ped(ii,4))/4.0*alpha*sol(ped(ii,3))
end if
if(ped(ii,3) /=0) then
! d,d
diag(ped(ii,3)) = diag(ped(ii,3)) + dd(ped(ii,4))/4.0*alpha
! d,i
rhs(ped(ii,3)) = rhs(ped(ii,3)) + dd(ped(ii,4))/2.0*alpha*sol(ii)
! d,s
if(ped(ii,2) /=0) &
rhs(ped(ii,3)) = rhs(ped(ii,3)) - dd(ped(ii,4))/4.0*alpha*sol(ped(ii,2))
end if
end do
else ! fixed effect
do ii = level_s(if), level_e(if)
sol(ii) = rhs(ii) / diag(ii)
end do
end if end do
print *, it, sol(1) ! too see convergence, print first solution
end do
write(11,*) 'Eqn sol'
do ii = 1, level_e(no_col-1)
write(11, *) ii, sol(ii)
end do
close(10)
close(11)
end program
위 프로그램을 gsiod3.f95 로 저장

 

위 자료, 혈통, 프로그램을 이용한 해

 

Eqn sol
1 -0.18691799
2 -0.11230095
3 -0.10126753
4 0.42026323
5 0.10768158
6 -0.35306439
7 0.36787596
8 -9.70994607E-02
9 1.35889370E-02
10 0.16486143
11 38.481461
12 31.308426
프로그램 컴파일 및 실행 예