Project 4

Download the file bodyfat.csv. This is a dataset of body fat, age, height, and weight for a set of participants in a study. BMI categories are as follows:

|Severely underweight | BMI < 16.0 | |Underweight | 16 <= BMI < 18.5 | |Normal | 18.5 <= BMI < 25 | |Overweight | 25 <= BMI < 30 | |Obese Class I | 30 <= BMI < 35 | |Obese Class II | 35 <= BMI < 40 | |Obese Class III | BMI > 40 |

Write a bmi_calculator module containing functions/subroutines for the following:

1. Convert pounds to kilograms. Use the actual conversion factor, not the approximate one. Look it up on Google.
2. Convert feet/inches to meters. Look up the conversion factor, do not guess at it.
3. Compute BMI.
4. Determine where the user falls in the table supplied and return that informationin an appropriate form.

Write a module stats that implements the following:

1. Mean of an array
2. Standard deviation of an array
3. Outlier rejection using Chauvenet’s criterion. Pseudocode given further down. Make as much use of Fortran intrinsics/array operations as you can.

Write a main program that implements the following:

2. Reads the input file into appropriate allocatable arrays (use one-dimensional arrays for this project). Don’t assume you know the length of the file (but you can assume the number of header lines is fixed).
3. Pass appropriate arrays to a subroutine that computes an array of BMI data based on height and weight and returns the BMI array.
4. Rejects the outlier(s). The function should return an array of logicals that you can apply to the original data using WHERE or similar. Create new arrays with the outlier(s) deleted.

Write a file that contains the corrected data for bodyfat and BMI. Use Excel or whatever you normally use to plot BMI as a function of percentage body fat. Be sure to plot it as a scatter plot (points only, no connecting lines).

Chauvenet’s criterion: It’s not the state of the art but works pretty well.

1. Compute the mean and standard deviations of the observations.
2. Compute the absolute values of the deviations, i.e. abs(A-mean(A))/std(A)
3. Use the tails devs=devs/sqrt(2.)
4. Compute the probabilities prob=erfc(devs) : erfc is an intrinsic in any fairly recent Fortran compiler.
5. The criterion is that we retain data with prob>=1./(2*N_obs) (number of observations).
Example solution

module stats
implicit none
!This module contains some simple functions to work with imperfect data
!  Author:    Katherine Holcomb
!  Changelog: Initial code 20120220
!             Minor improvements 20160307

integer, parameter  :: rk=kind(1.0)
real(rk), parameter :: MISSING=-999._rk

contains

function mean(A)
! Computes the mean of an array A
real(rk)                           :: mean
real(rk), intent(in), dimension(:) :: A
integer                            :: lenA

lenA=size(A)
mean=sum(A)/size(A)

end function mean

function std(A)
! Computes the standard deviation of an array A
real(rk)                           :: std
real(rk), intent(in), dimension(:) :: A
integer                            :: lenA

lenA=size(A)
std =sqrt((sum(A**2)/lenA)-mean(A)**2)

end function std

function reject_outliers(A)  result(valid)
!Returns the indices of A with true for valid, false for outlier
! A must be inout if it is passed to something else as an argument
! (since the compiler doesn't try to trace the whole call tree to
! figure out whether it's changed or not).
real(rk),    dimension(:),       intent(inout) :: A
logical,     dimension(size(A))                :: valid
real(rk),    dimension(size(A))                :: devs, prob
real(rk)                                       :: mean_A, stdv_A
real(rk)                                       :: criterion
integer                                        :: lenA

mean_A=mean(A)
stdv_A=std(A)
lenA = size(A)
criterion = 1.0_rk/(2*lenA)     ! Chauvenet's criterion
devs=abs(A-mean_A)/stdv_A
devs = devs/sqrt(2.0_rk)        ! The left and right tail threshold values
prob = erfc(devs)
valid=prob>=criterion

end function reject_outliers

!Returns A with actual data substituted for missing data
!Relying on a side effect here, hence we use a subroutine.
real,    dimension(:),       intent(inout) :: A
integer                                    :: lenA
integer                                    :: i

lenA=size(A)
do i=1,lenA
if (A(i)==MISSING) then
if (i==0) then
A(i)=A(i+1)
else
A(i)=A(i-1)
endif
endif
enddo

end subroutine fix_missing

end module stats


module bmi_calculator
implicit none

!This module computes BMI from height and weight.
!Units are converted from metric (kg, cm) or Imperial (ft/in, lb) to
!to (kg,m).

!Globals

integer, parameter               :: rk=kind(1.0)
real(rk)                         :: invalid= 999._rk
real(rk)                         :: signal =-999._rk

contains

!--------------------------------------------------------------------------
!

elemental function cm_to_m(l)
real(rk)             :: cm_to_m
real(rk), intent(in) :: l
cm_to_m=0.01*l
end function cm_to_m

elemental function inch_to_m(l)
real(rk)             :: inch_to_m
real(rk), intent(in) :: l
real(rk), parameter  :: inch2m=0.0254_rk
inch_to_m=inch2m*l
end function inch_to_m

elemental function lb_to_kg(l)
real(rk)             :: lb_to_kg
real(rk), intent(in) :: l
real(rk), parameter  :: lb2kg=0.45359237_rk
lb_to_kg=lb2kg*l
end function lb_to_kg

subroutine imperial_to_metric(l_inch,wt_lb,l_m,wt_kg)
real(rk), intent(in) :: l_inch, wt_lb
real(rk), intent(out):: l_m, wt_kg
l_m   = inch_to_m(l_inch)
wt_kg = lb_to_kg(wt_lb)
end subroutine imperial_to_metric

elemental function BMI(ht,wt)
real(rk)             :: BMI
real(rk), intent(in) :: ht,wt
BMI=wt/ht**2
end function BMI

function BMI_table(bmi_val)
integer          :: BMI_table
real(rk), intent(in) :: bmi_val

if (bmi_val <  16.0 )                      BMI_table=1
if ( bmi_val>= 16.0 .and. bmi_val < 18.5 ) BMI_table=2
if ( bmi_val>= 18.5 .and. bmi_val < 25.0 ) BMI_table=3
if ( bmi_val>= 25.0 .and. bmi_val < 30.0 ) BMI_table=4
if ( bmi_val>= 30.0 .and. bmi_val < 35.0 ) BMI_table=5
if ( bmi_val>= 35.0 .and. bmi_val < 40.0 ) BMI_table=6
if (bmi_val >= 40.0 )                      BMI_table=7
end function BMI_table

function is_inrange(inval,bounds)
logical                       :: is_inrange
real(rk),              intent(in) :: inval
real(rk),dimension(2), intent(in) :: bounds
if ( bounds(1) <= inval .and. inval <= bounds(2) ) then
is_inrange=.true.
else
is_inrange=.false.
endif
end function is_inrange

end module bmi_calculator


program bmi_data
use bmi_calculator
use stats

! This program accepts data on weight, height, and body fat, computes the
!   corresponding BMIs, tabulates the BMI, and outputs the data
!   (excluding outliers).

real, dimension(:), allocatable     :: percent_bf, age, weight, height
real, dimension(:), allocatable     :: cleaned_bf, cleaned_bmi
real, dimension(:), allocatable     :: wt_kg, ht_m, bmi_values
logical, dimension(:), allocatable  :: valid
integer                                 :: nargs, inunit
integer                                 :: nlines, nobs
logical                                 :: file_exists
character(len=40)                       :: infile
interface
function get_unit() result(unit_number)
implicit none
integer                           :: unit_number
end function get_unit

function line_count(unit_number)
implicit none
integer                           :: line_count
integer, intent(in)               :: unit_number
end function line_count

subroutine convert_imp_metr(nobs,height,weight,ht_m,wt_kg)
use bmi_calculator
implicit none
integer,                   intent(in)    :: nobs
real, dimension(nobs), intent(in)    :: height, weight
real, dimension(nobs), intent(inout) :: ht_m,   wt_kg
end subroutine convert_imp_metr
end interface

!Read file name from the command line
nargs=command_argument_count()
if (nargs .ne. 1) then
stop "No data file specified"
else
call get_command_argument(1,infile)
endif

inquire(file=infile,exist=file_exists)

if (file_exists) then
inunit=get_unit()
open(inunit,file=infile,iostat=ios)
if ( ios .ne. 0 ) then
stop "Unable to open data file"
endif
else
endif

nlines=line_count(inunit)

nobs=nlines-1

allocate(percent_bf(nobs), age(nobs), weight(nobs), height(nobs))
allocate(bmi_values(nobs), wt_kg(nobs),   ht_m(nobs))
allocate(valid(nobs))

do n=1,nobs
if ( ios .ne. 0 ) stop "Error reading file"
enddo

!Convert units
call convert_imp_metr(nobs,height,weight,ht_m,wt_kg)

!Compute BMI array
bmi_values=BMI(ht_m,wt_kg)

!Reject outliers
valid=reject_outliers(bmi_values)
nvalid=count(valid)
allocate(cleaned_bf(nvalid),cleaned_bmi(nvalid))
cleaned_bf =pack(percent_bf,valid)
cleaned_bmi=pack(bmi_values,valid)

!Write data to CSV file for plotting
call dump_results(infile,nvalid,cleaned_bf,cleaned_bmi)

end program bmi_data

function get_unit() result(unit_number)
implicit none
integer               :: unit_number
integer               :: base_unit=10
logical               :: first=.true.

if ( first ) then
unit_number=base_unit
first=.false.
else
unit_number=unit_number+1
endif
end function get_unit

function line_count(inunit)
implicit none
integer             :: line_count
integer, intent(in) :: inunit
integer             :: nlines
nlines=0
do
nlines=nlines+1
enddo
10 continue

rewind(inunit)
line_count=nlines
end function line_count

subroutine convert_imp_metr(nobs,height,weight,ht_m,wt_kg)
use bmi_calculator
implicit none
integer,               intent(in)    :: nobs
real, dimension(nobs), intent(in)    :: height, weight
real, dimension(nobs), intent(inout) :: ht_m,   wt_kg
integer                              :: i

do i=1,nobs
call imperial_to_metric(height(i),weight(i),ht_m(i),wt_kg(i))
enddo

end subroutine convert_imp_metr

subroutine dump_results(infile,nobs,percent_bf,bmi_values)
implicit none
integer,               intent(in)  :: nobs
real, dimension(nobs), intent(in)  :: percent_bf, bmi_values
integer                            :: period
integer                            :: outunit, ios
integer                            :: n
character(len=40)                  :: infile,outfile
interface
function get_unit() result(unit_number)
implicit none
integer                           :: unit_number
end function get_unit
end interface
!Prepare file to plot with plotting software.
!Using the same name as infile but a different suffix.
outunit=get_unit()
period=index(infile,'.')
if ( period > 0 ) then
outfile=infile(1:period-1)//"_plt.csv"
else
outfile=infile//"_plt.csv"
endif
open(outunit,file=outfile,iostat=ios)
if ( ios .ne. 0 ) then
stop "Unable to open output file"
endif

do n=1,nobs
write(outunit,'(f8.2,a,f8.2)') percent_bf(n),",",bmi_values(n)
enddo
end subroutine dump_results



Previous
Next