cov.py

import casadi as ca
import pylab as pl

# [...]

# A is of type SX, column vector, shape (788, 1), dense
# B is of type sx, column vector, shape (984, 1), dense

# V is of type struct_symSX, shape (1594, 1), dense
# W is a numpy array, matrix with shape (788, 788), almost diagonal

J1 = ca.mul(W, ca.jacobian(A, V))
J2 = ca.jacobian(B, V)

h = ca.blockcat([[ca.mul([J1.T, W, J1]), J2.T], [J2, ca.DMatrix(B.size, B.size)]])
rhs = ca.vertcat((J1.T, ca.DMatrix(B.size, A.size)))
Jp = ca.solve(h, rhs)[:V.size, :] # Note: first V.size rows

fcov = ca.SXFunction([V], [ca.mul(Jp, Jp.T)])
fcov.init()

# [...] followed by an evaluation of fcov