def primal_dual_solve(func, x0, Dfunc, A, b, G, h, args=(), tol=1e-10,
                      maxiter=100, show_progress=True, full_output=False):
    """Wrap calls to function and Jacobian"""
    fcalls, func = wrap_function(func, args)
    Dfcalls, Dfunc = wrap_function(Dfunc, args)

    M, N = G.shape
    P, N = A.shape

    """Total number of inequality constraints"""
    m = M

    def gap(z):
        r"""Gap-function"""
        l = z[N+P:N+P+M]
        s = z[N+P+M:]
        return mydot(l, s)/m

    def centrality(z):
        r"""Centrality function"""
        l = z[N+P:N+P+M]
        s = z[N+P+M:]
        return np.min(l*s)

    def KKT(z, sigma=0.0):
        r"""KKT system (possible perturbed)."""

        """Primal variable"""
        x = z[0:N]

        """Multiplier for equality constraints"""
        nu = z[N:N+P]

        """Multiplier for inequality constraints"""
        l = z[N+P:N+P+M]

        """Slacks"""
        s = z[N+P+M:]

        """Evaluate objective function"""
        F = func(x)

        """Dual infeasibility"""
        rdual = F+mydot(A.transpose(), nu)+mydot(G.transpose(), l)

        """Primal infeasibilities"""
        rprim1 = mydot(A, x)-b
        rprim2 = mydot(G, x)-h+s

        """Complementary slackness (perturbed)"""
        mu = gap(z)
        rslack = l*s-sigma*mu

        return np.hstack((rdual, rprim1, rprim2, rslack))

    def step_fast(z, KKTval, LU, G, A, mu, beta, gamma, alpha0):
        r"""Affine scaling step."""
        dz = solve_factorized(z, KKTval, LU, G, A)

        """Reduce step length until slacks s and multipliers l are positive"""
        alpha = 1.0*alpha_0
        k = 0
        for k in range(10):
            z_new = z + alpha*dz
            if np.all( z_new[N+P:] > 0.0 ):
                break
            alpha *= 0.5
            k += 1
        if k == 10 - 1:
            raise RuntimeError("Maximum steplength reduction reached")

        """Reduce step length until iterates lie in correct neighborhood"""
        for k in range(10):
            z_new = z + alpha*dz
            KKTval_new = KKT(z_new)
            dual = mynorm(KKTval_new[0:N])
            prim = mynorm(KKTval_new[N:N+P+M])
            mu_new = gap(z_new)
            cent_new = centrality(z_new)

            if (dual <= beta*mu_new and prim <= beta*mu_new and
                cent_new >= gamma*mu_new):
                break
            alpha *= 0.5
            # alpha *= 0.95
        if k == 10 - 1:
            raise RuntimeError("Maximum steplength reduction reached")
        return z_new, mu_new

    def step_safe(z, KKTval, LU, G, A, mu, beta, gamma):
        r"""Centering step."""
        dz = solve_factorized(z, KKTval, LU, G, A)

        """Reduce step length until slacks s and multipliers l are positive"""
        alpha = 1.0
        k = 0
        for k in range(10):
            z_new = z + alpha*dz
            if np.all( z_new[N+P:] > 0.0 ):
                break
            alpha *= 0.5
            k += 1
        if k == 10 - 1:
            raise RuntimeError("Maximum steplength reduction (pos.) reached")

        """Reduce step length until iterates lie in correct neighborhood
        and mu fulfills Armijo condition"""
        k = 0
        for k in range(10):
            z_new = z+alpha*dz
            KKTval_new = KKT(z_new, sigma=SIGMA)
            dual = mynorm(KKTval_new[0:N])
            prim = mynorm(KKTval_new[N:N+P+M])
            mu_new = gap(z_new)
            cent_new = centrality(z_new)

            if (dual <= beta*mu_new and prim <= beta*mu_new and
                cent_new >= gamma*mu_new and
                mu_new<=(1.0-KAPPA*alpha*(1.0-SIGMA))*mu):
                break
            alpha *= 0.5
        if k == 10 - 1:
            raise RuntimeError("Maximum steplength reduction reached")
        return z_new, mu_new

    """INITIALIZATION"""

    """Initial Slacks for inequality constraints"""
    s0 = -1.0*(mydot(G, x0)-h)

    """Initial multipliers for inequality constraints"""
    l0 = 1.0*np.ones(M)

    """Initial multipliers for equality constraints"""
    nu0 = np.zeros(P)

    """Initial point"""
    z0 = np.hstack((x0, nu0, l0, s0))

    """Initial KKT-values"""
    KKTval0 = KKT(z0, sigma=0.0)
    mu0 = gap(z0)
    dual0 = mynorm(KKTval0[0:N])
    prim0 = mynorm(KKTval0[N:N+P+M])

    """Initial neighborhood"""
    beta = BETA * np.sqrt(dual0**2 + prim0**2)/mu0
    gamma = 1.0 * GAMMA_MAX

    """Number of fast steps"""
    t = 0

    """Number of iterations"""
    n = 0

    """Dummy variable for step type"""
    step_type = " "

    if show_progress:
        print("%s %s %s %s %s %s %s" %("iter", "gap", "dual", "primal",
                                       "min", "max", "step"))
    """MAIN LOOP"""
    z = z0
    x = z0[0:N]
    KKTval = KKTval0
    dual = dual0
    prim = prim0
    mu = mu0
    Dfunc_val = Dfunc(x)
    LU = factor(z, Dfunc_val, G, A)

    if full_output:
        info = {'z': []}
        info['z'].append(z)

    for n in range(maxiter):
        if show_progress:
            l=z[N+P:N+P+M]
            s=z[N+P+M:]
            print("%i %.6e %.6e %.6e %.6e %.6e %s" %(n+1, mu, dual, prim,
                                                     np.min(l*s), np.max(l*s),
                                                     step_type))
        """Attempt fast step"""
        beta_new = (1.0 + GAMMA_BAR**(t+1)) * beta
        gamma_new = GAMMA_MIN + GAMMA_BAR**(t+1)*(GAMMA_MAX - GAMMA_MIN)
        alpha_0 = 1.0 - np.sqrt(mu)/GAMMA_BAR**t

        if alpha_0 > 0.0:
            z_new, mu_new = step_fast(z, KKTval, LU, G, A, mu,
                                      beta_new, gamma_new, alpha_0)
            if mu_new < RHO * mu:
                """Fast successful"""
                z = z_new
                mu = mu_new
                beta = beta_new
                gamma = gamma_new
                t += 1
                step_type = "f"
            else:
                """Perturbed right-had side"""
                KKTval_pert = 1.0*KKTval
                KKTval_pert[N+P+M:] -= SIGMA * mu
                z, mu = step_safe(z, KKTval_pert, LU, G, A, mu, beta, gamma)
                step_type = "s"
        else:
            """Perturbed right-hand side"""
            KKTval_pert = 1.0*KKTval
            KKTval_pert[N+P+M:] -= SIGMA * mu
            z, mu = step_safe(z, KKTval_pert, LU, G, A, mu, beta, gamma)
            step_type = "s"

        """Compute new iterates"""
        KKTval = KKT(z, sigma=0.0)
        dual = mynorm(KKTval[0:N])
        prim = mynorm(KKTval[N:N+P+M])
        x = z[0:N]
        Dfunc_val = Dfunc(x)
        LU = factor(z, Dfunc_val, G, A)
        if full_output:
            info['z'].append(z)
        if mu < tol and dual < tol and prim < tol:
            break
    if n == maxiter - 1:
        raise RuntimeError("Maximum number of iterations reached")

    if show_progress:
        l=z[N+P:N+P+M]
        s=z[N+P+M:]
        print("%i %.6e %.6e %.6e %.6e %.6e %s" %(n+1, mu, dual, prim,
                                                 np.min(l*s), np.max(l*s),
                                                 step_type))
    if full_output:
        return z[0:N], info
    else:
        return z[0:N]
def primal_dual_solve ( func , x0 , Dfunc , A , b , G , h , args = ( ) , tol = 1e-10 , maxiter = 100 , show_progress = True , full_output = False ) : fcalls , func = wrap_function ( func , args ) Dfcalls , Dfunc = wrap_function ( Dfunc , args ) M , N = G . shape P , N = A . shape """Total number of inequality constraints""" m = M def gap ( z ) : l = z [ N + P : N + P + M ] s = z [ N + P + M : ] return mydot ( l , s ) / m def centrality ( z ) : l = z [ N + P : N + P + M ] s = z [ N + P + M : ] return np . min ( l * s ) def KKT ( z , sigma = 0.0 ) : """Primal variable""" x = z [ 0 : N ] """Multiplier for equality constraints""" nu = z [ N : N + P ] """Multiplier for inequality constraints""" l = z [ N + P : N + P + M ] """Slacks""" s = z [ N + P + M : ] """Evaluate objective function""" F = func ( x ) """Dual infeasibility""" rdual = F + mydot ( A . transpose ( ) , nu ) + mydot ( G . transpose ( ) , l ) """Primal infeasibilities""" rprim1 = mydot ( A , x ) - b rprim2 = mydot ( G , x ) - h + s """Complementary slackness (perturbed)""" mu = gap ( z ) rslack = l * s - sigma * mu return np . hstack ( ( rdual , rprim1 , rprim2 , rslack ) ) def step_fast ( z , KKTval , LU , G , A , mu , beta , gamma , alpha0 ) : dz = solve_factorized ( z , KKTval , LU , G , A ) """Reduce step length until slacks s and multipliers l are positive""" alpha = 1.0 * alpha_0 k = 0 for k in range ( 10 ) : z_new = z + alpha * dz if np . all ( z_new [ N + P : ] > 0.0 ) : break alpha *= 0.5 k += 1 if k == 10 - 1 : raise RuntimeError ( "Maximum steplength reduction reached" ) """Reduce step length until iterates lie in correct neighborhood""" for k in range ( 10 ) : z_new = z + alpha * dz KKTval_new = KKT ( z_new ) dual = mynorm ( KKTval_new [ 0 : N ] ) prim = mynorm ( KKTval_new [ N : N + P + M ] ) mu_new = gap ( z_new ) cent_new = centrality ( z_new ) if ( dual <= beta * mu_new and prim <= beta * mu_new and cent_new >= gamma * mu_new ) : break alpha *= 0.5 if k == 10 - 1 : raise RuntimeError ( "Maximum steplength reduction reached" ) return z_new , mu_new def step_safe ( z , KKTval , LU , G , A , mu , beta , gamma ) : dz = solve_factorized ( z , KKTval , LU , G , A ) """Reduce step length until slacks s and multipliers l are positive""" alpha = 1.0 k = 0 for k in range ( 10 ) : z_new = z + alpha * dz if np . all ( z_new [ N + P : ] > 0.0 ) : break alpha *= 0.5 k += 1 if k == 10 - 1 : raise RuntimeError ( "Maximum steplength reduction (pos.) reached" ) """Reduce step length until iterates lie in correct neighborhood
        and mu fulfills Armijo condition""" k = 0 for k in range ( 10 ) : z_new = z + alpha * dz KKTval_new = KKT ( z_new , sigma = SIGMA ) dual = mynorm ( KKTval_new [ 0 : N ] ) prim = mynorm ( KKTval_new [ N : N + P + M ] ) mu_new = gap ( z_new ) cent_new = centrality ( z_new ) if ( dual <= beta * mu_new and prim <= beta * mu_new and cent_new >= gamma * mu_new and mu_new <= ( 1.0 - KAPPA * alpha * ( 1.0 - SIGMA ) ) * mu ) : break alpha *= 0.5 if k == 10 - 1 : raise RuntimeError ( "Maximum steplength reduction reached" ) return z_new , mu_new """INITIALIZATION""" """Initial Slacks for inequality constraints""" s0 = - 1.0 * ( mydot ( G , x0 ) - h ) """Initial multipliers for inequality constraints""" l0 = 1.0 * np . ones ( M ) """Initial multipliers for equality constraints""" nu0 = np . zeros ( P ) """Initial point""" z0 = np . hstack ( ( x0 , nu0 , l0 , s0 ) ) """Initial KKT-values""" KKTval0 = KKT ( z0 , sigma = 0.0 ) mu0 = gap ( z0 ) dual0 = mynorm ( KKTval0 [ 0 : N ] ) prim0 = mynorm ( KKTval0 [ N : N + P + M ] ) """Initial neighborhood""" beta = BETA * np . sqrt ( dual0 ** 2 + prim0 ** 2 ) / mu0 gamma = 1.0 * GAMMA_MAX """Number of fast steps""" t = 0 """Number of iterations""" n = 0 """Dummy variable for step type""" step_type = " " if show_progress : print ( "%s %s %s %s %s %s %s" % ( "iter" , "gap" , "dual" , "primal" , "min" , "max" , "step" ) ) """MAIN LOOP""" z = z0 x = z0 [ 0 : N ] KKTval = KKTval0 dual = dual0 prim = prim0 mu = mu0 Dfunc_val = Dfunc ( x ) LU = factor ( z , Dfunc_val , G , A ) if full_output : info = { 'z' : [ ] } info [ 'z' ] . append ( z ) for n in range ( maxiter ) : if show_progress : l = z [ N + P : N + P + M ] s = z [ N + P + M : ] print ( "%i %.6e %.6e %.6e %.6e %.6e %s" % ( n + 1 , mu , dual , prim , np . min ( l * s ) , np . max ( l * s ) , step_type ) ) """Attempt fast step""" beta_new = ( 1.0 + GAMMA_BAR ** ( t + 1 ) ) * beta gamma_new = GAMMA_MIN + GAMMA_BAR ** ( t + 1 ) * ( GAMMA_MAX - GAMMA_MIN ) alpha_0 = 1.0 - np . sqrt ( mu ) / GAMMA_BAR ** t if alpha_0 > 0.0 : z_new , mu_new = step_fast ( z , KKTval , LU , G , A , mu , beta_new , gamma_new , alpha_0 ) if mu_new < RHO * mu : z = z_new mu = mu_new beta = beta_new gamma = gamma_new t += 1 step_type = "f" else : KKTval_pert = 1.0 * KKTval KKTval_pert [ N + P + M : ] -= SIGMA * mu z , mu = step_safe ( z , KKTval_pert , LU , G , A , mu , beta , gamma ) step_type = "s" else : KKTval_pert = 1.0 * KKTval KKTval_pert [ N + P + M : ] -= SIGMA * mu z , mu = step_safe ( z , KKTval_pert , LU , G , A , mu , beta , gamma ) step_type = "s" """Compute new iterates""" KKTval = KKT ( z , sigma = 0.0 ) dual = mynorm ( KKTval [ 0 : N ] ) prim = mynorm ( KKTval [ N : N + P + M ] ) x = z [ 0 : N ] Dfunc_val = Dfunc ( x ) LU = factor ( z , Dfunc_val , G , A ) if full_output : info [ 'z' ] . append ( z ) if mu < tol and dual < tol and prim < tol : break if n == maxiter - 1 : raise RuntimeError ( "Maximum number of iterations reached" ) if show_progress : l = z [ N + P : N + P + M ] s = z [ N + P + M : ] print ( "%i %.6e %.6e %.6e %.6e %.6e %s" % ( n + 1 , mu , dual , prim , np . min ( l * s ) , np . max ( l * s ) , step_type ) ) if full_output : return z [ 0 : N ] , info else : return z [ 0 : N ]
Wrap calls to function and Jacobian