Feature enhancement : subsampling the answer space

Dear Rick, dear community. I come here with a proposal of feature / enhancement. Please note that there might already be a way (maybe with your lib or with "plain" prolog), that I may have missed, because I am fairly new to Prolog.

Here is a small example concerning the link between cartesian 2d coordinates, and its representation in circular coordinates.

%SWI-Prolog (threaded, 64 bits, version 9.1.10)
:- use_module(library(clpBNR)). %documentation at : https://ridgeworks.github.io/clpBNR/CLP_BNR_Guide/CLP_BNR_Guide.html

cartesian_circular(X, Y, R, Theta):-
    [X, Y, R]::real(-1.0e+16, 1.0e+16),
    [Theta]::real(-pi, pi),
    {
        X == R * cos(Theta),
        Y == R * sin(Theta)
    }.

test1_(R, Theta):-
    X::real(-50.0, -5.0),
    Y::real(1.0, 20.0),
    cartesian_circular(X, Y, R, Theta), 
    solve([R, Theta]).

    %answers : 
    %R::real(-53.851650867757314, -5.099019332794144),
    %Theta::real(-1.3258177199522887, -0.019997279208893784) ;
    %R::real(5.099019290218574, 53.85164889904919),
    %Theta::real(1.8157749415368631, 3.121595456378679) ;

This is nice, but, how to I get a list of valid answers ? My wish / idea would be ask for "automated" subsampling of the answer space (subsambling of the "inputs", "outputs" and both ).

Bellow, some examples with comments :


expected_feature_subsampling_1_1(R, Theta):-
    %Here, the expected feature would be to subssample the answer set based on a subsampling of the "INPUTS".
    %Practically, in this example, this would be "get the R,T coords of a grid of points in a rectangle between xmin, xmax, ymin, ymax"

    %idea 1 : using an absolute subsampling step
    X::real(-50.0, -5.0),
    Y::real(1.0, 20.0),
    cartesian_circular(X, Y, R, Theta), 
    %subsample_absolute_step(X, 5.0), % this would be like [-50.0, -45.0, -40.0, .. -10.0, -5.0]
    %subsample_absolute_step(Y, 2.0), % this would be like [1.0, 3.0, 5.0, .. 15.0, 17.0, 19.0, 20.0] <- yes, I think this is good to keep the extrem values, even if this does not match the abs step.
    solve([R, Theta]).

expected_feature_subsampling_1_2(R, Theta):-
    %idea 2 : using a relative subsampling step (a ratio)
    X::real(-50.0, -5.0),
    Y::real(1.0, 20.0),
    cartesian_circular(X, Y, R, Theta), 
    %subsample_relative_step(X, 0.1), %i.e. 10% relative : (-5.0 - (-50.0)) * 0.1 = 45.0 * 0.1 = 4.5 of absolute step (see subsample_absolute_step)
    %subsample_relative_step(Y, 0.1), %i.e. 10% relative : (20.0 - 1.0) * 0.1 = 19.0  * 0.1 = 1.9 of absolute step (see subsample_absolute_step)
    solve([R, Theta]).

expected_feature_subsampling_1_3(R, Theta):-
    %idea 3 : give the number of samples wanted
    X::real(-50.0, -5.0),
    Y::real(1.0, 20.0),
    cartesian_circular(X, Y, R, Theta), 
    %subsample_nb_elems(X, 15), %15 subsampled elements, and we keep the extremums (min and max of interval)
    %subsample_nb_elems(Y, 10), %10 subsampled elements, and we keep the extremums (min and max of interval)
    solve([R, Theta]).

expected_feature_subsampling_2_1(R, Theta):-
    %Here, the expected feature would be to subssample the answer set based on a subsampling of the "OUTPUTS".
    %Practically, in this example, this would be "get the R,T coords in a rectangle between xmin, xmax, ymin, ymax, and subsample the output space R,T"

    %idea 1 : using an absolute subsampling step
    X::real(-50.0, -5.0),
    Y::real(1.0, 20.0),
    cartesian_circular(X, Y, R, Theta), 
    %subsample_absolute_step(R, 5.0), % same idea as examples expected_feature_subsampling_1 (but on "outputs")
    %subsample_absolute_step(Theta, pi/10.0), % same idea as examples expected_feature_subsampling_1 (but on "outputs")
    solve([R, Theta]).

expected_feature_subsampling_2_2(R, Theta):-
    %idea 2 : using a relative subsampling step (a ratio)
    X::real(-50.0, -5.0),
    Y::real(1.0, 20.0),
    cartesian_circular(X, Y, R, Theta), 
    %subsample_relative_step(X, 0.1), % same idea as examples expected_feature_subsampling_1 (but on "outputs")
    %subsample_relative_step(Theta, 0.1), % same idea as examples expected_feature_subsampling_1 (but on "outputs")
    solve([R, Theta]).

expected_feature_subsampling_2_3(R, Theta):-
    %idea 3 : give the number of samples wanted
    X::real(-50.0, -5.0),
    Y::real(1.0, 20.0),
    cartesian_circular(X, Y, R, Theta), 
    %subsample_nb_elems(R, 15), % same idea as examples expected_feature_subsampling_1 (but on "outputs")
    %subsample_nb_elems(Theta, 10), % same idea as examples expected_feature_subsampling_1 (but on "outputs")
    solve([R, Theta]).

expected_feature_subsampling_3_1(R, Theta):-
    %Here, the expected feature would be to subssample the answer set based on a subsampling of variables, whatever if they are "inputs" or "ouputs"
    %Here, this would be very nice that the lib "tries its best" when the subsampling is no longer feasible, because previous subsampling already reduced the space.

    X::real(-50.0, -5.0),
    Y::real(1.0, 20.0),
    cartesian_circular(X, Y, R, Theta),
    %subsample_nb_elems(X, 15), %15 subsampled elements, and we keep the extremums (min and max of interval)
    %subsample_nb_elems(Theta, 10), %10 subsampled elements, and we keep the extremums (min and max of interval)
    solve([R, Theta]).

expected_feature_subsampling_3_2(R, Theta):-
    %Here, this would be very nice that the lib "tries its best" when the subsampling is no longer feasible, because previous subsampling already reduced the space.

    X::real(-50.0, -5.0),
    Y::real(1.0, 20.0),
    cartesian_circular(X, Y, R, Theta),
    %subsample_nb_elems(X, 15), %15 subsampled elements, and we keep the extremums (min and max of interval)
    %subsample_absolute_step(Theta, pi/10.0), % Here, this would be very nice that the lib tries its best
    solve([R, Theta]).

expected_feature_subsampling_3_3(R, Theta):-
    %Here, this would be very nice that the lib "tries its best" when the subsampling is no longer feasible, because previous subsampling already reduced the space.

    X::real(-50.0, -5.0),
    Y::real(1.0, 20.0),
    cartesian_circular(X, Y, R, Theta),
    %subsample_nb_elems(X, 15), %15 subsampled elements, and we keep the extremums (min and max of interval)
    %subsample_absolute_step(Theta, pi/10.0), % Here, this would be very nice that the lib tries its best
    %subsample_absolute_step(Y, 2.0), %here, I am not even sure this is feasible / compatible with the other subsamplings. But this would be nice if the lib tries to do it (if feasible) or simply don't subsample any more.

    solve([R, Theta]).

Thanks in advance.

In general terms constraints aren't answers, they just constrain answers. So I think you pretty much just have to write the sampling predicates shown in your examples. Here's one for subsample_absolute_step:

subsample_absolute_step(X,X,_U,_Step).    % X=L
subsample_absolute_step(X,L,U,Step) :-    % or bump L and retry
    NxtL is L+Step,
    NxtL =< U,                            % fail if NxtL > U
    subsample_absolute_step(X,NxtL,U,Step).

Note this this works even if X isn't a constrained variable. Also note that if the unification in the head of the first clause fails due to a constraint violation, no answer is produced and the alternative clause (bumping L) is tried. A test predicate:

test2_(X, Y, R, Theta) :-
    cartesian_circular(X, Y, R, Theta),
    subsample_absolute_step(X,-50,-5,5),
    subsample_absolute_step(Y,1,20,2),
    solve([R,Theta]).

I also think R should be constrained to be positive, so:

cartesian_circular(X, Y, R, Theta):-
    [X, Y]::real,            % finite X,Y
    R::real(0, _),           % positive, finite R
    [Theta]::real(-pi, pi),
    {
        X == R * cos(Theta),
        Y == R * sin(Theta)
    }.

Now:

?- test2_(X,Y,R,Theta).
X = -50,
Y = 1,
R:: 50.0099990001999...,
Theta:: 3.12159531961664... ;
X = -50,
Y = 3,
R:: 50.0899191454727...,
Theta:: 3.08166449846858... ;
X = -50,
Y = 5,
R:: 50.2493781056044...,
Theta:: 3.04192400109863... ;
X = -50,
Y = 7,
R:: 50.4876222454573...,
Theta:: 3.00249671210772... .

% manually terminated

I didn't look closely at the other sampling strategies, but they can probably be implemented in a similar fashion.

The general idea of CLP's is to apply constraints and then use normal Prolog unification to generate answers subject to those constraints. A variation would be to add further constraints in a systematic way to solve sub-problems; that's essentially what solve does.

Does this address your question?

Dear Rick, thank you for your answer. Thanks for the clarification "and then use normal Prolog unification to generate answers". I was testing your proposal but then I realized something weird that may require help first, before going back to the subject of subsampling.

You are correct concerning the range of R : it should be positive (see polar coordinates on wikipedia).

But suprisingly, I got weird output range for R when using your implementation of cartesian_circular.

I now have :

?- test1_(R, Theta).
R::real(5.000999886333218, 206.15533777544633),
Theta::real(1.815774921412059, 3.121595456457498) ;
false.

The range for R looks unrealistic to me (R max is too high).

Whereas before, with a constraint of R::real(-1.0e+16, 1.0e+16) , I had the following :

?- test1_(R, Theta).
R::real(-53.851650867757314, -5.099019332794144),
Theta::real(-1.3258177199522887, -0.019997279208893784) ;
R::real(5.099019290218574, 53.85164889904919),
Theta::real(1.8157749415368631, 3.121595456378679) ;
false.

And here the range for R looks OK (eventhough I can have a negative R, which is like "back drawing" the radius ;) )

Do you have the same ? Can you explain me why ?

But surpisingly, I got weird output range for R when using your implementation of cartesian_circular.

Good catch and a little subtle. First a recommended solution:

Without constraining R to be positive, you get two solutions as you observed:

% R::real
?- test1_(R,Theta).
R::real(-53.851650867757314, -5.099019332794144),
Theta::real(-1.3258177199522887, -0.019997279208893784) ;
R::real(5.099019290218574, 53.85164889904919),
Theta::real(1.8157749415368631, 3.121595456378679) ;
false.

which I think is undesirable, but

% R::real(0,_)
?- test1_(R,Theta).
R::real(5.000999886333218, 206.15533777544633),
Theta::real(1.815774921412059, 3.121595456457498) ;
false.

The reason that R doesn't narrow sufficiently is that the current constraints in cartesian_circular don't enable much interval narrowing to be done, e.g.,

?- Theta::real(-pi,pi), R::real(0,_), {X == R*cos(Theta)}, X::real(-50,-5).
Theta::real(-3.1415926535897936, 3.1415926535897936),
R::real(5, 1.0e+16),
X::real(-50, -5).

?- Theta::real(-pi,pi), R::real(0,_), {Y == R*sin(Theta)}, Y::real(1,20).
Theta::real(-3.1415926535897936, 3.1415926535897936),
R::real(1, 1.0e+16),
Y::real(1, 20).

solve manages to forcibly narrow the intervals; it does a good job on Theta for falls short for R (more on this below).

A good way to address this problem is to add the "redundant" constraint { X**2 + Y**2 == R**2 }. It's redundant because it can be mathematically inferred from the other constraints (since sin(A)**2 + cos(A)**2 ==1, but the clpBNR constraint system isn't clever enough to sort that out. Conventional wisdom might say you should avoid redundant constraints (unnecessary overhead and an added source of rounding error) but here it helps; given the bounds of X and Y, the bounds of R can be calculated directly. This is more efficient (less work for solve); in this particular case, means you only have to solve for Theta.

Now:

?- test1_(R,Theta).
R::real(5.099019513592782, 53.85164807134507),
Theta::real(1.8157748914243488, 3.1215953577549764) ;
false.

with:

cartesian_circular(X, Y, R, Theta):-
    [X, Y]::real,            % finite X,Y
    R::real(0,_),            % positive, finite R
    [Theta]::real(-pi, pi),
    {
        X == R * cos(Theta),
        Y == R * sin(Theta),
        X**2 + Y**2 == R**2
    }.

I won't go into details, but the existing version of cartesian_circular works better if you use solve([Theta, R]), i.e., reverse the order of the R and Theta in the list. Although I don't think that's the right solution of this particular example (use redundant constraint instead), it's not very user friendly so I'm going to see if that can be improved.

Thank you Rick, the last problem is corrected adding the additionnal constraint. Now lets come back to the subsampling.

I have taken your subsample_absolute_step proposal, and modified it a bit to include the Upper bound (I am sad that it took me 30 minutes to find the proper way to get this... Prolog is far from other languages).

Now I am using it to sample the answer space on R. And the weird thing is that I don't get the full range I was expecting.

Bellow, the full code :

%Example for the issue : 
%https://github.com/ridgeworks/clpBNR/issues/19

%SWI-Prolog (threaded, 64 bits, version 9.1.10)
:- use_module(library(clpBNR)). %documentation at : https://ridgeworks.github.io/clpBNR/CLP_BNR_Guide/CLP_BNR_Guide.html

cartesian_circular(X, Y, R, Theta):-
    [X, Y]::real,            % finite X,Y
    R::real(0, _),           % positive, finite R
    [Theta]::real(-pi, pi),
    {
        X == R * cos(Theta),
        Y == R * sin(Theta),
        X**2 + Y**2 == R**2 %https://github.com/ridgeworks/clpBNR/issues/19#issuecomment-1605572520
    }.

subsample_absolute_step(X,X,_U,_Step).    % X=L

subsample_absolute_step(X,L,U,Step) :-    % or bump L and retry
    NxtL is L+Step,
    NxtL < U -> subsample_absolute_step(X,NxtL,U,Step) ; X is U.

test2_2_is_it_a_good_idea(X, Y, R, Theta, [StartR, EndR], [StartT, EndT]) :-
    X::real(-50.0, -5.0),
    Y::real(1.0, 20.0),    
    cartesian_circular(X, Y, R, Theta),
    range(R, [StartR, EndR]),
    range(Theta, [StartT, EndT]),
    subsample_absolute_step(R, StartR, EndR, 5.0).

Now testing query :

?-test2_2_is_it_a_good_idea(X, Y, R, Theta, [StartR, EndR], [StartT, EndT]).

Got result :

R = 10.099019513592783,
Theta = 2.1991148575128547,
StartR = 5.099019513592783,
EndR = 53.85164807134507,
StartT = -3.1415926535897936,
EndT = 3.1415926535897936,
X:: -5.93605473270373...,
Y:: 8.17027841302078... .

Pressing Tab, I never see Rbellow 10.099 eventhough StartR = 5.099019513592783

I decide to concat the results with query :

?-findall(R, test2_2_is_it_a_good_idea(X, Y, R, Theta, [StartR, EndR], [StartT, EndT]), L).

Answer is :

L = [10.099019513592783, 15.099019513592783, 20.099019513592783, 25.099019513592783, 30.099019513592783, 35.09901951359278, 40.09901951359278, 45.09901951359278, 50.09901951359278].

So no R values starting at StartR.

Am I missing something ?

Thanks in advance

Profuse apologies - the email from GitHub ended up in my Junk and I just discovered it. I'm a bit busy today but will get back to you by the end of the week (if this is still and issue for you).

Short answer: interval arithmetic is conservative with respect to floating point rounding errors so there may not be a solution at the exact lower bound value of R (StartR). If you start with a slightly higher value of StartR (or a smaller step size), you may get closer to the answers you're looking for, e.g., using a step size of 0.1 rather than 5.0:

?- test2_2_is_it_a_good_idea(X, Y, R, Theta, [StartR, EndR], [StartT, EndT]).
R = 5.1990195135927815,
StartR = 5.099019513592782,
EndR = 53.85164807134506,
StartT = -3.1415926535897936,
EndT = 3.1415926535897936,
X::real(-5.101941189656986, -4.999999999999999),
Y::real(0.9999999999999999, 1.4247118665605811),
Theta::real(-3.1415926535897936, 3.1415926535897936) ;
R = 5.299019513592781,
StartR = 5.099019513592782,
EndR = 53.85164807134506,
StartT = -3.1415926535897936,
EndT = 3.1415926535897936,
X::real(-5.203807049212825, -4.999999999999999),
Y::real(0.9999999999999999, 1.7548811371249877),
Theta::real(-3.1415926535897936, 3.1415926535897936) 
%%% manually terminated

To get a better idea of what's happening, consider the following query which calculates R two ways; using cartesion_circular/4 to calculate R0 and by direct computation (sqrt(X**2+Y**2)) to calculate R1:

?- X::real(-50.0, -5.0), Y::real(1.0, 20.0), cartesian_circular(X,Y,R0,Theta),{R1 is sqrt(X**2+Y**2)}.
X::real(-50.00000000000001, -4.999999999999999),
Y::real(0.9999999999999999, 20.000000000000004),
R0::real(5.099019513592782, 53.85164807134506),
Theta::real(-3.1415926535897936, 3.1415926535897936),
R1::real(5.099019513592783, 53.85164807134505).

Note that the lower bound of R0 is slightly out of range (by 1 floating point distance?) of R1 (similarly for the upper bound). So unifying R1 with the exact bound(s) will fail, although it will probably succeed for almost any other value in R0's range. (The root cause of this is probably the use of incorrectly rounded libc floating point arithmetic requiring extra outward rounding on trig functions but I need to dig into this further.)

All of which begs the question of what the intervals bounds mean? An interval's width starts out being fairly wide (possibly infinite) and then gets narrowed by the application of constraints. Such narrowing implies there is no solution consistent with the constraints outside the interval. But there are no guarantees about what's inside the interval. There may be 0, 1, or many (infinite?) solutions contained between (and including) the bounds. The same applies to bounds that are only 1 floating point distance apart (the domain of reals is continuous).

This would seem to be a bit worrisome but, for all practical purposes, problem knowledge can be used to anticipate the meaning. In this particular example, we know the functions are continuous in the range of initial range of values, so there nominally an infinite set of (real) solutions. But the vagueries of floating point arithmetic and the supporting libraries sometimes prevent the bounds from being as tight as possible. (Using the actual boundary values in interval unification, as in this example, requires particular care.)

So my next question is: Do you see a fundamental requirement for "sampling" in your application in your application, or is it more of an interesting experiment to come to grips with how to use CLP over intervals? Generally, it's preferable to maintain intervals in their current state while applying additional constraints, e.g., minimize travel distance, before converting them to a final point value for "output".

 ... Prolog is far from other languages).

Certainly true and Prolog literacy will be needed to progress. Thinking in terms of an imperative language implementation will often be counter-productive. Fortunately there are lots of resources (textbooks, consulting expertise, mailing lists) to help on that front.

And in case you haven't seen this already, here's another paper on using CLP over intervals for robotics type applications. https://arxiv.org/pdf/cs/0007002.pdf The authors list might be another source of expertise if you're looking for outside help.

Short answer: interval arithmetic is conservative with respect to floating point rounding errors so there may not be a solution at the exact lower bound value of R (StartR). If you start with a slightly higher value of StartR (or a smaller step size), you may get closer to the answers you're looking for, e.g., using a step size of 0.1 rather than 5.0:

Okay, indeed, starting a bit away from StartR returns the expected result.

test2_2_is_it_a_good_idea(X, Y, R, Theta, [StartR, EndR], [StartT, EndT]) :-
    X::real(-50.0, -5.0),
    Y::real(1.0, 20.0),    
    cartesian_circular(X, Y, R, Theta),
    range(R, [StartR, EndR]),
    range(Theta, [StartT, EndT]),
    S2 is StartR + 0.1, %Interval arithmetic is conservative with respect to floating point rounding errors so there may not be a solution at the exact lower bound value . See : https://github.com/ridgeworks/clpBNR/issues/19#issuecomment-1642201945
    subsample_absolute_step(R, S2, EndR, 5.0).

So my next question is: Do you see a fundamental requirement for "sampling" in your application in your application, or is it more of an interesting experiment to come to grips with how to use CLP over intervals?

I tried to solve some 6 degree of freedom (3 translations, 3 rotations) problems (such as the "change of basis" problem in the SWI prolog discourse), where there are several compositions of transforms. This leads to a probably big search space. Your lib may return intervals in a nice iteractive time. But when I called the solve , this was too long for my patience. I my case, I don't need optimization. Or sometime, I cannot even think of an objective function that ensure me to set all variables to a minimum (or maximum). So I though : how about sampling within the intervals of some of my variables (the ones I can practically subsample, such as the Translations X and Y), and then solve the other variables ? But maybe this was not a good idea.

But more wisely, I concluded that : this is not because you can write it (without simplification hypothesis) with CLP(R) and/or CLP(BNR) that you should / that it will solve your problems in a reasonable execution time. On my former problem, we found an alternative way to formulate, with a lot of hypothesis (and drawbacks), but now we can solve it with linear programming, and this is fast as light.

(if this is still and issue for you)

Well I need to train more. Currently I still feel sad that I can't fully achieve my goals. But I am no longer in a hurry. So let's take it gently with pleasure when I find some time. And this is not a matter if I cannot succeed.

I tried to solve some 6 degree of freedom (3 translations, 3 rotations) problems (such as the "change of basis" problem in the SWI prolog discourse), where there are several compositions of transforms.

This sounds familiar to Benhamou et al.'s paper, but I could be wrong.

But when I called the solve , this was too long for my patience. I my case, I don't need optimization.

It may not help, but solve has an optional precision parameter which roughly corresponds to the number of digits of precision in the solution. Default is 6 but sounds like 2 or 3 may be sufficient for your purposes and should be faster.

Also solve tries to find bifurcating (split) points which aren't solutions, which is useful for separating solutions, e.g., roots of a polynomial. But, if your example is typical, this will be difficult since pretty much all such split points are solutions. This can result in much fruitless searching with the conclusion that no splitting is possible. splitsolve (with precision option) may be a better choice for "labelling" in these kinds of situations.

but now we can solve it with linear programming, and this is fast as light.

Yes, if the problem can be reduced to a linear model there are lots of efficient, well known techniques for tackling it.

So let's take it gently with pleasure when I find some time. And this is not a matter if I cannot succeed.

That's probably wise - it's hard to learn under pressure.

I find the topic intriguing so if you can find the time to compose a small relevant sub-problem (6 degrees of freedom?), I would be interested in looking at it. Or if you have any experiments that failed to meet requirements/expectations that you could share, I'd be happy to look at those as well. My motivation is to determine whether there are fundamental issues with how things work, or if the documentation can be improved to avoid common pitfalls, or ??.

ridgeworks / clpBNR

Feature enhancement : subsampling the answer space #19