자료
animal_model.dat
1 1 8 28
1 2 9 19
1 3 10 25
2 1 11 30
2 3 12 21
2 3 13 30
3 2 14 21
3 1 15 27
3 3 16 25
3 2 17 19
1st - birth year2nd - dam of age3rd -animal (random)4th - observation
혈통자료
animal_pedi.dat
1 0 0
2 0 0
3 0 0
4 0 0
5 0 0
6 0 0
7 0 0
8 1 6
9 2 6
10 3 7
11 1 6
12 2 7
13 4 7
14 1 6
15 2 6
16 3 7
17 5 7
lhs와 rhs를 만드는 프로그램
animal_model_setup_1.f95
PROGRAM animal_model_setup
! program name - 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.1.
! update - none
IMPLICIT NONE
! data dictionary
CHARACTER(LEN = 256) :: data_filename ! data file name
CHARACTER(LEN = 256) :: pedi_filename ! pedigree 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, ALLOCATABLE :: effect_data(:) ! fixed and animal effects for each line of data file
REAL :: value ! value of trait for each line of data file
INTEGER :: animal, sire, dam ! animal, sire, dam name for pedigree file
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 pedigree file name
WRITE (*,*) "Enter pedigree file name : "
READ (*,*) pedi_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 ! initialize xpy - right-hand side
No_Of_Equations = SUM(no_of_levels)
ALLOCATE(xpy(No_Of_Equations))
xpy = 0 ! initialize
level_of_fixed_effect = SUM(no_of_levels(1:no_of_effects-1))
! 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' ! 1 means fixed effect
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
value = 0.0 ! 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 ! find location of equations
animal = animal + level_of_fixed_effect
IF (sire /= 0) sire = sire + level_of_fixed_effect ! if sire is 0, sire is unknown
IF (dam /= 0) dam = dam + level_of_fixed_effect ! if dam is 0, dam is unknown ! write animal effect
IF ( sire /= 0 .AND. dam /= 0) THEN ! both parents are known
WRITE(UNIT = 12, FMT = *) animal, sire, dam, value, '2'
WRITE(UNIT = 12, FMT = *) sire, animal, dam, value, '3'
WRITE(UNIT = 12, FMT = *) dam, animal, sire, value, '3'
ELSE IF (sire /= 0 .AND. dam == 0) THEN ! only sire is known
WRITE(UNIT = 12, FMT = *) animal, sire, dam, value, '4'
WRITE(UNIT = 12, FMT = *) sire, animal, dam, value, '5'
ELSE IF (sire == 0 .AND. dam /= 0) THEN ! only dam is known
WRITE(UNIT = 12, FMT = *) animal, sire, dam, value, '4'
WRITE(UNIT = 12, FMT = *) dam, animal, sire, value, '5'
ELSE ! no parents are known
WRITE(UNIT = 12, FMT = *) animal, sire, dam, value, '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 animal_model_setup
프로그램 실행을 편하게 하기 위하여 만든 파라미터 파일
animal_model_setup_1.par
animal_model.dat
animal_pedi.dat
3
3 3 17
위 프로그램을 이용하여 만든 lhs.dat
lhs.dat
1 4 14 28.000000 1
4 1 14 28.000000 1
14 1 4 28.000000 1
1 5 15 19.000000 1
5 1 15 19.000000 1
15 1 5 19.000000 1
1 6 16 25.000000 1
6 1 16 25.000000 1
16 1 6 25.000000 1
2 4 17 30.000000 1
4 2 17 30.000000 1
17 2 4 30.000000 1
2 6 18 21.000000 1
6 2 18 21.000000 1
18 2 6 21.000000 1
2 6 19 30.000000 1
6 2 19 30.000000 1
19 2 6 30.000000 1
3 5 20 21.000000 1
5 3 20 21.000000 1
20 3 5 21.000000 1
3 4 21 27.000000 1
4 3 21 27.000000 1
21 3 4 27.000000 1
3 6 22 25.000000 1
6 3 22 25.000000 1
22 3 6 25.000000 1
3 5 23 19.000000 1
5 3 23 19.000000 1
23 3 5 19.000000 1
7 0 0 0.0000000 6
8 0 0 0.0000000 6
9 0 0 0.0000000 6
10 0 0 0.0000000 6
11 0 0 0.0000000 6
12 0 0 0.0000000 6
13 0 0 0.0000000 6
14 7 12 0.0000000 2
7 14 12 0.0000000 3
12 14 7 0.0000000 3
15 8 12 0.0000000 2
8 15 12 0.0000000 3
12 15 8 0.0000000 3
16 9 13 0.0000000 2
9 16 13 0.0000000 3
13 16 9 0.0000000 3
17 7 12 0.0000000 2
7 17 12 0.0000000 3
12 17 7 0.0000000 3
18 8 13 0.0000000 2
8 18 13 0.0000000 3
13 18 8 0.0000000 3
19 10 13 0.0000000 2
10 19 13 0.0000000 3
13 19 10 0.0000000 3
20 7 12 0.0000000 2
7 20 12 0.0000000 3
12 20 7 0.0000000 3
21 8 12 0.0000000 2
8 21 12 0.0000000 3
12 21 8 0.0000000 3
22 9 13 0.0000000 2
9 22 13 0.0000000 3
13 22 9 0.0000000 3
23 11 13 0.0000000 2
11 23 13 0.0000000 3
13 23 11 0.0000000 3
위 프로그램을 이용하여 만든 rhs.dat
72.000000
81.000000
92.000000
85.000000
59.000000
101.00000
0.0000000
0.0000000
0.0000000
0.0000000
0.0000000
0.0000000
0.0000000
28.000000
19.000000
25.000000
30.000000
21.000000
30.000000
21.000000
27.000000
25.000000
19.000000
lhs의 정렬
명령어 : sort -n -o sorted_lhs.dat lhs.dat
정렬된 파일 : sorted_lhs.dat
1 4 14 28.000000 1
1 5 15 19.000000 1
1 6 16 25.000000 1
2 4 17 30.000000 1
2 6 18 21.000000 1
2 6 19 30.000000 1
3 4 21 27.000000 1
3 5 20 21.000000 1
3 5 23 19.000000 1
3 6 22 25.000000 1
4 1 14 28.000000 1
4 2 17 30.000000 1
4 3 21 27.000000 1
5 1 15 19.000000 1
5 3 20 21.000000 1
5 3 23 19.000000 1
6 1 16 25.000000 1
6 2 18 21.000000 1
6 2 19 30.000000 1
6 3 22 25.000000 1
7 0 0 0.0000000 6
7 14 12 0.0000000 3
7 17 12 0.0000000 3
7 20 12 0.0000000 3
8 0 0 0.0000000 6
8 15 12 0.0000000 3
8 18 13 0.0000000 3
8 21 12 0.0000000 3
9 0 0 0.0000000 6
9 16 13 0.0000000 3
9 22 13 0.0000000 3
10 0 0 0.0000000 6
10 19 13 0.0000000 3
11 0 0 0.0000000 6
11 23 13 0.0000000 3
12 0 0 0.0000000 6
12 14 7 0.0000000 3
12 15 8 0.0000000 3
12 17 7 0.0000000 3
12 20 7 0.0000000 3
12 21 8 0.0000000 3
13 0 0 0.0000000 6
13 16 9 0.0000000 3
13 18 8 0.0000000 3
13 19 10 0.0000000 3
13 22 9 0.0000000 3
13 23 11 0.0000000 3
14 1 4 28.000000 1
14 7 12 0.0000000 2
15 1 5 19.000000 1
15 8 12 0.0000000 2
16 1 6 25.000000 1
16 9 13 0.0000000 2
17 2 4 30.000000 1
17 7 12 0.0000000 2
18 2 6 21.000000 1
18 8 13 0.0000000 2
19 2 6 30.000000 1
19 10 13 0.0000000 2
20 3 5 21.000000 1
20 7 12 0.0000000 2
21 3 4 27.000000 1
21 8 12 0.0000000 2
22 3 6 25.000000 1
22 9 13 0.0000000 2
23 3 5 19.000000 1
23 11 13 0.0000000 2
위 sorted_lhs.dat와 rhs.dat를 이용하여 해를 구하는 프로그램
animal_model_solve_1.f95
PROGRAM animal_model_solve ! program name - animal_model_solve
! programmer - Park Byoungho
! usage - animal_model_solve < animal_model_solve.par
! 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.5.1.
! update - none IMPLICIT NONE
! data dictionary
INTEGER :: no_of_effects ! number of effects
INTEGER, ALLOCATABLE :: level_of_effects(:) ! level of each fixed and animal effect
INTEGER :: no_of_eq ! number of equations
INTEGER :: loc_of_animal ! starting location of animal effect 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 :: data_type ! data type : 1:fixed, 2-6:animal REAL :: diag_ele ! diagonal element of lhs
REAL :: var_ratio ! variance ratio between error variance and animal variance INTEGER :: no_of_lhs ! number of lhs lines
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 squares 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 ! 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 animal effect : "
READ (*,*) level_of_effects
no_of_eq = SUM(level_of_effects) ! total number of equations
loc_of_animal = 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 animal 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 = 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
READ(UNIT = 12,FMT = *) rhs
close(12)
! iteration
DO iteration = 1, MAX_ITER pre_eq_no = 1
diag_ele = 0.0
rhs_element = rhs(1)
epsilon = 0.0 ! read and process each line of 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))
cur_eq_no = lhs(i,1) ! nrew equation, or calculate solution
IF (pre_eq_no /= cur_eq_no) THEN
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 diag_ele = 0.0
rhs_element = rhs(INT(lhs(i,1)))
END IF
! process off-diagnoal
SELECT CASE(data_type)
CASE(1) ! fixed effect
rhs_element = rhs_element - SUM(solutions(loc_in_eq))
diag_ele = diag_ele + 1
pre_eq_no = lhs(i,1)
CASE(2) ! animal effect, both parents are known, animal is diagonal
rhs_element = rhs_element - SUM(solutions(loc_in_eq)*(/-0.5,-0.5/) * 2.0 * var_ratio)
diag_ele = diag_ele + 2.0 * var_ratio
pre_eq_no = lhs(i,1)
CASE(3) ! animal effect, both parents are known, sire or dam is diagonal
rhs_element = rhs_element - SUM(solutions(loc_in_eq)*(/-0.5,0.25/) * 2.0 * var_ratio)
diag_ele = diag_ele + 0.25 * 2.0 * var_ratio
pre_eq_no = lhs(i,1)
CASE(4) ! animal effect, one parent is known, animal is diagonal
rhs_element = rhs_element - SUM(solutions(loc_in_eq(1)) * -0.5 * (4.0/3.0) * var_ratio)
diag_ele = diag_ele + (4.0/3.0) * var_ratio
pre_eq_no = lhs(i,1)
CASE(5) ! animal effect, one parents is known, parent is diagonal
rhs_element = rhs_element - SUM(solutions(loc_in_eq(1)) * -0.5 * (4.0/3.0) * var_ratio)
diag_ele = diag_ele + 0.25 * (4.0/3.0) * var_ratio
pre_eq_no = lhs(i,1)
CASE(6) ! animal effect, no parents is known
diag_ele = diag_ele + 1.0 * var_ratio
pre_eq_no = lhs(i,1)
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
epsilon = epsilon + (pre_sol - solutions(pre_eq_no)) ** 2 ! calculate sum of squares of difference between old solution and new solutions ! write iteration number 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 = *) i, solutions(i)
END DO
CLOSE(13)END PROGRAM animal_model_solve
실행을 편하게 위한 파라미터 파일
animal_model_solve_1.par
3
3 3 17
2.33
위 프로그램을 이용하여 구한 해
sol.dat
1 24.127184
2 25.431318
3 24.230770
4 3.7396011
5 -4.4202266
6 0.51047230
7 0.27496710
8 -0.93301475
9 0.11148594
10 0.66864324
11 -0.12174667
12 0.12215573
13 -0.12114576
14 0.18701611
15 -0.45870286
16 6.00416027E-02
17 0.30996057
18 -1.3070647
19 0.94238120
20 0.37363133
21 -0.50524253
22 4.17404920E-02
23 -0.24319477
컴파일 및 실행 화면
1241166711_animal_model.dat
1241166711_animal_pedi.dat
1241166711_animal_model_setup_1.f95
1241166711_animal_model_setup_1.par
1241166711_lhs.dat
1241166711_rhs.dat
1241166711_sorted_lhs.dat
1244217206_animal_model_solve_1.f95
1241166711_animal_model_solve_1.par
1241166711_sol.dat
'Animal Breeding > Fortran program' 카테고리의 다른 글
| 혈통을 세대별로 나누기 (0) | 2009.05.25 |
|---|---|
| multiple trait sire model - lhs와 rhs를 setup하기 (0) | 2009.05.21 |
| sire model - gauss seidel 방법을 이용한 해 찾기 2 (0) | 2009.04.26 |
| animal model - gauss seidel, iteration on data (0) | 2009.04.22 |
| 근교계수(inbreeding coefficient)와 Dii 구하기 (0) | 2009.04.21 |