def jacvec_product(self, inputs, d_inputs, d_outputs, mode):
    if (self.metadata['jacobian_type'] == 'matvec'):
        k = self.metadata['num_comp']
        num_var = self.metadata['num_var']
        theta_out = self._cycle_names['theta_out']
        theta = self._cycle_names['theta']
        psi = self._cycle_names['psi']
        if (mode == 'fwd'):
            if (theta_out in d_outputs):
                if (theta in d_inputs):
                    d_outputs[theta_out] += (0.5 * d_inputs[theta])
                if (psi in d_inputs):
                    d_outputs[theta_out] += ((- d_inputs[psi]) / ((2 * k) - 2))
            for i in range(num_var):
                in_var = self._cycle_names['x'].format(i)
                if ((in_var in d_inputs) and ('x_norm2' in d_outputs)):
                    d_outputs['x_norm2'] += np.dot(inputs[in_var].flat, d_inputs[in_var].flat)
        elif (mode == 'rev'):
            if ('x_norm2' in d_outputs):
                dxnorm = d_outputs['x_norm2']
                for i in range(num_var):
                    x_i_name = self._cycle_names['x'].format(i)
                    if (x_i_name in d_inputs):
                        d_inputs[x_i_name] += (inputs[x_i_name] * dxnorm)
            if (theta_out in d_outputs):
                dtheta_out = d_outputs[theta_out]
                if (theta in d_inputs):
                    d_inputs[theta] += (0.5 * dtheta_out)
                if (psi in d_inputs):
                    d_inputs[psi] += ((- dtheta_out) / ((2 * k) - 2))
def jacvec_product(self, inputs, d_inputs, d_outputs, mode):
    if (self.metadata['jacobian_type'] == 'matvec'):
        k = self.metadata['num_comp']
        num_var = self.metadata['num_var']
        theta_out = self._cycle_names['theta_out']
        theta = self._cycle_names['theta']
        psi = self._cycle_names['psi']
        if (mode == 'fwd'):
            if (theta_out in d_outputs):
                if (theta in d_inputs):
                    d_outputs[theta_out] += (0.5 * d_inputs[theta])
                if (psi in d_inputs):
                    d_outputs[theta_out] += ((- d_inputs[psi]) / ((2 * k) - 2))
            for i in range(num_var):
                in_var = self._cycle_names['x'].format(i)
                if ((in_var in d_inputs) and ('x_norm2' in d_outputs)):
                    d_outputs['x_norm2'] += np.dot(inputs[in_var].flat, d_inputs[in_var].flat)
        elif (mode == 'rev'):
            if ('x_norm2' in d_outputs):
                dxnorm = d_outputs['x_norm2']
                for i in range(num_var):
                    x_i_name = self._cycle_names['x'].format(i)
                    if (x_i_name in d_inputs):
                        d_inputs[x_i_name] += (inputs[x_i_name] * dxnorm)
            if (theta_out in d_outputs):
                dtheta_out = d_outputs[theta_out]
                if (theta in d_inputs):
                    d_inputs[theta] += (0.5 * dtheta_out)
                if (psi in d_inputs):
                    d_inputs[psi] += ((- dtheta_out) / ((2 * k) - 2))