11

I want to generate several random arrays {x1,x2,x3,x4,x5} that satisfy the following complex constraints:

(2.3 x1 + 4.1 x2 + 2.2 x3 + 3.5 x4 + 1.8 x5  <= 50) && (0 <= x1 <= 5 && 
   0 <= x2 && 1 <= x3 <= 5 && 1 <= x4 && 0 <= x5 ) && (2 <= x3 + x4 <= 5)

How I can do it.

In addition, I want to generate integer points in the region and improve the speed of generating these points

2 Answers2

16
constraints = (2.3 x1 + 4.1 x2 + 2.2 x3 + 3.5 x4 + 1.8 x5 <= 
     50) && (0 <= x1 <= 5 && 0 <= x2 && 1 <= x3 <= 5 && 1 <= x4 && 
     0 <= x5) && (2 <= x3 + x4 <= 5);

vars = {x1, x2, x3, x4, x5};

Use constraints and vars to define an ImplicitRegion:

ir = ImplicitRegion[constraints, Evaluate@vars];

Use RandomPoint[region, n] to draw n random points from region:

sample = RandomPoint[ir, 3]

{{0.569525,0.655619,1.19525,2.5036,15.6858},{0.666914,5.79318,1.71149,2.1047,7.46705},{3.56459,3.98713,3.0594,1.16909,7.3919}}

Verify that each point satisfies the constraints:

constraints /. Thread[vars -> #] & /@ sample

{True, True, True}

Update: "to generate integer points in the region"

Two approaches:

  1. Generate all integer tuples inside the cuboid that contains ir and select the ones that satisfy the constraints; then sample from the resulting list of tuples.

feasibleQ = Function[Evaluate@vars, Evaluate@constraints];

integertuples = Tuples @ MapThread[Range[Ceiling@#, Floor@#2] &, 
    Transpose @ RegionBounds[ir]];

integerfeasible = Select[Apply[feasibleQ]]@integertuples

Length@integerfeasible

5319

RandomChoice[integerfeasible, 3]

{{4, 1, 1, 1, 12}, {0, 4, 2, 3, 9}, {0, 0, 3, 1, 3}}

feasibleQ @@@ %

{True, True, True}

  1. Rationalize the constraints and use Solve:

solutions = vars /. Solve[Rationalize@constraints, vars, Integers];

Length@solutions

5319

feasibleQ = Function[Evaluate@vars, Evaluate@constraints];

And @@ (feasibleQ @@@ solutions)

True

RandomChoice[solutions, 3]

{{0, 5, 4, 1, 8}, {3, 1, 2, 1, 4}, {2, 0, 1, 2, 1}}

kglr
  • 394,356
  • 18
  • 477
  • 896
  • How to deal with the condition that these points are all integer points and how to improve the speed of generating these points. –  Feb 06 '20 at 08:04
  • 2
    @bmc013 you should update your question with this additional set of questions. In the future it might be better to start with this sort of information/requirement in your original post. While kglr is super cool and will likely update with these additional questions being answered, it is more helpful to future users if you update your question with these requirements. – CA Trevillian Feb 06 '20 at 08:25
  • 1
    Thank you @CATrevillian for the kind words:) – kglr Feb 06 '20 at 08:51
  • @bmc013, don't know how to add the integer constraint (adding the domain Element[vars, Integers] into the constraints does not work). – kglr Feb 06 '20 at 08:58
  • 2
    @CATrevillian Thank you for your advice, I have updated the question. –  Feb 06 '20 at 09:00
  • @kglr dude nice edit! I found a way to use FindInstance with your definition of ImplicitRegion and posted it as an answer, not sure what way is better, but your answer surely is (& more complete, too, what with a defined testQ n what-not) :D – CA Trevillian Feb 06 '20 at 09:30
4

Jumping off of @kglr's answer (pre-edit of the integer method), one can use FindInstance with RandomSeeding->Automatic to accomplish exactly the OP's updated task. I am not sure if you can get more efficient than this, as it takes negligible time on my PC for the production of a list of 10 Integer sets using RandomSeeding->Automatic. I suppose it would depend on your precise use case if this is efficient enough.

vars /. FindInstance[vars ∈ ir, vars, Integers, 10, RandomSeeding -> Automatic]

{{1, 4, 1, 1, 9}, {0, 3, 2, 1, 8}, {3, 0, 3, 1, 0}, {4, 1, 1, 3, 9}, {5, 0, 2, 3, 1}, {4, 1, 1, 1, 4}, {0, 0, 2, 3, 0}, {5, 3, 2, 3, 0}, {2, 1, 1, 2, 15}, {5, 0, 1, 2, 8}}

Here, I use 10 for the number of produced sets as not defining it, i.e., leaving it as 1, only produces the same solution each time.

vars /. FindInstance[vars ∈ ir, vars, Integers, RandomSeeding -> Automatic]

{{0, 0, 1, 1, 0}}

kglr's method is faster:

(vars /. FindInstance[vars ∈ ir, vars, Integers, 10, 
 RandomSeeding -> Automatic]) // feasibleQ @@@ # & // AbsoluteTiming

(integertuples = 
Tuples@MapThread[{Ceiling@#, Floor@#2} &, 
  Transpose@RegionBounds[ir]];
feasibleQ = Function[Evaluate@vars, Evaluate@constraints];
integerfeasible = Select[Apply[feasibleQ]]@integertuples;
RandomChoice[integerfeasible, 10]) // feasibleQ @@@ # & // AbsoluteTiming

{0.350308, {True, True, True, True, True, True, True, True, True, True}}
{0.0725979, {True, True, True, True, True, True, True, True, True, True}}

Both methods yield appropriate results, however.

With an updated definition of integertuples

integertuples = Tuples@MapThread[Range[Ceiling@#, Floor@#2] &, Transpose@RegionBounds[ir]];

Timing increases slightly

{0.215071, {True, True, True, True, True, True, True, True, True, True}}

However, if one is timing these appropriately by not computing integertuples again each time upon evaluation, one gets

RandomChoice[oldintegerfeasible, 10] // feasibleQ @@@ # & // AbsoluteTiming
RandomChoice[integerfeasible, 10] // feasibleQ @@@ # & // AbsoluteTiming

{0.0000818, {True, True, True, True, True, True, True, True, True, True}}
{0.0000852, {True, True, True, True, True, True, True, True, True, True}}

Illustrating comparable timings between versions of kglr's integertuples (definition shown above is oldintegertuples with oldintegerfeasible in comparison to kglr's correction using Range, which gives all solutions held within the cuboid.)

MikeY
  • 7,153
  • 18
  • 27
CA Trevillian
  • 3,342
  • 2
  • 8
  • 26
  • 1
    please check the corrected version of integertuples in my updated post. – kglr Feb 06 '20 at 09:51
  • @kglr if it is the edit using Range, I find comparable timings, the whole of the method takes a bit longer, but I realized I was timing inappropriately, and updated with that. Seems like your old method was slightly quicker? Maybe? What is the reasoning behind the corrected version? – CA Trevillian Feb 06 '20 at 10:11
  • 1
    the older (incorrect) version takes only the corners of the bounding cuboid (it contains only 2^5 = 32 tuples) it ignores the integer tuples inside the cuboid. – kglr Feb 06 '20 at 10:21
  • I think FindInstance doesn't guarantee lack of bias on instances it returns. I would be wary of it for this purpose if there are requirements to randomness of results. – kirma Feb 29 '20 at 10:55
  • 1
    @kirma agreed. When you only use 1 you always get the same output. – CA Trevillian Feb 29 '20 at 18:54
  • @CATrevillian That is one thing. Another is that even if you would supply different RandomSeeding values to it ad infinitum, nothing really guarantees that the distribution would approach uniform distribution over the region. RandomPoint is really designed for that purpose, while FindInstance... I'm not quite certain what one should say about it. It's definitely practical for many tasks, but maybe random sampling with it comes with caveats. – kirma Mar 01 '20 at 07:02