Single Trait Random Regression Model - 305 육종가 계산

 

Single Trait Random Regression Model - 305일 육종가 계산

 

이전 프로그램은 개체 육종가의 회귀계수만 표시했는데 그걸 이용해서 6일부터 310일까지의 매일 매일의 육종가를 더해서 305일 육종가를 계산하는 프로그램 로직을 더했음

 

매일 매일의 육종가를 계산하여 더하는 것이 아니라, 6일부터 310일까지의 르장드르 폴리노미얼을 더해서 한 방에 개체 육종가의 회귀계수와 곱하여 305일 육종가 계산

 

개체 육종가 파일 - sol_bv.dat

 

해를 구하는 프로그램

 

PROGRAM strrm_solve

! program name
! : iteration on data of single trait random regression model
! : solve program
! prgrammer
! : Park Byoungho
! Update
! : 2010.2.11. - 1st edition
! : 2010.2.22. - calculate 305-day breeding value

USE gi

IMPLICIT NONE

! data dictionary

INTEGER, PARAMETER :: no_of_cv = 1 ! number of class variables
INTEGER, PARAMETER :: no_of_rc = 5 ! number of regression coefficient
INTEGER, PARAMETER :: no_of_fe = no_of_cv + no_of_rc ! number of fixed effect: class variable + no of fixed regression coefficient
INTEGER, PARAMETER :: no_of_effect = no_of_cv + no_of_rc + 2 ! number of effects

REAL(KIND = 8), PARAMETER :: animal_vr(3,3) = & ! variance ratio of animal genetic variance and residual variance
RESHAPE((/ 3.297, 0.594, -1.381, &
0.594, 0.921, -0.289, &
-1.381, -0.289, 1.005/),(/3,3/))

REAL(KIND = 8), PARAMETER :: pe_vr(3,3) = & ! variance ratio of permanent environmental efeect variance and residual variance
RESHAPE((/ 6.872, -0.254, -1.101, &
-0.254, 3.171, 0.167, &
-1.101, 0.167, 2.457/),(/3,3/))

REAL(KIND = 8), PARAMETER :: residual_vr = 3.710

REAL(KIND = 8) :: animal_vr_i(3,3), pe_vr_i(3,3) ! (inverse of animal and permanent environment variance-covariance) * residual variance

INTEGER, PARAMETER :: levels_of_cv = 10 ! levels of each calss variables
INTEGER, PARAMETER :: no_of_animal = 8
INTEGER, PARAMETER :: no_of_eq = levels_of_cv + no_of_rc + no_of_animal * 2

INTEGER :: no_of_lhs ! number of left hand side lines
REAL(KIND = 8), ALLOCATABLE :: lhs(:,:), rhs(:,:) ! left-hand side, right-hand side, x'y
REAL(KIND = 8), ALLOCATABLE :: solutions(:,:) ! solutions
INTEGER :: data_type ! fixed, animal, permanent environment
INTEGER :: dim ! days in milk

INTEGER :: pre_eq_no ! previous equation number
INTEGER :: cur_eq_no ! current equation number

REAL(KIND = 8) :: diag_ele(3,3) ! diagonal element
REAL(KIND = 8) :: diag_ele_i(3,3) ! diagonal element
REAL(KIND = 8) :: temp_rhs(3) ! temporary rihgt hand side
REAL(KIND = 8) :: pre_sol(3) ! previous solution

! Legendre polynomial
REAL(KIND = 8), PARAMETER :: lambda(3,3) = &
RESHAPE((/SQRT(.5) , 0.0 , 0.0 , &
0.0 , SQRT(1.5) , 0.0 , &
-SQRT(2.5)*.5 , 0.0 , SQRT(2.5)*1.5 /),(/3,3/))

INTEGER :: dim_min, dim_max, dim_range ! minimun, maximum and range of days in milk
REAL(KIND = 8) :: sut ! standardized unit of time
REAL(KIND = 8) :: sut_poly(3) ! polynomial of standardized unit of time

REAL(KIND = 8) :: dim_lp(500,3) ! Legendre polynomial of days in milk
REAL(KIND = 8) :: dim_lp_m(500,3,3) ! outer product of Legendre polynomial of days in milk

! update 2010.02.22.
REAL(KIND = 8) :: t_for305(3) ! t vector for 305-day breeding value
! sum 6 to 310 elements of dim_lp by column

INTEGER :: iteration ! iteration
REAL(KIND = 8) :: epsilon ! sum of squares of difference between old and new solutions
INTEGER :: i,j ! loop
INTEGER :: status ! i/o status
CHARACTER(LEN = 40) :: error_msg ! error message
INTEGER, PARAMETER :: MAX_ITER = 1000 ! maximum number of iteration
REAL(KIND = 8), PARAMETER :: criteria = 1.E-12 ! criteria for stopping

! calculate the inverse of animal and permanent environment variance-covariance
CALL geninv(animal_vr, animal_vr_i)
animal_vr_i = animal_vr_i * residual_vr

CALL geninv(pe_vr, pe_vr_i)
pe_vr_i = pe_vr_i * residual_vr

! allocate rhs, solution
ALLOCATE(rhs(no_of_eq, 3))
ALLOCATE(solutions(no_of_eq, 3))

! initialize
solutions = 0.0

!open left-hand side file
OPEN(UNIT = 10, 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 = 10, FMT = *, IOSTAT = status)
IF (status /= 0) EXIT ! reach end of file
no_of_lhs = no_of_lhs + 1
END DO

! decide array size of lhs
ALLOCATE(lhs(no_of_lhs, 2 + (no_of_cv + no_of_rc + 2) * 2))

!write(*,*) 'number of lines of left hand side = ', no_of_lhs

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

!DO i = 1, no_of_lhs
! write(*,*) i, lhs(i,:)
!END DO

!open right-hand side file
OPEN(UNIT = 20, 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 = 20,FMT = *) rhs(i,:)
END DO
close(20)

! write(*,*) 'right hand side ...'
! DO i = 1, levels_of_cv + no_of_rc
! write(*,*) i, rhs(i)
! END DO

! prepare Legendre polynomial of standardized unit of time(days in milk)
! open setup result file
OPEN(UNIT = 25, FILE = 'strrm_setup.rst', STATUS = 'OLD', ACTION = 'READ', IOSTAT = status, IOMSG = error_msg)

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

READ(UNIT = 25, FMT ='(/ I12 / I12 / I12)', IOSTAT = status) dim_min, dim_max, dim_range

!initialize
dim_lp = 0.0
dim_lp_m = 0.0

! prepare Legendre polynomial of standardized unit of time(days in milk)
DO i = dim_min, dim_max
sut = -1. + 2. * (i - dim_min) / dim_range ! standardized unit of time
sut_poly(1) = 1.0 ! polynomial of standardized unit of time
sut_poly(2) = sut
sut_poly(3) = sut_poly(2) * sut
dim_lp(i,:) = MATMUL(sut_poly, lambda)

DO j = 1, 3
dim_lp_m(i,j,:) = dim_lp(i,j) * dim_lp(i,:)
END DO
END DO

! update 2010.02.22.
t_for305 = SUM(dim_lp(6:310,:),1) ! t vector for 305-day breeding value calculation

! WRITE(*,*) (t_for305(i), i = 1, 3)

! solve MME
DO iteration = 1, MAX_ITER

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

! start reading and processing each line of lhs
DO i = 1, no_of_lhs

data_type = INT(lhs(i,1))
dim = INT(lhs(i,2))
cur_eq_no = INT(lhs(i,3))

IF(pre_eq_no /= cur_eq_no) THEN
pre_sol = solutions(pre_eq_no,:) ! store old solution
CALL geninv(diag_ele, diag_ele_i)
solutions(pre_eq_no,:) = MATMUL(diag_ele_i, 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
diag_ele = 0.0
temp_rhs = rhs(cur_eq_no,:)
END IF

! process off-diagnoal
SELECT CASE(data_type)
CASE(1) ! fixed effect
! adjust right-hand side
! temp_rhs - independend variable * solutions
DO j = 2, no_of_fe
temp_rhs(1) = temp_rhs(1) - lhs(i, 2*j+2) * solutions(INT(lhs(i,2*j+1)),1)
END DO

temp_rhs(1) = temp_rhs(1) - lhs(i, 16) * SUM(dim_lp(dim,:) * solutions(INT(lhs(i,15)),:)) &
- lhs(i, 18) * SUM(dim_lp(dim,:) * solutions(INT(lhs(i,17)),:))

! add diagonal element
diag_ele(1,1) = diag_ele(1,1) + lhs(i,4)
pre_eq_no = cur_eq_no

CASE(2) ! fixed part of animal effect and permanent environmental effect
! adjust right-hand side
DO j = 2, (no_of_effect - 1)
temp_rhs = temp_rhs - dim_lp(dim,:) * lhs(i,2*j+2) * solutions(INT(lhs(i,2*j+1)),1)
END DO

temp_rhs = temp_rhs - MATMUL(dim_lp_m(dim,:,:) , solutions(INT(lhs(i,17)),:))

! add diagonal element
diag_ele = diag_ele + dim_lp_m(dim,:,:)
pre_eq_no = cur_eq_no
CASE(1001) ! animal effect, both parents are known, animal is diagonal
temp_rhs = temp_rhs - (-0.5) * 2.0 * MATMUL(animal_vr_i , solutions(INT(lhs(i,4)),:))
temp_rhs = temp_rhs - (-0.5) * 2.0 * MATMUL(animal_vr_i , solutions(INT(lhs(i,5)),:))
diag_ele = diag_ele + 2.0 * animal_vr_i
pre_eq_no = cur_eq_no
CASE(1002) ! animal effect, both parents are known, sire or dam is diagonal
temp_rhs = temp_rhs - (-0.5) * 2.0 * MATMUL(animal_vr_i , solutions(INT(lhs(i,4)),:))
temp_rhs = temp_rhs - 0.25 * 2.0 * MATMUL(animal_vr_i , solutions(INT(lhs(i,5)),:))
diag_ele = diag_ele + 0.25 * 2.0 * animal_vr_i
pre_eq_no = cur_eq_no
CASE(1003) ! animal effect, one parent is known, animal is diagonal
temp_rhs = temp_rhs - (-0.5) * (4.0/3.0) * MATMUL(animal_vr_i , solutions(INT(lhs(i,3)),:))
diag_ele = diag_ele + (4.0/3.0) * animal_vr_i
pre_eq_no = cur_eq_no
CASE(1004) ! animal effect, one parents is known, parent is diagonal
temp_rhs = temp_rhs - (-0.5) * (4.0/3.0) * MATMUL(animal_vr_i , solutions(INT(lhs(i,3)),:))
diag_ele = diag_ele + 0.25 * (4.0/3.0) * animal_vr_i
pre_eq_no = cur_eq_no
CASE(1005) ! animal effect, no parents is known
diag_ele = diag_ele + animal_vr_i
pre_eq_no = cur_eq_no
CASE(2001) ! permanent environmental effect
diag_ele = diag_ele + pe_vr_i
pre_eq_no = cur_eq_no
END SELECT

END DO
! end reading lhs

! calculate last solution
pre_sol = solutions(pre_eq_no,:) ! store old solution
CALL geninv(diag_ele, diag_ele_i)
solutions(pre_eq_no,:) = MATMUL(diag_ele_i, 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,'th iteration''s solutions'
! DO i = 1, no_of_rc
! WRITE(*,*) i, solutions(i)
! END DO

epsilon = epsilon / (no_of_eq)

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

END DO
! end iteration

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

! update 2010.02.22.
! open file for animal's breeding value
OPEN(UNIT=40, FILE='sol_bv.dat', STATUS='REPLACE', ACTION='WRITE', IOSTAT=status)
DO i = 1, no_of_animal
WRITE(40,*) i, SUM(solutions(i + levels_of_cv + no_of_rc , :) * t_for305)
END DO
CLOSE(40)

END PROGRAM strrm_solve

 

sol_bv.dat 내용

 

1 -12.374094083249899
2 -15.734206104887509
3 28.105126488412509
4 74.813685726560394
5 -98.416456204591583
6 -118.42790491986700
7 184.16653015093632
8 47.687014160733646

1st col : animal number

2nd col : 305-day breeding value

 

프로그램 컴파일 및 실행 화면

 

 

관련 파일

1266887541_strrm2.zip

다운로드

1266887541_strrm2.zip