Skip to content

Commit d4536b5

Browse files
jacobwilliamsjacobwilliams
authored
Merge pull request #11 from jacobwilliams/develop
added a new test case with pyplot-fortran dependency
2 parents a80d35a + 10fd6be commit d4536b5

4 files changed

Lines changed: 190 additions & 1 deletion

File tree

‎.gitignore‎

Lines changed: 5 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -41,4 +41,8 @@ Debug
4141
x64
4242

4343
# mac
44-
.DS_Store
44+
.DS_Store
45+
46+
# test artifacts:
47+
/*.pdf
48+
/*.py

‎README.md‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -89,6 +89,8 @@ For example, to build a single precision version of the library, use:
8989
fpm build --profile release --flag "-DREAL32"
9090
```
9191

92+
Note that the [pyplot-fortran](https://github.com/jacobwilliams/pyplot-fortran) library is a dependency for one of the test cases. FPM will automatically download the right version.
93+
9294
## Documentation
9395

9496
The latest API documentation can be found [here](https://jacobwilliams.github.io/NumDiff/). This was generated from the source code using [FORD](https://github.com/Fortran-FOSS-Programmers/ford) (note that the included `build.sh` script will also generate these files).

‎fpm.toml‎

Lines changed: 8 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -18,6 +18,9 @@ auto-executables = false
1818
auto-examples = false
1919
auto-tests = false
2020

21+
[dev-dependencies]
22+
pyplot-fortran = { git="https://github.com/jacobwilliams/pyplot-fortran.git", tag = "3.2.1" }
23+
2124
[[test]]
2225
name = "dsm_test"
2326
source-dir = "tests"
@@ -27,3 +30,8 @@ main = "dsm_test.f90"
2730
name = "test1"
2831
source-dir = "tests"
2932
main = "test1.f90"
33+
34+
[[test]]
35+
name = "test2"
36+
source-dir = "tests"
37+
main = "test2.f90"

‎tests/test2.f90‎

Lines changed: 175 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,175 @@
1+
!*******************************************************************************
2+
!> author: Jacob Williams
3+
! date: October 1, 2022
4+
!
5+
! Test2 for the numerical differentiation module.
6+
! Makes a plot of the errors using different methods and step sizes.
7+
8+
program test2
9+
10+
use iso_fortran_env
11+
use numerical_differentiation_module
12+
use numdiff_kinds_module, only: wp
13+
use pyplot_module
14+
15+
implicit none
16+
17+
integer,parameter :: n = 1 !! number of variables
18+
integer,parameter :: m = 1 !! number of functions
19+
real(wp),dimension(n),parameter :: x = 1.0_wp !! point at which to compute the derivative
20+
real(wp),dimension(n),parameter :: xlow = -100000.0_wp !! bounds not really needed
21+
real(wp),dimension(n),parameter :: xhigh = 100000.0_wp !!
22+
integer,parameter :: perturb_mode = 1 !! absolute step
23+
integer,parameter :: cache_size = 0 !! `0` indicates not to use cache
24+
integer,parameter :: sparsity_mode = 1 !! assume dense
25+
integer,dimension(*),parameter :: methods = [1,3,10,21,36,500,600,700,800]
26+
!! array of method IDs:
27+
!! [1,3,10,21,36,500,600,700,800]
28+
!! forward + all the central diff ones
29+
integer,parameter :: exp_star = -ceiling(log10(epsilon(1.0_wp))) !! exponent to start with
30+
integer,parameter :: exp_stop = -2 !! exponent to end with
31+
integer,parameter :: exp_step = -1 !! ecponent step
32+
integer,parameter :: exp_scale = 5 !! number of substeps from one to the next
33+
34+
type(numdiff_type) :: my_prob !! main class to compute the derivatives
35+
integer :: i !! method counter
36+
integer :: j !! counter
37+
real(wp),dimension(:),allocatable :: jac !! jacobian
38+
integer :: func_evals !! function evaluation counter
39+
character(len=:),allocatable :: error_msg !! error message string
40+
real(wp),dimension(n) :: dpert !! perturbation step size
41+
integer :: ipert !! perturbation step size counter
42+
real(wp) :: error !! diff from true derivative
43+
integer :: num_dperts !! number of dperts to test
44+
integer :: num_methods !! number of methods to test
45+
real(wp),dimension(:),allocatable :: results_dpert !! results array - dpert
46+
real(wp),dimension(:),allocatable :: results_errors !! results array - errors
47+
type(pyplot) :: plt !! for plotting the results
48+
character(len=:),allocatable :: formula !! finite diff forumla for the plot legend
49+
integer :: idx !! index in results arrays
50+
character(len=:),allocatable :: real_kind_str !! real kind for the plot title
51+
real(wp),dimension(3) :: color !! line color array
52+
real(wp),dimension(2) :: ylim !! plot y limit array
53+
54+
! size the arrays:
55+
num_methods = size(methods)
56+
num_dperts = size([(i, i = exp_star*exp_scale, exp_stop*exp_scale, exp_step)])
57+
allocate(results_errors(num_dperts)); results_errors = -huge(1.0_wp)
58+
59+
! for plot title:
60+
select case(wp)
61+
case(REAL32); real_kind_str = '[Single Precision]'
62+
case(REAL64); real_kind_str = '[Double Precision]'
63+
case(REAL128); real_kind_str = '[Quad Precision]'
64+
case default; error stop 'Invalid real kind'
65+
end select
66+
67+
! initialize the plot:
68+
call plt%initialize(grid = .true.,&
69+
figsize = [20,10],&
70+
axes_labelsize = 30, &
71+
xtick_labelsize = 30, &
72+
ytick_labelsize = 30, &
73+
font_size = 30, &
74+
xlabel = 'Finite Difference Perturbation Step Size $h$',&
75+
ylabel = 'Finite Difference Derivative Error',&
76+
title = 'Derivative of $x + \sin(x)$ at $x=1$ '//real_kind_str,&
77+
legend = .true., &
78+
legend_fontsize = 10,&
79+
usetex = .true.)
80+
81+
! try different finite diff methods
82+
do j=1,size(methods)
83+
84+
idx = 0
85+
i = methods(j) ! method id
86+
call get_finite_diff_formula(i,formula)
87+
88+
! cycle through perturbation step sizes:
89+
do ipert = exp_star*exp_scale, exp_stop*exp_scale, exp_step
90+
91+
idx = idx + 1 ! index for arrays
92+
dpert = 10.0_wp ** (-ipert/real(exp_scale,wp)) ! compute perturbation step size
93+
if (j==1) then
94+
if (.not. allocated(results_dpert)) allocate(results_dpert(0))
95+
results_dpert = [results_dpert, dpert(1)] ! save dpert
96+
end if
97+
98+
func_evals = 0
99+
call my_prob%destroy()
100+
call my_prob%initialize(n,m,xlow,xhigh,perturb_mode,dpert,&
101+
problem_func=my_func,&
102+
sparsity_mode=sparsity_mode,&
103+
jacobian_method=i,&
104+
partition_sparsity_pattern=.false.,&
105+
cache_size=cache_size)
106+
if (my_prob%failed()) then
107+
call my_prob%get_error_status(error_msg=error_msg)
108+
error stop error_msg
109+
end if
110+
111+
call my_prob%compute_jacobian(x,jac)
112+
if (my_prob%failed()) then
113+
call my_prob%get_error_status(error_msg=error_msg)
114+
error stop error_msg
115+
end if
116+
error = abs(deriv(x(1)) - jac(1))
117+
!write(output_unit,'(I5,1X,*(F27.16))') i , dpert, error
118+
if (error<=epsilon(1.0_wp)) error = 0.0_wp
119+
results_errors(idx) = error ! save result
120+
121+
end do
122+
123+
! line color for the plot:
124+
select case (j)
125+
case(1); color = [1.0_wp, 0.0_wp, 0.0_wp] ! red
126+
case(2); color = [0.0_wp, 1.0_wp, 0.0_wp] ! green
127+
case default
128+
! blue gradient for the others
129+
color = [real(j-2,wp)/num_methods, &
130+
real(j-2,wp)/num_methods, &
131+
0.9_wp]
132+
end select
133+
ylim = [10.0_wp**(ceiling(log10(epsilon(1.0_wp)))), 1.0_wp]
134+
135+
! plot for this method:
136+
call plt%add_plot(results_dpert,results_errors,&
137+
xscale='log', yscale='log',&
138+
label=formula,linestyle='.-',markersize=5,linewidth=2, &
139+
color = color,&
140+
xlim = ylim,&
141+
ylim = ylim)
142+
143+
end do
144+
145+
! save plot:
146+
call plt%savefig('results '//real_kind_str//'.pdf', pyfile='results.py')
147+
148+
contains
149+
150+
function func(x)
151+
!! Problem function
152+
real(wp), intent(in) :: x
153+
real(wp) :: func
154+
func = x + sin(x)
155+
end function func
156+
function deriv(x)
157+
!! Problem function true derivative
158+
real(wp), intent(in) :: x
159+
real(wp) :: deriv
160+
deriv = 1.0_wp + cos(x)
161+
end function deriv
162+
163+
subroutine my_func(me,x,f,funcs_to_compute)
164+
!! Problem function interface for numdiff
165+
class(numdiff_type),intent(inout) :: me
166+
real(wp),dimension(:),intent(in) :: x
167+
real(wp),dimension(:),intent(out) :: f
168+
integer,dimension(:),intent(in) :: funcs_to_compute
169+
if (any(funcs_to_compute==1)) f(1) = func(x(1))
170+
func_evals = func_evals + 1
171+
end subroutine my_func
172+
173+
!*******************************************************************************
174+
end program test2
175+
!*******************************************************************************

0 commit comments

Comments
 (0)