multiple trait animal model setup
DATA
1 1 8 28 82 570
1 2 9 19 72 530
1 3 10 25 69 490
2 1 11 30 90 600
2 3 12 21 83 559
2 3 13 30 79 530
3 2 14 21 76 560
3 1 15 27 88 602
3 3 16 25 91 0
3 2 17 19 83 0
1st col : year of birth2nd col : age of dam3dr col : animal4th col : birth weight5th col : weaning weight6th col : yearling weight
위 자료를 mt_animal_model.dat로 저장
혈통자료
1 0 0
2 0 0
3 0 0
4 0 0
5 0 0
6 0 0
7 0 0
8 6 1
9 6 2
10 7 3
11 6 1
12 7 2
13 7 4
14 6 1
15 6 2
16 7 3
17 7 5
1st : animal2nd : sire3rd : dam
위 자료를 mt_animal_pedi.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
1st row : generic variance-covariance2nd row : residual variance-covariancerest : trait combination of 3 triat
위 자료를 mt_animal_model_par.dat 로 저장
프로그램
PROGRAM mt_animal_model_setup
! program name - multiple trait animal model setup
! programmer - Park Byoungho
! usage - animal_model_setup < parameter_file
! purpose : read data file and write left-hand side and right-hand side
! Date - 2009.5.31.
! update - none
USE gi
IMPLICIT NONE
! data dictionary
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_Effects ! number of effects
INTEGER, ALLOCATABLE :: No_Of_Levels(:) ! number of levels of each fixed and animal effect
INTEGER :: level_of_fixed_effect ! total level of fixed effect INTEGER :: no_of_trait ! number of trait
REAL, 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_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 :: effect_data(:) ! fixed and animal effects for each line
INTEGER :: animal, sire, dam ! animal, sire, dam name for pedigree file
REAL(KIND = 8), ALLOCATABLE :: value(:) ! value of trait for each line
REAL(KIND = 8), ALLOCATABLE :: xpy(:,:) ! right-hand side
INTEGER :: No_Of_Equations ! number of equation
INTEGER :: i, j, k ! 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
! get data file name
WRITE (*,*) "Enter data file name : "
READ (*,*) data_filename ! get pedigree file name
WRITE (*,*) "Enter pedigree file name : "
READ (*,*) pedi_filename
! get parameter file name
WRITE (*,*) "Enter parameter file name : "
READ (*,*) par_filename
! get number of effects including fixed and animal effect
WRITE (*,*) "Enter the number of effects including fixed and animal effects: "
READ (*,*) No_Of_Effects ! allocate dimensions of no_of_levels and fixed_data
ALLOCATE(No_Of_Levels(No_Of_Effects), effect_data(No_Of_Effects))
WRITE (*,*) "Enter the number of levels for each fixed and animal effect : "
READ (*,*) No_Of_Levels WRITE (*,*) "Enter the number of traits : "
READ (*,*) 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(value(no_of_trait)) ! observation
ALLOCATE(gvcv(no_of_ud)) ! genetic variance-covariance
ALLOCATE(rvcv(no_of_ud)) ! residual variance-covariance
ALLOCATE(rvcv_tc_m(no_of_tc,no_of_trait,no_of_trait))
ALLOCATE(temp_rvcv(no_of_ud)) ! temporary residual variance-covariance
ALLOCATE(trait_combi(no_of_ud)) ! trait combination
! initialize xpy - right-hand side
No_Of_Equations = SUM(No_Of_Levels)
ALLOCATE(xpy(No_Of_Equations, no_of_trait))
xpy = 0 ! initialize
level_of_fixed_effect = SUM(no_of_levels(1:no_of_effects-1)) ! 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
! genetic variance-covariance
READ(UNIT = 21, FMT = *) gvcv
WRITE(31,*) 'genetic variance-covariance = ', gvcv
! residual variance-covariance
READ(UNIT = 21, FMT = *) rvcv
WRITE(31,*) 'residual variance-covariance = ', rvcv
! 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)
WRITE(31,*) 'inverse of trait combination rvcv', temp_rvcv
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
WRITE(31,*) 'matrix of rvcv', rvcv_tc_m(k,:,:)
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) effect_data, value
IF (status /= 0) EXIT ! end of file
! find trait combination number
tc_no = 0
DO i = 1, no_of_trait
IF (value(i) /= 0.0) tc_no = tc_no + 2 ** (i - 1)
END DO ! find location of equations
DO i = 2, No_Of_Effects
effect_data(i) = effect_data(i) + SUM(No_Of_Levels(1:i - 1))
END DO ! write location of non-zero elements :
DO i = 1, No_Of_Effects
! swap
temp_value = effect_data(i)
effect_data(i) = effect_data(1)
effect_data(1) = temp_value
WRITE(UNIT = 12, FMT = *) effect_data, tc_no, '1'
END DO
! add value to xpy
DO i = 1, No_Of_Effects
xpy(effect_data(i),:) = xpy(effect_data(i),:) + MATMUL(rvcv_tc_m(tc_no,:,:),value)
END DO
END DO ! end of reading data 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 tc_no = 0! read each line
DO
READ(UNIT = 13, FMT = *, IOSTAT = status) animal, sire, dam
IF (status /= 0) EXIT ! end of file ! find location of equations
animal = animal + level_of_fixed_effect
IF (sire /= 0) sire = sire + level_of_fixed_effect
IF (dam /= 0) dam = dam + level_of_fixed_effect ! write animal effect
IF ( sire /= 0 .AND. dam /= 0) THEN
WRITE(UNIT = 12, FMT = *) animal, sire, dam, tc_no, '2'
WRITE(UNIT = 12, FMT = *) sire, animal, dam, tc_no, '3'
WRITE(UNIT = 12, FMT = *) dam, animal, sire, tc_no, '3'
ELSE IF (sire /= 0 .AND. dam == 0) THEN
WRITE(UNIT = 12, FMT = *) animal, sire, dam, tc_no, '4'
WRITE(UNIT = 12, FMT = *) sire, animal, dam, tc_no, '5'
ELSE IF (sire == 0 .AND. dam /= 0) THEN
WRITE(UNIT = 12, FMT = *) animal, sire, dam, tc_no, '4'
WRITE(UNIT = 12, FMT = *) dam, animal, sire, tc_no, '5'
ELSE
WRITE(UNIT = 12, FMT = *) animal, sire, dam, tc_no, '6'
END IF
END DO ! end of reading data 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_Equations
WRITE(UNIT = 14, FMT = *) xpy(i,:)
END DO
CLOSE(14) ! close file for right-hand sideEND PROGRAM mt_animal_model_setup
위 프로그램을 mt_animal_model_setup.f95 로 저장
위에서 inverse 하는 프로그램 module
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로 저장
위 프로그램을 실행하기 편하게 만든 파라미터 파일
mt_animal_model.dat
mt_animal_pedi.dat
mt_animal_model_par.dat
3
3 3 17
3
위 자료를 mt_animal_model_setup.par로 저장
프로그램을 실행하여나온 lhs.dat
1 4 14 7 1
4 1 14 7 1
14 1 4 7 1
1 5 15 7 1
5 1 15 7 1
15 1 5 7 1
1 6 16 7 1
6 1 16 7 1
16 1 6 7 1
2 4 17 7 1
4 2 17 7 1
17 2 4 7 1
2 6 18 7 1
6 2 18 7 1
18 2 6 7 1
2 6 19 7 1
6 2 19 7 1
19 2 6 7 1
3 5 20 7 1
5 3 20 7 1
20 3 5 7 1
3 4 21 7 1
4 3 21 7 1
21 3 4 7 1
3 6 22 3 1
6 3 22 3 1
22 3 6 3 1
3 5 23 3 1
5 3 23 3 1
23 3 5 3 1
7 0 0 0 6
8 0 0 0 6
9 0 0 0 6
10 0 0 0 6
11 0 0 0 6
12 0 0 0 6
13 0 0 0 6
14 7 12 0 2
7 14 12 0 3
12 14 7 0 3
15 8 12 0 2
8 15 12 0 3
12 15 8 0 3
16 9 13 0 2
9 16 13 0 3
13 16 9 0 3
17 7 12 0 2
7 17 12 0 3
12 17 7 0 3
18 8 13 0 2
8 18 13 0 3
13 18 8 0 3
19 10 13 0 2
10 19 13 0 3
13 19 10 0 3
20 7 12 0 2
7 20 12 0 3
12 20 7 0 3
21 8 12 0 2
8 21 12 0 3
12 21 8 0 3
22 9 13 0 2
9 22 13 0 3
13 22 9 0 3
23 11 13 0 2
11 23 13 0 3
13 23 11 0 3
프로그램을 실행하여 나온 rhs.dat
1.9571226080793762 1.7957476966690291 0.23437987243090011
2.3159107016300489 2.0812473423104185 0.23651523742026931
2.5084927007139184 3.2069154257848993 0.17655138199858258
2.4685802031656037 2.1218426647767545 0.25010205527994334
1.1732186978527308 2.1414437904528079 0.17449326718639263
3.1397271094050092 2.8206240095347845 0.22285116938341604
0.0000000000000000 0.0000000000000000 0.0000000000000000
0.0000000000000000 0.0000000000000000 0.0000000000000000
0.0000000000000000 0.0000000000000000 0.0000000000000000
0.0000000000000000 0.0000000000000000 0.0000000000000000
0.0000000000000000 0.0000000000000000 0.0000000000000000
0.0000000000000000 0.0000000000000000 0.0000000000000000
0.0000000000000000 0.0000000000000000 0.0000000000000000
0.86120954405858741 0.65400425230333104 8.03756201275691168E-002
0.26765887077722650 0.60836286321757616 8.53862508858965336E-002
0.82825419324356231 0.53338058114812192 6.86180014174344599E-002
0.91778880226789483 0.73564847625797325 8.22820694542877457E-002
0.29996456413890848 0.73255846917080092 8.56647767540751337E-002
1.0981573352232457 0.61304039688164425 6.85683912119064570E-002
0.35825655563430187 0.63387668320340185 8.91070163004961113E-002
0.68958185683912121 0.73218993621545003 8.74443656980864736E-002
0.91335101679929265 0.94164456233421745 0.0000000000000000
0.54730327144120250 0.89920424403183019 0.0000000000000000
프로그램 컴파일 및 실행 화면 캡쳐
1243780584_mt_animal_model.dat
1244216071_mt_animal_pedi.dat
1243780584_mt_animal_model_par.dat
1243780584_mt_animal_model_setup.f95
1243780584_inverse.f95
1243780584_mt_animal_model_setup.par
1243780584_lhs.dat
1243780584_rhs.dat
'Animal Breeding > Fortran program' 카테고리의 다른 글
| same model, different level, multiple trait animal model 1 (0) | 2009.06.16 |
|---|---|
| multiple trait animal model - solve (0) | 2009.06.06 |
| multiple trait sire model - Gauss seidel solver (0) | 2009.05.25 |
| 혈통을 세대별로 나누기 (0) | 2009.05.25 |
| multiple trait sire model - lhs와 rhs를 setup하기 (0) | 2009.05.21 |