In [1]:
%matplotlib inline
import numpy as np
import matplotlib as mpl

import matplotlib.pyplot as plt



A = np.mat('[1, 2; 3, 4; 5 0; -1, -5; -1 1; -1 0; 0.1 1]')
b = np.mat('[1;2;1; 0; 2; 1.35; 0.8]')
#print(A)

def oracle(x, A, b):
    y = A*x - b
    return np.all(y <= 0)

def binsearch(x, d, A, b):
    t = d
    while oracle(x+t, A, b):
        t = 2*t    
    l = x
    r = x+t 
    for i in range(20):
        m = (l+r)/2
        if oracle(m, A, b):
            l = m
        else:
            r = m   
    return m

for i in range(1000):
    x0 = np.random.randn(2,1)
    if oracle(x0, A, b):
        break

#print(y.T)
#plt.plot(x[0], x[1] ,'.')
x = x0
plt.figure(figsize=(5,5))
plt.gca().set_xlim((-1.5,0.3))
plt.gca().set_ylim((-0.2,1))
for i in range(500):
    d = np.random.randn(2,1)
    d = d/np.linalg.norm(d)
    
    left = binsearch(x, d, A, b)
    right = binsearch(x, -d, A, b)
    
    lam = np.random.rand(1,1)
    x = (lam*left+(1-lam)*right)
    # plt.plot([x[0],y[0]],[x[1],y[1]])
    plt.plot(x[0],x[1],'.')


Sample from $\exp(-\frac{1}{T} c^\top x)$

where $x = \lambda x_0 + (1-\lambda) x_1$

$c_0 = c^\top x_0$ and $c_1 = c^\top x_1$

$\exp(-\frac{1}{T} ( \lambda (c_0-c_1) + c_1)) \propto \exp(-\frac{\lambda}{T} (c_0-c_1)) $

\begin{eqnarray} \int \exp(-\frac{(c_0-c_1)}{T} \tau ) d \tau & = & -\frac{T}{(c_0-c_1)} \exp(-\frac{(c_0-c_1)}{T} \tau ) |_{0}^\lambda \\ & = & -\frac{T}{(c_0-c_1)} (1-\exp(-\frac{(c_0-c_1)}{T} \lambda) ) \\ Z & = & \frac{1}{\beta} (\exp(-\beta)-1) \end{eqnarray}\begin{eqnarray} u & = & \frac{1}{Z \beta} (\exp(-\beta \lambda) -1 ) \end{eqnarray}

draw $u \sim \mathcal{U}[0, 1] $

Compute

$$ -\frac{1}{\beta} \log(Z \beta u + 1) = \lambda $$$$ -\frac{1}{\beta} \log(exp(-\beta)u - u + 1) = \lambda $$

In [2]:
c0 = 1000
c1 = 1000
T = 0.001;

def draw_lambda(c0, c1, T):
    beta = (c0-c1)/T
    if beta== 0:
        return np.random.rand(1)
    else:
        u = np.random.rand(1)
        r = np.max([np.log(1-u),-beta+np.log(u)])
        lam = -1.0/beta * (r + np.log(np.exp(-beta + np.log(u) - r) + np.exp(np.log(1-u)-r)))
        return lam

In [3]:
%matplotlib inline
import time
import pylab as plt
from IPython import display

x = x0
x = np.mat('[-0.5;0.4]')
#c = np.mat('[-1;-0.5]')
c = np.mat('[0;0]')
T = 0.1

fig = plt.figure(figsize=(5,5))
plt.gca().set_xlim((-1.5,0.3))
plt.gca().set_ylim((-0.2,1))
ax = fig.gca()
ln = plt.Line2D(xdata=(x0[0], x0[0]), ydata=(x0[1], x0[1]), marker='o', linestyle='-',linewidth=1)
ax.add_line(ln)

for i in range(500):
    d = np.random.randn(2,1)
    d = d/np.linalg.norm(d)
    
    left  = binsearch(x, d, A, b)
    right = binsearch(x,-d, A, b)
    
    ln.set_xdata((left[0],right[0]))
    ln.set_ydata((left[1],right[1]))
    
    c0 = c.T*left
    c1 = c.T*right
    
    lam = draw_lambda(c0, c1, T=1/(i+1))
    #x = (lam[0,0]*left+(1-lam[0,0])*right)
    x = (lam[0]*left+(1-lam[0])*right)
    # plt.plot([x[0],y[0]],[x[1],y[1]])

    display.clear_output(wait=True)
    display.display(plt.gcf())
    time.sleep(0.05)
    plt.plot(x[0],x[1],'bo', alpha=0.5, markersize=5)



In [ ]:
%matplotlib inline
import time
import pylab as plt
from IPython import display

x = x0
x = np.mat('[-0.5;0.4]')
#c = np.mat('[-1;-0.5]')
c = np.mat('[0;0]')
T = 0.1

fig = plt.figure(figsize=(5,5))
plt.gca().set_xlim((-1.5,0.3))
plt.gca().set_ylim((-0.2,1))
ax = fig.gca()
ln = plt.Line2D(xdata=(x0[0], x0[0]), ydata=(x0[1], x0[1]), marker='o', linestyle='-',linewidth=1)
ax.add_line(ln)

for i in range(500):
    d = np.random.randn(2,1)
    d = d/np.linalg.norm(d)
    
    left  = binsearch(x, d, A, b)
    right = binsearch(x,-d, A, b)
    
#    ln.set_xdata((left[0],right[0]))
#    ln.set_ydata((left[1],right[1]))
    
    c0 = c.T*left
    c1 = c.T*right
    
    lam = draw_lambda(c0, c1, T=1/(i+1))
    #x = (lam[0,0]*left+(1-lam[0,0])*right)
    xx = (lam[0]*left+(1-lam[0])*right)
    # plt.plot([x[0],y[0]],[x[1],y[1]])

    display.clear_output(wait=True)
    display.display(plt.gcf())
    time.sleep(0.05)
#    plt.plot(xx[0],xx[1],'bo', alpha=0.5, markersize=5)   
    plt.plot(np.hstack((left[0],right[0])), np.hstack((left[1],right[1])) , '-', alpha=0.5)

In [ ]:
np.hstack((left[0],right[0]))

In [ ]:
%matplotlib inline
fig = plt.figure(figsize=(5,5))
plt.gca().set_xlim((-1,1))
plt.gca().set_ylim((-1,1))
ax = fig.gca()
theta = np.pi/128
T = np.mat([[np.cos(theta), np.sin(theta)],[-np.sin(theta), np.cos(theta)]])

r = np.mat('[1, -1;0, 0]')

ln = plt.Line2D(xdata=r[0,:], ydata=r[1,:], marker='o', linestyle='-',linewidth=1)
ax.add_line(ln)

for i in range(20):
    ln.set_xdata(r[0,:])
    ln.set_ydata(r[1,:])
    #plt.plot()
    
    r = T*r
    
    display.clear_output(wait=True)
    display.display(plt.gcf())
    time.sleep(0.1)
    #plt.show()

In [ ]:
a = np.mat('[1,2]')
b = np.mat('[1;2]')

print(a*b)

a = np.array(a)
b = np.array(b)

a+b
a*b

In [3]:
import pylab as plt

ax = plt.gca()
ax.set_aspect

In [14]:
np.arctan(1/-0.0001)


Out[14]:
-1.5706963267952299