same model, different level, animal model 3
(same model, different level, animal model 2는 실패)
same model, different level, animal model 1에서는CASE가 너무 많았습니다.이것이 형질이 늘어나면 좀 심각해 집니다.형질이 3개이면 30개, 5개이면 100개 이상으로 늘어납니다.그래서 좀 다르게 생각을 해서 CASE 문을 줄여 봤습니다.rhs를 adjust할 때, 빼기만 했는데 빼고 더하는 것으로 생각을 바꾸었습니다.예를 들어 2 fixed effects, 2 trait2 일 경우LHS의 한 라인이 다음과 같습니다.1 5 8 2 6 81과 2는 첫째 fixed effect의 서로 다른 level이고5와 6은 두째 fixed effect의 서로 다른 level이고8은 animal입니다.지금까지는 위를1 5 2 6 8로 만들고, 1은 diagonal, 나머지는 off-diagonal로서 솔루션을 이용하여 rhs를 adjust했습니다.
이것을 다음과 같이 바꾸었습니다.1 1 5 8 2 6 8첫째 1은 diagonal이고, 나머지는 none-zero 원소들입니다.즉 1 5 8 2 6 8에 대한 솔루션을 rhs에서 빼주고, 첫째 1의 solution을 더해 주면 위와 같은 것이 됩니다.이걸 기초로 해서 방정식을 풀었습니다.기본적으로 equation 별로 그러니까, diag_ele의 차수는 no_of_trait X no_of_trait가 되도록 풀었습니다.즉 block diagonal이 아니지만 그런식으로 풀었습니다. 왜냐하면 animal쪽으로 가면 block diagonal이 되니까요.
자 보시죠.
사용한 자료
8 0 0 1 6 0 9570
9 1 5 1 6 9219 10530
10 2 6 4 7 10025 11490
11 2 6 3 7 7230 8600
12 1 7 0 0 8121 0
1st col : cowid2nd col : age 1 - fixed3rd col : contemporary group - fixed4th col : age 2 - fixed5th col : contemporary group - fixedmilk1 : observation 1milk2 : observation 2
위 데이터를 dairy2.dat로 저장
혈통 자료
8 13 16
9 14 16
10 15 17
11 13 16
12 14 17
13 0 0
14 0 0
15 0 0
16 0 0
17 0 0
animal, dam, sire
위 데이터를 pedi2.dat로 저장
파라미터 파일
30000 1150 31000
8000 4000 8300
1 0 0
0 0 1
1 1 1
1st row : residual variance - covariance 2nd row : genetic variance - covariance others : trait combination
위 데이터를 vcvtrt.par로 저장
프로그램
PROGRAM mt_dl_animal_model3_setup
! program name - multiple trait, different fixed level, animal model setup
! same model but different fixed level
! programmer - Park Byoungho
! usage - mt_df_animal_model3_setup
! purpose : read data file and write left-hand side and right-hand side
! Date - 2009.7.7.
! update - none
USE gi
IMPLICIT NONE
! data dictionary
INTEGER, PARAMETER :: no_of_trait = 2 ! number of trait
INTEGER, PARAMETER :: no_of_fixed = 2 ! number of fixed effects
INTEGER, PARAMETER :: no_of_eq = 17 ! number of equation CHARACTER(LEN = 256) :: data_filename ! data file name
CHARACTER(LEN = 256) :: pedi_filename ! pedigree file name
CHARACTER(LEN = 256) :: par_filename ! parameter file name INTEGER :: no_of_ud ! number of upper digonal
INTEGER :: no_of_tc ! number of trait combination
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 :: 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(:,:) ! rvcv according to trait combination and then inverse
REAL(KIND = 8), ALLOCATABLE :: rvcv_tc_m(:,:,:) ! rvcv matrix according to trait combination and then inverse
INTEGER, ALLOCATABLE :: trait_combi(:) ! trait combination
INTEGER :: tc_no ! trait combination number INTEGER, ALLOCATABLE :: temp_effects(:) ! temporary animal and fixed effect
REAL(KIND = 8) :: milk1, milk2, r1, r2 ! data of trait for each line
INTEGER, ALLOCATABLE :: effects(:,:) ! array for fixed and animal effects
REAL(KINd = 8), ALLOCATABLE :: observations(:) ! array for observations
REAL(KINd = 8), ALLOCATABLE :: observations_mbr(:) ! array for observations multiplied by rvcv
INTEGER :: animal, sire, dam ! animal, sire, dam name for pedigree file
REAL(KIND = 8), ALLOCATABLE :: xpy(:,:) ! right-hand side
INTEGER, ALLOCATABLE :: nobs(:,:) ! number of observations
INTEGER :: i, j, k, l ! do loop
REAL :: temp_value ! temporary storate for swap
INTEGER, DIMENSION(99) :: iw ! for inverse
REAL(KIND = 8) :: z ! for inverse
INTEGER :: mr ! for inverse
INTEGER :: status ! I/O status
CHARACTER(LEN = 40) :: error_msg ! error message
INTEGER :: zero = 0
! get data file name
data_filename = 'dairy2.dat'
! get pedigree file name
pedi_filename = 'pedi2.dat'
! get parameter file name
par_filename = 'vcvtrt.par' 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(no_of_tc, no_of_ud)) ! residual variance-covariance according to trait combination
ALLOCATE(rvcv_tc_m(no_of_tc, no_of_trait, no_of_trait)) ! residual variance-covariance matrix according to trait combination
ALLOCATE(trait_combi(no_of_ud)) ! trait combination
ALLOCATE(xpy(no_of_eq, no_of_trait)) ! right hand side
ALLOCATE(nobs(no_of_eq, no_of_trait)) ! number of observations
ALLOCATE(temp_effects(no_of_trait * no_of_fixed + 1)) ! temporary array for fixed and animal effects
ALLOCATE(effects(no_of_trait, no_of_fixed + 1)) ! array for processing fixed and animal effects
ALLOCATE(observations(no_of_trait)) ! array for observations
ALLOCATE(observations_mbr(no_of_trait)) ! array for observations multiplied by rvcv
xpy = 0
nobs = 0 ! open file for data processing procedure
OPEN(UNIT = 31, FILE = 'ongoing.dat', STATUS = 'REPLACE', ACTION = 'WRITE')
! open parameter file
OPEN(UNIT = 21, FILE = par_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 ! residual variance-covariance
READ(UNIT = 21, FMT = *) rvcv
WRITE(31,*) 'residual variance-covariance = ', rvcv
! genetic variance-covariance
READ(UNIT = 21, FMT = *) gvcv
WRITE(31,*) 'genetic variance-covariance = ', gvcv
! trait combination
DO k = 1, no_of_tc
READ(UNIT = 21, FMT = *) trait_combi
WRITE(31,*) 'trait combination ', k , '=', trait_combi
temp_rvcv = rvcv * trait_combi
WRITE(31,*) 'trati combination residual vcv', temp_rvcv
z = 0.0
CALL DJNVHF(temp_rvcv,no_of_trait,iw,z,mr) ! inverse of upper diagonal matrix(residual)
WRITE(31,*) 'inverse of trait combination rvcv', temp_rvcv
rvcv_tc(k,:) = temp_rvcv ! sotre the inverse of upper diagonal matrix(residual)
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 data file
OPEN(UNIT = 11, FILE = data_filename, STATUS = 'OLD', ACTION = 'READ', IOSTAT = status, IOMSG = error_msg)
IF (status /= 0) THEN ! file open failed
WRITE (*,'(1X, A, A)') 'Data file open failed -- error message : ', error_msg
STOP
END IF ! open data file for writing left-hand side
OPEN(UNIT = 12, FILE = 'lhs.dat', STATUS = 'REPLACE', ACTION = 'WRITE') ! read each line
DO
READ(UNIT = 11, FMT = *, IOSTAT = status) temp_effects, observations
IF (status /= 0) EXIT ! end of file ! copy temp_effects(fixed) to effects
DO i = 1, no_of_trait
effects(i,1:no_of_fixed) = temp_effects((i-1) * no_of_fixed + 1 : i * no_of_fixed)
END DO
! copy temp_effects(animal) to effects
DO i = 1, no_of_trait
effects(i,no_of_fixed + 1) = temp_effects(no_of_trait * no_of_fixed + 1)
END DO
! find trait combination number
tc_no = 0
DO i = 1, no_of_trait
IF (observations(i) /= 0.0) tc_no = tc_no + 2 ** (i - 1)
END DO observations_mbr = MATMUL(rvcv_tc_m(tc_no,:,:), observations) ! left hand side
DO i = 1, no_of_trait
IF (effects(i,1) /= 0) THEN
DO j = 1, no_of_fixed + 1
IF (j < no_of_fixed + 1) THEN
WRITE(UNIT = 12, FMT = *) effects(i,j), i, ((effects(k,l),l = 1, no_of_fixed + 1), k = 1, no_of_trait), '1', tc_no
ELSE
WRITE(UNIT = 12, FMT = *) effects(i,j), i, ((effects(k,l),l = 1, no_of_fixed + 1), k = 1, no_of_trait), '2', tc_no
END IF
END DO
END IF
END DO
! right hand side
DO i = 1, no_of_trait
IF (effects(i,1) /= 0) THEN
DO j = 1, no_of_fixed + 1
xpy(effects(i,j),i) = xpy(effects(i,j),i) + observations_mbr(i)
nobs(effects(i,j),i) = nobs(effects(i,j),i) + 1
END DO
END IF
END DO END DO CLOSE(11) ! close data file
! open pedigree file
OPEN(UNIT = 13, FILE = pedi_filename, STATUS = 'OLD', ACTION = 'READ', IOSTAT = status, IOMSG = error_msg)
IF (status /= 0) THEN ! file open failed
WRITE (*,'(1X, A, A)') 'pedigree file open failed -- error message : ', error_msg
STOP
END IF ! read each line
DO
READ(UNIT = 13, FMT = *, IOSTAT = status) animal, sire, dam
IF (status /= 0) EXIT ! end of file ! write animal effect
! DO i = 1, no_of_trait
IF ( sire /= 0 .AND. dam /= 0) THEN
WRITE(UNIT = 12, FMT = *) animal, zero, animal, sire, dam, animal, sire, dam, '1001', zero
WRITE(UNIT = 12, FMT = *) sire, zero, animal, dam, sire, animal, sire, dam, '1002', zero
WRITE(UNIT = 12, FMT = *) dam, zero, animal, sire, dam, animal, sire, dam, '1002', zero
ELSE IF (sire /= 0 .AND. dam == 0) THEN
WRITE(UNIT = 12, FMT = *) animal, zero, animal, sire, dam, animal, sire, dam, '1003', zero
WRITE(UNIT = 12, FMT = *) sire, zero, animal, sire, dam, animal, sire, dam, '1004', zero
ELSE IF (sire == 0 .AND. dam /= 0) THEN
WRITE(UNIT = 12, FMT = *) animal, zero, animal, dam, sire, animal, dam, sire, '1003', zero
WRITE(UNIT = 12, FMT = *) dam, zero, animal, dam, sire, animal, dam, sire, '1004', zero
ELSE
WRITE(UNIT = 12, FMT = *) animal, zero, animal, sire, dam, animal, sire, dam, '1005', zero
END IF
! END DO
END DO ! end of reading datas CLOSE(13) ! close animal file
CLOSE(12) ! close file for left-hand side ! open file for right hand side
OPEN(UNIT = 14, FILE = 'rhs.dat', STATUS = 'REPLACE', ACTION = 'WRITE') ! write xpy
DO i = 1, no_of_eq
WRITE(UNIT = 14, FMT = *) xpy(i,:), nobs(i,:)
END DO
CLOSE(14) ! close file for right-hand sideEND PROGRAM mt_dl_animal_model3_setup
컴파일 : gfortran inverse.f95 mt_dl_animal_model3_setup.f95 -o mt_dl_animal_model3_setup
위 프로그램을 컴파일할 때 필요한 inverse.f95
MODULE gi CONTAINS
SUBROUTINE geninv(x,gix)
! subroutine name - genive, generalized inverse
! programmer - Park Byoungho
! usage - CALL genive(x, gix)
! purpose - to find out a generalized inverse of matrix using ihmssf and djnvhf
! ihmssf = matrix <-> vector using upper diagonal element
! djnvhf = solve generalized inverse of vector
! Date - 2009.7.7.
! last update - 2009.7.10. : precision of x and gix changes single precision to double precision REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: x ! square matrix to inverse
REAL(KIND=8), 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로 저장
프로그램으로 구한 LHS
1 2 0 0 8 1 6 8 1 2
6 2 0 0 8 1 6 8 1 2
8 2 0 0 8 1 6 8 2 2
1 1 1 5 9 1 6 9 1 3
5 1 1 5 9 1 6 9 1 3
9 1 1 5 9 1 6 9 2 3
1 2 1 5 9 1 6 9 1 3
6 2 1 5 9 1 6 9 1 3
9 2 1 5 9 1 6 9 2 3
2 1 2 6 10 4 7 10 1 3
6 1 2 6 10 4 7 10 1 3
10 1 2 6 10 4 7 10 2 3
4 2 2 6 10 4 7 10 1 3
7 2 2 6 10 4 7 10 1 3
10 2 2 6 10 4 7 10 2 3
2 1 2 6 11 3 7 11 1 3
6 1 2 6 11 3 7 11 1 3
11 1 2 6 11 3 7 11 2 3
3 2 2 6 11 3 7 11 1 3
7 2 2 6 11 3 7 11 1 3
11 2 2 6 11 3 7 11 2 3
1 1 1 7 12 0 0 12 1 1
7 1 1 7 12 0 0 12 1 1
12 1 1 7 12 0 0 12 2 1
8 0 8 13 16 8 13 16 1001 0
13 0 8 16 13 8 13 16 1002 0
16 0 8 13 16 8 13 16 1002 0
9 0 9 14 16 9 14 16 1001 0
14 0 9 16 14 9 14 16 1002 0
16 0 9 14 16 9 14 16 1002 0
10 0 10 15 17 10 15 17 1001 0
15 0 10 17 15 10 15 17 1002 0
17 0 10 15 17 10 15 17 1002 0
11 0 11 13 16 11 13 16 1001 0
13 0 11 16 13 11 13 16 1002 0
16 0 11 13 16 11 13 16 1002 0
12 0 12 14 17 12 14 17 1001 0
14 0 12 17 14 12 14 17 1002 0
17 0 12 14 17 12 14 17 1002 0
13 0 13 0 0 13 0 0 1005 0
14 0 14 0 0 14 0 0 1005 0
15 0 15 0 0 15 0 0 1005 0
16 0 16 0 0 16 0 0 1005 0
17 0 17 0 0 17 0 0 1005 0
프로그램으로 구한 RHS
0.56539810563947124 0.63745474769401966 2 2
0.55110789267533677 0.0000000000000000 2 0
0.0000000000000000 0.26886136468257282 0 1
0.0000000000000000 0.35875882639560025 0 1
0.29469810563947119 0.0000000000000000 1 0
0.55110789267533677 0.63745474769401966 2 2
0.27070000000000000 0.62762019107817313 1 2
0.0000000000000000 0.30870967741935484 0 1
0.29469810563947119 0.32874507027466482 1 1
0.32041424498816867 0.35875882639560025 1 1
0.23069364768716805 0.26886136468257282 1 1
0.27070000000000000 0.0000000000000000 1 0
0.0000000000000000 0.0000000000000000 0 0
0.0000000000000000 0.0000000000000000 0 0
0.0000000000000000 0.0000000000000000 0 0
0.0000000000000000 0.0000000000000000 0 0
0.0000000000000000 0.0000000000000000 0 0
sort -n -o sorted_lhs.dat lhs.dat
정렬된 sorted_lhs.dat
1 1 1 5 9 1 6 9 1 3
1 1 1 7 12 0 0 12 1 1
1 2 0 0 8 1 6 8 1 2
1 2 1 5 9 1 6 9 1 3
2 1 2 6 10 4 7 10 1 3
2 1 2 6 11 3 7 11 1 3
3 2 2 6 11 3 7 11 1 3
4 2 2 6 10 4 7 10 1 3
5 1 1 5 9 1 6 9 1 3
6 1 2 6 10 4 7 10 1 3
6 1 2 6 11 3 7 11 1 3
6 2 0 0 8 1 6 8 1 2
6 2 1 5 9 1 6 9 1 3
7 1 1 7 12 0 0 12 1 1
7 2 2 6 10 4 7 10 1 3
7 2 2 6 11 3 7 11 1 3
8 0 8 13 16 8 13 16 1001 0
8 2 0 0 8 1 6 8 2 2
9 0 9 14 16 9 14 16 1001 0
9 1 1 5 9 1 6 9 2 3
9 2 1 5 9 1 6 9 2 3
10 0 10 15 17 10 15 17 1001 0
10 1 2 6 10 4 7 10 2 3
10 2 2 6 10 4 7 10 2 3
11 0 11 13 16 11 13 16 1001 0
11 1 2 6 11 3 7 11 2 3
11 2 2 6 11 3 7 11 2 3
12 0 12 14 17 12 14 17 1001 0
12 1 1 7 12 0 0 12 2 1
13 0 8 16 13 8 13 16 1002 0
13 0 11 16 13 11 13 16 1002 0
13 0 13 0 0 13 0 0 1005 0
14 0 9 16 14 9 14 16 1002 0
14 0 12 17 14 12 14 17 1002 0
14 0 14 0 0 14 0 0 1005 0
15 0 10 17 15 10 15 17 1002 0
15 0 15 0 0 15 0 0 1005 0
16 0 8 13 16 8 13 16 1002 0
16 0 9 14 16 9 14 16 1002 0
16 0 11 13 16 11 13 16 1002 0
16 0 16 0 0 16 0 0 1005 0
17 0 10 15 17 10 15 17 1002 0
17 0 12 14 17 12 14 17 1002 0
17 0 17 0 0 17 0 0 1005 0
정렬된 sorted_lhs.dat와 rhs.dat를 이용하여 해를 구하는 프로그램
same model, different level, multiple trait animal model3 - slove
PROGRAM mt_dl_animal_model3_solve ! program name - multiple trait, different fixed level, animal_model_solve
! programmer - Park Byoungho
! usage - mt_dl_animal_model3_solve
! purpose : read sorted_lhs.dat and rhs.dat which are made in animal_model_setup.exe and are sorted
! find a solution using Gauss-Seidel iteration
! Date - 2009.7.7
! update - none USE gi IMPLICIT NONE
! data dictionary
INTEGER, PARAMETER :: no_of_trait = 2 ! number of trait
INTEGER, PARAMETER :: no_of_fixed = 2 ! number of fixed effects
INTEGER, PARAMETER :: no_of_eq = 17 ! number of equation INTEGER :: no_of_ud ! number of upper digonal
INTEGER :: no_of_tc ! number of trait combination 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(:,:)
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(:,:) ! rvcv according to trait combination and then inverse
REAL(KIND = 8), ALLOCATABLE :: rvcv_tc_m(:,:,:) ! rvcv matrix according to trait combination and then inverse
INTEGER, ALLOCATABLE :: trait_combi(:) ! trait combination
INTEGER :: tc_no ! trait combination number INTEGER, DIMENSION(99) :: iw ! for inverse
REAL(KIND = 8) :: z ! for inverse
INTEGER :: mr ! for inverse
INTEGER :: pre_eq_no ! previous equation number
INTEGER :: cur_eq_no ! current equation number
INTEGER :: pre_trait_no ! previous trait number
INTEGER :: cur_trait_no ! current trati number
INTEGER, ALLOCATABLE :: loc_of_nz(:,:) ! location of none-zero element
INTEGER, ALLOCATABLE :: loc_of_nz_a(:,:) ! location of none-zero element
INTEGER :: data_type ! data type : 1:fixed, 2-6:animal REAL(KIND = 8), ALLOCATABLE :: diag_ele(:,:) ! diagonal element of lhs
REAL(KIND = 8) :: inv_diag_ele(2,2) ! inverse of diagonal INTEGER :: no_of_lhs ! number of lhs lines
INTEGER, ALLOCATABLE :: lhs(:,:) ! left-hand side
REAL(KIND = 8), ALLOCATABLE :: rhs(:,:) ! right-hand side, x'y
REAL(KIND = 8), ALLOCATABLE :: temp_rhs(:) ! right-hand side one element
REAL(KIND = 8), ALLOCATABLE :: solutions(:,:) ! solutions
REAL(KIND = 8), ALLOCATABLE :: pre_sol(:) ! previous solution
REAL(KIND = 8), ALLOCATABLE :: sum_sol(:) ! sum solution INTEGER :: iteration ! iteration
REAL(KIND = 8) :: epsilon ! sum of squares 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(KIND = 8), PARAMETER :: criteria = 1.E-12 ! criteria for stopping 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(gvcv_m(no_of_trait, no_of_trait)) ! genetic variance-covariance
ALLOCATE(rvcv(no_of_ud)) ! residual variance-covariance
ALLOCATE(temp_rvcv(no_of_ud)) ! temporary residual variance-covariance
ALLOCATE(rvcv_tc(no_of_tc, no_of_ud)) ! residual variance-covariance according to trait combination
ALLOCATE(rvcv_tc_m(no_of_tc, no_of_trait, no_of_trait)) ! residual variance-covariance matrix according to trait combination
ALLOCATE(trait_combi(no_of_ud)) ! trait combination
ALLOCATE(rhs(no_of_eq, no_of_trait * 2)) ! right hand side
ALLOCATE(loc_of_nz(no_of_trait,no_of_fixed + 1)) ! location of none-zero element
ALLOCATE(loc_of_nz_a(no_of_trait,no_of_fixed))
ALLOCATE(sum_sol(no_of_trait)) ! temporary solution
ALLOCATE(solutions(no_of_eq, no_of_trait)) ! solution
ALLOCATE(pre_sol(no_of_trait))
ALLOCATE(diag_ele(no_of_trait,no_of_trait))
ALLOCATE(temp_rhs(no_of_trait)) ! initialiaze solutions to zero
solutions = 0
! open file for data processing procedure
OPEN(UNIT = 31, FILE = 'ongoing_solve.dat', STATUS = 'REPLACE', ACTION = 'WRITE')
! open parameter file
OPEN(UNIT = 21, FILE = 'vcvtrt.par', 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 ! residual variance-covariance
READ(UNIT = 21, FMT = *) rvcv
write(31,*) 'residual variance-covariance = ', rvcv
! genetic variance-covariance
READ(UNIT = 21, FMT = *) gvcv
write(31,*) 'genetic variance-covariance = ', 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
write(31,*) 'inverse of genetic variance-covariance', ((gvcv_m(i,j),i = 1, no_of_trait),j = 1, no_of_trait)
! trait combination
DO k = 1, no_of_tc
READ(UNIT = 21, FMT = *) trait_combi
write(31,*) 'trait combination ', k , '=', trait_combi
temp_rvcv = rvcv * trait_combi
write(31,*) 'trati combination residual vcv', temp_rvcv
z = 0.0
CALL DJNVHF(temp_rvcv,no_of_trait,iw,z,mr) ! inverse of upper diagonal matrix(residual)
write(31,*) 'inverse of trait combination rvcv', temp_rvcv
rvcv_tc(k,:) = temp_rvcv ! sotre the inverse of upper diagonal matrix(residual)
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-han 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_trait * (no_of_fixed + 1) + 4))
write(31,*) 'number of lines of left hand side = ', no_of_lhs ! 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) write(31,*) 'left hand side, lines = ', no_of_lhs
DO i = 1, no_of_lhs
write(31,*) lhs(i,:)
END DO !open right-hand side file
OPEN(UNIT = 12, 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 = 12,FMT = *) rhs(i,:)
END DO
close(12) write(31,*) 'right hand side, number of equations = ', no_of_eq
DO i = 1, no_of_eq
write(31,*) i, rhs(i,:)
END DO ! start iteration
DO iteration = 1, MAX_ITER
pre_eq_no = 1
diag_ele = 0.0
temp_rhs = rhs(lhs(1,1),1:no_of_trait)
epsilon = 0.0 ! start reading and processing each line of lhs
DO i = 1, no_of_lhs cur_eq_no = lhs(i,1)
cur_trait_no = lhs(i,2)
data_type = lhs(i,9)
tc_no = lhs(i,10) ! if new equation, calculate solution
IF (pre_eq_no /= cur_eq_no) THEN
pre_sol = solutions(pre_eq_no, :)
CALL geninv(diag_ele, inv_diag_ele)
solutions(pre_eq_no, :) = MATMUL(inv_diag_ele, temp_rhs)
epsilon = epsilon + SUM((pre_sol - solutions(pre_eq_no, :))**2) diag_ele = 0.0
temp_rhs = rhs(lhs(i,1), 1:no_of_trait) pre_eq_no = cur_eq_no
END IF ! adjust rhs and add diag_ele
SELECT CASE(data_type) CASE(1) ! fixed effect
DO k = 1, no_of_trait
loc_of_nz(k,:) = lhs(i,(k-1)*(no_of_fixed + 1)+ 3 : k * (no_of_fixed + 1) + 2)
END DO sum_sol = 0
DO k = 1, no_of_trait
If (loc_of_nz(k,1) /= 0) THEN
sum_sol(k) = SUM(solutions(loc_of_nz(k,:),k))
END IF
END DO temp_rhs(cur_trait_no) = temp_rhs(cur_trait_no) &
- SUM(rvcv_tc_m(tc_no,cur_trait_no,:) * sum_sol) &
+ rvcv_tc_m(tc_no,cur_trait_no,cur_trait_no) * solutions(cur_eq_no, cur_trait_no)
diag_ele(cur_trait_no,cur_trait_no) = diag_ele(cur_trait_no,cur_trait_no) + rvcv_tc_m(tc_no, cur_trait_no, cur_trait_no) CASE(2) ! animal as a fixed effect
DO k = 1, no_of_trait
loc_of_nz_a(k,:) = lhs(i,(k-1)*(no_of_fixed + 1)+ 3 : k * (no_of_fixed + 1) + 1)
END DO sum_sol = 0
DO k = 1, no_of_trait
If (loc_of_nz_a(k,1) /= 0) THEN
sum_sol(k) = SUM(solutions(loc_of_nz_a(k,:),k))
END IF
END DO temp_rhs(cur_trait_no) = temp_rhs(cur_trait_no) &
- SUM(rvcv_tc_m(tc_no,cur_trait_no,:) * sum_sol) diag_ele(cur_trait_no,:) = diag_ele(cur_trait_no,:) + rvcv_tc_m(tc_no, cur_trait_no,:)
CASE(1001) ! animal, both parents are known, animal is diagonal temp_rhs = temp_rhs - (-1.) * MATMUL(gvcv_m,solutions(lhs(i,4),:)) &
- (-1.) * MATMUL(gvcv_m,solutions(lhs(i,5),:)) diag_ele = diag_ele + 2.0 * gvcv_m CASE(1002) ! animal, both parents are known, sire or dam is diagonal
temp_rhs = temp_rhs - (-1.0) * MATMUL(gvcv_m,solutions(lhs(i,3),:)) &
- ( 0.5) * MATMUL(gvcv_m,solutions(lhs(i,4),:)) diag_ele = diag_ele + ( 0.5) * gvcv_m CASE(1003) ! animal, one parent is known, animal is diagonal
temp_rhs = temp_rhs - (-2.0/3.0) * MATMUL(gvcv_m,solutions(lhs(i,4),:)) diag_ele = diag_ele + (4.0/3.0) * gvcv_m CASE(1004) ! animal effect, one parents is known, parent is diagonal
temp_rhs = temp_rhs - (-2.0/3.0) * MATMUL(gvcv_m,solutions(lhs(i,3),:)) diag_ele = diag_ele + (1.0/3.0) * gvcv_m CASE(1005) ! animal effect, no parents is known
diag_ele = diag_ele + 1.0 * gvcv_m CASE DEFAULT
WRITE(31,*) cur_eq_no, 'equation number has a ', data_type,' Invalid data_type. ', i,'th line of lhs'
END SELECT END DO
! end reading lhs ! calculate last solution
pre_sol = solutions(pre_eq_no, :)
CALL geninv(diag_ele, inv_diag_ele)
solutions(pre_eq_no, :) = MATMUL(inv_diag_ele, temp_rhs)
epsilon = (epsilon + SUM((pre_sol - solutions(pre_eq_no, :))**2))/(no_of_eq * no_of_trait) WRITE(31,*) iteration,'th iteration''s solutions'
DO k = 1, no_of_eq
WRITE(31,*) k, solutions(k,:)
END DO ! write iteration number and epsilon
write(*,*) 'iteration = ', iteration , ', epsilon = ', epsilon
IF (epsilon < criteria) THEN
EXIT
END IF END DO
! end iteration
CLOSE(31) ! open file for writing solution
OPEN(UNIT=13, FILE='sol.dat', STATUS='REPLACE', ACTION='WRITE', IOSTAT=status)
DO i = 1, no_of_eq
WRITE(13,*) i, solutions(i,:)
END DO
CLOSE(13)
END PROGRAM mt_dl_animal_model3_solve
컴파일 : gfortran inverse.f95 mt_dl_animal_model3_solve.f95 -o mt_dl_animal_model3_solve
위 프로그램을 이용하여 구한 해
1 8454.2182152398800 9927.6791905101018
2 8261.1167651059532 0.0000000000000000
3 0.0000000000000000 8800.0109257839267
4 0.0000000000000000 11287.160013707018
5 786.53649245948100 0.0000000000000000
6 372.58462226476178 177.24048702883138
7 -418.84723805070314 14.282410520178871
8 -183.66081971739675 -150.43353735092694
9 -36.017904252047643 40.594179480067496
10 292.90497151403378 146.45248593441889
11 -305.30774751739148 -172.18824815595565
12 85.629022375033102 62.348888048670780
13 -171.25804250755650 -124.69777383415753
14 24.805558020331283 51.471532363231937
15 146.45248554791553 73.226242808573829
16 -146.45248319577800 -73.226239349709402
17 146.45248605247176 73.226243198593409
파일
1247413635_dairy2.dat
1247413635_pedi2.dat
1247413635_vcvtrt.par
1247413635_mt_dl_animal_model3_setup.f95
1247413635_inverse.f95
1247413635_lhs.dat
1247413635_rhs.dat
1247413635_sorted_lhs.dat
1247413635_mt_dl_animal_model3_solve.f95
1247413635_sol.dat
한 방에 묶은 파일
1247413635_mt_dl_animal_model3.zip
'Animal Breeding > Fortran program' 카테고리의 다른 글
| 5형질 multple trait animal model - same model but different leve. (0) | 2009.10.09 |
|---|---|
| same model, different level, multiple trait animal model 3 - 3형질 (0) | 2009.07.14 |
| same model, different level, multiple trait animal model 1 (0) | 2009.06.16 |
| multiple trait animal model - solve (0) | 2009.06.06 |
| multiple trait animal model setup (0) | 2009.05.31 |