minimal sphere enclosing points with MOSEK fusion

gistfile1.py

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)

plot.py

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')

sphere_enclosing_points_api.py

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)

sphere_enclosing_points_fusion.py

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)
#------------------------------------------------------------