Skip to content

Instantly share code, notes, and snippets.

@andreacassioli
Last active August 29, 2015 13:56
Embed
What would you like to do?
minimal sphere enclosing points with MOSEK fusion
import sys
import mosek
from mosek.fusion import *
import numpy as np
import time
import plot
def sphere_enclosing_primal(n,k,p):
print('Creating the Fusion optimization model...')
start=time.clock()
M = Model("minimal sphere enclosing a set of points primal")
print('done!\n')
print('Declaring the variables...')
r0 = M.variable("r_0", Domain.unbounded())
p0 = M.variable("p_0", NDSet(1,n), Domain.unbounded())
print('done!\n')
Expr.sub( Variable.repeat(p0,k), DenseMatrix(p) )
print('Defining the constraints...')
M.constraint( Expr.hstack( Variable.repeat(r0,k),Expr.sub( Variable.repeat(p0,k), DenseMatrix(p) ) ), Domain.inQCone())
print('done!\n')
print('Defining the objective function...')
M.objective('min_radius', ObjectiveSense.Minimize, r0)
print('done!\n')
print('Model building time: '+ str(time.clock()-start))
M.setLogHandler(sys.stdout)
M.solve()
print("\nSolution status= "+ str(M.getPrimalSolutionStatus()))
rstar= r0.level()[0]
pstar= p0.level()
print("r*= %f" % rstar)
print("p*= ")
print(pstar)
#comment out this block if you haven't matplotlib installed
if n==2:
plot.plot_points(p,pstar,rstar)
#------------------------------------------------------------
def sphere_enclosing_dual(n,k,p):
print('Creating the Fusion optimization model...')
start=time.clock()
M = Model("minimal sphere enclosing a set of points dual")
print('Declaring the variables...')
y= M.variable('y',NDSet(k,n+1),Domain.inQCone(k,n+1))
print('done!\n')
c=[0. for i in range(n+1)]
c[0]=1.
print('Defining the constraints...')
M.constraint('lin', Expr.mul(DenseMatrix(1,k,1.0), y), Domain.equalsTo(c) )
print('done!\n')
print('Defining the objective function...')
M.objective('dual',ObjectiveSense.Maximize, Expr.sum(Expr.mulDiag(\
DenseMatrix(p.tolist()),\
y.slice([0,1],[k,n+1]).transpose())))
print('done!\n')
print('Model building time: '+ str(time.clock()-start))
M.setLogHandler(sys.stdout)
M.solve()
print("\nSolution status= "+ str(M.getPrimalSolutionStatus()))
k=int(sys.argv[1]) #number of points
n=int(sys.argv[2]) #dimension
if n<2:
print("the dimension must be greater than one!")
exit(1)
print('Generate '+str(k)+' points in R^'+str(n)+' normally distributed around the origin...')
p= np.random.normal(size=(n,k))
print('done!\n')
print('call the primal..')
sphere_enclosing_primal(n,k,p)
print('----------------------------------------------------------------------\n\n')
print('call the dual..')
sphere_enclosing_dual(n,k,p)
import numpy as np
import numpy.linalg
import scipy.spatial
import math
#comment out this block if you haven't matplotlib installed
import matplotlib
import matplotlib.pyplot as plt
#--------------------------------------------------------
def plot_points(p,pstar,rstar):
(k,n)=p.shape
print('Plotting...')
fig, ax = plt.subplots()
ax.set_aspect('equal')
import matplotlib.patches as mpatches
c = mpatches.Circle( pstar, rstar , fc="w", ec="r", lw=1.5)
ax.add_patch(c)
ax.plot( p[:,0], p[:,1], 'b.')
ax.plot( pstar[0], pstar[1], 'ro' )
hull= scipy.spatial.ConvexHull(p)
for simplex in hull.simplices:
plt.plot(p[simplex,0], p[simplex,1], 'g*--', lw=1.5)
boundary= [ p[i,:] for i in range(n) if math.fabs(np.linalg.norm(pstar -p[i,:] )-rstar )<1e-06]
print(boundary)
for b in boundary:
ax.plot( b[0], b[1], 'go')
plt.show()
print('done!\n')
import sys
import mosek
import numpy as np
import plot
#--------------------------------------------------------
def streamprinter(text):
sys.stdout.write(text)
sys.stdout.flush()
k=int(sys.argv[1]) #number of points
n=int(sys.argv[2]) #dimension
p= np.random.normal(size=(k,n))
env = mosek.Env()
task= env.Task(0,0)
numvar= 1 + n + k*n + k
offset_r= 0
offset_p= 1
offset_w= 1+n
offset_t= 1+n+n*k
task.appendvars(numvar)
task.putvarboundslice(0,numvar,[mosek.boundkey.fr for i in range(numvar)], np.zeros(numvar), np.zeros(numvar))
numcon= 2*k + k
task.appendcons(numcon)
for i in range(k):#for each point
task.putarow(i, [0,i+offset_t], [1., -1,]) #delta coord#radius aux
task.putconbound(i, mosek.boundkey.fx, 0. ,0.)
for j in range(n):
task.putconbound(k+ i*n + j, mosek.boundkey.fx, p[i][j] , p[i][j])
task.putarow(k+i*n + j,[offset_p + j, offset_w+i*n + j],[1.0,-1.0])
task.appendcone(mosek.conetype.quad,0., np.hstack([ offset_t+i, [ offset_w+i*n + j for j in range(n)]] ))
task.putcj(0,1.0)
task.putobjsense(mosek.objsense.minimize)
task.set_Stream (mosek.streamtype.log, streamprinter)
task.writedata("dump.opf")
task.optimize()
xx = np.zeros(numvar, float)
task.getxx(mosek.soltype.itr,xx)
r0= xx[0]
p0= xx[1:n+1]
print(r0,p0)
#comment out this block if you haven't matplotlib installed
if n==2:
plot.plot_points(p,p0,r0)
import sys
import mosek
from mosek.fusion import *
import numpy as np
import plot
k=int(sys.argv[1]) #number of points
n=int(sys.argv[2]) #dimension
if n<2:
print("the dimension must be greater than one!")
exit(1)
print('Generate '+str(k)+' points in R^'+str(n)+' normally distributed around the origin...')
p= np.random.normal(size=(k,n))
print('done!\n')
print('Creating the Fusion optimization model...')
M = Model("minimal sphere enclosing a set of points")
print('done!\n')
print('Declaring the variables...')
r = M.variable("r", Domain.unbounded())
p0 = M.variable("p_0", NDSet(1,n), Domain.unbounded())
print('done!\n')
print('Defining the constraints...')
M.constraint('norm', Expr.hstack( [ Expr.mul(np.ones(k),r), Expr.sub(Expr.mul(np.ones( (k,1) ), p0 ), DenseMatrix(p) ) ] ), Domain.inQCone(k,n+1))
print('done!\n')
print('Defining the objective function...')
M.objective('min_radius', ObjectiveSense.Minimize, r)
print('done!\n')
M.setLogHandler(sys.stdout)
M.solve()
print("\nSolution status= "+ str(M.getPrimalSolutionStatus()))
rstar= r.level()[0]
pstar= p0.level()
print("r*= %f" % rstar)
print("p*= ")
print(pstar)
#comment out this block if you haven't matplotlib installed
if n==2:
plot.plot_points(p,pstar,rstar)
#------------------------------------------------------------
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment