multiple trait sire model - lhs와 rhs를 setup하기

 

multiple trait sire model - lhs와 rhs를 setup하기

 

 

자료

 

1 1 1 28 82 570
1 2 1 19 72 530
1 3 2 25 69 490
2 1 1 30 90 600
2 3 2 21 83 559
2 3 2 30 79 530
3 2 1 21 76 560
3 1 1 27 88 602
3 3 2 25 91 0
3 2 2 19 83 0
위 자료를 mt_sire_model.dat로 저장

 

1st - fixed, birth-year2nd - fixed, age of dam3rd - random, sire4th - observation, birth weight5th - observation,weaning weight6th - observation,yearling weight

 

sire list

 

1
2
위 sire list를 sire_list.dat로 저장

 

 

유전 분산-공분산, 잔차 분산-공분산, trait combination을 담고 있는 파일

 

1.25 1.75 2.5 7.5 62.5 500
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
위 자료를 mt_sire_model_par.dat로 저장

 

 

위 자료를 이용하여 multiple trait sire model setup하는 프로그램

 

PROGRAM mt_sire_model_setup


! program name - multiple trait sire model setup
! programmer - Park Byoungho
! usage - sire_model_setup < parameter_file
! purpose : read data file and write left-hand side and right-hand side
! Date - 2009.5.20.
! update - none

USE gi

IMPLICIT NONE

! data dictionary
CHARACTER(LEN = 256) :: data_filename ! data file name
CHARACTER(LEN = 256) :: sire_filename ! sire 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 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
INTEGER, ALLOCATABLE :: trait_combi(:,:) ! trait combination
REAL(KIND = 8), ALLOCATABLE :: rvcv_tc(:,:) ! rvcv according to trait combination and then inverse
REAL(KIND = 8), ALLOCATABLE :: rmat_temp(:,:) ! temorary residual variance-covariance matirx according to 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 sire effects for each line
INTEGER :: sire ! sire name of sire file
REAL, ALLOCATABLE :: value(:) ! value of trait for each line

REAL, ALLOCATABLE :: xpy(:,:) ! right-hand side

INTEGER :: No_Of_Equations ! number of equation
INTEGER :: i, j ! 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 :: dummy = 0

! get data file name
WRITE (*,*) "Enter data file name : "
READ (*,*) data_filename

! get sire file name
WRITE (*,*) "Enter sire file name : "
READ (*,*) sire_filename

! get parameter file name
WRITE (*,*) "Enter parameter file name : "
READ (*,*) par_filename

! get number of effects including fixed and sire effect
WRITE (*,*) "Enter the number of effects including fixed and sire 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 sire 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(trait_combi(no_of_tc, no_of_ud)) ! trait combination
ALLOCATE(rvcv_tc(no_of_tc,no_of_ud)) !

ALLOCATE(rmat_temp(no_of_trait,no_of_trait))

! initialize xpy - right-hand side
No_Of_Equations = SUM(No_Of_Levels)
ALLOCATE(xpy(No_Of_Equations,no_of_trait))
xpy = 0 ! initialize

! 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

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

! trait combination
DO i = 1, no_of_tc
READ(UNIT = 21, FMT = *) trait_combi(i,:)
END DO

! making rvcv_tc
DO i = 1, no_of_tc
rvcv_tc(i,:) = rvcv * trait_combi(i,:)
END DO

! inverse of each residual variance-covariance according to trait combination
DO i = 1, no_of_tc z = 0.0
CALL DJNVHF(rvcv_tc(i,:),no_of_trait,iw,z,mr)
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

DO i = 1, no_of_trait
DO j = 1, no_of_trait
rmat_temp(i,j) = rvcv_tc(tc_no,IHMSSF(i,j,no_of_trait))
END DO
END DO

! add value to xpy
DO i = 1, No_Of_Effects
xpy(effect_data(i),:) = xpy(effect_data(i),:) + MATMUL(rmat_temp,value)
END DO

END DO ! end of reading data

CLOSE(11) ! close data file


! open sire file
OPEN(UNIT = 13, FILE = sire_filename, STATUS = 'OLD', ACTION = 'READ', IOSTAT = status, IOMSG = error_msg)

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

! read each line
DO
READ(UNIT = 13, FMT = *, IOSTAT = status) sire
IF (status /= 0) EXIT ! end of file

! write sire effect
WRITE(UNIT = 12, FMT = *) sire + SUM(no_of_levels(1:no_of_effects - 1)) , dummy, dummy, dummy, '2'
END DO ! end of reading data

CLOSE(13) ! close sire 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 side

END PROGRAM mt_sire_model_setup

위 프로그램 소스를 mt_sire_model_setup.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 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_sire_model.dat
sire_list.dat
mt_sire_model_par.dat
3
3 3 2
3

 

위를 mt_sire_model_setup.par로 저장

 

 

위 프로그램을 이용하여 나온 left hand side

 

1 4 7 7 1
4 1 7 7 1
7 1 4 7 1
1 5 7 7 1
5 1 7 7 1
7 1 5 7 1
1 6 8 7 1
6 1 8 7 1
8 1 6 7 1
2 4 7 7 1
4 2 7 7 1
7 2 4 7 1
2 6 8 7 1
6 2 8 7 1
8 2 6 7 1
2 6 8 7 1
6 2 8 7 1
8 2 6 7 1
3 5 7 7 1
5 3 7 7 1
7 3 5 7 1
3 4 7 7 1
4 3 7 7 1
7 3 4 7 1
3 6 8 3 1
6 3 8 3 1
8 3 6 3 1
3 5 8 3 1
5 3 8 3 1
8 3 5 3 1
7 0 0 0 2
8 0 0 0 2
위 프로그램을 이용하여 구한 right hand side

 

1.9571227 1.7957478 0.23437987
2.3159108 2.0812473 0.23651524
2.5084927 3.2069154 0.17655139
2.4685802 2.1218426 0.25010207
1.1732187 2.1414437 0.17449327
3.1397271 2.8206241 0.22285117
3.0944958 3.3640823 0.42459533
3.6870303 3.7198284 0.22285117

 

컴파일 명령문gfortran inverse.f95 mt_sire_model_setup.f95 -o mt_sire_model_setup.exe

 

실행 명령문mt_sire_model_setup < mt_sire_model_setup.par

 

 

컴파일 및 실행 화면

 

1242833081_mt_sire_model.dat

다운로드

1242833081_mt_sire_model.dat

1242833081_sire_list.dat

다운로드

1242833081_sire_list.dat

1242833081_mt_sire_model_par.dat

다운로드

1242833081_mt_sire_model_par.dat

1242833081_mt_sire_model_setup.f95

다운로드

1242833081_mt_sire_model_setup.f95

1242833081_inverse.f95

다운로드

1242833081_inverse.f95

1242833081_mt_sire_model_setup.par

다운로드

1242833081_mt_sire_model_setup.par

1242833081_lhs.dat

다운로드

1242833081_lhs.dat

1242833081_rhs.dat

다운로드

1242833081_rhs.dat