johanlofberg/moresos.m

## moresos.m
%% Naive sos model
x = sdpvar(1,1);
y = sdpvar(1,1);
z = sdpvar(1,1);
p = 12+y^2-2*x^3*y+2*y*z^2+x^6-2*x^3*z^2+z^4+x^2*y^2;
options = sdpsettings('sos.newton',0,'sos.congruence',0);
[sol,v,Q] = solvesos(sos(p),[],options);

%% Newton polytope
options = sdpsettings('sos.newton',1,'sos.congruence',0);
[sol,v,Q] = solvesos(sos(p),[],options);

%% Visualize Newton polytope
clf
N = getexponentbase(newtonmonoms(p),[x y z]);
s = sdpvar(length(N),1);q = sdpvar(3,1);
plot([q == N'*s, sum(s)==1,s>=0],q,'b',[],sdpsettings('plot.shade',0.2))
hold on
m = plot3(N(:,1),N(:,2),N(:,3),'yo');
set(m,'Markersize',10);
set(m,'Markerface','yellow');

%% Symmetry reduction
options = sdpsettings('sos.newton',1,'sos.congruence',1);
[sol,v,Q] = solvesos(sos(p),[],options);


%% Post-processing
options = sdpsettings('sos.numblkdg',1e-4);
[sol,v,Q] = solvesos(sos(p),[],options);
Q{1}

x = sdpvar(1,1);
y = sdpvar(1,1);
z = sdpvar(1,1);
w = sdpvar(1,1);
p = 8*x^6*w^6*y^2+4*z^4*y^8*x^2+4*z^4*y^4*x^4+8*z^4*y^6*x^2+...
4*x^8*w^6*y^2+4*z^4*y^4*x^2+4*z^6*y^8*x^2+8*y^2*x^6*w^4+...
4*x^4*y^4*w^4-2*x^6*y^6-16*z^4*y^4*x^4*w^2-16*z^4*y^4*x^2*w^4-...
20*z^4*y^4*x^6*w^2-24*z^4*y^4*x^6*w^4+12*z^4*y^6*x^2*w^2-...
8*z^4*y^6*x^2*w^4+4*z^4*y^6*x^4*w^2-24*z^4*y^6*x^4*w^4-...
16*z^4*y^6*x^6*w^2-40*z^4*y^4*x^4*w^4-4*z^4*y^2*x^2*w^2-...
8*z^4*y^2*x^2*w^4-12*z^4*y^2*x^4*w^2-16*z^4*y^2*x^4*w^4-...
12*z^4*y^2*x^6*w^2-8*z^4*y^2*x^6*w^4-4*z^4*y^2*x^8*w^2-...
16*z^4*y^6*x^6*w^4+8*z^4*y^8*x^2*w^2+12*z^4*y^8*x^4*w^2-...
6*y^2*x^6*z^2*w^2-2*y^2*x^8*z^2*w^2-14*x^4*y^4*z^2*w^2-...
6*x^2*y^4*z^2*w^2-6*x^4*y^2*z^2*w^2-6*x^2*y^6*z^2*w^2-...
10*x^4*y^6*z^2*w^2-10*y^4*x^6*z^2*w^2-2*x^2*y^8*z^2*w^2-...
20*x^6*y^6*z^2*w^2+6*y^4*x^8*z^2*w^2+6*x^4*y^8*z^2*w^2+...
12*z^2*w^4*y^2*x^6-12*z^2*w^4*x^2*y^4-16*z^2*w^4*x^4*y^4+...
4*z^2*w^4*x^6*y^4-12*z^2*w^4*x^2*y^6-20*z^2*w^4*x^4*y^6+...
8*z^2*w^4*x^8*y^2+4*x^4*y^2*w^4+x^4*w^4+2*x^6*w^4+x^8*w^4+...
z^4*y^4+2*z^4*y^6+z^4*y^8+x^8*y^4+x^4*y^8+8*x^6*y^4*w^4+4*x^8*y^2*w^4+...
8*z^4*y^6*x^4+6*z^4*y^8*x^4-4*z^4*y^6*x^6+2*z^4*y^4*x^8+...
3*z^4*y^4*w^2+2*z^4*y^4*w^4+6*z^4*y^6*w^2+4*z^4*y^6*w^4+...
3*z^4*y^8*w^2+2*z^4*y^8*w^4-6*x^6*y^6*z^2+3*y^4*x^8*z^2+...
3*x^4*y^8*z^2-6*x^6*y^6*w^2+3*x^8*y^4*w^2+3*x^4*y^8*w^2+...
3*z^2*w^4*x^4+6*z^2*w^4*x^6+3*z^2*w^4*x^8+2*z^4*w^4*x^4+...
4*z^4*w^4*x^6+2*z^4*w^4*x^8-16*z^2*w^4*x^6*y^6-4*z^2*w^4*x^2*y^8+...
12*z^2*w^4*x^8*y^4-4*x^2*w^4*y^2*z^2-2*x^2*y^2*z^2*w^2+...
8*x^8*y^4*z^2*w^6-4*x^6*w^4*y^6+6*x^8*w^4*y^4+2*x^4*w^4*y^8+...
2*x^4*z^2*w^6+2*y^4*z^6*w^2+y^4*z^6+x^4*w^6+2*y^6*z^6+4*x^4*y^8*z^6+...
4*x^8*y^4*w^6+4*x^2*y^4*z^6+8*x^2*y^6*z^6+4*x^4*y^4*z^6+...
2*y^8*z^6*w^2+8*x^6*w^6*y^4+4*x^4*w^6*y^2+4*x^4*w^6*y^4+...
2*x^8*w^6*z^2+8*x^4*y^8*z^6*w^2+8*x^2*y^4*z^6*w^2+16*x^2*y^6*z^6*w^2+...
8*x^4*y^4*z^6*w^2+16*x^6*w^6*y^4*z^2+8*x^8*w^6*y^2*z^2+...
8*x^4*w^6*y^4*z^2+8*x^4*w^6*y^2*z^2+16*x^6*w^6*y^2*z^2+...
y^8*z^6+2*x^6*w^6+x^8*w^6+4*x^6*w^6*z^2+16*x^4*y^6*z^6*w^2+...
8*z^6*y^6*x^4+4*z^6*y^6*w^2+8*x^2*y^8*z^6*w^2
[sol,v,Q,res] = solvesos(sos(p))
[sol,v,Q,res] = solvesos(sos(p),[],sdpsettings('sos.numblk',1e-6))
spy(Q{1})
	%% Naive sos model
	x = sdpvar(1,1);
	y = sdpvar(1,1);
	z = sdpvar(1,1);
	p = 12+y^2-2x^3y+2yz^2+x^6-2x^3z^2+z^4+x^2*y^2;
	options = sdpsettings('sos.newton',0,'sos.congruence',0);
	[sol,v,Q] = solvesos(sos(p),[],options);

	%% Newton polytope
	options = sdpsettings('sos.newton',1,'sos.congruence',0);
	[sol,v,Q] = solvesos(sos(p),[],options);

	%% Visualize Newton polytope
	clf
	N = getexponentbase(newtonmonoms(p),[x y z]);
	s = sdpvar(length(N),1);q = sdpvar(3,1);
	plot([q == N'*s, sum(s)==1,s>=0],q,'b',[],sdpsettings('plot.shade',0.2))
	hold on
	m = plot3(N(:,1),N(:,2),N(:,3),'yo');
	set(m,'Markersize',10);
	set(m,'Markerface','yellow');

	%% Symmetry reduction
	options = sdpsettings('sos.newton',1,'sos.congruence',1);
	[sol,v,Q] = solvesos(sos(p),[],options);


	%% Post-processing
	options = sdpsettings('sos.numblkdg',1e-4);
	[sol,v,Q] = solvesos(sos(p),[],options);
	Q{1}

	x = sdpvar(1,1);
	y = sdpvar(1,1);
	z = sdpvar(1,1);
	w = sdpvar(1,1);
	p = 8x^6w^6y^2+4z^4y^8x^2+4z^4y^4x^4+8z^4y^6x^2+...
	4x^8w^6y^2+4z^4y^4x^2+4z^6y^8x^2+8y^2x^6w^4+...
	4x^4y^4w^4-2x^6y^6-16z^4y^4x^4w^2-16z^4y^4x^2*w^4-...
	20z^4y^4x^6w^2-24z^4y^4x^6w^4+12z^4y^6x^2w^2-...
	8z^4y^6x^2w^4+4z^4y^6x^4w^2-24z^4y^6x^4w^4-...
	16z^4y^6x^6w^2-40z^4y^4x^4w^4-4z^4y^2x^2w^2-...
	8z^4y^2x^2w^4-12z^4y^2x^4w^2-16z^4y^2x^4w^4-...
	12z^4y^2x^6w^2-8z^4y^2x^6w^4-4z^4y^2x^8w^2-...
	16z^4y^6x^6w^4+8z^4y^8x^2w^2+12z^4y^8x^4w^2-...
	6y^2x^6z^2w^2-2y^2x^8z^2w^2-14x^4y^4z^2w^2-...
	6x^2y^4z^2w^2-6x^4y^2z^2w^2-6x^2y^6z^2w^2-...
	10x^4y^6z^2w^2-10y^4x^6z^2w^2-2x^2y^8z^2w^2-...
	20x^6y^6z^2w^2+6y^4x^8z^2w^2+6x^4y^8z^2w^2+...
	12z^2w^4y^2x^6-12z^2w^4x^2y^4-16z^2w^4x^4y^4+...
	4z^2w^4x^6y^4-12z^2w^4x^2y^6-20z^2w^4x^4y^6+...
	8z^2w^4x^8y^2+4x^4y^2w^4+x^4w^4+2x^6w^4+x^8*w^4+...
	z^4y^4+2z^4y^6+z^4y^8+x^8y^4+x^4y^8+8x^6y^4w^4+4x^8y^2w^4+...
	8z^4y^6x^4+6z^4y^8x^4-4z^4y^6x^6+2z^4y^4x^8+...
	3z^4y^4w^2+2z^4y^4w^4+6z^4y^6w^2+4z^4y^6w^4+...
	3z^4y^8w^2+2z^4y^8w^4-6x^6y^6z^2+3y^4x^8z^2+...
	3x^4y^8z^2-6x^6y^6w^2+3x^8y^4w^2+3x^4y^8w^2+...
	3z^2w^4x^4+6z^2w^4x^6+3z^2w^4x^8+2z^4w^4x^4+...
	4z^4w^4x^6+2z^4w^4x^8-16z^2w^4x^6y^6-4z^2w^4x^2y^8+...
	12z^2w^4x^8y^4-4x^2w^4y^2z^2-2x^2y^2z^2w^2+...
	8x^8y^4z^2w^6-4x^6w^4y^6+6x^8w^4y^4+2x^4w^4*y^8+...
	2x^4z^2w^6+2y^4z^6w^2+y^4z^6+x^4w^6+2y^6z^6+4x^4y^8*z^6+...
	4x^8y^4w^6+4x^2y^4z^6+8x^2y^6z^6+4x^4y^4z^6+...
	2y^8z^6w^2+8x^6w^6y^4+4x^4w^6y^2+4x^4w^6y^4+...
	2x^8w^6z^2+8x^4y^8z^6w^2+8x^2y^4z^6w^2+16x^2y^6z^6*w^2+...
	8x^4y^4z^6w^2+16x^6w^6y^4z^2+8x^8w^6y^2z^2+...
	8x^4w^6y^4z^2+8x^4w^6y^2z^2+16x^6w^6y^2z^2+...
	y^8z^6+2x^6w^6+x^8w^6+4x^6w^6z^2+16x^4y^6z^6*w^2+...
	8z^6y^6x^4+4z^6y^6w^2+8x^2y^8z^6w^2
	[sol,v,Q,res] = solvesos(sos(p))
	[sol,v,Q,res] = solvesos(sos(p),[],sdpsettings('sos.numblk',1e-6))
	spy(Q{1})