clnlbeam.jl

using Calculus
# make derivative-rules Dict with Function keys and anonymous-function output
deriv_dict = Dict{Function, Function}()
for (funsym, exp) in Calculus.derivative_rules
    @eval deriv_dict[$(funsym)] =
    let xp = 1
        x -> $exp
    end
end

function main(n)
h = 1/n
# nonlinear constraint functions f(x) are encoded
# in two sparse matrices, P and K, as follows:
# f(x) = K * prod(V,2)
# V[j,k] = P[j,k](x[k]) if P[j,k] is a function
# V[j,k] = x[k]^P[j,k]  otherwise
Prows = 1:4n+4
Pcols = [n+2:2n+2; 1:n+1; 1:n+1; 2n+3:3n+3]
Pvals = Array(Any, 4n+4)
Pvals[1:n+1] = 1
Pvals[n+2:2n+2] = sin
Pvals[2n+3:4n+4] = 1
P = sparse(Prows, Pcols, Pvals)
Pt = P.' # save transpose since row-major is more useful here
Krows = [1:n;       1:n;       1:n;       1:n;
         n+1:2n;    n+1:2n;    n+1:2n;    n+1:2n]
Kcols = [1:n;       2:n+1;     n+2:2n+1;  n+3:2n+2;
         2n+3:3n+2; 2n+4:3n+3; 3n+4:4n+3; 3n+5:4n+4]
Kvals = [-ones(n);  ones(n);   -0.5h*ones(2n);
         -ones(n);  ones(n);   -0.5h*ones(2n)]
K = sparse(Krows, Kcols, Kvals)

# evaluation operator for an element of P
function apply_op(x, op::Integer)
    if op >= 0
        return Base.power_by_squaring(x, op)
    else
        return Base.power_by_squaring(1.0/x, -op)
    end
end
apply_op(x, op::Function) = op(x)
apply_op(x, op) = x^op

# derivative operator for an element of P
function apply_der(x, op::Integer)
    if op >= 0
        return op*Base.power_by_squaring(x, op-1)
    else
        return Base.power_by_squaring(1.0/x, 1-op)
    end
end
function apply_der(x, op::Function)
    if haskey(deriv_dict, op)
        return deriv_dict[op](x)
    else
        return derivative(op)(x)
    end
end
function apply_der(x, op)
    if op == one(op) || op == zero(op)
        return op
    else
        return op*(x^(op-1))
    end
end

# for compatibility with Julia from earlier than commit
# c93ed16a093e33fe9e0de153120c821959641bd2
!isdefined(:nfilled) && (nfilled = nnz)

# evaluation of nonlinear function at vector x, given P' and K
function eval_fcn(x, Pt::SparseMatrixCSC, K::SparseMatrixCSC)
    prods = ones(typeof(x[1]), size(Pt,2))
    for i=1:size(Pt,2), j=Pt.colptr[i]:Pt.colptr[i+1]-1
        prods[i] *= apply_op(x[Pt.rowval[j]], Pt.nzval[j])
    end
    return K * prods
end

# evaluation of jacobian at vector x, given P' and K
function eval_jac(x, Pt::SparseMatrixCSC, K::SparseMatrixCSC)
    prodjacvals = ones(Float64, nfilled(Pt))
    for i=1:size(Pt,2), j=Pt.colptr[i]:Pt.colptr[i+1]-1
        # derivative of product term wrt this variable
        prodjacvals[j] *= apply_der(x[Pt.rowval[j]], Pt.nzval[j])
        for k = Pt.colptr[i]:j-1
            # product term as function of other variables
            prodjacvals[k] *= apply_op(x[Pt.rowval[j]], Pt.nzval[j])
        end
        for k = j+1:Pt.colptr[i+1]-1
            prodjacvals[k] *= apply_op(x[Pt.rowval[j]], Pt.nzval[j])
        end
    end
    return K * SparseMatrixCSC(Pt.m, Pt.n, Pt.colptr, Pt.rowval, prodjacvals).'
end

xv = randn(3n+3)
J = eval_jac(xv, Pt, K)

@time for i=1:100; J = eval_jac(xv, Pt, K); end
println("for 100 Jacobian evaluations")

#Profile.print()

end