FINDSENSORLAGPARAMS Sensor lag parameter estimation for profiles.
Syntax:
PARAMS = FINDSENSORLAGPARAMS(TIME1, DEPTH1, DATA1, TIME2, DEPTH2, DATA2)
PARAMS = FINDSENSORLAGPARAMS(TIME1, DEPTH1, DATA1, FLOW1, TIME2, DEPTH2, DATA2, FLOW2)
PARAMS = FINDSENSORLAGPARAMS(..., OPTIONS)
PARAMS = FINDSENSORLAGPARAMS(..., OPT1, VAL1, ...)
[PARAMS, EXITFLAG, RESIDUAL] = FINDSENSORLAGPARAMS(...)
Description:
PARAMS = FINDSENSORLAGPARAMS(TIME1, DEPTH1, DATA1, TIME2, DEPTH2, DATA2)
finds the sensor lag parameter for two profiles with constant flow speed
given by vectors TIME1, DEPTH1 and DATA1, and TIME2, DEPTH2 and DATA2.
The profiles are supposed to measure the same column of water in opposite
directions.
Under the assumption that both profiles should be as similar as possible,
the function solves the minimization problem of finding the time constant
such that the area between corrected profiles returned by CORRECTSENSORLAG
is minimal. This problem is solved by function FMINCON using default
values for the parameter bounds. See OPTIONS description below.
PARAMS = FINDSENSORLAGPARAMS(TIME1, DEPTH1, DATA1, FLOW1, TIME2, DEPTH2, DATA2, FLOW2)
performs the same estimation but for a pair profiles with variable flow
speed given by respective vectors FLOW1 and FLOW2. The estimated parameters
are returned in a two element vector PARAMS with the offset and the slope
of the sensor lag parameter with respect to the inverse flow speed.
PARAMS = FINDSENSORLAGPARAMS(..., OPTIONS) and
PARAMS = FINDSENSORLAGPARAMS(..., OPT1, VAL1, ...) allow passing extra
options given in key-value pairs OPT1, VAL1... or in a struct OPTIONS
with field names as option keys and field values as option values.
Recognized options are:
GRAPHICS: whether graphic output should be produced.
A boolean setting whether a nice figure showing the evolution of the
minimization process will be displayed. It includes plots for the
parameter value, the objective function value, a value-time diagram of
both sequences, and a depth-value diagram.
Default value: false.
GUESS: initial guess for minimization function FMINCON.
A one or two element vector with the initial guess for each parameter.
Default value:
For constant flow speed: 0.5
For variable flow speed: [0.3568 0.07]
LOWER: lower bounds of parameters for minimization function FMINCON.
A one or two element vector with the lower bound for each parameter.
Default value:
For constant flow speed: 0 (no correction)
For variable flow speed: [0 0] (no correction)
UPPER: upper bounds of parameters for minimization function FMINCON.
A one or two element vector with the upper bound for each parameter.
Default value:
For constant flow speed: 16
For variable flow speed: [16 7.5]
OPTIMOPTS: extra options for the minimization function FMINCON.
An option struct as needed by the function FMINCON.
Default value: default options for FMINCON, except for:
'Algorithm' : 'interior-point'
'FinDiffType': 'central'
'TolFun' : 1e-4
'TolCon' : 1e-5
'TolX' : 1e-5
'Display' : 'off'
[PARAMS, EXITFLAG, RESIDUAL] = FINDSENSORLAGPARAMS(...) also returns the
exit code of the minimization function FMINCON in EXITFLAG, and the
residual area in RESIDUAL. EXITFLAG is positive when minimization succeeds.
Notes:
This function is an improved version of a previous function by Tomeu Garau,
called ADJUSTTIMECONSTANT. He is the true glider man. Main changes are:
- Support for dynamic senosr lag parameter for sequences with variable
flow speed.
- Different minimization algorithm (interior-point instead of active-set).
- Support for custom minimization options.
- Optional predefined graphical output.
It remains to be assessed that the time constant found minimizes the area
between the corrected profiles globally, because the solver looks for local
minimizers only.
Examples:
% Constant flow speed profiles (pumped CTD):
data = load('private/test/ctd_sharp_pumped.dat');
profile = data(:, 10);
time1 = data(profile==3, 1);
temp1 = data(profile==3, 3);
pres1 = data(profile==3, 4);
time2 = data(profile==4, 1);
temp2 = data(profile==4, 3);
pres2 = data(profile==4, 4);
params = findSensorLagParams(time1, pres1, temp1, time2, pres2, temp2)
% Variable flow speed profiles (unpumped CTD):
data = load('private/test/ctd_sharp_unpumped.dat');
profile = data(:, 10);
time1 = data(profile==17, 1);
temp1 = data(profile==17, 3);
pres1 = data(profile==17, 4);
ptch1 = data(profile==17, 5);
lati1 = data(profile==17, 8);
dpth1 = sw_dpth(pres1, lati1);
flow1 = computeCTDFlowSpeed(time1, dpth1, ptch1, 'minpitch', deg2rad(11));
time2 = data(profile==18, 1);
temp2 = data(profile==18, 3);
pres2 = data(profile==18, 4);
ptch2 = data(profile==18, 5);
lati2 = data(profile==18, 8);
dpth2 = sw_dpth(pres2, lati2);
flow2 = computeCTDFlowSpeed(time2, dpth2, ptch2, 'minpitch', deg2rad(11));
params = findSensorLagParams(time1, pres1, temp1, flow1, ...
time2, pres2, temp2, flow2)
% Variable flow speed profiles with exit code, residual and extra options.
[params, exitflag, residual] = ...
findSensorLagParams(time1, pres1, temp1, flow1, ...
time2, pres2, temp2, flow2, ...
'graphics', 'true')
See also:
FMINCON
OPTIMSET
CORRECTSENSORLAG
COMPUTECTDFLOWSPEED
Authors:
Joan Pau Beltran <joanpau.beltran@socib.cat>0001 function [params, exitflag, residual] = findSensorLagParams(varargin) 0002 %FINDSENSORLAGPARAMS Sensor lag parameter estimation for profiles. 0003 % 0004 % Syntax: 0005 % PARAMS = FINDSENSORLAGPARAMS(TIME1, DEPTH1, DATA1, TIME2, DEPTH2, DATA2) 0006 % PARAMS = FINDSENSORLAGPARAMS(TIME1, DEPTH1, DATA1, FLOW1, TIME2, DEPTH2, DATA2, FLOW2) 0007 % PARAMS = FINDSENSORLAGPARAMS(..., OPTIONS) 0008 % PARAMS = FINDSENSORLAGPARAMS(..., OPT1, VAL1, ...) 0009 % [PARAMS, EXITFLAG, RESIDUAL] = FINDSENSORLAGPARAMS(...) 0010 % 0011 % Description: 0012 % PARAMS = FINDSENSORLAGPARAMS(TIME1, DEPTH1, DATA1, TIME2, DEPTH2, DATA2) 0013 % finds the sensor lag parameter for two profiles with constant flow speed 0014 % given by vectors TIME1, DEPTH1 and DATA1, and TIME2, DEPTH2 and DATA2. 0015 % The profiles are supposed to measure the same column of water in opposite 0016 % directions. 0017 % 0018 % Under the assumption that both profiles should be as similar as possible, 0019 % the function solves the minimization problem of finding the time constant 0020 % such that the area between corrected profiles returned by CORRECTSENSORLAG 0021 % is minimal. This problem is solved by function FMINCON using default 0022 % values for the parameter bounds. See OPTIONS description below. 0023 % 0024 % PARAMS = FINDSENSORLAGPARAMS(TIME1, DEPTH1, DATA1, FLOW1, TIME2, DEPTH2, DATA2, FLOW2) 0025 % performs the same estimation but for a pair profiles with variable flow 0026 % speed given by respective vectors FLOW1 and FLOW2. The estimated parameters 0027 % are returned in a two element vector PARAMS with the offset and the slope 0028 % of the sensor lag parameter with respect to the inverse flow speed. 0029 % 0030 % PARAMS = FINDSENSORLAGPARAMS(..., OPTIONS) and 0031 % PARAMS = FINDSENSORLAGPARAMS(..., OPT1, VAL1, ...) allow passing extra 0032 % options given in key-value pairs OPT1, VAL1... or in a struct OPTIONS 0033 % with field names as option keys and field values as option values. 0034 % Recognized options are: 0035 % GRAPHICS: whether graphic output should be produced. 0036 % A boolean setting whether a nice figure showing the evolution of the 0037 % minimization process will be displayed. It includes plots for the 0038 % parameter value, the objective function value, a value-time diagram of 0039 % both sequences, and a depth-value diagram. 0040 % Default value: false. 0041 % GUESS: initial guess for minimization function FMINCON. 0042 % A one or two element vector with the initial guess for each parameter. 0043 % Default value: 0044 % For constant flow speed: 0.5 0045 % For variable flow speed: [0.3568 0.07] 0046 % LOWER: lower bounds of parameters for minimization function FMINCON. 0047 % A one or two element vector with the lower bound for each parameter. 0048 % Default value: 0049 % For constant flow speed: 0 (no correction) 0050 % For variable flow speed: [0 0] (no correction) 0051 % UPPER: upper bounds of parameters for minimization function FMINCON. 0052 % A one or two element vector with the upper bound for each parameter. 0053 % Default value: 0054 % For constant flow speed: 16 0055 % For variable flow speed: [16 7.5] 0056 % OPTIMOPTS: extra options for the minimization function FMINCON. 0057 % An option struct as needed by the function FMINCON. 0058 % Default value: default options for FMINCON, except for: 0059 % 'Algorithm' : 'interior-point' 0060 % 'FinDiffType': 'central' 0061 % 'TolFun' : 1e-4 0062 % 'TolCon' : 1e-5 0063 % 'TolX' : 1e-5 0064 % 'Display' : 'off' 0065 % 0066 % [PARAMS, EXITFLAG, RESIDUAL] = FINDSENSORLAGPARAMS(...) also returns the 0067 % exit code of the minimization function FMINCON in EXITFLAG, and the 0068 % residual area in RESIDUAL. EXITFLAG is positive when minimization succeeds. 0069 % 0070 % Notes: 0071 % This function is an improved version of a previous function by Tomeu Garau, 0072 % called ADJUSTTIMECONSTANT. He is the true glider man. Main changes are: 0073 % - Support for dynamic senosr lag parameter for sequences with variable 0074 % flow speed. 0075 % - Different minimization algorithm (interior-point instead of active-set). 0076 % - Support for custom minimization options. 0077 % - Optional predefined graphical output. 0078 % 0079 % It remains to be assessed that the time constant found minimizes the area 0080 % between the corrected profiles globally, because the solver looks for local 0081 % minimizers only. 0082 % 0083 % Examples: 0084 % % Constant flow speed profiles (pumped CTD): 0085 % data = load('private/test/ctd_sharp_pumped.dat'); 0086 % profile = data(:, 10); 0087 % time1 = data(profile==3, 1); 0088 % temp1 = data(profile==3, 3); 0089 % pres1 = data(profile==3, 4); 0090 % time2 = data(profile==4, 1); 0091 % temp2 = data(profile==4, 3); 0092 % pres2 = data(profile==4, 4); 0093 % params = findSensorLagParams(time1, pres1, temp1, time2, pres2, temp2) 0094 % 0095 % % Variable flow speed profiles (unpumped CTD): 0096 % data = load('private/test/ctd_sharp_unpumped.dat'); 0097 % profile = data(:, 10); 0098 % time1 = data(profile==17, 1); 0099 % temp1 = data(profile==17, 3); 0100 % pres1 = data(profile==17, 4); 0101 % ptch1 = data(profile==17, 5); 0102 % lati1 = data(profile==17, 8); 0103 % dpth1 = sw_dpth(pres1, lati1); 0104 % flow1 = computeCTDFlowSpeed(time1, dpth1, ptch1, 'minpitch', deg2rad(11)); 0105 % time2 = data(profile==18, 1); 0106 % temp2 = data(profile==18, 3); 0107 % pres2 = data(profile==18, 4); 0108 % ptch2 = data(profile==18, 5); 0109 % lati2 = data(profile==18, 8); 0110 % dpth2 = sw_dpth(pres2, lati2); 0111 % flow2 = computeCTDFlowSpeed(time2, dpth2, ptch2, 'minpitch', deg2rad(11)); 0112 % params = findSensorLagParams(time1, pres1, temp1, flow1, ... 0113 % time2, pres2, temp2, flow2) 0114 % 0115 % % Variable flow speed profiles with exit code, residual and extra options. 0116 % [params, exitflag, residual] = ... 0117 % findSensorLagParams(time1, pres1, temp1, flow1, ... 0118 % time2, pres2, temp2, flow2, ... 0119 % 'graphics', 'true') 0120 % 0121 % See also: 0122 % FMINCON 0123 % OPTIMSET 0124 % CORRECTSENSORLAG 0125 % COMPUTECTDFLOWSPEED 0126 % 0127 % Authors: 0128 % Joan Pau Beltran <joanpau.beltran@socib.cat> 0129 0130 % Copyright (C) 2013-2016 0131 % ICTS SOCIB - Servei d'observacio i prediccio costaner de les Illes Balears 0132 % <http://www.socib.es> 0133 % 0134 % This program is free software: you can redistribute it and/or modify 0135 % it under the terms of the GNU General Public License as published by 0136 % the Free Software Foundation, either version 3 of the License, or 0137 % (at your option) any later version. 0138 % 0139 % This program is distributed in the hope that it will be useful, 0140 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0141 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0142 % GNU General Public License for more details. 0143 % 0144 % You should have received a copy of the GNU General Public License 0145 % along with this program. If not, see <http://www.gnu.org/licenses/>. 0146 0147 error(nargchk(6, 18, nargin, 'struct')); 0148 0149 0150 %% Parse basic input arguments. 0151 % Get numeric (non-option) arguments. 0152 nargnum = find(~cellfun(@isnumeric, varargin), 1, 'first') - 1; 0153 if isempty(nargnum) 0154 nargnum = nargin; 0155 end 0156 switch(nargnum) 0157 case 6 0158 % Constant flow speed (pumped CTD). 0159 [time1, depth1, data1] = varargin{1:3}; 0160 [time2, depth2, data2] = varargin{4:6}; 0161 constant_flow = true; 0162 case 8 0163 % Variable flow speed (unpumped CTD). 0164 [time1, depth1, data1, flow1] = varargin{1:4}; 0165 [time2, depth2, data2, flow2] = varargin{5:8}; 0166 constant_flow = false; 0167 end 0168 0169 0170 %% Configure default options. 0171 options.graphics = false; 0172 if constant_flow 0173 options.guess = 0.5; % from original implementation. 0174 options.lower = 0; % no correction. 0175 options.upper = 16; % from original implementation. 0176 else 0177 options.guess = [0.3568 0.07]; 0178 options.lower = [0 0]; 0179 options.upper = [16 7.5]; 0180 end 0181 options.optimopts = optimset(optimset('fmincon'), ... 0182 'Algorithm', 'interior-point', ... 0183 'TolFun', 1e-4, 'TolCon', 1e-5, 'TolX', 1e-5, ... 0184 'FinDiffType', 'central', ... 0185 'Display', 'off'); 0186 0187 0188 %% Parse extra arguments. 0189 % Get option key-value pairs in any accepted call signature. 0190 argopts = varargin(nargnum+1:end); 0191 if isscalar(argopts) && isstruct(argopts{1}) 0192 % Options passed as a single option struct argument: 0193 % field names are option keys and field values are option values. 0194 opt_key_list = fieldnames(argopts{1}); 0195 opt_val_list = struct2cell(argopts{1}); 0196 elseif mod(numel(argopts), 2) == 0 0197 % Options passed as key-value argument pairs. 0198 opt_key_list = argopts(1:2:end); 0199 opt_val_list = argopts(2:2:end); 0200 else 0201 error('glider_toolbox:findSensorLagParams:InvalidOptions', ... 0202 'Invalid optional arguments (neither key-value pairs nor struct).'); 0203 end 0204 % Overwrite default options with values given in extra arguments. 0205 for opt_idx = 1:numel(opt_key_list) 0206 opt = lower(opt_key_list{opt_idx}); 0207 val = opt_val_list{opt_idx}; 0208 if isfield(options, opt) 0209 options.(opt) = val; 0210 else 0211 error('glider_toolbox:findSensorLagParams:InvalidOption', ... 0212 'Invalid option: %s.', opt); 0213 end 0214 end 0215 0216 0217 %% Enable graphical output, if needed. 0218 % The functions are defined below. 0219 if options.graphics 0220 if constant_flow 0221 data_cor1 = correctSensorLag(time1, data1, options.guess); 0222 data_cor2 = correctSensorLag(time2, data2, options.guess); 0223 else 0224 data_cor1 = correctSensorLag(time1, data1, flow1, options.guess); 0225 data_cor2 = correctSensorLag(time2, data2, flow2, options.guess); 0226 end 0227 options.optimopts.OutputFcn = @optimoutUpdateCorrectedData; 0228 options.optimopts.PlotFcns = ... 0229 {@optimplotx @optimplotfval @optimplotDataTime @optimplotDepthData}; 0230 end 0231 0232 0233 %% Perform estimation through minimization. 0234 % Definition of minimization objective function. 0235 objective_function = @optimobjDepthValueArea; 0236 0237 % Run minimization procedure. 0238 [params, residual, exitflag] = ... 0239 fmincon(objective_function, options.guess, ... 0240 [], [], [], [], options.lower, options.upper, [], ... 0241 options.optimopts); 0242 0243 % This is the equivalent version of the original code by Tomeu Garau, 0244 % with external constant functions inlined and coding style adaptions: 0245 %{ 0246 options = optimset(optimset('fmincon'), ... 0247 'LargeScale', 'off', ... 0248 'Algorithm', 'active-set', ... 0249 'TolFun', 1e-4, ... 0250 'TolCon', 1e-5, ... 0251 'TolX', 1e-5, ... 0252 ... % 'Plotfcns', [], % {@optimplotfval, @optimplotfirstorderopt, @optimplotx}); 0253 'Display', 'off'); 0254 first_guess = 0.5; 0255 upper_bound = 16; 0256 lower_bound = eps; 0257 constant = fmincon(@objective_function, first_guess, [], [], [], [], ... 0258 lower_bound, upper_bound, [], options); 0259 %} 0260 0261 0262 %% Definition of auxiliary objective and plotting functions. 0263 % They should be nested to access cast data. 0264 function area = optimobjDepthValueArea(params) 0265 %OPTIMOBJDEPTHVALUEAREA Compute area enclosed by profiles in depth-value diagram. 0266 if constant_flow 0267 data_cor1 = correctSensorLag(time1, data1, params); 0268 data_cor2 = correctSensorLag(time2, data2, params); 0269 else 0270 data_cor1 = correctSensorLag(time1, data1, flow1, params); 0271 data_cor2 = correctSensorLag(time2, data2, flow2, params); 0272 end 0273 depth_min = max(min(depth1), min(depth2)); 0274 depth_max = min(max(depth1), max(depth2)); 0275 depth_mask1 = (depth_min <= depth1) & (depth1 <= depth_max); 0276 depth_mask2 = (depth_min <= depth2) & (depth2 <= depth_max); 0277 min_idx1 = find(depth_mask1, 1, 'first'); 0278 min_idx2 = find(depth_mask2, 1, 'first'); 0279 max_idx1 = find(depth_mask1, 1, 'last'); 0280 max_idx2 = find(depth_mask2, 1, 'last'); 0281 area = ... 0282 profileArea(data_cor1(min_idx1:max_idx1), depth1(min_idx1:max_idx1), ... 0283 data_cor2(min_idx2:max_idx2), depth2(min_idx2:max_idx2)); 0284 end 0285 0286 function stop = optimoutUpdateCorrectedData(params, ~, state) 0287 %OPTIMOUTUPDATEPLOTDATA Update corrected data sequences. 0288 switch state 0289 case 'init' 0290 case 'iter' 0291 if constant_flow 0292 data_cor1 = correctSensorLag(time1, data1, params); 0293 data_cor2 = correctSensorLag(time2, data2, params); 0294 else 0295 data_cor1 = correctSensorLag(time1, data1, flow1, params); 0296 data_cor2 = correctSensorLag(time2, data2, flow2, params); 0297 end 0298 case 'interrupt' 0299 case 'done' 0300 end 0301 stop = false; 0302 end 0303 0304 % Depth-value diagram. 0305 function stop = optimplotDepthData(~, ~, state) 0306 %OPTIMPLOTDEPTHDATA Depth-value diagram for a pair of casts. 0307 switch state 0308 case 'init' 0309 case 'iter' 0310 valid_cor1 = ~(isnan(data_cor1) | isnan(depth1)); 0311 valid_cor2 = ~(isnan(data_cor2) | isnan(depth2)); 0312 valid1 = ~(isnan(data1) | isnan(depth1)); 0313 valid2 = ~(isnan(data2) | isnan(depth2)); 0314 plot(data_cor1(valid_cor1), depth1(valid_cor1), '-r', ... 0315 data_cor2(valid_cor2), depth2(valid_cor2), '-b', ... 0316 data1(valid1), depth1(valid1), ':r', ... 0317 data2(valid2), depth2(valid2), ':b'); 0318 title('Depth-Data diagram'); 0319 xlabel('Data'); 0320 ylabel('Depth'); 0321 set(gca, 'YDir', 'reverse') 0322 case 'interrupt' 0323 case 'done' 0324 end 0325 stop = false; 0326 end 0327 0328 % Value-time diagram. 0329 function stop = optimplotDataTime(~, ~, state) 0330 %OPTIMPLOTDATATIME Value-time diagram for a pair of casts. 0331 switch state 0332 case 'init' 0333 case 'iter' 0334 valid_cor1 = (time1 > 0) & ~isnan(data_cor1); 0335 valid_cor2 = (time2 > 0) & ~isnan(data_cor2); 0336 valid1 = (time1 > 0) & ~isnan(data1); 0337 valid2 = (time2 > 0) & ~isnan(data2); 0338 time_offset = min(min(time1(valid1)), min(time2(valid2))); 0339 plot(time1(valid_cor1) - time_offset, data_cor1(valid_cor1), '-r', ... 0340 time2(valid_cor2) - time_offset, data_cor2(valid_cor2), '-b', ... 0341 time1(valid1) - time_offset, data1(valid1), ':r', ... 0342 time2(valid2) - time_offset, data2(valid2), ':b'); 0343 title('Data-Time diagram'); 0344 xlabel('Time'); 0345 ylabel('Data'); 0346 case 'interrupt' 0347 case 'done' 0348 end 0349 stop = false; 0350 end 0351 0352 end