pkg load interval

## Datasets (observations or known error bounds)
global X
X = vertcat (infsup (0.1, 0.2), ...
             infsup (0.2, 0.3), ...
             infsup (0.8, 0.9));
global Y
Y = vertcat (infsup (-0.2, 1.5), ...
             infsup (-0.3, 1.9), ...
             infsup (1.2, 2.5));

## This function is to be inverted
function [fval, a1, a2, a3, b1, b2] = F (y, a1, a2, a3, b1, b2)
    global X
    global Y
    # Vectorized forward evaluation
    A = [a1; a2; a3];
    lnX = log (X);
    A_lnX = A .* lnX;
    b2_X = b2 .* X;
    b1_b2_X = b1 + b2_X;
    fval = A_lnX .* b1_b2_X;
    
    # Backward evaluation and refinement of the parameters
    y = intersect (y, fval);
    # Contract the factors of the multiplication
    A_lnX = mulrev (b1_b2_X, y, A_lnX);
    b1_b2_X = mulrev (A_lnX, y, b1_b2_X);
    A = mulrev (lnX, A_lnX, A);
    a1 = A(1, :); a2 = A(2, :); a3 = A(3, :);
    # Contract the sum
    b2_X = intersect (b2_X, b1_b2_X - b1);
    b1 = intersect (b1, union (b1_b2_X - b2_X, [], 1));
    # Contract b2
    b2 = union (mulrev (X, b2_X, b2), [], 1);
endfunction


## Initial search ranges
A = vertcat (infsup (3.1, 3.2), ...
             infsup (3.3, 3.4), ...
             infsup (3.5, 3.6));
B = vertcat (infsup (-10, 10), ...
             infsup (-10, 10));

[parameters, paving] = fsolve(@F, [A; B], Y, struct ("Contract", true));
a1 = parameters(1) # [3.1, 3.2001]
a2 = parameters(2) # [3.2999, 3.4]
a3 = parameters(3) # [3.5, 3.6001]
b1 = parameters(4) # [0.55553, 1.8562]
b2 = parameters(5) # [-9.3644, -4.5616]

plot (paving(4, :), paving(5, :))
xlabel ("b1")
ylabel ("b2")