Find all partitions of a plane created by a set of lines

Finding the number of partitions follows this formula (when no 3 or more lines intersect at the same point - this assumption is carried throughout this post):

num_partitions = (n_lines * (n_lines + 1) / 2 ) + 1

An explanation can be found here, and another one there.

The desired output would be a complete list of such sets: {(x,y)|y <= -30x+28 && 5x+3 <= y <= 60x+2}, {(x,y)| -30x+28<= y && 60x+2 <= y <= 90x+7}... etc... They don't have to be given in this notation, of course.

The lack of a precise notation is a handicap here.

Below is my back of the envelop attempt. As you see, it is possible to identify the numbered areas based on their relative positioning to each line.
There are 5 empty sets, which won't be the same if the lines ordering changes.
It would probably be easier to partition a set of points on the plane with lines, trying to determine which point belongs to which set; in such a case, exploring the 2^n potential partitions, and returning their content would be easy. (easier than trying to find a good notation for identifying an abstract set)

This doesn't fully answer your question, but may be a good strating point for someone able/willing to push this further.
Here is a note on partitioning a set of points with two lines in the plane. it is a different problem, but some of its approach could be useful.

other approaches:

Identify the polygons formed by the line segments, calculate a convex hull, determine if a point is in that hull.

Here is a solution in Mathematica. The method involves finding the line intersection points, the line fragments and then the partitions, while keeping track of which lines the points are connected to.

y[1] = 5 x + 3;
y[2] = 90 x + 7;
y[3] = -30 x + 28;
y[4] = 60 x + 2;

findpoints[n_] := Module[{},
  xp = DeleteCases[{{First[#1], First[#2]},
       Solve[Last[#1] == Last[#2], x]} & @@@
     DeleteCases[
      Tuples[Array[{#, y[#]} &, n], 2],
      {{x_, _}, {x_, _}}], {_, {}}];
  yp = y[#[[1, 1]]] /. #[[2, 1]] & /@ xp;
  MapThread[{#1, {#2, #3}} &,
   {xp[[All, 1]], xp[[All, 2, 1, 1, 2]], yp}]]

xyp = findpoints[4];

{xmin, xmax} = Through[{Min, Max}@
    xyp[[All, 2, 1]]] + {-0.7, 0.7};

outers = Flatten[Array[Function[n,
     MapThread[List[{{n, 0}, {##}}] &,
      {{xmin, xmax}, y[n] /.
    List /@ Thread[x -> {xmin, xmax}]}]], 4], 2];

xyp = Join[outers, xyp];

findlines[p_] := Module[{},
  pt = DeleteCases[
    Cases[xyp, {{p[[1, 1]], _}, _}], p];
  n = First@Cases[pt,
     {_, First@Nearest[Last /@ pt, Last[p]]}];
  {{First[p], First[n]}, {Last[p], Last[n]}}]

lines = Map[findlines, xyp];

(* boundary lines *)
{ymin, ymax} = Through[{Min, Max}@outers[[All, 2, 2]]];
{lbtm, rbtm, ltop, rtop} = {{xmin, ymin},
   {xmax, ymin}, {xmin, ymax}, {xmax, ymax}};
xminlines = Partition[Union@Join[{ymin, ymax},
      Cases[xyp, {_, {xmin, _}}][[All, 2, 2]]], 2, 1] /.
   x_Real :> {xmin, x};
xmaxlines = Partition[Union@Join[{ymin, ymax},
      Cases[xyp, {_, {xmax, _}}][[All, 2, 2]]], 2, 1] /. 
   x_Real :> {xmax, x};
lines2 = Join[Last /@ lines, xminlines, xmaxlines,
   {{lbtm, rbtm}}, {{ltop, rtop}}];

ListLinePlot[lines2]

(* add vertex points *)
xyp2 = Join[xyp, {
   {{SortBy[Cases[outers, {_, {xmin, _}}],
       Last][[-1, 1, 1]], -1}, ltop},
   {{SortBy[Cases[outers, {_, {xmax, _}}],
       Last][[-1, 1, 1]], -1}, rtop},
   {{SortBy[Cases[outers, {_, {xmin, _}}],
       Last][[1, 1, 1]], -1}, lbtm},
   {{SortBy[Cases[outers, {_, {xmax, _}}],
       Last][[1, 1, 1]], -1}, rbtm}}];

anglecalc[u_, v_] := Mod[(ArcTan @@ u) - (ArcTan @@ v), 2 π]

getlineangles[] := Module[{},
  (* find the angles from current line
     to all the linked lines *)
  angle = Map[
    anglecalc[{c, d} - {g, h}, # - {g, h}] &,
    union = DeleteCases[Union@Join[
        Last /@ Cases[lines2, {{g, h}, _}],
        First /@ Cases[lines2, {_, {g, h}}]],
      {c, d}]];
  Sort[Transpose[{N@angle, union}]]]

getpolygon[pt_, dir_] := Module[{},
  Clear[p];
  p[n = 1] = {{a, b}, {c, d}} = pt;
  (* find the angles from vector (0, -1) or (0, 1)
     to all the linked lines *)

  angle = Map[anglecalc[If[dir == 1, {0, -1}, {0, 1}], # - {c, d}] &,
    union = Union@Join[
       Last /@ Cases[lines2, {{c, d}, _}],
       First /@ Cases[lines2, {_, {c, d}}]]];
  lineangles = Sort[Transpose[{N@angle, union}]];
  (* next point *)
  p[++n] = {{e, f}, {g, h}} = First@
     Cases[xyp2, {_, lineangles[[1, 2]]}];

  While[Last[p[n]] != Last[p[1]],
   lineangles = getlineangles[];
   (* reset first point *)
   {{a, b}, {c, d}} = {{e, f}, {g, h}};
   (* next point *)
   p[++n] = {{e, f}, {g, h}} = First@
      Cases[xyp2, {_, lineangles[[1, 2]]}]];

  Array[p, n]]

len = Length[xyp];

polygons = Join[Array[(poly[#] = getpolygon[xyp[[#]], 1]) &, len],
   Array[(poly[# + len] = getpolygon[xyp[[#]], 2]) &, len]];

graphics = DeleteDuplicates /@ Array[Last /@ poly[#] &, 2 len];

sortedgraphics = Sort /@ graphics;

positions = Map[Position[sortedgraphics, #] &,
    DeleteDuplicates[sortedgraphics]][[All, 1, 1]];

unique = poly /@ positions;

poly2 = unique[[All, All, 2]];

poly2 = Delete[poly2,
   Array[If[Length[Intersection[poly2[[#]],
         Last /@ Take[xyp2, -4]]] == 4, {#}, Nothing] &,
    Length[poly2]]];

len2 = Length[poly2];

poly3 = Polygon /@ Rest /@ poly2;

Array[(centroid[#] = RegionCentroid[poly3[[#]]]) &, len2];

Show[Graphics[Array[{ColorData[24][#],
     poly3[[#]]} &, len2], AspectRatio -> 1/GoldenRatio],
 Graphics[Array[Text[#, centroid[#]] &, len2]]]

unique2 = Extract[unique,
   Flatten[Position[unique[[All, All, 2]], #] & /@ poly2, 1]];

makerelations[oneconnection_, areanumber_] := Module[{i},
  If[Intersection @@ oneconnection == {} ||
    (i = First[Intersection @@ oneconnection]) < 1,
   Nothing,
   centroidx = First[centroid[areanumber]];
   linepos = y[i] /. x -> centroidx;
   relation = If[linepos < Last[centroid[areanumber]],
     " >= ", " < "];
   string = StringJoin["y", relation, ToString[y[i]]]]]

findrelations[n_] := Module[{},
  areanumber = n;
  onearea = unique2[[areanumber]];
  connections = Partition[First /@ onearea, 2, 1];
  strings = DeleteDuplicates@
    Map[makerelations[#, areanumber] &, connections];
  StringJoin["Area ", ToString[areanumber],
   If[areanumber > 9, ": ", ":  "],
   StringRiffle[strings, " &&\n         "]]]

Show[Plot[Evaluate@Array[y, 4], {x, -1, 1.5},
  PlotLegends -> "Expressions", Axes -> None],
 Graphics[Array[Text[#, centroid[#]] &, len2]]]

Column@Array[findrelations, len2]

Area 1:  y >= 28 - 30 x &&
         y < 3 + 5 x
Area 2:  y >= 2 + 60 x &&
         y >= 28 - 30 x &&
         y < 7 + 90 x
Area 3:  y < 28 - 30 x &&
         y < 7 + 90 x &&
         y < 2 + 60 x &&
         y < 3 + 5 x
Area 4:  y >= 3 + 5 x &&
         y >= 28 - 30 x &&
         y < 2 + 60 x
Area 5:  y < 3 + 5 x &&
         y >= 2 + 60 x &&
         y < 7 + 90 x
Area 6:  y < 28 - 30 x &&
         y >= 2 + 60 x &&
         y >= 3 + 5 x &&
         y < 7 + 90 x
Area 7:  y < 28 - 30 x &&
         y >= 3 + 5 x &&
         y < 2 + 60 x
Area 8:  y < 28 - 30 x &&
         y >= 7 + 90 x &&
         y >= 3 + 5 x
Area 9:  y < 2 + 60 x &&
         y >= 7 + 90 x
Area 10: y >= 28 - 30 x &&
         y >= 7 + 90 x
Area 11: y < 3 + 5 x &&
         y >= 7 + 90 x &&
         y >= 2 + 60 x

Full Matlab solution. 120 lines 7021 partitions, in 0.4 seconds raw calculation (+2 seconds for drawing).

Final.

clear
tic
%rng(15) % fix rng for debug
NLines=9;
DRAW_FINAL_POLY= true==true ; %false true
DRAW_FINALTEXT= true==true ; %false true
DRAW_LINES= false==true; %false true
DRAW_DEBUG= false==true; %false true

% part a - generate lines
NORM_EPS=1e-10;
Lines=10*(rand(NLines,4)-0.5); %(x1,y1,x2,y2)
Color=rand(NLines,3);
Lines(1,:)=[-10,0,+10,0];% x axis if we want to add asix as lines
Lines(2,:)=[0,-10,0,+10];% y axis
Color(1,:)=[1,0,0];% x axis red
Color(2,:)=[0,1,0];% y axis green
Color(3,:)=[0,0,1];% third blue

AllPairs=sortrows(combnk(1:NLines,2));
NPairs=size(AllPairs,1);
[px,py,isok,d] = LineIntersection(Lines(AllPairs(:,1),:) , Lines(AllPairs(:,2),:));
% draw lines and intersections
figure(7); cla; ax=gca;axis(ax,'auto');
if DRAW_LINES
    for iline=1:size(Lines,1)
        line(ax,Lines(iline,[1 3]),Lines(iline,[2 4]),'LineWidth',2','Color',Color(iline,:)); %draw partial line defined by points x1y1x2y2
    end
    for i=1:NPairs %extraploate line to intersection point
        iline1=AllPairs(i,1);
        iline2=AllPairs(i,2);
        line(ax,[Lines(iline1,[1]),px(i)],[Lines(iline1,[2]),py(i)],'LineWidth',2','Color',Color(iline1,:));
        line(ax,[Lines(iline2,[1]),px(i)],[Lines(iline2,[2]),py(i)],'LineWidth',2','Color',Color(iline2,:));
    end
    line(ax,px,py,'LineStyle','none','Marker','*','MarkerSize',8,'Color','g');%draw intersection points
    for i=1:NPairs
        text(px(i)+0.2,py(i)+0.2,sprintf('%d',i),'FontSize',14,'Color','c')
    end
end

% part b - find regions
% 1 for each line sort all its intersections
SortIntr=cell(1,NLines); %(if no parallel lines than NintrsctnsPerline = NLines-1)
IpairiIdx1=cell(1,NLines);
IpairiIdx2=cell(1,NLines);
IpairiIdxI=cell(1,NLines);
for iline=1:NLines
    [idx1,idx2]=find(AllPairs==iline);
    intr=[px(idx1),py(idx1)];
    [~, Isortedintr]=sortrows(intr.*[1,1]); %sort by increaing x y. (prepare clockwise travers)
    SortIntr{iline}=intr(Isortedintr,:);
    IpairiIdx1{iline}=idx1(Isortedintr);
    IpairiIdx2{iline}=3-idx2(Isortedintr); %keep indexes of second line
    IpairiIdxI{iline}=Isortedintr;
end
% 2 traverse from every point and find next closest intersctionn
% we go clockwise x+x-y
PointsInPartition={};
SolIndexInPartition={};
LineIndexInPartition={};
Line4PointsInPartition={};
VisitedSequence=false(NPairs,NPairs); %skip same sequence
count_added=0; count_skipped=0;

for iline=1:NLines
    for ipoint_idx=1:length(SortIntr{iline})-1 %cant start from last point in line
        ipoint=SortIntr{iline}(ipoint_idx,:);
        if DRAW_DEBUG
            delete(findall(ax,'Tag','tmppoint'));
            line(ax,ipoint(1),ipoint(2),'LineStyle','none','Marker','O','MarkerSize',12,'Color','r','tag','tmppoint');%draw intersection points
        end
        current_line=iline;
        isol_idx=IpairiIdx1{current_line}(ipoint_idx);
        current_p_idx=ipoint_idx;
        current_l_next_p_idx=current_p_idx+1;
        next_line=AllPairs(IpairiIdx1{current_line}(current_l_next_p_idx), IpairiIdx2{current_line}(current_l_next_p_idx));
        % next_point_idx = find(IpairiIdx1{current_line}(next_point_idx)==IpairiIdx1{next_line});
        sol_idx_list=[isol_idx]; 
%       if ismember(isol_idx,[7,8,12]),keyboard;end
        point_list=[ipoint];
        line_list=[iline];
        while next_line~=iline
            if DRAW_DEBUG
                delete(findall(ax,'Tag','tmpline'));
                line(ax,Lines(current_line,[1 3]),Lines(current_line,[2 4]),'LineWidth',4','Color',[ 0,0,0 ],'Tag','tmpline');
                line(ax,Lines(next_line,[1 3]),Lines(next_line,[2 4]),'LineWidth',4','Color',[ 1,1,1 ],'Tag','tmpline');
            end
            current_sol_idx=IpairiIdx1{current_line}(current_l_next_p_idx);
            current_p_idx = find(IpairiIdx1{next_line}==current_sol_idx);
            current_line=next_line;
            current_point=SortIntr{current_line}(current_p_idx,:);
            current_nrm=norm(current_point-ipoint);
            current_o=atan2d(-current_point(2)+ipoint(2),current_point(1)-ipoint(1));

            sol_idx_list(end+1)=current_sol_idx; %#ok<SAGROW>
            point_list(end+1,:)=current_point; %#ok<SAGROW>
            line_list(end+1)=current_line;
            if DRAW_DEBUG,line(ax,current_point(1),current_point(2),'LineStyle','none','Marker','O','MarkerSize',12,'Color','m','tag','tmppoint');end %draw intersection points
            %select between two options. next clockwise point is the one with higher angle
            if current_p_idx+1<=length(SortIntr{current_line}) && current_p_idx-1>0
                next_point_1=SortIntr{current_line}(current_p_idx+1,:);
                next_point_2=SortIntr{current_line}(current_p_idx-1,:);
                if norm(next_point_1-ipoint)<NORM_EPS
                    current_l_next_p_idx=current_p_idx+1;
                elseif norm(next_point_2-ipoint)<NORM_EPS
                    current_l_next_p_idx=current_p_idx-1;
                else
                    o1=atan2d(-next_point_1(2)+ipoint(2),next_point_1(1)-ipoint(1));
                    o2=atan2d(-next_point_2(2)+ipoint(2),next_point_2(1)-ipoint(1));
                    if o1>o2
                        current_l_next_p_idx=current_p_idx+1;
                    else
                        current_l_next_p_idx=current_p_idx-1;
                    end
                end
            elseif current_p_idx-1>0
                current_l_next_p_idx=current_p_idx-1;
            else
                current_l_next_p_idx=current_p_idx+1;
            end

            next_p=SortIntr{current_line}(current_l_next_p_idx,:);
            next_o=atan2d(-next_p(2)+ipoint(2),next_p(1)-ipoint(1));
            next_nrm=norm(next_p-ipoint);
            if DRAW_DEBUG,disp([current_nrm,next_nrm,current_o,next_o]);end
            next_line=AllPairs(IpairiIdx1{current_line}(current_l_next_p_idx), IpairiIdx2{current_line}(current_l_next_p_idx));
            next_sol_idx=IpairiIdx1{current_line}(current_l_next_p_idx);
            if VisitedSequence(current_sol_idx,next_sol_idx)
                next_line=-2;
                if DRAW_DEBUG,disp('seq visited');end
                break;
            end
            if next_nrm>NORM_EPS && next_o<current_o || next_o-current_o>180
                next_line=-2;
                if DRAW_DEBUG,disp('next_o<current_o');end
                break; %non clockwise
            end
            assert(next_nrm<NORM_EPS && next_line==iline || next_nrm>=NORM_EPS && next_line~=iline);
        end

        if next_line==iline
            sol_idx_list(end+1)=next_sol_idx; %#ok<SAGROW>
            point_list(end+1,:)=next_p; %#ok<SAGROW>
            PointsInPartition{end+1}=point_list; %#ok<SAGROW>
            SolIndexInPartition{end+1}=sol_idx_list; %#ok<SAGROW>
            %Line4PointsInPartition{end+1}=(Lines(AllPairs(sol_idx_list,1),:), [1; 0])+kron(Lines(AllPairs(sol_idx_list,2),:), [0; 1]);%#ok<SAGROW>
            Line4PointsInPartition{end+1}=Lines(line_list,:);%#ok<SAGROW>

            for i=1:length(sol_idx_list)-1
                VisitedSequence(sol_idx_list(i),sol_idx_list(i+1))=true;
            end
            count_added=count_added+1;
        else
            count_skipped=count_skipped+1;
        end
        if DRAW_DEBUG, disp([next_line==iline, count_added,count_skipped]);end
    end
end

% draw all segments
if DRAW_DEBUG %clear debug
    delete(findall(ax,'Tag','tmppoint'));
    delete(findall(ax,'Tag','tmpline'));
end
NPartition=length(PointsInPartition);
s=sprintf('Lines=%d, Segments=%d, RunTime=%1.2fsec',NLines,NPartition,toc);
title(ax,s);
fprintf([s,newline]);
if DRAW_FINAL_POLY
    hold(ax,'on');
    for i=1:NPartition
        plist=PointsInPartition{i};
        patch(plist(:,1),plist(:,2),i,'FaceAlpha',.3)
        [cx,cy]=ploygon_centroid(plist(:,1),plist(:,2));
        if(DRAW_FINALTEXT),text(cx,cy,sprintf('%d',i),'FontSize',12,'Color','k');end
    end
end

function [px,py,isok,d] = LineIntersection(Line1, Line2)
    %Line1=[x1,y1,x2,y2] Line2=[x3,y3,x4,y4]
    x1=Line1(:,1); y1=Line1(:,2); x2=Line1(:,3); y2=Line1(:,4);
    x3=Line2(:,1); y3=Line2(:,2); x4=Line2(:,3); y4=Line2(:,4);
    d=(x1-x2).*(y3-y4)-(y1-y2).*(x3-x4);%determinant
    px0=(x1.*y2-y1.*x2).*(x3-x4)-(x1-x2).*(x3.*y4-y3.*x4);
    py0=(x1.*y2-y1.*x2).*(y3-y4)-(y1-y2).*(x3.*y4-y3.*x4);
    isok=abs(d)>1e-6;
    px=px0./d;
    py=py0./d;
end

function [xc,yc] = ploygon_centroid(x,y)
    xs = circshift(x,-1);
    ys = circshift(y,-1);
    area = 0.5*sum (x.*ys-xs.*y);
    xc = sum((x.*ys-xs.*y).*(x+xs))/(6*area);
    yc = sum((x.*ys-xs.*y).*(y+ys))/(6*area);
end

EDIT New solution, using itertools. Old solution below. Using the python package itertools, I believe this solution is faster than the original. However, it is still extremely slow and is not feasible beyond about 20-ish lines. For 120 lines, none of these will terminate.

import itertools

strings = ["<=", ">="]

fxs = ["5x + 3", "-60x + 7", ...]

parts = []

if __name__ == "__main__":
    print(len(fxs))
    out = itertools.product(strings, repeat=len(fxs))
    for i in out:
        curpart = ""
        for j in range(len(i)):
            print(curpart)
            if j != len(i):
                curpart = curpart + "y " + i[j] + fxs[j] + " && "
            else:
                curpart = curpart + "y " + i[j] + fxs[j]
        parts.append(curpart)
    print(parts)

OLD SOLUTION

Here's what ended up working for me, in python:

fxs = [list of functions as strings, ex: "5x+3","6x+4"... etc.]

parts = []


if __name__ == "__main__":
    for f in fxs:
        fhalf = "y >= " + f
        shalf = "y <= " + f
        if len(parts) == 0:
            parts.append(fhalf)
            parts.append(shalf)
        else:
            parts1 = [s + " && " + fhalf for s in parts]
            parts2 = [s + " && " + shalf for s in parts]
            parts = parts1 + parts2
    print(parts)

For the above example, this outputs:

['y >= 5x+3 && y >= 90x+7 && y >= -30x+28 && y >= 60x+2', 'y <= 5x+3 && y >= 90x+7 && y >= -30x+28 && y >= 60x+2', 'y >= 5x+3 && y <= 90x+7 && y >= -30x+28 && y >= 60x+2', 'y <= 5x+3 && y <= 90x+7 && y >= -30x+28 && y >= 60x+2', 'y >= 5x+3 && y >= 90x+7 && y <= -30x+28 && y >= 60x+2', 'y <= 5x+3 && y >= 90x+7 && y <= -30x+28 && y >= 60x+2', 'y >= 5x+3 && y <= 90x+7 && y <= -30x+28 && y >= 60x+2', 'y <= 5x+3 && y <= 90x+7 && y <= -30x+28 && y >= 60x+2', 'y >= 5x+3 && y >= 90x+7 && y >= -30x+28 && y <= 60x+2', 'y <= 5x+3 && y >= 90x+7 && y >= -30x+28 && y <= 60x+2', 'y >= 5x+3 && y <= 90x+7 && y >= -30x+28 && y <= 60x+2', 'y <= 5x+3 && y <= 90x+7 && y >= -30x+28 && y <= 60x+2', 'y >= 5x+3 && y >= 90x+7 && y <= -30x+28 && y <= 60x+2', 'y <= 5x+3 && y >= 90x+7 && y <= -30x+28 && y <= 60x+2', 'y >= 5x+3 && y <= 90x+7 && y <= -30x+28 && y <= 60x+2', 'y <= 5x+3 && y <= 90x+7 && y <= -30x+28 && y <= 60x+2']

The code is extremely simple and the output is exactly as desired.

The idea is basically that each line divides the plane into two halves, and so the loop (which iterates over the set of lines) intersects each partition already found with each half and adds the two new partitions to the set, while removing the original, in-intersected partition. The result doesn't always give the simplest possible condition (some of the conditions may be redundant), but they do give a complete description of each partition.

While this solution works, with 120 lines, it is pretty slow. I'd be interested in seeing if there are more efficient ways to accomplish this, using this method or otherwise.