pkg load optim

% Initialize some other values
BBR_S = [ ...
    120.272     93.782     69.465     49.861     34.546     22.970; ...
    276.088    230.073    185.954    145.986    111.780     83.573; ...
    550.313    483.285    414.183    347.850    286.163    228.998; ...
    720.538    671.688    605.032    547.171    480.848    415.814; ...
   1006.866    947.639    886.043    822.754    739.724    667.029];

BBR_D = 1./BBR_S;

BBR_tr = [ ...
434631.32360 814933.73175 1629867.46349 3259734.92699 6519469.85398 13038939.70796; ...
 29101.54173  54565.39075  109130.78150  218261.56301  436523.12602   873046.25204; ...
  1565.26858   2934.87860    5869.75719   11739.51438   23479.02877    46958.05753; ...
   149.31155    279.95915     559.91830    1119.83659    2239.67319     4479.34638; ...
     8.00000     15.00000      30.00000      60.00000     120.00000      240.00000];

% initialise the GPL parameters
GPL_x = [0; 0; 0];  % in order D0 D1 m

% the first estimate of model parameters
% in the GPL model D(t) = D0 + D1*(t)^m
GPL_D0 = 0.0006; %0.000626249601623718;
GPL_D1 = 0.0001; %0.000107587010307332;
GPL_m  = 0.3       %0.312356329396664;

% D_test = fn_GPL_model (GPL_D0, GPL_D1, GPL_m, BBR_tr)

% this should not be necessary, but for simplicity let us put data in one linear array
GPL_D  = BBR_D(1,:);
GPL_tr = BBR_tr(1,:);
for ii = 2 : rows(BBR_D)
  GPL_D  = [GPL_D, BBR_D(ii,:)];
  GPL_tr = [GPL_tr, BBR_tr(ii,:)];
end

fun_GPL = @(GPL_x)fn_BBR_GPL_RMSE (GPL_x, GPL_D, GPL_tr);
%GPL_x0_lb = [1000.0; 1.0e6; 0.1; 0.0];
%GPL_x0_ub = [3000.0; 4.0e6; 0.9; 1.0];
GPL_x0 = [GPL_D0; GPL_D1; GPL_m];
[GPL_x, GPL_RMSE, cvg, outp] =...
   nonlin_min(fun_GPL, GPL_x0);

disp(" ");
disp ("Final values")
GPL_D0 = GPL_x(1)
GPL_D1 = GPL_x(2)
GPL_m  = GPL_x(3)