Is this enough?

# p a polynom

type(p, polynom(anything, x)) and `not`(has(p, I))

returns true if p is a polynom in x and no coefficient is complex

I have often noticed that integration of a piecewise function can be difficult  (or even impossible), but that once the function is expressed using Heavisides,  the integration is done well.

uv := convert(u*v, Heaviside):
int(uv, x=-1..1, y=-1..1);
                               1 
                               --
                               12

I think it only be done by hand:

1. below the plot insert this

   P := [
   
   ];
   
   MyPoints := Matrix(nops(P)/2, 2, P);
   
   alignment := align=right:   # to be set in an ad hoc way
   plots:-display(
     p,
     seq(plots:-textplot(evalf[4]([MyPoints[n, 1], (MyPoints[n, 2])$2]), alignment), n=1..nops(P)/2),
     seq(plots:-pointplot([MyPoints[n, 1], MyPoints[n, 2]], symbol=cross, symbolsize=15, color=blue), n=1..nops(P)/2),
     seq(plot([[0, MyPoints[n, 2]], [MyPoints[n, 1], MyPoints[n, 2]]], linestyle=3, color=blue), n=1..nops(P)/2),
     seq(plot([[MyPoints[n, 1], 0], [MyPoints[n, 1], MyPoints[n, 2]]], linestyle=3, color=blue), n=1..nops(P)/2)
   ); 

2. right-click (control+click with Mac OSX) on the plot to select Prob Info / Nearest point on line
3. move the probe with the mouse to select a point
4. right-click (control+click) and select Probe info / Copy Data  (one single point at the times)
5. In the P:=[..] structure : paste the data (ctrl+v or command+v on Mac OSX)
6. go to step 3 to add a new point or to the next step
7. add a comma after each line but the last one in the P:=[..] structure
8. replace the decimal separator (comma) by a point
9. execute the block 

probe.mw
 

What is your definition of "minimal set"?

If you mean "the set which contains the minimum number of sets", then your second example produces a minimal set among others.

For instance the set 

{{}, {5}, {1, 6}, {4, 6}, {2, 3, 6}}

has also 5 elements whose unions generate all intersections of the elements the set of your second example

{{1, 2, 3, 6}, {1, 4, 5, 6}, {2, 3, 4, 6}}

contains.
 

Here is an idea to obtain "a set which contains the minimum number of sets".

S := {{1, 2, 3}, {1, 4, 5}, {2, 3, 4}};
N := nops(S);

eS   := op~(S);
T__0 := {{}, (`{}`~(eS))[]};

cross := table([]):
for m from 1 to N-1 do
  for n from m+1 to N do
    mn := [m, n]:
    cross[mn] := S[m] intersect S[n];
    print(mn, %);
  end do:
end do:

T__1 := {{}, entries(cross, nolist)};
T__2 := T__1 union {(eS minus op~(T__1))};

minimal_sets.mw

REMARK 1
Replacing the line

T__1 union {(eS minus op~(T__1))}

by

T__1 union {(eS minus op~(T__1))} minus {{}}

gives an "even more minimal set" which has the good properties... at least if no intersection between two elements of the initial set is the empty set.


REMARK 2
My method doesn't return a correct answer if the initial set contains an element U whose intersection with any other element is empty.
But improvements could be done

As @dharr mentioned it, using useassumptions is extremely slow (I stopped the computation even before it completed).

Here is a workaround.

STEP 1: discard obviously wrong solutions

nops(sol4cons):
                               18

# A simple way to discard solutions with b=c
sol4cons_cleaned := select(`not`@has, sol4cons, b=c):
nops(%);
                               17

# Unfortunately this works only if the solution contains the relation b=c.
# It doesn't discard sol4cons[6] which contains b=0 and c=0.
#
# A more general method
sol4cons_cleaned := map(s -> if evalb(eval(b<>c, s)) then s end if, sol4cons):
nops(%);

# numbers discarded are 

sol4cons_cleaned := map(n -> if `not`(evalb(eval(b<>c, sol4cons[n]))) then n end if, [$1..nops(sol4cons)]);
                               16
                             [4, 6]

The problem is that the remaining solutions are not correct, for instance 

sol4cons[1]
 /                              1           2                  
{ a = a, b = b, c = c, lambda = - alpha A[0] , mu = mu, p = p, 
 \                              3                              

                                  \ 
  A[-1] = 0, A[0] = A[0], A[1] = 0 }
                                  / 

and nothing insures that a^2-b^2+c^2 > 0.


STEP 2: do a second pass on sol4cons_cleaned
So there is another workaround:
For each remaining solution in sol4cons_cleaned : "re-solve" the solutions with the constraints a^2-b^2+c^2 > 0 and b<>c:

goodsols := map(s -> solve(s union {b<>c, a^2-b^2-c^2>0}), sol4cons_cleaned):

There remain two problems

• Two  solve fail:

• n := 14: # and n:=15
  sol4cons_cleaned[n]:
  solve(sol4cons_cleaned[n] union {b<>c, a^2-b^2-c^2>0})
  Error, (in unknown) numeric exception: division by zero
  

I guess this mean that "preliminary" solutions number 14 and 15 cannot be solutions that respect the constraints, but I'm not at all sure of this interpretation
 

• the 5th solution takes a long time to be re-solved (I stopped the process)
   

Final remark
I guess useassumption leads to extremely time consuming computations because there is a huge number of solutions to compute.
For instance, re-solving sol4cons_cleaned[1] generates 12 solutions:

n := 1:
sol4cons_cleaned[n]:
[solve(sol4cons_cleaned[n] union {b<>c, a^2-b^2-c^2>0})]:
nops(%);
                           12
print~(%%)

You were very close to the solution

n := `2`:
p := ListTools:-Search(n, lhs~(S));
[n, rhs(S[p]):-mu, rhs(S[p]):-sigma]
                               3
             [2, 1208.001018040, 141.512246146314]

A solution which reuses the most part of what you did (this was my my constrraint)
 

restart:
with(ListTools):
C := x*u[]*u[1,2,2]^3*u[1,1,2];
indx_set := [ (`[]`@ op)~(indets(C, indexed))[] ];

deg      := degree~(C, [indets(C, indexed)[]]);
Indx_set := zip((u,v) -> u$v, indx_set, deg);

Flatten(convert(Indx_set,list));

ListTools:-Occurrences(1, %);
                        5

I'm not sure at all I understood correctly, but maybe this could help?
if I'm out of line just forget it
 

# If only the sums of terms with ranks 5, 10, 15, ...matter, what is the need to 
# compute U and V?

restart:
u := (a, b) -> exp(-(a-Lmu)^2)*exp(-(b-Wmu)^2);
llim := 1:
ulim := 20000:
step := 5:
nok := iquo(ulim-llim+1, 5);
        4000
W := CodeTools:-Usage( Vector(nok, i -> add(u(5*i, x__2)*u(x__4, 5*i))) );

memory used=19.89MiB, alloc change=36.00MiB, cpu time=313.00ms, real time=314.00ms, gc time=21.15ms

Download one_way.mw

This was intended to answer your last question, but it seems to have disappeared.
The procedure Circle in the attached file finds, if it does exist, the "tritangent" circle for arbitrary rational fraction P(x)/Q(x) where:

• P(x) is of degree at most 2
• Q(x) has degree 2 and two real roots

The procedure Circle solves the equation 

(x-p)^2+(f(x)-q)^2-r^2=0

whose numerator is a polynom of degree 6 in x.
The circle of centre [p, q] and radius r has 3 tangency point if this numerator is of the form 

(x-A[1])^2 * (x-A[2])^2 * (x-A[3])^2

After having equalized the coefficients of these two equations the procedure solves a system of 6 equations (p, q, r, plus  the abscissas A[1], A[2], A[3] of the three tangency points).


The file contains two "pathological" examples:

• In the first one, even if it seems visually that a solution exists, procedure Circle fails to find one.
  In procedure Circle change the line

  sol := {map(s -> if keep(s)  then s end if, sols)[]};

  by the line

  sol := sols;

  to see what happens: 

  • the function keep is aimed to discard solutions with nerative r, identical contact points or imaginary contact points

  • writting sol := sols instead will return the "unfiltered" solutions

• The second one is very interesting because it shows that for some functions f(x) with three branches, there can exist a circle which tangents a branch once and a second branch twice.
   

Another_way.mw

Kitonum's variant.
Use tubeplot instead of spacecurve:

• pros: doesn't require a fictitious shift of the ellipse (with is to be adjusted in an ad hoc way)
• cons: for scaling=unconstrained (the default value) the tube doesn't have a circular cross section but an elliptical one: if you increase the value of radius the ellipse will ressemble more to a short thin sylibder than to a line.

So make your choice

tubeplot([cos(t), sin(t), 1 - 12*cos(t) + 15*sin(t)], t = 0 .. 2*Pi, radius=0.05, color=red, orientation = [-15, 68, 5])

Seems there is no solution.

Solution without range constraints

N := 4:
sys := NULL:
for k to N while k <> j do 
sys := sys, 
       sinh(x[k]+I*a)^N*(product(sinh(1/2*(x[k]-x[j]-(2*I)*a)), j = 1 .. N))
       /
       (sinh(x[k]-I*a)^N*(product(sinh(1/2*(x[k]-x[j]+(2*I)*a)), j = 1 .. N))) 
end do:
sys := [sys]:

# check
# print~(sys):

fsolve(sys=~[-1$4], {x[1], x[2], x[3], x[4]});
                               
            {x[1] = -3.799588601 - 1.570796327 I, 

              x[2] = -3.799588601 - 1.570796327 I, 

              x[3] = 0.1196098914 - 1.570796327 I, 

              x[4] = 1.110549099 - 1.570796327 I}

# check solution (t[1]..t[4] are your original variables) 
# eval([seq(t[i], i=1..4)], %);

You cannot impose ranges the way you did.
The correct way is to write something like a1+b1*I .. a2+b2*I
For instance, with a1=-10 and a2=+10, b1=-1, b2=+1 one gets

fsolve(sys=~[-1$4], {x[1]=-10-1*I..10+1*I, x[2]=-10-1*I..10+1*I, x[3]=-10-1*I..+10+1*I, x[4]=-10-1*I..+10+1*I});
            {x[1] = 0.2343030315 - 0.9904555257 I, 

              x[2] = 0.9120742051 - 0.5349352395 I, 

              x[3] = 0.2227793496 + 0.3919592237 I, 

              x[4] = 0.2227793496 + 0.3919592237 I}

For the imaginary range you want there is unfortunately no solutions found (mayba you will find for other values of a1 and a2???)

fsolve(sys=~[-1$4], {x[1]=-10-I..10-0.5*I, x[2]=-10-I..10-0.5*I, x[3]=-10+0.5*I..+10+I, x[4]=-10+0.5*I..+10+I});

Remark
Depending on the way the system is written  fsolve will find different solutions.
Example:

sys_nd := map(simplify, numer~(sys)+~denom~(sys)):
fsolve(sys_nd);

            {x[1] = 0.7772186938 + 0.1314607556 I, 

              x[2] = 0.7772186938 + 0.1314607556 I, 

              x[3] = -0.7600984061 - 1.168722678 I, 

              x[4] = 0.7772186938 + 0.1314607556 I}

#check the solution
#eval([seq(t[i], i=1..4)], %);

But even formulated this way there is still nosolution found for the imaginary ranges you want:

fsolve(sys_nd, {x[1]=-10-I..10-0.5*I, x[2]=-10-I..10-0.5*I, x[3]=-10+0.5*I..+10+I, x[4]=-10+0.5*I..+10+I});

Shorter but less elegant than @MapleMathMatt 's, it requires a few knowledge on MathML syntax
 

sinx := `#mrow(mo("sin",mathcolor="red"),mo("(",mathcolor="red"),mo("x",mathcolor="red"),mo(")",mathcolor="red"))`:
plot( sin(x), x = 0 .. 2 * Pi, title = typeset( "Plot of %1", sinx ) );

If the space between "(" and "x") doesn't suit you, replace the definition of sinx by
 

sinx := `#mrow(mo("sin",mathcolor="red"),mo("(",mathcolor="red"),mo(" "),mo("x",mathcolor="red"),mo(")",mathcolor="red"))`:

A more formal way
 

restart;
with(codegen):

B := p -> x^4 - 4*x^3/(2 + p) - 6*(p - 1)*x^2/((3 + p)*(2 + p)) - 4*p*(p - 5)*x/((4 + p)*(3 + p)*(2 + p)):

# formal roots of B(p)

s := solve(B(p), x):
a  := [allvalues(s)]:

# The problem could end here, but plotting the solution a for different values of 
# p would require a lot of time for a has a very lengthy expression as cost 
# (codegen:-cost) demonstrates

cost(a[1]);
1246*additions + 6298*multiplications + 136*divisions + 162*functions

# the idea is to use the package codegen to rewrite a into a form which 
# minimizes the number of operations it contains.
#
# first step : transform each a[k] into a procedure with argument p

for k from 1 to 4 do
  F||k := makeproc(evalf(a[k]),p):
end do:

# second step : join the 4 procedures F1..F4

G := joinprocs([seq(F||k, k=1..4)]):

# third step : optimize the joint procedure G

H := optimize(G):

# cost of H (must be compared to 4 times the cost of a[1])

cost(H)
41 storage + 41 assignments + 107 multiplications + 90 additions + 15 divisions + 10 functions

# plot

CodeTools:-Usage(
  plots:-complexplot(
    [seq(H(q), q=0..5000)], 
    style = point, symbol=solidcircle, symbolsize=6
  )
);
memory used=170.05MiB, alloc change=32.00MiB, cpu time=1.97s, real time=1.80s, gc time=258.09ms

# for comparison

CodeTools:-Usage(
  plots:-complexplot(
    [seq(op(evalf(eval(a, p=q))), q=0..5000)], 
    style = point, symbol=solidcircle, symbolsize=6
  )
);
memory used=0.72GiB, alloc change=32.00MiB, cpu time=6.72s, real time=6.29s, gc time=624.36ms

As you can see in the attached file
exact_plus_codegen.mw
the formal approach is about 3 times faster than the numerical one and uses about one third of the memory this latter uses.
 

If you have written something like

a := 0.4:
:
:
:
plot(a, a=0..1)

and if you want ti unassign a before the plot command, just do

a := 0.4:
:
:
:
a := 'a':
plot(a, a=0..1)

This answer is based upon Maple 2015.2
If you have more recent version of Maple (which I don't at home) , DrawGraph proposes a lot feautures to customize the graph.
In any case, whatever the version, you can always modify the plot the way you want.

Here is an example.
Note that at some point I use GetVertexPosition to catch the position of the vertices: you can modify them as you want (SetVertexPosition).
As shown in tha attached file, you can also dilate a graph in a simple way by using  plottools:-transform.
 

small_vertex_names.mw