Below there is a sample code where the BSplineComp is combined either with an ExplicitComp or ExternalCodeComp. Both of these two do the same calculation and both of the components' gradients are calculated using finite difference.
If I run the version Bspline+ExplicitComp the result is achieved within 2,3 iterations. If I run the version Bspline+ExternalCodeComp I have to wait a lot. In this case it is trying to find the gradient of the output with respect to each input. So for example there are 9 control points that are interpolated to 70 points in the bspline component. Then the external component has to be evaluated as many as the interpolated points (70 times)
So in a case where the bspline is combined with an expensive external code the finite difference requires as much as the number of points it is interpolated to which becomes the bottleneck of the computation.
Based on this input I have two questions
1- If external code component is based on the explicit component what is the major difference that causes the behaviour difference? (considering both have an input of shape=70)
2- In the previously mentioned scneario where the bspline is combined with an expensive external code would there be a more efficient way of combining them apart from the way it is shown here.
MAIN CODE: 'external' variable is the flag for toggling external/explicit code comp. set that true/false for running the two cases explained above.
from openmdao.components.bsplines_comp import BsplinesComp
from openmdao.api import IndepVarComp, Problem, ExplicitComponent,ExecComp,ExternalCodeComp
from openmdao.api import ScipyOptimizeDriver, SqliteRecorder, CaseReader
import matplotlib.pyplot as plt
import numpy as np
external=True # change this to true for the case with external code comp. or false for the case with explicit comp.
rr=np.arange(0,70,1)
"Explicit component for the area under the line calculation"
class AreaComp(ExplicitComponent):
def initialize(self):
self.options.declare('lenrr', int)
self.options.declare('rr', types=np.ndarray)
def setup(self):
self.add_input('h', shape=lenrr)
self.add_output('area')
self.declare_partials(of='area', wrt='h', method='fd')
def compute(self, inputs, outputs):
rr = self.options['rr']
outputs['area'] = np.trapz(rr,inputs['h'])
class ExternalAreaComp(ExternalCodeComp):
def setup(self):
self.add_input('h', shape=70)
self.add_output('area')
self.input_file = 'paraboloid_input.dat'
self.output_file = 'paraboloid_output.dat'
# providing these is optional; the component will verify that any input
# files exist before execution and that the output files exist after.
self.options['external_input_files'] = [self.input_file]
self.options['external_output_files'] = [self.output_file]
self.options['command'] = [
'python', 'extcode_paraboloid.py', self.input_file, self.output_file
]
# this external code does not provide derivatives, use finite difference
self.declare_partials(of='*', wrt='*', method='fd')
def compute(self, inputs, outputs):
h = inputs['h']
# generate the input file for the paraboloid external code
np.savetxt(self.input_file,h)
# the parent compute function actually runs the external code
super(ExternalAreaComp, self).compute(inputs, outputs)
# parse the output file from the external code and set the value of f_xy
f_xy=np.load('a.npy')
outputs['area'] = f_xy
prob = Problem()
model = prob.model
n_cp = 9
lenrr = len(rr)
"Initialize the design variables"
x = np.random.rand(n_cp)
model.add_subsystem('px', IndepVarComp('x', val=x))
model.add_subsystem('interp', BsplinesComp(num_control_points=n_cp,
num_points=lenrr,
in_name='h_cp',
out_name='h'))
if external:
comp=ExternalAreaComp()
model.add_subsystem('AreaComp', comp)
else:
comp = AreaComp(lenrr=lenrr, rr=rr)
model.add_subsystem('AreaComp', comp)
case_recorder_filename2 = 'cases4.sql'
recorder2 = SqliteRecorder(case_recorder_filename2)
comp.add_recorder(recorder2)
comp.recording_options['record_outputs']=True
comp.recording_options['record_inputs']=True
model.connect('px.x', 'interp.h_cp')
model.connect('interp.h', 'AreaComp.h')
model.add_constraint('interp.h', lower=0.9, upper=1, indices=[0])
prob.driver = ScipyOptimizeDriver()
prob.driver.options['optimizer'] = 'SLSQP'
prob.driver.options['disp'] = True
#prob.driver.options['optimizer'] = 'COBYLA'
#prob.driver.options['disp'] = True
prob.driver.options['tol'] = 1e-9
model.add_design_var('px.x', lower=1,upper=10)
model.add_objective('AreaComp.area',scaler=1)
prob.setup(check=True)
#prob.run_model()
prob.run_driver()
cr = CaseReader(case_recorder_filename2)
case_keys = cr.system_cases.list_cases()
cou=-1
for case_key in case_keys:
cou=cou+1
case = cr.system_cases.get_case(case_key)
plt.plot(rr,case.inputs['h'],'-*')
The external code extcode_paraboloid.py is below
import numpy as np
if __name__ == '__main__':
import sys
input_filename = sys.argv[1]
output_filename = sys.argv[2]
h=np.loadtxt(input_filename)
rr=np.arange(0,70,1)
rk= np.trapz(rr,h)
np.save('a',np.array(rk))
