multiple trait animal model - solve

multiple trait animal mdoel setup에서 lhs.dat와 rhs.dat 준비

 

lhs.dat의 sort

 

명령어 : sort -n -o sorted_lhs.dat lhs.dat

 

다음은 sorted_lhs.dat

 

1 4 14 7 1
1 5 15 7 1
1 6 16 7 1
2 4 17 7 1
2 6 18 7 1
2 6 19 7 1
3 4 21 7 1
3 5 20 7 1
3 5 23 3 1
3 6 22 3 1
4 1 14 7 1
4 2 17 7 1
4 3 21 7 1
5 1 15 7 1
5 3 20 7 1
5 3 23 3 1
6 1 16 7 1
6 2 18 7 1
6 2 19 7 1
6 3 22 3 1
7 0 0 0 6
7 14 12 0 3
7 17 12 0 3
7 20 12 0 3
8 0 0 0 6
8 15 12 0 3
8 18 13 0 3
8 21 12 0 3
9 0 0 0 6
9 16 13 0 3
9 22 13 0 3
10 0 0 0 6
10 19 13 0 3
11 0 0 0 6
11 23 13 0 3
12 0 0 0 6
12 14 7 0 3
12 15 8 0 3
12 17 7 0 3
12 20 7 0 3
12 21 8 0 3
13 0 0 0 6
13 16 9 0 3
13 18 8 0 3
13 19 10 0 3
13 22 9 0 3
13 23 11 0 3
14 1 4 7 1
14 12 7 0 2
15 1 5 7 1
15 12 8 0 2
16 1 6 7 1
16 13 9 0 2
17 2 4 7 1
17 12 7 0 2
18 2 6 7 1
18 13 8 0 2
19 2 6 7 1
19 13 10 0 2
20 3 5 7 1
20 12 7 0 2
21 3 4 7 1
21 12 8 0 2
22 3 6 3 1
22 13 9 0 2
23 3 5 3 1
23 13 11 0 2

 

 

해를 구하는 프로그램

 

PROGRAM mt_animal_model_solve

! program name - multiple trait animal_model_solve
! programmer - Park Byoungho
! usage - mt_animal_model_solve < mt_animal_model_solve.par
! purpose : read sorted_lhs.dat and rhs.dat which are made in mt_animal_model_setup.exe and are sorted
! find a solution using Gauss-Seidel iteration
! Date - 2009.6.6
! update - none

USE gi

IMPLICIT NONE

! data dictionary
INTEGER :: no_of_effects ! number of effects
INTEGER :: pre_eq_no ! previous equation number
INTEGER :: cur_eq_no ! current equation number
INTEGER, ALLOCATABLE :: loc_of_off(:) ! location of off-diagonal in the equation
INTEGER :: class_code ! to identify fixed and random : 1-fixed, 2-random
INTEGER, ALLOCATABLE :: level_of_effects(:) ! level of each fixed and animal effect
REAL, ALLOCATABLE :: diagonal(:,:) ! diagonal element of lhs
REAL, ALLOCATABLE :: i_diagonal(:,:) ! inverse of diagonal element of lhs
INTEGER :: no_of_eq ! number of equations
INTEGER :: loc_of_animal ! starting location of animal effect
INTEGER :: no_of_trait ! number of trait

CHARACTER(LEN = 256) :: vcv_filename ! variance-covariance parameter file name

REAL(KIND = 8), ALLOCATABLE :: gvcv(:) ! upper diagonal part of genetic variance-covariance, (no_of_trait)*(no_of_trait + 1) / 2 dimension
REAL(KIND = 8), ALLOCATABLE :: gvcv_m(:,:) ! genetic variance-covariance matrix
REAL(KIND = 8), ALLOCATABLE :: rvcv(:) ! upper diagonal part of residual variance-covariance, (no_of_trait)*(no_of_trait + 1) / 2 dimension
REAL(KIND = 8), ALLOCATABLE :: temp_rvcv(:) ! temporary rvcv
REAL(KIND = 8), ALLOCATABLE :: rvcv_tc_m(:,:,:) ! rvcv according to trait combination and then inverse
INTEGER, ALLOCATABLE :: trait_combi(:) ! trait combination
INTEGER :: no_of_ud ! number of upper digonal
INTEGER :: no_of_tc ! number of trait combination
INTEGER :: tc_no ! trait combination number

INTEGER, ALLOCATABLE :: lhs(:,:) ! left-hand side
REAL(KIND = 8), ALLOCATABLE :: rhs(:,:)! left-hand side and right-hand side, x'y
REAL(KIND = 8), ALLOCATABLE :: temp_rhs(:) ! right-hand side one line
REAL(KIND = 8), ALLOCATABLE :: solutions(:,:) ! solutions
REAL(KIND = 8), ALLOCATABLE :: pre_sol(:) ! previous solution
INTEGER :: iteration ! iteration
REAL :: epsilon ! sum of square between old and new solutions
INTEGER :: i,j,k ! loop
INTEGER :: status ! i/o status
CHARACTER(LEN = 40) :: error_msg ! error message
INTEGER, PARAMETER :: MAX_ITER = 50 ! maximum number of iteration
REAL, PARAMETER :: criteria = 1.E-8 ! criteria for stopping
INTEGER :: no_of_lhs ! number of lhs lines
INTEGER :: data_type ! data type : 1-fixed 2-6:animal

INTEGER, DIMENSION(99) :: iw ! for inverse
REAL(KIND = 8) :: z ! for inverse
INTEGER :: mr ! for inverse

! get the number of effectss
WRITE (*,*) "Enter the number of effects including fixed and animal effects: "
READ (*,*) no_of_effects

ALLOCATE(loc_of_off(no_of_effects - 1), level_of_effects(no_of_effects))

! get the number of equations
WRITE (*,*) "Enter the level of each fixed and animal effect : "
READ (*,*) level_of_effects

! get the number of traits
WRITE (*,*) "Enter the number of traits : "
READ (*,*) no_of_trait

ALLOCATE(gvcv_m(no_of_trait,no_of_trait))
ALLOCATE(diagonal(no_of_trait,no_of_trait))
ALLOCATE(i_diagonal(no_of_trait,no_of_trait))
ALLOCATE(temp_rhs(no_of_trait))
ALLOCATE(pre_sol(no_of_trait)) ! previous solution

no_of_eq = SUM(level_of_effects) ! total number of equations
loc_of_animal = SUM(level_of_effects(1:no_of_effects-1)) + 1 ! (sum of level of fixed effects) + 1
ALLOCATE(rhs(no_of_eq ,no_of_trait), solutions(no_of_eq ,no_of_trait))

no_of_ud = (no_of_trait)*(no_of_trait + 1) / 2 ! number of upper diagonal of variance-covariance
no_of_tc = 0 ! number of trait combination
DO i = 1, no_of_trait
no_of_tc = no_of_tc + 2 ** (i - 1)
END DO

ALLOCATE(gvcv(no_of_ud)) ! genetic variance-covariance
ALLOCATE(rvcv(no_of_ud)) ! residual variance-covariance
ALLOCATE(temp_rvcv(no_of_ud)) ! temporary residual variance-covariance
ALLOCATE(rvcv_tc_m(no_of_tc,no_of_trait,no_of_trait))
ALLOCATE(trait_combi(no_of_ud))

! initialiaze solutions to zero
solutions = 0

! get the genetic and residual variance-covariance file name
WRITE (*,*) "Enter the genetic and residual variance-covariance file name "
READ (*,*) vcv_filename

! open parameter file
OPEN(UNIT = 21, FILE = vcv_filename, STATUS = 'OLD', ACTION = 'READ', IOSTAT = status, IOMSG = error_msg)

IF (status /= 0) THEN ! file open failed
WRITE (*,'(1X, A, A)') 'parameter file open failed -- error message : ', error_msg
STOP
END IF

! genetic variance-covariance
READ(UNIT = 21, FMT = *) gvcv

! inverse of genetic variance-covariance
z = 0.0
CALL DJNVHF(gvcv,no_of_trait,iw,z,mr)

DO i = 1, no_of_trait
DO j = 1, no_of_trait
gvcv_m(i,j) = gvcv(IHMSSF(i,j,no_of_trait))
END DO
END DO

! residual variance-covariance
READ(UNIT = 21, FMT = *) rvcv

! read trait combination and inverse rvcv
DO k = 1, no_of_tc
READ(UNIT = 21, FMT = *) trait_combi
temp_rvcv = rvcv * trait_combi
z = 0.0
CALL DJNVHF(temp_rvcv,no_of_trait,iw,z,mr)
DO i = 1, no_of_trait
DO j = 1, no_of_trait
rvcv_tc_m(k,i,j) = temp_rvcv(IHMSSF(i,j,no_of_trait))
END DO
END DO
END DO

!open left-hand side file
OPEN(UNIT = 11, FILE = 'sorted_lhs.dat', STATUS = 'OLD', ACTION = 'READ', IOSTAT = status, IOMSG = error_msg)
IF (status /= 0) THEN ! file open failed
WRITE (*,'(1X, A, A)') 'Sorted lhs file open failed -- error message : ', error_msg
STOP
END IF

! count the lines of lhs
no_of_lhs = 0
DO
READ(UNIT = 11, FMT = *, IOSTAT = status)
IF (status /= 0) EXIT ! reach end of file
no_of_lhs = no_of_lhs + 1
END DO

ALLOCATE(lhs(no_of_lhs,no_of_effects + 2))

! store lhs to array
REWIND(11)
DO i = 1, no_of_lhs
READ(UNIT = 11, FMT = *, IOSTAT = status) lhs(i,:)
IF (status /= 0) EXIT ! reach end of file
END DO
CLOSE(11)

!open right-hand side file
OPEN(UNIT = 11, FILE = 'rhs.dat', STATUS = 'OLD', ACTION = 'READ', IOSTAT = status, IOMSG = error_msg)

IF (status /= 0) THEN ! file open failed
WRITE (*,'(1X, A, A)') 'RHS file open failed -- error message : ', error_msg
STOP
END IF
! read and store right-hand side
DO i = 1, no_of_eq
READ(UNIT = 11,FMT = *, IOSTAT = status) rhs(i,:)
IF (status /= 0) EXIT ! reach end of file
END DO
close(11)
! iteration
DO iteration = 1, MAX_ITER

pre_eq_no = 1
diagonal = 0.0
temp_rhs = rhs(1,:)
epsilon = 0.0

! lhs
DO i = 1, no_of_lhs
cur_eq_no = lhs(i,1)
loc_of_off = lhs(i,2:no_of_effects) ! location of off diagonal
tc_no = lhs(i,no_of_effects + 1) ! trati combination number
data_type = lhs(i,no_of_effects + 2)

! nrew equation, or calculate solution
IF (pre_eq_no /= cur_eq_no) THEN
pre_sol = solutions(pre_eq_no,:) ! store old solution
CALL geninv(diagonal, i_diagonal)
solutions(pre_eq_no,:) = MATMUL(i_diagonal, temp_rhs) ! get a new solution
epsilon = epsilon + SUM((pre_sol - solutions(pre_eq_no,:)) * (pre_sol - solutions(pre_eq_no,:))) ! calculate sum of squares of difference between old solution and new solution

diagonal = 0.0
temp_rhs = rhs(cur_eq_no,:)
END IF

! process off-diagnoal
SELECT CASE(data_type)
CASE(1) ! fixed effect
temp_rhs = temp_rhs - MATMUL(rvcv_tc_m(tc_no,:,:), SUM(solutions(loc_of_off,:),1))
diagonal = diagonal + rvcv_tc_m(tc_no,:,:)
pre_eq_no = cur_eq_no
CASE(2) ! animal effect, both parents are known, animal is diagonal
temp_rhs = temp_rhs - (-0.5) * 2.0 * MATMUL(gvcv_m, solutions(loc_of_off(1),:)) &
- (-0.5) * 2.0 * MATMUL(gvcv_m, solutions(loc_of_off(2),:))
diagonal = diagonal + 2.0 * gvcv_m
pre_eq_no = cur_eq_no
CASE(3) ! animal effect, both parents are known, sire or dam is diagonal
temp_rhs = temp_rhs - (-0.5) * 2.0 * MATMUL(gvcv_m, solutions(loc_of_off(1),:)) &
- 0.25 * 2.0 * MATMUL(gvcv_m, solutions(loc_of_off(2),:))
diagonal = diagonal + 0.25 * 2.0 * gvcv_m
pre_eq_no = cur_eq_no
CASE(4) ! animal effect, one parent is known, animal is diagonal
temp_rhs = temp_rhs - (-0.5) * (4.0/3.0) * MATMUL(gvcv_m, solutions(loc_of_off(1),:))
diagonal = diagonal + (4.0/3.0) * gvcv_m
pre_eq_no = cur_eq_no
CASE(5) ! animal effect, one parents is known, parent is diagonal
temp_rhs = temp_rhs - (-0.5) * (4.0/3.0) * MATMUL(gvcv_m, solutions(loc_of_off(1),:))
diagonal = diagonal + 0.25 * (4.0/3.0) * gvcv_m
pre_eq_no = cur_eq_no
CASE(6) ! animal effect, no parents is known
diagonal = diagonal + 1.0 * gvcv_m
pre_eq_no = cur_eq_no
END SELECT
END DO ! one iteration ends - end of array

! calculate last solution
pre_sol = solutions(pre_eq_no,:) ! store old solution
CALL geninv(diagonal, i_diagonal)
solutions(pre_eq_no,:) = MATMUL(i_diagonal, temp_rhs) ! get a new solution
epsilon = epsilon + SUM((pre_sol - solutions(pre_eq_no,:)) * (pre_sol - solutions(pre_eq_no,:))) ! calculate sum of squares of difference between old solution and new solution

! write iteration number and epsilon
epsilon = epsilon / (no_of_eq * no_of_trait)
WRITE(*,*) 'iteration = ', iteration , ', epsilon = ', epsilon
IF (epsilon < criteria) THEN
EXIT
END IF

END DO ! iteration

! open file for writing solutions
OPEN(UNIT = 13, FILE = 'sol.dat', STATUS = 'REPLACE', ACTION = 'WRITE', IOSTAT = status, IOMSG = error_msg)
DO i = 1, no_of_eq
WRITE(UNIT = 13, FMT = *) solutions(i,:)
END DO
CLOSE(13)

END PROGRAM mt_animal_model_solve
위 소스를 mt_animal_model_solve.f95로 저장

 

위에서 역행렬을 구하는데 필요한 소스

 

MODULE gi CONTAINS

SUBROUTINE geninv(x,gix) REAL, DIMENSION(:,:), INTENT(IN) :: x ! square matrix to inverse
REAL, DIMENSION(:,:), INTENT(OUT) :: gix ! square matirx after inverse
INTEGER :: i, j, no_of_eq, size_x, mr
REAL(KIND = 8), ALLOCATABLE :: upperd(:)
INTEGER, DIMENSION(99) :: iw
REAL(KIND = 8) :: z no_of_eq = SIZE(x,1) ALLOCATE(upperd((no_of_eq*(no_of_eq +1)/2)))

upperd = 0

DO i = 1, no_of_eq
DO j = 1, no_of_eq
upperd(ihmssf(i,j,no_of_eq))=x(i,j)
END DO
END DO !WRITE(*,*) 'upperd = ', upperd z = 0.0
CALL DJNVHF(upperd,no_of_eq,iw,z,mr) DO i = 1, no_of_eq
DO j = 1, no_of_eq
gix(i,j) = upperd(ihmssf(i,j,no_of_eq))
END DO
END DO END SUBROUTINE geninv

 

SUBROUTINE DJNVHF(A,N,IS,Z,IR)
!
! MATRIX INVERSION SUBROUTINE
!
! 'A' IS A HALF-STORED MATRIX TO BE INVERTED. THE INVERTED
! RESULT IS RETURNED IN 'A'.
! 'N' IS THE ORDER OF THE MATRIX TO BE INVERTED.
! 'IS' IS A WORK VECTOR.
! 'Z' IS SET TO ZERO ON CALLS TO THIS SUBROUTINE FROM THE BEEF
! SIRE MONITORING SYSTEM.
! 'IR' RETURNS THE RANK OF THE MATRIX.
!
! THIS SUBROUTINE CALLS THE SUBROUTINE DREARN.
!
DIMENSION A(1),IS(1)
REAL*8 A,BIG,Z,RECIP,R
201 NP=(N*(N+1))/2
IF(N-1) 21,23,96
21 WRITE(6,22)
22 FORMAT(' N.LT.1')
23 IF(A(1)) 25,24,25
24 IR=0
26 RETURN
25 IR=1
A(1)=1.D0/A(1)
RETURN
96 DO 50 L=1,N
BIG=0.D1
DO 2 I=L,N
II=-(I*(I-3))/2+N*(I-1)
IF(DABS(A(II))-DABS(BIG)) 2,2,1
1 IS(L)=I
BIG=A(II)
2 CONTINUE
IF(DABS(BIG)-Z) 63,63,62
63 IR=L-1
IF(IR) 26,26,98
98 DO 64 I=1,L
LI=-(I*(I-1))/2+N*(I-1)
DO 64 J=L,N
LIJ=LI+J
64 A(LIJ)=-0.D1
IF(L-N) 66,69,67
67 STOP 67
66 LP1=L+1
DO 68 I=LP1,N
LI=-(I*(I-1))/2+N*(I-1)
DO 68 J=I,N
LIJ=LI+J
68 A(LIJ)=-0.D1
69 IF(IR-2) 56,27,27
56 CALL DREARN(A,N,IS,1)
DO 57 I=1,NP
57 A(I)=-A(I)
RETURN
27 LP1=L+1
DO 65 I=2,L
J=LP1-I
65 CALL DREARN(A,N,IS,J)
DO 101 I=1,NP
101 A(I)=-A(I)
RETURN
62 CALL DREARN(A,N,IS,L)
28 K3=N*(L-1)-(L*(L-3))/2
RECIP=1./A(K3)
A(K3)=-RECIP
DO 50 I=1,N
IF(I-L) 6,50,7
6 K11=N*(I-1)-(I*(I-3)/2)
K1=K11+L-I
GO TO 8
7 K11=N*(I-1)-(I*(I-3))/2
K1=K3+I-L
8 R=A(K1)*RECIP
301 DO 12 J=I,N
K4=K11+J-I
IF(J-L) 10,12,11
10 K5=N*(J-1)-(J*(J-3)/2)+L-J
GO TO 14
11 K5=K3+J-L
14 A(K4)=A(K4)-R*A(K5)
12 CONTINUE
A(K1)=R
50 CONTINUE
NP=(N*(N+1))/2
DO 100 I=1,NP
100 A(I)=-A(I)
NP1=N+1
DO 61 I=2,N
L=NP1-I
61 CALL DREARN(A,N,IS,L)
IR=N
RETURN
END SUBROUTINE

SUBROUTINE DREARN(A,N,IS,L)
!
! THIS IS A SUPPORT SUBROUTINE CALLED BY THE MATRIX INVERSION
! SUBROUTINE DJNVHF.
!
DIMENSION A(1),IS(1)
DOUBLE PRECISION A,SAVE
ISL=IS(L)
IF(ISL-L) 3,28,4
3 STOP 3
4 LM1=L-1
IF(LM1) 22,22,5
5 DO 21 I=1,LM1
IL=-(I*(I-1))/2+N*(I-1)+L
IISL=IL-L+ISL
SAVE=A(IL)
A(IL)=A(IISL)
21 A(IISL)=SAVE
22 LP1=L+1
ISLM1=ISL-1
IF(LP1-ISLM1) 23,23,25
23 DO 24 I=LP1,ISLM1
LI=-(L*(L-1))/2+N*(L-1)+I
IISL=-(I*(I-1))/2+N*(I-1)+ISL
SAVE=A(LI)
A(LI)=A(IISL)
24 A(IISL)=SAVE
25 ISLP1=ISL+1
IF(ISLP1-N) 26,26,38
26 DO 27 I=ISLP1,N
ISLI=-(ISL*ISLM1)/2+N*ISLM1+I
LI=-(L*LM1)/2+N*LM1+I
SAVE=A(ISLI)
A(ISLI)=A(LI)
27 A(LI)=SAVE
38 LL=-(L*(L-3))/2+N*(L-1)
ISLISL=-(ISL*(ISL-3))/2+ISLM1*N
SAVE=A(LL)
A(LL)=A(ISLISL)
A(ISLISL)=SAVE
28 RETURN
END SUBROUTINE FUNCTION IHMSSF(IR,IC,N)

JR=IR
JC=IC
IF(IR-IC)40,40,41
41 JC=IR
JR=IC
40 NR=JR-1
KK=(JR*NR)/2
KK=KK+NR*(N-NR)
KK=KK+JC-NR
IHMSSF=KK
RETURN
END FUNCTION

END MODULE gi

 

위 소스를 inverse.f95로 저장

 

 

컴파일 명령 : gfortran inverse.f95 mt_animal_model_solve.f95 -o mt_animal_model_solve.exe

 

 

실행을 편하게 하기 위한 파라미터 파일

 

3
3 3 17
3
mt_animal_model_par.dat

 

위 파일을 mt_animal_model_solve.par로 저장

 

 

위의 mt_animal_model_par.dat 파일 내용

 

6 7 10 20 250 2000
15 12 90 85 200 4500
1 0 0 0 0 0
0 0 0 1 0 0
1 1 0 1 0 0
0 0 0 0 0 1
1 0 1 0 0 1
0 0 0 1 1 1
1 1 1 1 1 1

 

 

 

구한 해(sol.dat)

 

24.125706755520547 74.268295670354391 527.08028644687056
25.441921398396907 82.593543141471272 566.55416969797432
24.220001408184167 85.631965785314890 578.04662150918193
3.7632704190436970 5.8866115785656232 35.057122362870764
-4.4138517712974776 -4.8152560853272135 -4.7788844942106072
0.48793465747628662 -0.80363360304065923 -22.617193956001003
0.22614643375046239 0.14233968977258890 -3.3275142494486567
-0.87479454301318160 -0.31552433737757035 8.3615355467185211
0.13023931959343374 0.12005625102548459 0.96781996466511400
0.60757599137083063 -3.91612018471323500E-002 -9.8066393081367575
-8.91080747632977399E-002 9.23043298191184025E-002 3.8044053103661692
8.91798161271940698E-002 -9.21895471626397456E-002 -3.8039798436513075
-8.90043748020188791E-002 9.22574144874573204E-002 3.8029771193220707
0.15130746719018148 0.29818052929107419 5.22709065466071196E-002
-0.40596247329896329 -4.52885376604869327E-002 6.3882628792527765
4.01384376810988222E-002 -0.32551976603110727 -5.3422089438072096
0.28788631511787954 0.14515758579948962 -2.5974880821666866
-1.2196884076215417 -0.19255006665593399 14.920597638143102
0.86685978756745230 -1.26122030640144749E-002 -12.808446372532977
0.25993496240125380 -0.22578722906801330 -11.479592485088428
-0.51666770389439098 -0.59703505012241087 -2.3074453032995699
0.13132829565375642 0.65789677148818060 11.080980775099574
-0.17816500450516137 0.18458251313510871 7.6080572099350157

 

컴파일 및 실행 화면

 

 

1244219139_mt_animal_model_solve.f95

다운로드

1244219139_mt_animal_model_solve.f95

1244219139_mt_animal_model_solve.par

다운로드

1244219139_mt_animal_model_solve.par

1244219139_inverse.f95

다운로드

1244219139_inverse.f95

 

1244219139_mt_animal.zip

다운로드

1244219139_mt_animal.zip