add Scipy's conjugate gradient algorithm

nifty4/minimization/scipy_minimizer.py

@@ -113,3 +113,46 @@ def L_BFGS_B(ftol, gtol, maxiter, maxcor=10, disp=False, bounds=None):
     options = {"ftol": ftol, "gtol": gtol, "maxiter": maxiter,
                "maxcor": maxcor, "disp": disp}
     return ScipyMinimizer("L-BFGS-B", options, False, bounds)
+
+
+class ScipyCG(Minimizer):
+    def __init__(self, tol, maxiter):
+        super(ScipyCG, self).__init__()
+        if not dobj.is_numpy():
+            raise NotImplementedError
+        self._tol = tol
+        self._maxiter = maxiter
+
+    def __call__(self, energy, preconditioner=None):
+        from scipy.sparse.linalg import LinearOperator as scipy_linop, cg
+        from .quadratic_energy import QuadraticEnergy
+        if not isinstance(energy, QuadraticEnergy):
+            raise ValueError("need a quadratic energy for CG")
+
+        def toNdarray(fld):
+            return fld.to_global_data().reshape(-1)
+
+        def toField(arr, dom):
+            return Field.from_global_data(dom, arr.reshape(dom.shape))
+
+        class mymatvec(object):
+            def __init__(self, op):
+                self._op = op
+
+            def __call__(self, inp):
+                return toNdarray(self._op(toField(inp, self._op.domain)))
+
+        op = energy._A
+        b = toNdarray(energy._b)
+        sx = toNdarray(energy.position)
+        sci_op = scipy_linop(shape=(op.domain.size, op.target.size),
+                             matvec=mymatvec(op))
+        prec_op = None
+        if preconditioner is not None:
+            prec_op = scipy_linop(shape=(op.domain.size, op.target.size),
+                                  matvec=mymatvec(preconditioner))
+        res, stat = cg(sci_op, b, x0=sx, tol=self._tol, M=prec_op,
+                       maxiter=self._maxiter)
+        stat = (IterationController.CONVERGED if stat >= 0 else
+                IterationController.ERROR)
+        return energy.at(toField(res, op.domain)), stat