Using external function in nlpsol codegen

[LinusTxtonomy]

Hi,
I am building an optimization problem that should use external code as part of the objective function. I have built a small example that reproduces the error:

from casadi import *

x = MX.sym('x', 2)
slack = MX.sym('slack')

radius = MX.sym('radius')

J = x[0] + x[1]
obj_function = Function('objective', [x], [J])
obj_function.generate('J.c')
external_obj_function = external('objective', 'J.so', {'enable_fd': True})
J = external_obj_function(x)

g = x[0]**2 + x[1]**2 - radius + slack

prob = {'f': J, 'x': vertcat(x, slack), 'g': g, 'p': radius}

solver = nlpsol('solver', 'ipopt', prob, {'ipopt.print_level': 5})

solver.generate_dependencies('nlp.c')

res = solver(lbg = 0, ubg = 0, lbx = [-inf, -inf, 0], ubx = [inf, inf, inf], p = [4])
print(res)

The objective function code generates and builds without issues but I get compiler errors when trying to compile the generated nlp.c file even though the optimization problem runs fine inside the python example:

nlp.c: In function ‘nlp_smoothing_diff_err’:
nlp.c:211:55: error: ‘static_cast’ undeclared (first use in this function)
  211 |     wk = casadi_smoothing_diff_weights(k, yb, yc, yf, static_cast<casadi_real*>(0));
      |                                                       ^~~~~~~~~~~
nlp.c:211:55: note: each undeclared identifier is reported only once for each function it appears in
nlp.c:24:21: error: expected expression before ‘double’
   24 | #define casadi_real double
      |                     ^~~~~~
nlp.c:211:67: note: in expansion of macro ‘casadi_real’
  211 |     wk = casadi_smoothing_diff_weights(k, yb, yc, yf, static_cast<casadi_real*>(0));
      |                                                                   ^~~~~~~~~~~
nlp.c: In function ‘nlp_f6’:
nlp.c:812:39: error: expected ‘)’ before ‘;’ token
  812 |         if (mid = casadi_f1_checkout();
      |            ~                          ^
      |                                       )

I could also manually code the derivatives but I am not sure what functions are needed to do this. In my project I need to have part of the objective function as an external function. I assume I need atleast the jacobian and hessian but not sure about the naming conventions for that.
Regards
Linus

[ernestds]

I don't know exactly what's up with the first errors, but the last one has happened to me and is reproducible as such (3.6.5):

import casadi as ca
import os
import pathlib
import numpy as np
TEMP_FOLDER = pathlib.Path(os.getcwd()) / 'temp'

codegen_folder = TEMP_FOLDER / 'codegen'
codegen_folder.mkdir(exist_ok=True,parents=True)

n = 40
map_size = 20
A = ca.SX.sym('A',n,n)
B = ca.SX.sym('B',n,2)
C = ca.sumsqr(ca.cos(A)@B)

fun = ca.Function('fun_inner',[A,B],[C])
cg = ca.CodeGenerator('fun_inner')
cg.add(fun)
cg.generate(str(codegen_folder) + "/")

os.system(f"""gcc -O3 -march=native -shared {str(codegen_folder / 'fun_inner.c')} -o {str(codegen_folder / ('fun_inner.so'))}""")


fun_inner_compiled = ca.external('fun_inner',str(codegen_folder / ("fun_inner.so")))

A_map = ca.MX.sym('A',n,n*map_size)
B_map = ca.MX.sym('B',n,2*map_size)
# asdfasdffsd
fun_map_MX = ca.Function('fun_inner_mapped_MX',[A_map,B_map],[fun_inner_compiled.map(map_size)(A_map,B_map)])
# vs
fun_map_MX = ca.Function('fun_inner_mapped_MX',[A_map,B_map],[fun_inner_compiled(A_map,B_map)])

cg = ca.CodeGenerator('fun_inner_mapped_MX')
cg.add(fun_map_MX)
cg.generate(str(codegen_folder) + "/")
"""
# bug, code is generated wrong if .map(map_size) is used:

/workspaces/misc/temp/codegen/fun_inner_mapped_MX.c: In function ‘fun_inner_mapped_MX_f1’:
/workspaces/misc/temp/codegen/fun_inner_mapped_MX.c:126:35: error: expected ‘)’ before ‘;’ token
  126 |     if (mid = casadi_f2_checkout();
      |        ~                          ^
      |                                   )
"""
# print size of generated file
print(f"Size of generated file: {os.path.getsize(codegen_folder / 'fun_inner_mapped_MX.c')/1024:.1f} KB")

os.system(f"""gcc -O3 -march=native -c {str(codegen_folder / ('fun_inner_mapped_MX.c'))} -o {str(codegen_folder / ("fun_inner_mapped_MX.o"))}""")
os.system(f"""gcc -shared -o {str(codegen_folder / ("fun_inner_mapped_MX.so"))} {str(codegen_folder / ('fun_inner_mapped_MX.o'))}  {str(codegen_folder / ('fun_inner.so'))}""")
    # os.system(f'gcc -O3 -c {str(codegen_folder / (this_name + ".c"))} -o {str(codegen_folder / (this_name + ".o"))}')
    # os.system(f'gcc -shared -o {str(codegen_folder / (this_name + ".so"))} {str(codegen_folder / (this_name + ".o"))} -lm {str(codegen_folder / (diff_co_single_name + ".so"))}')

Comment one of the fun_map_MX to check it out. I've had this problem in the above case and also in a similar case as OP's.

Quick fix hack is to use this (works in my code, IDK if there's other edge cases):

def fix_code(path):
    import re
    # Define the pattern to find and the replacement pattern
    pattern = re.compile(r'if \((mid = (.*?))flag\) return 1;', re.MULTILINE | re.DOTALL)
    # pattern = re.compile(r'if \((mid(.*?))flag\) return 1;', re.MULTILINE | re.DOTALL)
    replacement = r'\1if (flag) return 1;'

    # Read the file
    with open(path, 'r') as file:
        file_contents = file.read()
    
    # Apply the correction
    corrected_contents = re.sub(pattern, replacement, file_contents)
    with open(path, 'w') as file:
        file.write(corrected_contents)
fix_code(nlp.c)

Also don't forget to link the external function to the final nlp.so.

In any case:

from casadi import *

x = MX.sym('x', 2)
slack = MX.sym('slack')

radius = MX.sym('radius')

J = x[0] + x[1]
obj_function = Function('objective', [x], [J])
cg = CodeGenerator('objective')
cg.add(obj_function)
cg.add(obj_function.jacobian())
cg.add(obj_function.jacobian().jacobian())
cg.generate()

os.system(f"""gcc -O1 -march=native -shared objective.c -o J.so""")
external_obj_function = external('objective', 'J.so', )
J = external_obj_function(x)

g = x[0]**2 + x[1]**2 - radius + slack

prob = {'f': J, 'x': vertcat(x, slack), 'g': g, 'p': radius}

solver = nlpsol('solver', 'ipopt', prob, {'ipopt.print_level': 5})

solver.generate_dependencies('nlp.c')

res = solver(lbg = 0, ubg = 0, lbx = [-inf, -inf, 0], ubx = [inf, inf, inf], p = [4])
print(res)
def fix_code(path):
    import re
    # Define the pattern to find and the replacement pattern
    pattern = re.compile(r'if \((mid = (.*?))flag\) return 1;', re.MULTILINE | re.DOTALL)
    # pattern = re.compile(r'if \((mid(.*?))flag\) return 1;', re.MULTILINE | re.DOTALL)
    replacement = r'\1if (flag) return 1;'

    # Read the file
    with open(path, 'r') as file:
        file_contents = file.read()
    
    # Apply the correction
    corrected_contents = re.sub(pattern, replacement, file_contents)
    with open(path, 'w') as file:
        file.write(corrected_contents)
        
        
fix_code('nlp.c')
os.system(f"""g++ -O1 -march=native -c nlp.c -o nlp.o""")
os.system(f"""g++ -shared -o nlp.so nlp.o  J.so""")

file = 'nlp.so'
oracle = external('nlp', file)

nlp_f = external('nlp_f', file,)
nlp_g = external('nlp_g', file,)
nlp_grad = external('nlp_grad', file,)
nlp_grad_f = external('nlp_grad_f', file,)
nlp_hess_l = external('nlp_hess_l', file,)
nlp_jac_g = external('nlp_jac_g', file,)
cache={'nlp_f':nlp_f,
       'nlp_grad':nlp_grad,
                                     'nlp_g':nlp_g,
                                    'nlp_grad_f':nlp_grad_f,
                                    'nlp_hess_l':nlp_hess_l,
                                    'nlp_jac_g':nlp_jac_g}
# cache = {}
oracle_in = oracle.mx_in()
oracle_out = oracle.call(oracle_in)
oracle_wrap = dict(zip(oracle.name_in() + oracle.name_out(), oracle_in + oracle_out))
solver = nlpsol('nlp', 'ipopt', oracle_wrap, dict(cache = cache))


res = solver(lbg = 0, ubg = 0, lbx = [-inf, -inf, 0], ubx = [inf, inf, inf], p = [4])
print(res)

Here's your code working + defining and loading external derivative information. Import thing is to name properly (jac_functionname for jacobian, jac_jac_functionname for the jacobian of the jacobian) and have the correct inputs and outputs. I don't remember if there's an easier way to load the NLP as an external function, but this works.

[ernestds]

Oh I also changed from gcc to g++, I don't know much about C but I guess the stuff in the other errors are C++ only?

[LinusTxtonomy]

Thanks a lot that worked. Do you know what part of the derivatives are actually used/needed by nlpsol? For example if I want to calculate the distance to a BSpline, I need to pass in all BSpline knots and x,y-coordinates. If I then use the created jacobians, casadi expects the jacobian of the distance with respect to the BSpline parameters. Is that actually required or can it be left out?

[ernestds]

You would need the jacobian and the hessian of the output of the function wrt. to the input decision variables. The Casadi jacobians have a specific format, you can call a function.jacobian() and see what it shows; they are always going to follow this format. If your function has as an input the BSpline knots, the BSpline control points and the x,y coordinates, then your jacobian would be a jacobian(knots,control_points,coord, some other inputs that are automatically created) -> dD/dknots, dD/dcontrol_points, dD/coord.
If your knots and coordinates are fixed, for example, then you might only need to specify dD/dcontrol_points because the other derivatives won't be used. You do need to specify the other ones, if I'm not mistaken, but you can just return a zero matrix of the appropriate size.
Also take a look at the "is_diff_in" : [True, False,...] option.