Sire Model - Gauss Seidel 방법을 이용

data

 

1 1 1 28
1 2 1 19
1 3 2 25
2 1 1 30
2 3 2 21
2 3 2 30
3 2 1 21
3 1 1 27
3 3 2 25
3 2 2 19
1st - birth year2nd - dam of age3rd - sire (random)4th - observation

 

위 자료를 sire_model.dat로 저장

 

sire list12

 

위 자료를 sire_list.dat로 저장

 

lhs.dat 및 rhs.dat를 만드는 프로그램

 

 

PROGRAM sire_model_setup
! program name - 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.4.20.
! update - none

IMPLICIT NONE

! data dictionary
CHARACTER(LEN = 40) :: data_filename ! data file name
CHARACTER(LEN = 40) :: sire_filename ! sire file name
INTEGER :: No_Of_Effects ! number of effects
INTEGER, ALLOCATABLE :: No_Of_Levels(:) ! number of levels of each fixed effect

INTEGER, ALLOCATABLE :: effect_data(:) ! fixed and sire effects for each line
INTEGER :: sire ! sire name of sire file
REAL :: 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 :: status ! I/O status
CHARACTER(LEN = 40) :: error_msg ! error message

! 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 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 ! initialize xpy - right-hand side
No_Of_Equations = SUM(No_Of_Levels)
ALLOCATE(xpy(No_Of_Equations))
xpy = 0 ! initialize

! 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 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, value, '1'
END DO

! add value to xpy
DO i = 1, No_Of_Effects
xpy(effect_data(i)) = xpy(effect_data(i)) + 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)) , '0 ', '0 ', '0 ', '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 sideEND PROGRAM sire_model_setup
위 프로그램을 sire_model_setup.f95 로 저장

 

위 프로그램에서 나온 lhs.dat

 

1 4 7 28.000000 1
4 1 7 28.000000 1
7 1 4 28.000000 1
1 5 7 19.000000 1
5 1 7 19.000000 1
7 1 5 19.000000 1
1 6 8 25.000000 1
6 1 8 25.000000 1
8 1 6 25.000000 1
2 4 7 30.000000 1
4 2 7 30.000000 1
7 2 4 30.000000 1
2 6 8 21.000000 1
6 2 8 21.000000 1
8 2 6 21.000000 1
2 6 8 30.000000 1
6 2 8 30.000000 1
8 2 6 30.000000 1
3 5 7 21.000000 1
5 3 7 21.000000 1
7 3 5 21.000000 1
3 4 7 27.000000 1
4 3 7 27.000000 1
7 3 4 27.000000 1
3 6 8 25.000000 1
6 3 8 25.000000 1
8 3 6 25.000000 1
3 5 8 19.000000 1
5 3 8 19.000000 1
8 3 5 19.000000 1
7 0 0 0 2
8 0 0 0 2

 

위의 lhs.dat를

 

sort -n -o sorted_lhs.dat lhs.dat

 

위의 명령어로 정렬(유닉스 sort 명령어를 사용할 것)

 

rhs.dat

 

72.000000
81.000000
92.000000
85.000000
59.000000
101.00000
125.00000
120.00000
위의 sorted_lhs.dat와 rhs.dat를 이용하여 해를 구하는 프로그램

 

PROGRAM sire_model_solve ! program name - sire_model_solve
! programmer - Park Byoungho
! usage - sire_model_solve < sire_model_solve.par
! purpose : read sorted_lhs.dat and rhs.dat which are made in sire_model_setup.exe and are sorted
! find a solution using Gauss-Seidel iteration
! Date - 2009.4.20
! update - none 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_in_eq(:) ! location of off-diagonal in the equation
REAL :: value ! observation
INTEGER :: class_code ! to identify fixed and random : 1-fixed, 2-random
INTEGER, ALLOCATABLE :: level_of_effects(:) ! level of each fixed and sire effect
REAL :: diag_ele ! diagonal element of lhs
INTEGER :: no_of_eq ! number of equations
INTEGER :: loc_of_sire ! starting location of sire effect
REAL :: var_ratio ! variance ratio between error variance and sire variance
REAL, ALLOCATABLE :: lhs(:,:), rhs(:)! left-hand side and right-hand side, x'y
REAL :: rhs_element ! right-hand side one element
REAL, ALLOCATABLE :: solutions(:) ! solutions
REAL :: pre_sol ! previous solution
INTEGER :: iteration ! iteration
REAL :: epsilon ! sum of square between old and new solutions
INTEGER :: i ! 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-sire ! get the number of effectss
WRITE (*,*) "Enter the number of effects : "
READ (*,*) no_of_effects ALLOCATE(loc_in_eq(no_of_effects - 1), level_of_effects(no_of_effects)) ! get the number of equations
WRITE (*,*) "Enter the level of each fixed and sire effect : "
READ (*,*) level_of_effects

no_of_eq = SUM(level_of_effects) ! total number of equations
loc_of_sire = SUM(level_of_effects(1:no_of_effects-1)) + 1 ! (sum of level of fixed effects) + 1

ALLOCATE(rhs(no_of_eq), solutions(no_of_eq))

! initialiaze solutions to zero
solutions = 0

! get the variance ratio
WRITE (*,*) "Enter the variance ratio between error variance and sire variance : "
READ (*,*) var_ratio

!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_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
READ(UNIT = 11,FMT = *) rhs
close(11)

 

! iteration
DO iteration = 1, MAX_ITER pre_eq_no = 0
diag_ele = 0.0
epsilon = 0.0 ! lhs
DO i = 1, no_of_lhs
loc_in_eq = INT(lhs(i,2:no_of_effects)) ! location of off diagonal
data_type = INT(lhs(i,no_of_effects + 2))

SELECT CASE(data_type)
CASE(1) ! fixed effect
IF (pre_eq_no == lhs(i,1)) THEN ! same equation number
rhs_element = rhs_element - sum(solutions(loc_in_eq))
diag_ele = diag_ele + 1
pre_eq_no = lhs(i,1)
ELSE ! new equation number
IF (lhs(i,1) > 1) THEN ! not first
pre_sol = solutions(pre_eq_no) ! store old solution
solutions(pre_eq_no) = rhs_element / diag_ele ! get a new solution
epsilon = epsilon + (pre_sol - solutions(pre_eq_no)) * (pre_sol - solutions(pre_eq_no)) ! calculate sum of squares of difference between old solution and new solution
END IF pre_eq_no = lhs(i,1)
rhs_element = rhs(INT(lhs(i,1))) - SUM(solutions(loc_in_eq))
diag_ele = 1 END IF
CASE(2) ! sire effect
diag_ele = diag_ele + var_ratio
END SELECT
END DO ! one iteration ends - end of array

! calculate last solution
pre_sol = solutions(pre_eq_no) ! store old solution
solutions(pre_eq_no) = rhs_element / diag_ele ! get a new solution
! calculate epsilon
epsilon = epsilon + (pre_sol - solutions(pre_eq_no)) ** 2 ! calculate sum of squares of difference between old solution and new solutions ! write iteration and epsilon
epsilon = epsilon / no_of_eq
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 sire_model_solve

 

위의 프로그램을 sire_model_solve.f95 로 저장

 

실행을 편하게 하기 위하여

 

sire_model.dat
sire_list.dat
3
3 3 2
위 자료를 sire_model_setup.par 로 저장

 

3
3 3 2
10
위 자료를 sire_model_solve.par 로 저장

 

구한 해

 

24.041245
25.403793
24.166222
3.7337582
-4.4787250
0.55872822
6.24965653E-02
-6.24977127E-02

 

프로그램 컴파일 및 실행