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:
- Convert pounds to kilograms. Use the actual conversion factor, not the approximate one. Look it up on Google.
- Convert feet/inches to meters. Look up the conversion factor, do not guess at it.
- Compute BMI.
- Determine where the user falls in the table supplied and return that informationin an appropriate form.
Write a module stats that implements the following:
- Mean of an array
- Standard deviation of an array
- 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:
- Uses your modules
- 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).
- Pass appropriate arrays to a subroutine that computes an array of BMI data based on height and weight and returns the BMI array.
- 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.
- Compute the mean and standard deviations of the observations.
- Compute the absolute values of the deviations, i.e. abs(A-mean(A))/std(A)
- Use the tails
devs=devs/sqrt(2.) - Compute the probabilities
prob=erfc(devs): erfc is an intrinsic in any fairly recent Fortran compiler. - 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
subroutine fix_missing(A,n_bad)
!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, intent(out) :: n_bad
integer :: lenA
integer :: i
lenA=size(A)
n_bad=0
do i=1,lenA
if (A(i)==MISSING) then
n_bad=n_bad+1
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
stop "Data file not found"
endif
nlines=line_count(inunit)
! Subtract header line
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))
!Read header
read(inunit,*)
!Read data
do n=1,nobs
read(inunit,*,iostat=ios) percent_bf(n),age(n),weight(n),height(n)
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
read(inunit,*,end=10)
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