# solve.py

# Imports
import argparse
import matplotlib.pyplot as plt
import numpy as np

from shard import Shard

# shard_gradient
def shard_gradient(I, J, alpha, X, y, k, outside_left=True, epsilon=1e-6):
    shape = I.shape[:2]
    domain = shape[::-1]
    shard = Shard(X, k, outside_left=outside_left)

    H, dX = shard(domain, return_dX=True, epsilon=epsilon)

    d = alpha * (J - y)
    R = (I - J) + d * H[..., np.newaxis]
    J = dX[..., np.newaxis] * d
    J = J.reshape(X.size, -1).transpose()
    dX = np.dot(R.ravel(), J)
    return dX.reshape(X.shape)

# main_test_shard_gradient
def main_test_shard_gradient():
    INPUT_PATH = 'data/180,180-40,40,90,127,140,40-1,0,0-2.png'
    I = plt.imread(INPUT_PATH).astype(np.float64)
    shape = I.shape[:2]
    domain = shape[::-1]

    X = np.r_[40.0, 50, 100, 127, 140, 40].reshape(-1, 2)
    outside_left = True
    y = np.r_[1.0, 0.0, 0.0]
    k = 1.0
    alpha = 0.5

    # `J` is the initial image which the shard is blended with
    J = np.zeros_like(I)
    dX = shard_gradient(I, J, alpha, X, y, k, outside_left=outside_left)

    def apply_gradient(eta):
        X_ = X + eta * dX
        shard = Shard(X_, k, outside_left=outside_left)

        H = shard(domain)
        d = alpha * (J - y)
        R = (I - J) + d * H[..., np.newaxis]
        r = R.ravel()
        E = np.dot(r, r)
        print '%.7g -> %.7g' % (eta, E)
        return H

    etas = np.arange(20) * -0.05
    n = etas.shape[0]
    rows = int(np.floor(np.sqrt(n)))
    cols = int(np.ceil(float(n) / rows))
    w = 1.0 / cols
    h = 1.0 / rows

    f = plt.figure()
    for i, eta in enumerate(etas):
        H = apply_gradient(eta)
        Ji = (J * (1.0 - alpha * H[..., np.newaxis]) 
              + alpha * y * H[..., np.newaxis])

        r = (rows - 1) - i / cols
        c = i % cols
        ax = f.add_axes([c * w, r * h, w, h], frameon=False)
        ax.imshow(Ji)
        ax.set_xticks([])
        ax.set_yticks([])
    plt.show()

if __name__ == '__main__':
    main_test_shard_gradient()
    

@@ -83,16 +83,51 @@
  as np%0A%0A
+from scipy.optimize import fmin_cg%0A
 from sha
@@ -597,16 +597,1108 @@
 shape)%0A%0A
+# fit_shard%0Adef fit_shard(I, J, alpha, X, y, k, outside_left=True, epsilon=1e-6, **kwargs):%0A    shape = I.shape%5B:2%5D%0A    domain = shape%5B::-1%5D%0A%0A    d = alpha * (J - y)%0A    R0 = (I - J)%0A%0A    def f(x):%0A        shard = Shard(x.reshape(X.shape), k, outside_left=outside_left)%0A        H = shard(domain)%0A        R = R0 + d * H%5B..., np.newaxis%5D%0A        r = R.ravel()%0A        return np.dot(r, r)%0A%0A    def fprime(x, return_energy=False):%0A        shard = Shard(x.reshape(X.shape), k, outside_left=outside_left)%0A        H, dX = shard(domain, return_dX=True, epsilon=epsilon)%0A        R = R0 + d * H%5B..., np.newaxis%5D%0A        J = dX%5B..., np.newaxis%5D * d%0A        J = J.reshape(X.size, -1).transpose()%0A        r = R.ravel()%0A        dx = np.dot(r, J)%0A%0A        if not return_energy:%0A            return dx%0A        else:%0A            return np.dot(r, r), dx%0A%0A    states = %5B%5D%0A    def save_state(xk):%0A        states.append(np.copy(xk).reshape(X.shape))%0A%0A    xk = X.ravel()%0A    kwargs%5B'disp'%5D = 0%0A    save_state(xk)%0A    x = fmin_cg(f, xk, fprime, callback=save_state, **kwargs)%0A%0A    return x.reshape(X.shape), states%0A%0A
 # main_t

Add `fit_shard` to solve.py