How to use the dymos.utils.lagrange.lagrange_matrices function in dymos

#
        # 1. Get the segment tau values of the given nodes.
        #
        i1, i2 = grid_data.subset_segment_indices['control_disc'][iseg, :]
        indices = grid_data.subset_node_indices['control_disc'][i1:i2]
        tau_s_given = grid_data.node_stau[indices]

        #
        # 2. Get the segment tau values of the evaluation nodes.
        #
        t_eval_iseg = t_eval_per_seg[iseg]
        t0_iseg, tf_iseg = time_seg_ends[iseg, :]
        tau_s_eval = 2.0 * (t_eval_iseg - t0_iseg) / (tf_iseg - t0_iseg) - 1

        L_block, D_block = lagrange_matrices(tau_s_given, tau_s_eval)
        _, Daa_block = lagrange_matrices(tau_s_given, tau_s_given)

        L_blocks.append(L_block)
        D_blocks.append(D_block)
        Daa_blocks.append(Daa_block)

    L_ae = block_diag(*L_blocks)
    D_ae = block_diag(*D_blocks)

    D_aa = block_diag(*Daa_blocks)
    D2_ae = np.dot(D_ae, D_aa)

    return L_ae, D_ae, D2_ae

'belongs')

        for name, options in self.options['polynomial_control_options'].items():
            disc_nodes, _ = lgl(options['order'] + 1)
            num_control_input_nodes = len(disc_nodes)
            shape = options['shape']
            size = np.prod(shape)
            units = options['units']
            rate_units = get_rate_units(units, self.options['time_units'], deriv=1)
            rate2_units = get_rate_units(units, self.options['time_units'], deriv=2)

            input_shape = (num_control_input_nodes,) + shape
            output_shape = (num_nodes,) + shape

            L_de, D_de = lagrange_matrices(disc_nodes, eval_nodes)
            _, D_dd = lagrange_matrices(disc_nodes, disc_nodes)
            D2_de = np.dot(D_de, D_dd)

            self._matrices[name] = L_de, D_de, D2_de

            self._input_names[name] = 'polynomial_controls:{0}'.format(name)
            self._output_val_names[name] = 'polynomial_control_values:{0}'.format(name)
            self._output_rate_names[name] = 'polynomial_control_rates:{0}_rate'.format(name)
            self._output_rate2_names[name] = 'polynomial_control_rates:{0}_rate2'.format(name)

            self.add_input(self._input_names[name], val=np.ones(input_shape), units=units)
            self.add_output(self._output_val_names[name], shape=output_shape, units=units)
            self.add_output(self._output_rate_names[name], shape=output_shape, units=rate_units)
            self.add_output(self._output_rate2_names[name], shape=output_shape, units=rate2_units)

            self.val_jacs[name] = np.zeros((num_nodes, size, num_control_input_nodes, size))
            self.rate_jacs[name] = np.zeros((num_nodes, size, num_control_input_nodes, size))

desc='Interpolated value of dynamic '
                         'entry {0} at timesteps'.format(dynamic_name))

        L_blocks = []
        segment_times = self.options['segment_times']
        for iseg in range(grid_data.num_segments):
            i1, i2 = grid_data.subset_segment_indices['all'][iseg, :]
            indices = grid_data.subset_node_indices['all'][i1:i2]
            nodes_given = grid_data.node_stau[indices]

            ts_start, ts_end = segment_times[iseg]
            segment_tsteps = normalized_times[ts_start:ts_end]
            timesteps_stau = (2.0 * segment_tsteps - (segment_tsteps[-1] + segment_tsteps[0])) \
                / (segment_tsteps[-1] - segment_tsteps[0])

            L_block, _ = lagrange_matrices(nodes_given, timesteps_stau)

            L_blocks.append(L_block)

        self._L = block_diag(*L_blocks)

        # Setup partials

        self.jacs = {'L': {}}
        self.matrices = {'L': self._L}

        for dynamic_name, opts in iteritems(control_options):
            if opts['dynamic']:
                shape = opts['shape']
                m = np.prod(shape)

                for key in self.jacs:

L: np.ndarray
        The lagrange interpolation matrix

    """
    L_blocks = []

    for iseg in range(old_grid.num_segments):
        i1, i2 = old_grid.subset_segment_indices['all'][iseg, :]
        indices = old_grid.subset_node_indices['all'][i1:i2]
        nodes_given = old_grid.node_stau[indices]

        i1, i2 = new_grid.subset_segment_indices['all'][iseg, :]
        indices = new_grid.subset_node_indices['all'][i1:i2]
        nodes_eval = new_grid.node_stau[indices]

        L_block, _ = lagrange_matrices(nodes_given, nodes_eval)

        L_blocks.append(L_block)

    L = block_diag(*L_blocks)

    return L

desc='duration of the phase to which this interpolated control group '
                            'belongs')

        for name, options in self.options['polynomial_control_options'].items():
            disc_nodes, _ = lgl(options['order'] + 1)
            num_control_input_nodes = len(disc_nodes)
            shape = options['shape']
            size = np.prod(shape)
            units = options['units']
            rate_units = get_rate_units(units, self.options['time_units'], deriv=1)
            rate2_units = get_rate_units(units, self.options['time_units'], deriv=2)

            input_shape = (num_control_input_nodes,) + shape
            output_shape = (num_output_nodes,) + shape

            L_do, D_do = lagrange_matrices(disc_nodes, output_nodes_ptau)
            _, D_dd = lagrange_matrices(disc_nodes, disc_nodes)
            D2_do = np.dot(D_do, D_dd)

            self._matrices[name] = L_do, D_do, D2_do

            self._input_names[name] = 'polynomial_controls:{0}'.format(name)
            self._output_val_names[name] = 'polynomial_control_values:{0}'.format(name)
            self._output_rate_names[name] = 'polynomial_control_rates:{0}_rate'.format(name)
            self._output_rate2_names[name] = 'polynomial_control_rates:{0}_rate2'.format(name)

            self.add_input(self._input_names[name], val=np.ones(input_shape), units=units)
            self.add_output(self._output_val_names[name], shape=output_shape, units=units)
            self.add_output(self._output_rate_names[name], shape=output_shape, units=rate_units)
            self.add_output(self._output_rate2_names[name], shape=output_shape, units=rate2_units)

            self.val_jacs[name] = np.zeros((num_output_nodes, size, num_control_input_nodes, size))

time_units = self.options['time_units']
        c = rk_methods[method]['c']

        stau_step = np.linspace(-1, 1, num_steps + 1)

        step_stau_span = 2.0 / num_steps

        # For each step, linear transform c (on [0, 1]) onto stau of the step boundaries
        # TODO: Change this to accommodate variable step sizes within the segment
        stau_stage = np.asarray([stau_step[i] + c * step_stau_span for i in range(num_steps)])
        stau_stage_flat = stau_stage.ravel()

        i1, i2 = gd.segment_indices[idx, :]
        stau_disc = gd.node_stau[i1:i2]

        self.L, self.D = lagrange_matrices(stau_disc, stau_stage_flat)
        _, D_dd = lagrange_matrices(stau_disc, stau_disc)
        self.D2 = np.dot(self.D, D_dd)

        self.add_input(name='dt_dstau', val=1.0, units=time_units,
                       desc='Ratio of segment time duration to segment tau duration (2).')

        self._setup_controls()

L_do_blocks = []
        D_do_blocks = []

        for iseg in range(num_seg):
            i1, i2 = gd.subset_segment_indices['control_disc'][iseg, :]
            indices = gd.subset_node_indices['control_disc'][i1:i2]
            nodes_given = gd.node_stau[indices]

            if output_nodes_per_seg is None:
                i1, i2 = gd.subset_segment_indices['all'][iseg, :]
                indices = gd.subset_node_indices['all'][i1:i2]
                nodes_eval = gd.node_stau[indices]
            else:
                nodes_eval = np.linspace(-1, 1, output_nodes_per_seg)

            L_block, D_block = lagrange_matrices(nodes_given, nodes_eval)

            L_do_blocks.append(L_block)
            D_do_blocks.append(D_block)

        L_do = block_diag(*L_do_blocks)
        D_do = block_diag(*D_do_blocks)

        self.L = np.dot(L_do, L_id)
        self.L_all = np.dot(L_da, L_id)
        self.D = np.dot(D_do, L_id)

        # Matrix D_dd interpolates rates at discretization nodes from values given at control
        # discretization nodes.
        _, D_dd = gd.phase_lagrange_matrices('control_disc', 'control_disc')

        # Matrix D2 provides second derivatives at output nodes given values at input nodes.

optau_segi = iptau_segi
            else:
                ptau_hi = igd.segment_ends[iseg+1]
                if iseg < igd.num_segments - 1:
                    idxs_in_iseg = np.where(output_nodes_ptau <= ptau_hi)[0]
                else:
                    idxs_in_iseg = np.arange(len(output_nodes_ptau))
                optau_segi = np.asarray(output_nodes_ptau)[idxs_in_iseg]
                # Remove the captured nodes so we don't accidentally include them again
                output_nodes_ptau = output_nodes_ptau[len(idxs_in_iseg):]

            # Now get the output nodes which fall in iseg in iseg's segment tau space.
            ostau_segi = 2.0 * (optau_segi - iptau_segi[0]) / (iptau_segi[-1] - iptau_segi[0]) - 1

            # Create the interpolation matrix and add it to the blocks
            L, _ = lagrange_matrices(istau_segi, ostau_segi)
            L_blocks.append(L)

        self.interpolation_matrix = block_diag(*L_blocks)
        r, c = np.nonzero(self.interpolation_matrix)

        output_num_nodes, input_num_nodes = self.interpolation_matrix.shape

        for (name, kwargs) in self._timeseries_outputs:

            input_kwargs = {k: kwargs[k] for k in ('units', 'desc')}
            input_name = 'input_values:{0}'.format(name)
            self.add_input(input_name,
                           shape=(input_num_nodes,) + kwargs['shape'],
                           **input_kwargs)

            output_name = name

#
        # 1. Get the segment tau values of the given nodes.
        #
        i1, i2 = grid_data.subset_segment_indices['control_disc'][iseg, :]
        indices = grid_data.subset_node_indices['control_disc'][i1:i2]
        tau_s_given = grid_data.node_stau[indices]

        #
        # 2. Get the segment tau values of the evaluation nodes.
        #
        t_eval_iseg = t_eval_per_seg[iseg]
        t0_iseg, tf_iseg = time_seg_ends[iseg, :]
        tau_s_eval = 2.0 * (t_eval_iseg - t0_iseg) / (tf_iseg - t0_iseg) - 1

        L_block, D_block = lagrange_matrices(tau_s_given, tau_s_eval)
        _, Daa_block = lagrange_matrices(tau_s_given, tau_s_given)

        L_blocks.append(L_block)
        D_blocks.append(D_block)
        Daa_blocks.append(Daa_block)

    L_ae = block_diag(*L_blocks)
    D_ae = block_diag(*D_blocks)

    D_aa = block_diag(*Daa_blocks)
    D2_ae = np.dot(D_ae, D_aa)

    return L_ae, D_ae, D2_ae

I: np.ndarray
        The integration matrix used to propagate initial states over segments

    """
    I_blocks = []

    for iseg in range(grid.num_segments):
        i1, i2 = grid.subset_segment_indices['all'][iseg, :]
        indices = grid.subset_node_indices['all'][i1:i2]
        nodes_given = grid.node_stau[indices]

        i1, i2 = grid.subset_segment_indices['all'][iseg, :]
        indices = grid.subset_node_indices['all'][i1:i2]
        nodes_eval = grid.node_stau[indices][1:]

        _, D_block = lagrange_matrices(nodes_given, nodes_eval)
        I_block = np.linalg.inv(D_block[:, 1:])
        I_blocks.append(I_block)

    I = block_diag(*I_blocks)

    return I