1 /* test_radiances.cpp - "main" function for Test 3
2
3 Copyright (C) 2012-2014 The University of Reading
4
5 Copying and distribution of this file, with or without modification,
6 are permitted in any medium without royalty provided the copyright
7 notice and this notice are preserved. This file is offered as-is,
8 without any warranty.
9 */
10
11 #include "adept.h"
12 #include "simulate_radiances.h"
13
14 using adept::Real;
15 using adept::aReal;
16
17 // This function provides an Adept interface to the simulate_radiances
18 // function
simulate_radiances_wrapper(int n,const aReal & surface_temperature,const aReal * temperature,aReal radiance[2])19 void simulate_radiances_wrapper(int n,
20 const aReal& surface_temperature,
21 const aReal* temperature,
22 aReal radiance[2]) {
23 // Create inactive (Real) versions of the active (aReal) inputs
24 Real st = value(surface_temperature);
25 std::vector<Real> t(n);
26 for (int i = 0; i < n; ++i) t[i] = value(temperature[i]);
27
28 // Declare variables to hold the inactive outputs and their Jacobians
29 Real r[2];
30 Real dr_dst[2];
31 std::vector<Real> dr_dt(2*n);
32
33 // Call the function with the non-Adept interface
34 simulate_radiances(n, st, &t[0], &r[0], dr_dst, &dr_dt[0]);
35
36 // Copy the results into the active variables, but use set_value in order
37 // not to write any equivalent derivative statement to the Adept stack
38 radiance[0].set_value(r[0]);
39 radiance[1].set_value(r[1]);
40
41 // Loop over the two radiances and add the derivative statements to
42 // the Adept stack
43 for (int i = 0; i < 2; ++i) {
44 // Add the first term on the right-hand-side of Equation 1 in the text
45 radiance[i].add_derivative_dependence(surface_temperature, dr_dst[i]);
46 // Now append the second term on the right-hand-side of Equation
47 // 1. The third argument "n" of the following function says that
48 // there are n terms to be summed, and the fourth argument "2"
49 // says to take only every second element of the Jacobian dr_dt,
50 // since the derivatives with respect to the two radiances have
51 // been interlaced. If the fourth argument is omitted then
52 // relevant Jacobian elements will be assumed to be contiguous in
53 // memory.
54 radiance[i].append_derivative_dependence(temperature, &dr_dt[i], n, 2);
55 }
56
57 for (int i = 0; i < 2; ++i) {
58 std::cout << "Channel " << i << "\n";
59 std::cout << "d[radiance]/d[surface_temperature] = " << dr_dst[i] << "\n";
60 std::cout << "d[radiance]/d[temperature] =";
61 for (int j = 0; j < n; ++j) {
62 std::cout << " " << dr_dt[i+j*2];
63 }
64 std::cout << "\n\n";
65 }
66
67 }
68
69
70 int
main(int argc,char ** argv)71 main(int argc, char** argv)
72 {
73 // Temperature (K) at 1000-m intervals from the mid-latitude summer
74 // standard atmosphere
75 static const int N_POINTS = 25;
76 static const Real temperature_profile[N_POINTS+1]
77 = {294.0, 290.0, 285.0, 279.0, 273.0, 267.0, 261.0, 255.0,
78 248.0, 242.0, 235.0, 229.0, 222.0, 216.0, 216.0, 216.0,
79 216.0, 216.0, 216.0, 217.0, 218.0, 219.0, 220.0, 222.0,
80 223.0, 224.0};
81
82 // Start the Adept stack
83 adept::Stack s;
84
85 // Copy the temperature profile information into active variables
86 aReal surface_temperature = temperature_profile[0];
87 aReal temperature[N_POINTS];
88 for (int i = 0; i < N_POINTS; i++) {
89 temperature[i] = temperature_profile[i+1];
90 }
91
92 // The simulated radiances will be put here...
93 aReal sim_radiance[2];
94
95 // ...and compared to the observed radiances here with their 1-sigma
96 // error
97 Real obs_radiance[2] = {0.00189, 0.00140};
98 Real radiance_error = 2.0e-5;
99
100 // Start recording derivative information
101 s.new_recording();
102
103 // Simulate the radiances for the input surface temperature and
104 // atmospheric temperature
105 simulate_radiances_wrapper(N_POINTS, surface_temperature,
106 temperature, sim_radiance);
107
108 std::cout << "Simulated radiances = "
109 << sim_radiance[0].value() << " "
110 << sim_radiance[1].value() << "\n";
111
112 // Compute a "cost function" (or "penalty function") expressing the
113 // sum of the squared number of error standard deviations the
114 // simulated radiances are from the observed radiances
115 aReal cost_function = 0.0;
116 for (int ichan = 0; ichan < 2; ichan++) {
117 cost_function
118 += (sim_radiance[ichan] - obs_radiance[ichan])
119 * (sim_radiance[ichan] - obs_radiance[ichan])
120 / (radiance_error*radiance_error);
121 }
122
123 std::cout << "Cost function = " << cost_function << "\n";
124
125 // We want the computed adjoints to be gradients of the cost
126 // function with respect to the surface temperature or atmospheric
127 // temperature
128 cost_function.set_gradient(1.0);
129
130 // Reverse-mode automatic differentiation
131 s.reverse();
132
133 // Extract the gradients
134 Real dcost_dsurface_temperature = 0;
135 Real dcost_dtemperature[N_POINTS];
136 surface_temperature.get_gradient(dcost_dsurface_temperature);
137 adept::get_gradients(temperature, N_POINTS, dcost_dtemperature);
138
139
140 std::cout << "d[cost_function]/d[surface_temperature] = "
141 << dcost_dsurface_temperature << "\n";
142 std::cout << "d[cost_function]/d[temperature] =";
143 for (int i = 0; i < N_POINTS; i++) {
144 std::cout << " " << dcost_dtemperature[i];
145 }
146 std::cout << "\n";
147
148
149 }
150