Gates_block.cpp

Go to the documentation of this file.

/*
Created on Fri Jun 26 14:13:26 2020
Copyright 2020 Peter Rakyta, Ph.D.

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

    http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.

@author: Peter Rakyta, Ph.D.
*/
#include <deque>
#include <set>
#include <tbb/enumerable_thread_specific.h>
#include "CU.h"
#include "CZ.h"
#include "CH.h"
#include "CNOT.h"
#include "CCX.h"
#include "SWAP.h"
#include "CSWAP.h"
#include "U1.h"
#include "U2.h"
#include "U3.h"
#include "RX.h"
#include "RY.h"
#include "R.h"
#include "RZ.h"
#include "CRY.h"
#include "CRX.h"
#include "CRZ.h"
#include "CP.h"
#include "CR.h"
#include "H.h"
#include "X.h"
#include "Y.h"
#include "Z.h"
#include "S.h"
#include "SDG.h"
#include "T.h"
#include "Tdg.h"
#include "SX.h"
#include "SXdg.h"
#include "SYC.h"
#include "CROT.h"
#include "RXX.h"
#include "RYY.h"
#include "RZZ.h"
#include "Adaptive.h"
#include "Gates_block.h"
#include "qgd_math.h"

#ifdef _WIN32
#include <cstdio>
#endif

#include "apply_large_kernel_to_input.h"

#ifdef USE_AVX 
#include "apply_large_kernel_to_input_AVX.h"
#endif

static tbb::enumerable_thread_specific<int> gates_block_derivative_depth([](){ return 0; });
static tbb::enumerable_thread_specific<int> gates_block_fusion_depth([](){ return 0; });


//static tbb::spin_mutex my_mutex;
Gates_block::Gates_block() : Gate() {

    // A string labeling the gate operation
    name = "Circuit";

    // A string describing the type of the operation
    type = BLOCK_OPERATION;
    // number of operation layers
    layer_num = 0;
    min_fusion = 14;
    involved_qubits_cache_valid = false;
    involved_target_qubits_cache_valid = false;
    involved_qubits_cache_mutex = std::make_shared<tbb::spin_mutex>();
}



Gates_block::Gates_block(int qbit_num_in) : Gate(qbit_num_in) {

    // A string labeling the gate operation
    name = "Circuit";

    // A string describing the type of the operation
    type = BLOCK_OPERATION;
    // number of operation layers
    layer_num = 0;
    
    min_fusion = 14;
    involved_qubits_cache_valid = false;
    involved_target_qubits_cache_valid = false;
    involved_qubits_cache_mutex = std::make_shared<tbb::spin_mutex>();
    
}

void
Gates_block::invalidate_structure_cache() {
    {
        tbb::spin_mutex::scoped_lock lock(*involved_qubits_cache_mutex);
        involved_qubits_cache_valid = false;
        involved_target_qubits_cache_valid = false;
        involved_qubits_cache.clear();
        involved_target_qubits_cache.clear();
    }
    fusion_block.update(std::shared_ptr<Gates_block>{});
}


Gates_block::~Gates_block() {

    release_gates();

}

void
Gates_block::release_gates() {

    //free the alloctaed memory of the stored gates
    for(std::vector<Gate*>::iterator it = gates.begin(); it != gates.end(); ++it) {

        Gate* gate = *it;
        delete gate;
        
    }
    
    gates.clear();
    layer_num = 0;
    parameter_num = 0;
    invalidate_structure_cache();

}


void
Gates_block::release_gate( int idx) {

    if ( idx>= (int)gates.size() ) return;

    // fist decrese the number of parameters
    Gate* gate = gates[idx];
    parameter_num -= gate->get_parameter_num();

    gates.erase( gates.begin() + idx );
    delete gate;
    
    // TODO: develop a more efficient method for large circuits. Now it is updating the whole circuit
    reset_parameter_start_indices();
    reset_dependency_graph();
    invalidate_structure_cache();


}

Matrix
Gates_block::get_matrix( Matrix_real& parameters ) {

    return get_matrix( parameters, false );


}



Matrix
Gates_block::get_matrix( Matrix_real& parameters, int parallel ) {

    //The stringstream input to store the output messages.
    std::stringstream sstream;

    // create matrix representation of the gate operations
    Matrix block_mtx = create_identity(matrix_size);

    apply_to(parameters, block_mtx, parallel);

#ifdef DEBUG
    if (block_mtx.isnan()) {
        std::string err( "Gates_block::get_matrix: block_mtx contains NaN." ) ;
         throw( err );       
    }
#endif

    return block_mtx;


}


Matrix_float
Gates_block::get_matrix( Matrix_real_float& parameters ) {
    return get_matrix( parameters, 0 );
}

Matrix_float
Gates_block::get_matrix( Matrix_real_float& parameters, int parallel ) {

    Matrix_float block_mtx = create_identity_float(matrix_size);
    apply_to(parameters, block_mtx, parallel);
    return block_mtx;

}




namespace {
template<typename Params, typename Input>
static void apply_to_list_impl(Gates_block* self, Params& parameters_mtx, std::vector<Input>& inputs, int parallel) {
    if (parallel == 0) {
        for (Input& input : inputs) {
            self->apply_to(parameters_mtx, input, parallel);
        }
        return;
    }

    int work_batch = 1;
    tbb::parallel_for( tbb::blocked_range<int>(0, static_cast<int>(inputs.size()), work_batch),
        [&](tbb::blocked_range<int> r) {
            for (int idx = r.begin(); idx < r.end(); ++idx) {
                self->apply_to(parameters_mtx, inputs[idx], parallel);
            }
        });
}

static Matrix_real precompute_block_sincos(const Matrix_real& parameters) {
    const int parameter_count = static_cast<int>(parameters.size());
    if (parameter_count <= 0) {
        return Matrix_real(0, 0);
    }

    Matrix_real sincos(parameter_count, 2);
    qgd_sincos_batch<double>(parameters.get_data(), sincos.get_data(), parameter_count, sincos.stride);

    return sincos;
}

static Matrix_real_float precompute_block_sincos(const Matrix_real_float& parameters) {
    const int parameter_count = static_cast<int>(parameters.size());
    if (parameter_count <= 0) {
        return Matrix_real_float(0, 0);
    }

    Matrix_real_float sincos(parameter_count, 2);
    qgd_sincos_batch<float>(parameters.get_data(), sincos.get_data(), parameter_count, sincos.stride);

    return sincos;
}

static void reset_identity(Matrix& mtx, int matrix_size) {
    if (mtx.rows != matrix_size || mtx.cols != matrix_size || mtx.stride != matrix_size) {
        mtx = Matrix(matrix_size, matrix_size);
    }
    memset(mtx.get_data(), 0, mtx.size()*sizeof(QGD_Complex16));
    for (int idx = 0; idx < matrix_size; ++idx) {
        mtx[idx*mtx.stride + idx].real = 1.0;
    }
}

static void reset_identity(Matrix_float& mtx, int matrix_size) {
    if (mtx.rows != matrix_size || mtx.cols != matrix_size || mtx.stride != matrix_size) {
        mtx = Matrix_float(matrix_size, matrix_size);
    }
    memset(mtx.get_data(), 0, mtx.size()*sizeof(QGD_Complex8));
    for (int idx = 0; idx < matrix_size; ++idx) {
        mtx[idx*mtx.stride + idx].real = 1.0f;
    }
}

struct BlockDerivativeWorkspace {
    std::vector<Gate> suffix_gates;
    std::vector<Matrix> suffix_matrices;
    std::vector<Matrix> forward_inputs;
    std::vector<std::vector<Matrix>> gate_gradients;
    std::vector<int> suffix_targets;
    Matrix suffix;
};

struct BlockDerivativeWorkspaceFloat {
    std::vector<Gate> suffix_gates;
    std::vector<Matrix_float> suffix_matrices;
    std::vector<Matrix_float> forward_inputs;
    std::vector<std::vector<Matrix_float>> gate_gradients;
    std::vector<int> suffix_targets;
    Matrix_float suffix;
};

static void prepare_suffix_gates(std::vector<Gate>& suffix_gates, size_t gate_count, int qbit_num) {
    if (suffix_gates.size() == gate_count
        && (suffix_gates.empty() || suffix_gates[0].get_qbit_num() == qbit_num)) {
        return;
    }

    suffix_gates.clear();
    suffix_gates.reserve(gate_count);
    for (size_t idx = 0; idx < gate_count; ++idx) {
        suffix_gates.emplace_back(qbit_num);
    }
}

static bool should_use_suffix(int qbit_num, size_t gate_count, int parameter_num) {
    if (qbit_num <= 4) return false;
    if (qbit_num >= 7) return true;
    double m = (double)(1 << qbit_num);
    return (double)parameter_num * (double)gate_count > ((double)gate_count + (double)parameter_num) * m;
}

static void prepare_all_qubits(std::vector<int>& qubits, int qbit_num) {
    if (qubits.size() != static_cast<size_t>(qbit_num)) {
        qubits.resize(static_cast<size_t>(qbit_num));
    }
    for (int idx = 0; idx < qbit_num; ++idx) {
        qubits[static_cast<size_t>(idx)] = idx;
    }
}

static void build_suffix_gates(
    const std::vector<Gate*>& gates,
    int qbit_num,
    int matrix_size,
    Matrix_real& parameters_mtx_in,
    const Matrix_real& precomputed_sincos,
    std::vector<Gate>& suffix_gates,
    std::vector<Matrix>& suffix_matrices,
    std::vector<int>& targets,
    Matrix& suffix) {

    prepare_suffix_gates(suffix_gates, gates.size(), qbit_num);
    suffix_matrices.resize(gates.size());

    if (gates.size() < 2) {
        return;
    }

    prepare_all_qubits(targets, qbit_num);
    reset_identity(suffix, matrix_size);

    for (int idx = static_cast<int>(gates.size()) - 1; idx >= 0; --idx) {
        if (static_cast<size_t>(idx) + 1 < gates.size()) {
            Matrix& suffix_after_current = suffix_matrices[static_cast<size_t>(idx)];
            suffix.copy_to(suffix_after_current);
            Gate& suffix_gate = suffix_gates[static_cast<size_t>(idx)];
            suffix_gate.set_target_qbits(targets);
            suffix_gate.set_matrix(suffix_after_current);
        }

        Gate* operation = gates[static_cast<size_t>(idx)];
        const int op_param_num = operation->get_parameter_num();
        const int op_param_start_idx = operation->get_parameter_start_idx();
        Matrix_real parameters_mtx(
            parameters_mtx_in.get_data() + op_param_start_idx,
            1,
            op_param_num
        );
        Matrix_real precomputed_sincos_loc(
            precomputed_sincos.get_data() + 2 * op_param_start_idx,
            op_param_num,
            2
        );

        if (op_param_num == 0 && operation->get_type() != BLOCK_OPERATION) {
            operation->apply_from_right(suffix);
        }
        else if (operation->get_type() == ADAPTIVE_OPERATION) {
            operation->apply_from_right(parameters_mtx, suffix);
        }
        else {
            operation->apply_from_right_inner(parameters_mtx, precomputed_sincos_loc, suffix);
        }
    }
}

static void build_forward_inputs(
    const std::vector<Gate*>& gates,
    Matrix& input,
    Matrix_real& parameters_mtx_in,
    const Matrix_real& precomputed_sincos,
    int parallel,
    std::vector<Matrix>& forward_inputs) {

    forward_inputs.resize(gates.size());
    if (gates.empty()) {
        return;
    }

    input.copy_to(forward_inputs[0]);
    for (size_t idx = 0; idx + 1 < gates.size(); ++idx) {
        forward_inputs[idx].copy_to(forward_inputs[idx + 1]);
        Gate* operation = gates[idx];
        const int op_param_num = operation->get_parameter_num();
        const int op_param_start_idx = operation->get_parameter_start_idx();
        Matrix_real parameters_mtx(
            parameters_mtx_in.get_data() + op_param_start_idx,
            1,
            op_param_num
        );
        Matrix_real precomputed_sincos_loc(
            precomputed_sincos.get_data() + 2 * op_param_start_idx,
            op_param_num,
            2
        );

        if (op_param_num == 0) {
            operation->apply_to(forward_inputs[idx + 1], parallel);
        }
        else if (operation->get_type() == ADAPTIVE_OPERATION) {
            operation->apply_to(parameters_mtx, forward_inputs[idx + 1], parallel);
        }
        else {
            operation->apply_to_inner(parameters_mtx, precomputed_sincos_loc, forward_inputs[idx + 1], parallel);
        }
    }
}

static void build_forward_inputs(
    const std::vector<Gate*>& gates,
    Matrix_float& input,
    Matrix_real_float& parameters_mtx_in,
    const Matrix_real_float& precomputed_sincos,
    int parallel,
    std::vector<Matrix_float>& forward_inputs) {

    forward_inputs.resize(gates.size());
    if (gates.empty()) {
        return;
    }

    input.copy_to(forward_inputs[0]);
    for (size_t idx = 0; idx + 1 < gates.size(); ++idx) {
        forward_inputs[idx].copy_to(forward_inputs[idx + 1]);
        Gate* operation = gates[idx];
        const int op_param_num = operation->get_parameter_num();
        const int op_param_start_idx = operation->get_parameter_start_idx();
        Matrix_real_float parameters_mtx(
            parameters_mtx_in.get_data() + op_param_start_idx,
            1,
            op_param_num
        );
        Matrix_real_float precomputed_sincos_loc(
            precomputed_sincos.get_data() + 2 * op_param_start_idx,
            op_param_num,
            2
        );

        if (op_param_num == 0) {
            operation->apply_to(forward_inputs[idx + 1], parallel);
        }
        else if (operation->get_type() == ADAPTIVE_OPERATION) {
            operation->apply_to(parameters_mtx, forward_inputs[idx + 1], parallel);
        }
        else {
            operation->apply_to_inner(parameters_mtx, precomputed_sincos_loc, forward_inputs[idx + 1], parallel);
        }
    }
}

static void build_suffix_gates(
    const std::vector<Gate*>& gates,
    int qbit_num,
    int matrix_size,
    Matrix_real_float& parameters_mtx_in,
    const Matrix_real_float& precomputed_sincos,
    std::vector<Gate>& suffix_gates,
    std::vector<Matrix_float>& suffix_matrices,
    std::vector<int>& targets,
    Matrix_float& suffix) {

    prepare_suffix_gates(suffix_gates, gates.size(), qbit_num);
    suffix_matrices.resize(gates.size());

    if (gates.size() < 2) {
        return;
    }

    prepare_all_qubits(targets, qbit_num);
    reset_identity(suffix, matrix_size);

    for (int idx = static_cast<int>(gates.size()) - 1; idx >= 0; --idx) {
        if (static_cast<size_t>(idx) + 1 < gates.size()) {
            Matrix_float& suffix_after_current = suffix_matrices[static_cast<size_t>(idx)];
            suffix.copy_to(suffix_after_current);
            Gate& suffix_gate = suffix_gates[static_cast<size_t>(idx)];
            suffix_gate.set_target_qbits(targets);
            suffix_gate.set_matrix(suffix_after_current);
        }

        Gate* operation = gates[static_cast<size_t>(idx)];
        const int op_param_num = operation->get_parameter_num();
        const int op_param_start_idx = operation->get_parameter_start_idx();
        Matrix_real_float parameters_mtx(
            parameters_mtx_in.get_data() + op_param_start_idx,
            1,
            op_param_num
        );
        Matrix_real_float precomputed_sincos_loc(
            precomputed_sincos.get_data() + 2 * op_param_start_idx,
            op_param_num,
            2
        );

        if (op_param_num == 0 && operation->get_type() != BLOCK_OPERATION) {
            operation->apply_from_right(suffix);
        }
        else if (operation->get_type() == ADAPTIVE_OPERATION) {
            operation->apply_from_right(parameters_mtx, suffix);
        }
        else {
            operation->apply_from_right_inner(parameters_mtx, precomputed_sincos_loc, suffix);
        }
    }
}
} // anonymous namespace


void 
Gates_block::apply_to_list( Matrix_real& parameters_mtx, std::vector<Matrix>& inputs, int parallel ) {

    apply_to_list_impl(this, parameters_mtx, inputs, parallel);

}


void
Gates_block::apply_to_list( Matrix_real_float& parameters_mtx, std::vector<Matrix_float>& inputs, int parallel ) {

    apply_to_list_impl(this, parameters_mtx, inputs, parallel);

}


void 
Gates_block::apply_to( Matrix_real& parameters_mtx_in, Matrix& input, int parallel ) {

    Matrix_real precomputed_sincos = precompute_block_sincos(parameters_mtx_in);
    apply_to_inner(parameters_mtx_in, precomputed_sincos, input, parallel);
}


void
Gates_block::apply_to_inner( Matrix_real& parameters_mtx_in, const Matrix_real& precomputed_sincos_in, Matrix& input, int parallel ) {



    std::vector<int> involved_qubits = get_involved_qubits();
    
    if (input.rows != matrix_size ) {
        std::string err("Gates_block::apply_to: Wrong input size in Gates_block gate apply.");
        throw err;    
    }

    if (qbit_num > 31) {
        std::string err("Gates_block::apply_to: Number of qubits supported up to 31"); 
        throw err;        
    }

    int size = static_cast<int>(involved_qubits.size());

    if (min_fusion != -1 && qbit_num >= min_fusion && size <= (input.cols == 1 ? 5 : 2) && qbit_num != size && gates.size() > 1) {
        auto fb = fusion_block.get();
        int& fusion_depth = gates_block_fusion_depth.local();
        static tbb::enumerable_thread_specific<std::deque<Matrix>> Umtx_mini_reusable_tls([](){ return std::deque<Matrix>{}; });
        std::deque<Matrix>& Umtx_mini_reusable = Umtx_mini_reusable_tls.local();
        if (Umtx_mini_reusable.size() <= static_cast<size_t>(fusion_depth)) {
            Umtx_mini_reusable.resize(fusion_depth + 1);
        }
        Matrix& Umtx_mini = Umtx_mini_reusable[fusion_depth];
        reset_identity(Umtx_mini, Power_of_2(size));
        if (fb == nullptr) {

            std::vector<int> old_to_new(qbit_num, -2);
            for (int i = 0; i < size; i++){
                old_to_new[i] = involved_qubits[i];
            }

            Gates_block* clone_block = clone();
            clone_block->reorder_qubits(old_to_new);
            clone_block->set_qbit_num(size);            
            fusion_block.update(clone_block);
            fb = fusion_block.get();
        }
        fusion_depth++;
        fb->apply_to_inner(parameters_mtx_in, precomputed_sincos_in, Umtx_mini, parallel);
        fusion_depth--;

        static tbb::enumerable_thread_specific<std::deque<Gate>> merged_gate_reusable_tls([](){ return std::deque<Gate>{}; });
        static tbb::enumerable_thread_specific<std::deque<int>> merged_gate_qbit_num_tls([](){ return std::deque<int>{}; });
        std::deque<Gate>& merged_gate_reusable = merged_gate_reusable_tls.local();
        std::deque<int>& merged_gate_qbit_num = merged_gate_qbit_num_tls.local();
        if (merged_gate_reusable.size() <= static_cast<size_t>(fusion_depth)) {
            merged_gate_reusable.resize(fusion_depth + 1);
            merged_gate_qbit_num.resize(fusion_depth + 1, -1);
        }
        Gate& merged_gate = merged_gate_reusable[fusion_depth];
        int& merged_gate_qbits = merged_gate_qbit_num[fusion_depth];
        if (merged_gate_qbits != qbit_num) {
            merged_gate = Gate(qbit_num);
            merged_gate_qbits = qbit_num;
        }
        merged_gate.set_target_qbits(involved_qubits);
        merged_gate.set_matrix(Umtx_mini);

        // Keep legacy high-qubit fused execution behavior, but route through Gate::apply_to.
        int fused_parallel = parallel;
        if (parallel != 0 && size <= 5 && qbit_num >= 14) {
            fused_parallel = 2;
        }
        merged_gate.apply_to(input, fused_parallel);
    } else {
        for( size_t idx=0; idx<gates.size(); idx++) {
            Gate* operation = gates[idx];
            const int op_param_num = operation->get_parameter_num();
            const int op_param_start_idx = operation->get_parameter_start_idx();
            
            Matrix_real parameters_mtx_loc(parameters_mtx_in.get_data() + op_param_start_idx, 1, op_param_num);
            Matrix_real precomputed_sincos_loc(precomputed_sincos_in.get_data() + 2 * op_param_start_idx, op_param_num, 2);
            
            if  ( parameters_mtx_loc.size() == 0 && operation->get_type() != BLOCK_OPERATION ) {
                operation->apply_to(input, parallel);            
            }
            else {
                if (operation->get_type() == ADAPTIVE_OPERATION) {
                    operation->apply_to(parameters_mtx_loc, input, parallel);
                }
                else {
                    operation->apply_to_inner(parameters_mtx_loc, precomputed_sincos_loc, input, parallel);
                }
            }
#ifdef DEBUG
            if (input.isnan()) {
                std::string err( "Gates_block::apply_to: transformed matrix contains NaN." );
                throw( err );
            }
#endif
        }
    }
}


void
Gates_block::apply_to( Matrix_float& input, int parallel ) {

    Matrix_real_float parameters_empty(1, 0);
    apply_to(parameters_empty, input, parallel);
}


void
Gates_block::apply_to( Matrix_real_float& parameters_mtx_in, Matrix_float& input, int parallel ) {

    Matrix_real_float precomputed_sincos = precompute_block_sincos(parameters_mtx_in);
    apply_to_inner(parameters_mtx_in, precomputed_sincos, input, parallel);
}


void
Gates_block::apply_to_inner( Matrix_real_float& parameters_mtx_in, const Matrix_real_float& precomputed_sincos_in, Matrix_float& input, int parallel ) {

    std::vector<int> involved_qubits = get_involved_qubits();

    if (input.rows != matrix_size ) {
        std::string err("Gates_block::apply_to(Matrix_real_float&, Matrix_float&): Wrong input size in Gates_block gate apply.");
        throw err;
    }

    if (qbit_num > 31) {
        std::string err("Gates_block::apply_to(Matrix_real_float&, Matrix_float&): Number of qubits supported up to 31");
        throw err;
    }

    int size = static_cast<int>(involved_qubits.size());

    if (min_fusion != -1 && qbit_num >= min_fusion && size <= (input.cols == 1 ? 5 : 2) && qbit_num != size && gates.size() > 1) {
        auto fb = fusion_block.get();
        int& fusion_depth = gates_block_fusion_depth.local();
        static tbb::enumerable_thread_specific<std::deque<Matrix_float>> Umtx_mini_reusable_tls([](){ return std::deque<Matrix_float>{}; });
        std::deque<Matrix_float>& Umtx_mini_reusable = Umtx_mini_reusable_tls.local();
        if (Umtx_mini_reusable.size() <= static_cast<size_t>(fusion_depth)) {
            Umtx_mini_reusable.resize(fusion_depth + 1);
        }
        Matrix_float& Umtx_mini = Umtx_mini_reusable[fusion_depth];
        reset_identity(Umtx_mini, Power_of_2(size));
        if (fb == nullptr) {

            std::vector<int> old_to_new(qbit_num, -2);
            for (int i = 0; i < size; i++) {
                old_to_new[i] = involved_qubits[i];
            }

            Gates_block* clone_block = clone();
            clone_block->reorder_qubits(old_to_new);
            clone_block->set_qbit_num(size);
            fusion_block.update(clone_block);
            fb = fusion_block.get();
        }
        fusion_depth++;
        fb->apply_to_inner(parameters_mtx_in, precomputed_sincos_in, Umtx_mini, parallel);
        fusion_depth--;

        static tbb::enumerable_thread_specific<std::deque<Gate>> merged_gate_reusable_tls([](){ return std::deque<Gate>{}; });
        static tbb::enumerable_thread_specific<std::deque<int>> merged_gate_qbit_num_tls([](){ return std::deque<int>{}; });
        std::deque<Gate>& merged_gate_reusable = merged_gate_reusable_tls.local();
        std::deque<int>& merged_gate_qbit_num = merged_gate_qbit_num_tls.local();
        if (merged_gate_reusable.size() <= static_cast<size_t>(fusion_depth)) {
            merged_gate_reusable.resize(fusion_depth + 1);
            merged_gate_qbit_num.resize(fusion_depth + 1, -1);
        }
        Gate& merged_gate = merged_gate_reusable[fusion_depth];
        int& merged_gate_qbits = merged_gate_qbit_num[fusion_depth];
        if (merged_gate_qbits != qbit_num) {
            merged_gate = Gate(qbit_num);
            merged_gate_qbits = qbit_num;
        }
        merged_gate.set_target_qbits(involved_qubits);
        merged_gate.set_matrix(Umtx_mini);

        // Keep legacy high-qubit fused execution behavior, but route through Gate::apply_to.
        int fused_parallel = parallel;
        if (parallel != 0 && size <= 5 && qbit_num >= 14) {
            fused_parallel = 2;
        }
        merged_gate.apply_to(input, fused_parallel);
    }
    else {
        for (size_t idx = 0; idx < gates.size(); idx++) {
            Gate* operation = gates[idx];
            const int op_param_num = operation->get_parameter_num();
            const int op_param_start_idx = operation->get_parameter_start_idx();

            Matrix_real_float parameters_mtx_loc(
                parameters_mtx_in.get_data() + op_param_start_idx,
                1,
                op_param_num
            );

            Matrix_real_float precomputed_sincos_loc(
                precomputed_sincos_in.get_data() + 2 * op_param_start_idx,
                op_param_num,
                2
            );

            if (parameters_mtx_loc.size() == 0 && operation->get_type() != BLOCK_OPERATION) {
                operation->apply_to(input, parallel);
            }
            else {
                if (operation->get_type() == ADAPTIVE_OPERATION) {
                    operation->apply_to(parameters_mtx_loc, input, parallel);
                }
                else {
                    operation->apply_to_inner(parameters_mtx_loc, precomputed_sincos_loc, input, parallel);
                }
            }
#ifdef DEBUG
            if (input.isnan()) {
                std::string err("Gates_block::apply_to(Matrix_real_float&, Matrix_float&): transformed matrix contains NaN.");
                throw(err);
            }
#endif
        }
    }
}


void
Gates_block::apply_to( Matrix_real_any& parameters_mtx_in, Matrix_any& input, int parallel ) {

    if (parameters_mtx_in.is_float64() && input.is_float64()) {
        apply_to(parameters_mtx_in.as_float64(), input.as_float64(), parallel);
        return;
    }

    if (parameters_mtx_in.is_float32() && input.is_float32()) {
        apply_to(parameters_mtx_in.as_float32(), input.as_float32(), parallel);
        return;
    }

    std::string err("Gates_block::apply_to(Matrix_real_any&, Matrix_any&): precision mismatch between parameters and input");
    throw err;
}




bool is_qbit_present(std::vector<int> involved_qubits, int new_qbit, int num_of_qbits){

    bool contained=false;
    
    for (int idx=0; idx<num_of_qbits; idx++) {
    
        if(involved_qubits[idx] == new_qbit) {
            contained=true;
        }
    }
    
    return contained;
    
}


void 
Gates_block::apply_from_right( Matrix_real& parameters_mtx, Matrix& input ) {

    Matrix_real precomputed_sincos = precompute_block_sincos(parameters_mtx);
    apply_from_right_inner(parameters_mtx, precomputed_sincos, input);
}


void
Gates_block::apply_from_right_inner( Matrix_real& parameters_mtx, const Matrix_real& precomputed_sincos, Matrix& input ) {

    for (int idx = static_cast<int>(gates.size()) - 1; idx >= 0; --idx) {

        Gate* operation = gates[idx];
        const int op_param_num = operation->get_parameter_num();
        const int op_param_start_idx = operation->get_parameter_start_idx();

        Matrix_real parameters_mtx_loc(
            parameters_mtx.get_data() + op_param_start_idx,
            1,
            op_param_num
        );
        Matrix_real precomputed_sincos_loc(
            precomputed_sincos.get_data() + 2 * op_param_start_idx,
            op_param_num,
            2
        );

        if (parameters_mtx_loc.size() == 0 && operation->get_type() != BLOCK_OPERATION) {
            operation->apply_from_right(input);
        }
        else {
            if (operation->get_type() == ADAPTIVE_OPERATION) {
                operation->apply_from_right(parameters_mtx_loc, input);
            }
            else {
                operation->apply_from_right_inner(parameters_mtx_loc, precomputed_sincos_loc, input);
            }
        }

#ifdef DEBUG
        if (input.isnan()) { 
            std::string err( "Gates_block::apply_from_right: transformed matrix contains NaN." );
             throw( err );            
        }
#endif


    }


}


void
Gates_block::apply_from_right( Matrix_real_float& parameters_mtx, Matrix_float& input ) {

    Matrix_real_float precomputed_sincos = precompute_block_sincos(parameters_mtx);
    apply_from_right_inner(parameters_mtx, precomputed_sincos, input);
}


void
Gates_block::apply_from_right_inner( Matrix_real_float& parameters_mtx, const Matrix_real_float& precomputed_sincos, Matrix_float& input ) {

    for (int idx = static_cast<int>(gates.size()) - 1; idx >= 0; --idx) {

        Gate* operation = gates[idx];
        const int op_param_num = operation->get_parameter_num();
        const int op_param_start_idx = operation->get_parameter_start_idx();

        Matrix_real_float parameters_mtx_loc(
            parameters_mtx.get_data() + op_param_start_idx,
            1,
            op_param_num
        );

        Matrix_real_float precomputed_sincos_loc(
            precomputed_sincos.get_data() + 2 * op_param_start_idx,
            op_param_num,
            2
        );

        if (parameters_mtx_loc.size() == 0 && operation->get_type() != BLOCK_OPERATION) {
            operation->apply_from_right(input);
        }
        else {
            if (operation->get_type() == ADAPTIVE_OPERATION) {
                operation->apply_from_right(parameters_mtx_loc, input);
            }
            else {
                operation->apply_from_right_inner(parameters_mtx_loc, precomputed_sincos_loc, input);
            }
        }

#ifdef DEBUG
        if (input.isnan()) {
            std::string err("Gates_block::apply_from_right(Matrix_real_float&, Matrix_float&): transformed matrix contains NaN.");
            throw(err);
        }
#endif
    }
}


std::vector<Matrix> 
Gates_block::apply_derivate_to( Matrix_real& parameters_mtx_in, Matrix& input, int parallel ) {

    std::vector<Matrix> grad;
    apply_derivate_to(parameters_mtx_in, input, parallel, grad);
    return grad;
}


void
Gates_block::apply_derivate_to( Matrix_real& parameters_mtx_in, Matrix& input, int parallel, std::vector<Matrix>& output ) {

    apply_derivate_to(parameters_mtx_in, input, parallel, output, 0, true);
}


void
Gates_block::apply_derivate_to( Matrix_real& parameters_mtx_in, Matrix& input, int parallel, std::vector<Matrix>& output, size_t output_offset, bool resize_output ) {

    const size_t output_size = output_offset + static_cast<size_t>(parameter_num);
    if (resize_output) {
        output.resize(output_size);
    }
    else if (output.size() < output_size) {
        std::string err("Gates_block::apply_derivate_to: output vector is too small.");
        throw err;
    }

    Matrix_real precomputed_sincos = precompute_block_sincos(parameters_mtx_in);
    static tbb::enumerable_thread_specific<std::deque<BlockDerivativeWorkspace>> workspace_tls;
    std::deque<BlockDerivativeWorkspace>& workspaces = workspace_tls.local();
    int& workspace_depth = gates_block_derivative_depth.local();
    if (workspaces.size() <= static_cast<size_t>(workspace_depth)) {
        workspaces.resize(static_cast<size_t>(workspace_depth) + 1);
    }
    BlockDerivativeWorkspace& workspace = workspaces[static_cast<size_t>(workspace_depth)];
    std::vector<Gate>& suffix_gates = workspace.suffix_gates;
    std::vector<Matrix>& suffix_matrices = workspace.suffix_matrices;
    std::vector<Matrix>& forward_inputs = workspace.forward_inputs;
    forward_inputs.resize(gates.size());
    workspace.gate_gradients.resize(gates.size());

    const bool use_suffix = should_use_suffix(qbit_num, gates.size(), parameter_num);

    if (parallel == 0 || gates.size() < 2) {
        if (use_suffix) {
            build_suffix_gates(gates, qbit_num, matrix_size, parameters_mtx_in, precomputed_sincos, suffix_gates, suffix_matrices, workspace.suffix_targets, workspace.suffix);
        }
        build_forward_inputs(gates, input, parameters_mtx_in, precomputed_sincos, parallel, forward_inputs);
    }
    else {
        if (use_suffix) {
            tbb::parallel_invoke(
                [&]() {
                    build_suffix_gates(gates, qbit_num, matrix_size, parameters_mtx_in, precomputed_sincos, suffix_gates, suffix_matrices, workspace.suffix_targets, workspace.suffix);
                },
                [&]() {
                    build_forward_inputs(gates, input, parameters_mtx_in, precomputed_sincos, parallel, forward_inputs);
                }
            );
        } else {
            build_forward_inputs(gates, input, parameters_mtx_in, precomputed_sincos, parallel, forward_inputs);
        }
    }

    auto apply_derivative_for_gate = [&](size_t deriv_idx) {
        Gate* gate_deriv = gates[deriv_idx];            

        // for constant gate no gradient component is calculated
        if ( gate_deriv->get_parameter_num() == 0 ) {
            return;
        }
            
        int deriv_parameter_idx = gate_deriv->get_parameter_start_idx();       

        int& derivative_depth = gates_block_derivative_depth.local();
        static tbb::enumerable_thread_specific<std::deque<Matrix>> input_loc_reusable_tls([](){ return std::deque<Matrix>{}; });
        std::deque<Matrix>& input_loc_reusable = input_loc_reusable_tls.local();
        if (input_loc_reusable.size() <= static_cast<size_t>(derivative_depth)) {
            input_loc_reusable.resize(derivative_depth + 1);
        }
        Matrix& input_loc = input_loc_reusable[derivative_depth];
        forward_inputs[deriv_idx].copy_to(input_loc);

        derivative_depth++;
        Gate* operation = gates[deriv_idx];
        const int op_param_num = operation->get_parameter_num();
        const int op_param_start_idx = operation->get_parameter_start_idx();
        Matrix_real parameters_mtx(parameters_mtx_in.get_data() + op_param_start_idx, 1, op_param_num);
        std::vector<Matrix>& grad_loc = workspace.gate_gradients[deriv_idx];

        if (operation->get_type() == ADAPTIVE_OPERATION) {
            operation->apply_derivate_to( parameters_mtx, input_loc, parallel, grad_loc );
        }
        else {
            operation->apply_derivate_to( parameters_mtx, input_loc, parallel, grad_loc );
        }

        if (deriv_idx + 1 < gates.size()) {
            if (use_suffix) {
                Gate& suffix_gate = suffix_gates[deriv_idx];
                if (parallel == 0) {
                    suffix_gate.apply_to_list(grad_loc, parallel);
                }
                else {
                    tbb::parallel_for(tbb::blocked_range<int>(0, static_cast<int>(grad_loc.size()), 1),
                        [&](tbb::blocked_range<int> r) {
                            for (int grad_idx = r.begin(); grad_idx < r.end(); ++grad_idx) {
                                suffix_gate.apply_to(grad_loc[static_cast<size_t>(grad_idx)], parallel);
                            }
                        });
                }
            } else {
                // Prefix-only: apply remaining gates sequentially using fast dedicated kernels.
                for (size_t g = deriv_idx + 1; g < gates.size(); ++g) {
                    Gate* op = gates[g];
                    const int pn = op->get_parameter_num();
                    const int ps = op->get_parameter_start_idx();
                    if (pn == 0) {
                        for (auto& m : grad_loc) op->apply_to(m, parallel);
                    } else {
                        Matrix_real pm(parameters_mtx_in.get_data() + ps, 1, pn);
                        Matrix_real psc(precomputed_sincos.get_data() + 2 * ps, pn, 2);
                        if (op->get_type() == ADAPTIVE_OPERATION) {
                            for (auto& m : grad_loc) op->apply_to(pm, m, parallel);
                        } else {
                            for (auto& m : grad_loc) op->apply_to_inner(pm, psc, m, parallel);
                        }
                    }
                }
            }
        }
        derivative_depth--;


        for ( int idx = 0; idx<(int)grad_loc.size(); idx++ ) {
            grad_loc[idx].copy_to(output[output_offset + static_cast<size_t>(deriv_parameter_idx + idx)]);
        }
    };

    if (parallel == 0) {
        for (size_t deriv_idx=0; deriv_idx<gates.size(); ++deriv_idx) {
            apply_derivative_for_gate(deriv_idx);
        }
        return;
    }

    int work_batch = 1;
    tbb::parallel_for( tbb::blocked_range<size_t>(0,gates.size(),work_batch), [&](tbb::blocked_range<size_t> r) {
        for (size_t deriv_idx=r.begin(); deriv_idx<r.end(); ++deriv_idx) {
            apply_derivative_for_gate(deriv_idx);
        }
    });
   

}


std::vector<Matrix_float>
Gates_block::apply_derivate_to( Matrix_real_float& parameters_mtx_in, Matrix_float& input, int parallel ) {

    std::vector<Matrix_float> grad;
    apply_derivate_to(parameters_mtx_in, input, parallel, grad);
    return grad;
}


void
Gates_block::apply_derivate_to( Matrix_real_float& parameters_mtx_in, Matrix_float& input, int parallel, std::vector<Matrix_float>& output ) {

    apply_derivate_to(parameters_mtx_in, input, parallel, output, 0, true);
}


void
Gates_block::apply_derivate_to( Matrix_real_float& parameters_mtx_in, Matrix_float& input, int parallel, std::vector<Matrix_float>& output, size_t output_offset, bool resize_output ) {

    const size_t output_size = output_offset + static_cast<size_t>(parameter_num);
    if (resize_output) {
        output.resize(output_size);
    }
    else if (output.size() < output_size) {
        std::string err("Gates_block::apply_derivate_to(Matrix_real_float): output vector is too small.");
        throw err;
    }

    Matrix_real_float precomputed_sincos = precompute_block_sincos(parameters_mtx_in);
    static tbb::enumerable_thread_specific<std::deque<BlockDerivativeWorkspaceFloat>> workspace_tls;
    std::deque<BlockDerivativeWorkspaceFloat>& workspaces = workspace_tls.local();
    int& workspace_depth = gates_block_derivative_depth.local();
    if (workspaces.size() <= static_cast<size_t>(workspace_depth)) {
        workspaces.resize(static_cast<size_t>(workspace_depth) + 1);
    }
    BlockDerivativeWorkspaceFloat& workspace = workspaces[static_cast<size_t>(workspace_depth)];
    std::vector<Gate>& suffix_gates = workspace.suffix_gates;
    std::vector<Matrix_float>& suffix_matrices = workspace.suffix_matrices;
    std::vector<Matrix_float>& forward_inputs = workspace.forward_inputs;
    forward_inputs.resize(gates.size());
    workspace.gate_gradients.resize(gates.size());

    const bool use_suffix = should_use_suffix(qbit_num, gates.size(), parameter_num);

    if (parallel == 0 || gates.size() < 2) {
        if (use_suffix) {
            build_suffix_gates(gates, qbit_num, matrix_size, parameters_mtx_in, precomputed_sincos, suffix_gates, suffix_matrices, workspace.suffix_targets, workspace.suffix);
        }
        build_forward_inputs(gates, input, parameters_mtx_in, precomputed_sincos, parallel, forward_inputs);
    }
    else {
        if (use_suffix) {
            tbb::parallel_invoke(
                [&]() {
                    build_suffix_gates(gates, qbit_num, matrix_size, parameters_mtx_in, precomputed_sincos, suffix_gates, suffix_matrices, workspace.suffix_targets, workspace.suffix);
                },
                [&]() {
                    build_forward_inputs(gates, input, parameters_mtx_in, precomputed_sincos, parallel, forward_inputs);
                }
            );
        } else {
            build_forward_inputs(gates, input, parameters_mtx_in, precomputed_sincos, parallel, forward_inputs);
        }
    }

    auto apply_derivative_for_gate = [&](size_t deriv_idx) {
        Gate* gate_deriv = gates[deriv_idx];            

        // for constant gate no gradient component is calculated
        if ( gate_deriv->get_parameter_num() == 0 ) {
            return;
        }
            
        int deriv_parameter_idx = gate_deriv->get_parameter_start_idx();       

        int& derivative_depth = gates_block_derivative_depth.local();
        static tbb::enumerable_thread_specific<std::deque<Matrix_float>> input_loc_reusable_tls([](){ return std::deque<Matrix_float>{}; });
        std::deque<Matrix_float>& input_loc_reusable = input_loc_reusable_tls.local();
        if (input_loc_reusable.size() <= static_cast<size_t>(derivative_depth)) {
            input_loc_reusable.resize(derivative_depth + 1);
        }
        Matrix_float& input_loc = input_loc_reusable[derivative_depth];
        forward_inputs[deriv_idx].copy_to(input_loc);

        derivative_depth++;
        Gate* operation = gates[deriv_idx];
        const int op_param_num = operation->get_parameter_num();
        const int op_param_start_idx = operation->get_parameter_start_idx();
        Matrix_real_float parameters_mtx(parameters_mtx_in.get_data() + op_param_start_idx, 1, op_param_num);
        std::vector<Matrix_float>& grad_loc = workspace.gate_gradients[deriv_idx];

        if (operation->get_type() == ADAPTIVE_OPERATION) {
            operation->apply_derivate_to( parameters_mtx, input_loc, parallel, grad_loc );
        }
        else {
            operation->apply_derivate_to( parameters_mtx, input_loc, parallel, grad_loc );
        }

        if (deriv_idx + 1 < gates.size()) {
            if (use_suffix) {
                Gate& suffix_gate = suffix_gates[deriv_idx];
                if (parallel == 0) {
                    suffix_gate.apply_to_list(grad_loc, parallel);
                }
                else {
                    tbb::parallel_for(tbb::blocked_range<int>(0, static_cast<int>(grad_loc.size()), 1),
                        [&](tbb::blocked_range<int> r) {
                            for (int grad_idx = r.begin(); grad_idx < r.end(); ++grad_idx) {
                                suffix_gate.apply_to(grad_loc[static_cast<size_t>(grad_idx)], parallel);
                            }
                        });
                }
            } else {
                // Prefix-only: apply remaining gates sequentially using fast dedicated kernels.
                for (size_t g = deriv_idx + 1; g < gates.size(); ++g) {
                    Gate* op = gates[g];
                    const int pn = op->get_parameter_num();
                    const int ps = op->get_parameter_start_idx();
                    if (pn == 0) {
                        for (auto& m : grad_loc) op->apply_to(m, parallel);
                    } else {
                        Matrix_real_float pm(parameters_mtx_in.get_data() + ps, 1, pn);
                        Matrix_real_float psc(precomputed_sincos.get_data() + 2 * ps, pn, 2);
                        if (op->get_type() == ADAPTIVE_OPERATION) {
                            for (auto& m : grad_loc) op->apply_to(pm, m, parallel);
                        } else {
                            for (auto& m : grad_loc) op->apply_to_inner(pm, psc, m, parallel);
                        }
                    }
                }
            }
        }
        derivative_depth--;


        for ( int idx = 0; idx<(int)grad_loc.size(); idx++ ) {
            grad_loc[idx].copy_to(output[output_offset + static_cast<size_t>(deriv_parameter_idx + idx)]);
        }
    };

    if (parallel == 0) {
        for (size_t deriv_idx=0; deriv_idx<gates.size(); ++deriv_idx) {
            apply_derivative_for_gate(deriv_idx);
        }
        return;
    }

    int work_batch = 1;
    tbb::parallel_for( tbb::blocked_range<size_t>(0,gates.size(),work_batch), [&](tbb::blocked_range<size_t> r) {
        for (size_t deriv_idx=r.begin(); deriv_idx<r.end(); ++deriv_idx) {
            apply_derivative_for_gate(deriv_idx);
        }
    });
   

}


std::vector<Matrix>
Gates_block::apply_to_combined( Matrix_real& parameters_mtx_in, Matrix& input, int parallel ) {

    Matrix_real precomputed_sincos = precompute_block_sincos(parameters_mtx_in);
    return apply_to_combined_inner(parameters_mtx_in, precomputed_sincos, input, parallel);
}


void
Gates_block::apply_to_combined( Matrix_real& parameters_mtx_in, Matrix& input, int parallel, std::vector<Matrix>& output ) {

    Matrix_real precomputed_sincos = precompute_block_sincos(parameters_mtx_in);
    apply_to_combined_inner(parameters_mtx_in, precomputed_sincos, input, parallel, output);
}


std::vector<Matrix>
Gates_block::apply_to_combined_inner( Matrix_real& parameters_mtx_in, const Matrix_real& precomputed_sincos, Matrix& input, int parallel ) {

    std::vector<Matrix> ret;
    apply_to_combined_inner(parameters_mtx_in, precomputed_sincos, input, parallel, ret);
    return ret;
}


void
Gates_block::apply_to_combined_inner( Matrix_real& parameters_mtx_in, const Matrix_real& precomputed_sincos, Matrix& input, int parallel, std::vector<Matrix>& output ) {

    output.resize(static_cast<size_t>(parameter_num) + 1);
    if (parallel == 0) {
        input.copy_to(output[0]);
        apply_to_inner(parameters_mtx_in, precomputed_sincos, output[0], parallel);
        apply_derivate_to(parameters_mtx_in, input, parallel, output, 1, false);
    }
    else {
        tbb::parallel_invoke(
            [&]() {
                input.copy_to(output[0]);
                apply_to_inner(parameters_mtx_in, precomputed_sincos, output[0], parallel);
            },
            [&]() {
                apply_derivate_to(parameters_mtx_in, input, parallel, output, 1, false);
            }
        );
    }
}


std::vector<Matrix_float>
Gates_block::apply_to_combined( Matrix_real_float& parameters_mtx_in, Matrix_float& input, int parallel ) {

    Matrix_real_float precomputed_sincos = precompute_block_sincos(parameters_mtx_in);
    return apply_to_combined_inner(parameters_mtx_in, precomputed_sincos, input, parallel);
}


void
Gates_block::apply_to_combined( Matrix_real_float& parameters_mtx_in, Matrix_float& input, int parallel, std::vector<Matrix_float>& output ) {

    Matrix_real_float precomputed_sincos = precompute_block_sincos(parameters_mtx_in);
    apply_to_combined_inner(parameters_mtx_in, precomputed_sincos, input, parallel, output);
}


std::vector<Matrix_float>
Gates_block::apply_to_combined_inner( Matrix_real_float& parameters_mtx_in, const Matrix_real_float& precomputed_sincos, Matrix_float& input, int parallel ) {

    std::vector<Matrix_float> ret;
    apply_to_combined_inner(parameters_mtx_in, precomputed_sincos, input, parallel, ret);
    return ret;
}


void
Gates_block::apply_to_combined_inner( Matrix_real_float& parameters_mtx_in, const Matrix_real_float& precomputed_sincos, Matrix_float& input, int parallel, std::vector<Matrix_float>& output ) {

    output.resize(static_cast<size_t>(parameter_num) + 1);
    if (parallel == 0) {
        input.copy_to(output[0]);
        apply_to_inner(parameters_mtx_in, precomputed_sincos, output[0], parallel);
        apply_derivate_to(parameters_mtx_in, input, parallel, output, 1, false);
    }
    else {
        tbb::parallel_invoke(
            [&]() {
                input.copy_to(output[0]);
                apply_to_inner(parameters_mtx_in, precomputed_sincos, output[0], parallel);
            },
            [&]() {
                apply_derivate_to(parameters_mtx_in, input, parallel, output, 1, false);
            }
        );
    }
}


void Gates_block::add_u1(int target_qbit) {
    
    // create the operation
    Gate* operation = static_cast<Gate*>(new U1( qbit_num, target_qbit ));
    
    // adding the operation to the end of the list of gates
    add_gate( operation );
}

void Gates_block::add_u1_to_front(int target_qbit ) {
    
    // create the operation
    Gate* gate = static_cast<Gate*>(new U1( qbit_num, target_qbit ));
    
    // adding the operation to the front of the list of gates
    add_gate_to_front( gate );
}

void Gates_block::add_u2(int target_qbit) {
    
    // create the operation
    Gate* operation = static_cast<Gate*>(new U2( qbit_num, target_qbit ));
    
    // adding the operation to the end of the list of gates
    add_gate( operation );
}

void Gates_block::add_u2_to_front(int target_qbit ) {
    
    // create the operation
    Gate* gate = static_cast<Gate*>(new U2( qbit_num, target_qbit ));
    
    // adding the operation to the front of the list of gates
    add_gate_to_front( gate );
}

void Gates_block::add_u3(int target_qbit) {

        // create the operation
        Gate* operation = static_cast<Gate*>(new U3( qbit_num, target_qbit ));

        // adding the operation to the end of the list of gates
        add_gate( operation );
}

void Gates_block::add_u3_to_front(int target_qbit) {

        // create the operation
        Gate* gate = static_cast<Gate*>(new U3( qbit_num, target_qbit ));

        // adding the operation to the front of the list of gates
        add_gate_to_front( gate );

}

void Gates_block::add_rx(int target_qbit) {

        // create the operation
        Gate* operation = static_cast<Gate*>(new RX( qbit_num, target_qbit));

        // adding the operation to the end of the list of gates
        add_gate( operation );
}

void Gates_block::add_rx_to_front(int target_qbit ) {

        // create the operation
        Gate* gate = static_cast<Gate*>(new RX( qbit_num, target_qbit ));

        // adding the operation to the front of the list of gates
        add_gate_to_front( gate );

}

void Gates_block::add_rxx(std::vector<int> target_qbits) {

        // create the operation
        Gate* operation = static_cast<Gate*>(new RXX( qbit_num, target_qbits));

        // adding the operation to the end of the list of gates
        add_gate( operation );
}

void Gates_block::add_rxx_to_front(std::vector<int> target_qbits ) {

        // create the operation
        Gate* gate = static_cast<Gate*>(new RXX( qbit_num, target_qbits ));

        // adding the operation to the front of the list of gates
        add_gate_to_front( gate );

}

void Gates_block::add_ryy(std::vector<int> target_qbits) {

        // create the operation
        Gate* operation = static_cast<Gate*>(new RYY( qbit_num, target_qbits));

        // adding the operation to the end of the list of gates
        add_gate( operation );
}

void Gates_block::add_ryy_to_front(std::vector<int> target_qbits ) {

        // create the operation
        Gate* gate = static_cast<Gate*>(new RYY( qbit_num, target_qbits ));

        // adding the operation to the front of the list of gates
        add_gate_to_front( gate );

}

void Gates_block::add_rzz(std::vector<int> target_qbits) {

        // create the operation
        Gate* operation = static_cast<Gate*>(new RZZ( qbit_num, target_qbits));

        // adding the operation to the end of the list of gates
        add_gate( operation );
}

void Gates_block::add_rzz_to_front(std::vector<int> target_qbits ) {

        // create the operation
        Gate* gate = static_cast<Gate*>(new RZZ( qbit_num, target_qbits ));

        // adding the operation to the front of the list of gates
        add_gate_to_front( gate );

}
void Gates_block::add_r(int target_qbit) {

        // create the operation
        Gate* operation = static_cast<Gate*>(new R( qbit_num, target_qbit));

        // adding the operation to the end of the list of gates
        add_gate( operation );
}

void Gates_block::add_r_to_front(int target_qbit ) {

        // create the operation
        Gate* gate = static_cast<Gate*>(new R( qbit_num, target_qbit ));

        // adding the operation to the front of the list of gates
        add_gate_to_front( gate );

}


void Gates_block::add_ry(int target_qbit) {

        // create the operation
        Gate* operation = static_cast<Gate*>(new RY( qbit_num, target_qbit));

        // adding the operation to the end of the list of gates
        add_gate( operation );
}


void Gates_block::add_ry_to_front(int target_qbit ) {

        // create the operation
        Gate* gate = static_cast<Gate*>(new RY( qbit_num, target_qbit ));

        // adding the operation to the front of the list of gates
        add_gate_to_front( gate );

}




void Gates_block::add_cry(int target_qbit, int control_qbit) {

        // create the operation
        Gate* operation = static_cast<Gate*>(new CRY( qbit_num, target_qbit, control_qbit));

        // adding the operation to the end of the list of gates
        add_gate( operation );
}



void Gates_block::add_cry_to_front(int target_qbit, int control_qbit ) {

        // create the operation
        Gate* gate = static_cast<Gate*>(new CRY( qbit_num, target_qbit, control_qbit ));

        // adding the operation to the front of the list of gates
        add_gate_to_front( gate );

}

void Gates_block::add_crx(int target_qbit, int control_qbit) {

        // create the operation
        Gate* operation = static_cast<Gate*>(new CRX( qbit_num, target_qbit, control_qbit));

        // adding the operation to the end of the list of gates
        add_gate( operation );
}



void Gates_block::add_crx_to_front(int target_qbit, int control_qbit ) {

        // create the operation
        Gate* gate = static_cast<Gate*>(new CRX( qbit_num, target_qbit, control_qbit ));

        // adding the operation to the front of the list of gates
        add_gate_to_front( gate );

}

void Gates_block::add_crz(int target_qbit, int control_qbit) {

        // create the operation
        Gate* operation = static_cast<Gate*>(new CRZ( qbit_num, target_qbit, control_qbit));

        // adding the operation to the end of the list of gates
        add_gate( operation );
}



void Gates_block::add_crz_to_front(int target_qbit, int control_qbit ) {

        // create the operation
        Gate* gate = static_cast<Gate*>(new CRZ( qbit_num, target_qbit, control_qbit ));

        // adding the operation to the front of the list of gates
        add_gate_to_front( gate );

}

void Gates_block::add_cp(int target_qbit, int control_qbit) {

        // create the operation
        Gate* operation = static_cast<Gate*>(new CP( qbit_num, target_qbit, control_qbit));

        // adding the operation to the end of the list of gates
        add_gate( operation );
}



void Gates_block::add_cp_to_front(int target_qbit, int control_qbit ) {

        // create the operation
        Gate* gate = static_cast<Gate*>(new CP( qbit_num, target_qbit, control_qbit ));

        // adding the operation to the front of the list of gates
        add_gate_to_front( gate );

}

void Gates_block::add_cr(int target_qbit, int control_qbit) {

        // create the operation
        Gate* operation = static_cast<Gate*>(new CR( qbit_num, target_qbit, control_qbit));

        // adding the operation to the end of the list of gates
        add_gate( operation );
}



void Gates_block::add_cr_to_front(int target_qbit, int control_qbit ) {

        // create the operation
        Gate* gate = static_cast<Gate*>(new CR( qbit_num, target_qbit, control_qbit ));

        // adding the operation to the front of the list of gates
        add_gate_to_front( gate );

}


void Gates_block::add_crot(int target_qbit, int control_qbit) {

        // create the operation
        Gate* operation = static_cast<Gate*>(new CROT( qbit_num, target_qbit, control_qbit));

        // adding the operation to the end of the list of gates
        add_gate( operation );
}



void Gates_block::add_crot_to_front(int target_qbit, int control_qbit ) {

        // create the operation
        Gate* gate = static_cast<Gate*>(new CROT( qbit_num, target_qbit, control_qbit ));

        // adding the operation to the front of the list of gates
        add_gate_to_front( gate );

}

void Gates_block::add_rz(int target_qbit) {

        // create the operation
        Gate* operation = static_cast<Gate*>(new RZ( qbit_num, target_qbit));

        // adding the operation to the end of the list of gates
        add_gate( operation );
}

void Gates_block::add_rz_to_front(int target_qbit ) {

        // create the operation
        Gate* gate = static_cast<Gate*>(new RZ( qbit_num, target_qbit ));

        // adding the operation to the front of the list of gates
        add_gate_to_front( gate );

}


void Gates_block::add_cnot(  int target_qbit, int control_qbit) {

        // new cnot operation
        Gate* gate = static_cast<Gate*>(new CNOT(qbit_num, target_qbit, control_qbit ));

        // append the operation to the list
        add_gate(gate);

}



void Gates_block::add_cnot_to_front( int target_qbit, int control_qbit) {

        // new cnot operation
        Gate* gate = static_cast<Gate*>(new CNOT(qbit_num, target_qbit, control_qbit ));

        // put the operation to tghe front of the list
        add_gate_to_front(gate);

}




void Gates_block::add_cz(  int target_qbit, int control_qbit) {

        // new cz operation
        Gate* gate = static_cast<Gate*>(new CZ(qbit_num, target_qbit, control_qbit ));

        // append the operation to the list
        add_gate(gate);

}



void Gates_block::add_cz_to_front(  int target_qbit, int control_qbit) {

        // new cz operation
        Gate* gate = static_cast<Gate*>(new CZ(qbit_num, target_qbit, control_qbit ));

        // put the operation to tghe front of the list
        add_gate_to_front(gate);

}

void Gates_block::add_h(int target_qbit) {

        // create the operation
        Gate* operation = static_cast<Gate*>(new H( qbit_num, target_qbit));

        // adding the operation to the end of the list of gates
        add_gate( operation );
}

void Gates_block::add_h_to_front(int target_qbit ) {

        // create the operation
        Gate* gate = static_cast<Gate*>(new H( qbit_num, target_qbit ));

        // adding the operation to the front of the list of gates
        add_gate_to_front( gate );

}





void Gates_block::add_x(int target_qbit) {

        // create the operation
        Gate* operation = static_cast<Gate*>(new X( qbit_num, target_qbit));

        // adding the operation to the end of the list of gates
        add_gate( operation );
}

void Gates_block::add_x_to_front(int target_qbit ) {

        // create the operation
        Gate* gate = static_cast<Gate*>(new X( qbit_num, target_qbit ));

        // adding the operation to the front of the list of gates
        add_gate_to_front( gate );

}



void Gates_block::add_y(int target_qbit) {

        // create the operation
        Gate* operation = static_cast<Gate*>(new Y( qbit_num, target_qbit));

        // adding the operation to the end of the list of gates
        add_gate( operation );
}

void Gates_block::add_y_to_front(int target_qbit ) {

        // create the operation
        Gate* gate = static_cast<Gate*>(new Y( qbit_num, target_qbit ));

        // adding the operation to the front of the list of gates
        add_gate_to_front( gate );

}



void Gates_block::add_z(int target_qbit) {

        // create the operation
        Gate* operation = static_cast<Gate*>(new Z( qbit_num, target_qbit));

        // adding the operation to the end of the list of gates
        add_gate( operation );
}

void Gates_block::add_z_to_front(int target_qbit ) {

        // create the operation
        Gate* gate = static_cast<Gate*>(new Z( qbit_num, target_qbit ));

        // adding the operation to the front of the list of gates
        add_gate_to_front( gate );

}




void Gates_block::add_s(int target_qbit) {

        // create the operation
        Gate* operation = static_cast<Gate*>(new S( qbit_num, target_qbit));

        // adding the operation to the end of the list of gates
        add_gate( operation );
}

void Gates_block::add_s_to_front(int target_qbit ) {

        // create the operation
        Gate* gate = static_cast<Gate*>(new S( qbit_num, target_qbit ));

        // adding the operation to the front of the list of gates
        add_gate_to_front( gate );

}



void Gates_block::add_sdg(int target_qbit) {

        // create the operation
        Gate* operation = static_cast<Gate*>(new SDG( qbit_num, target_qbit));

        // adding the operation to the end of the list of gates
        add_gate( operation );
}

void Gates_block::add_sdg_to_front(int target_qbit ) {

        // create the operation
        Gate* gate = static_cast<Gate*>(new SDG( qbit_num, target_qbit ));

        // adding the operation to the front of the list of gates
        add_gate_to_front( gate );

}




void Gates_block::add_t(int target_qbit) {

        // create the operation
        Gate* operation = static_cast<Gate*>(new T( qbit_num, target_qbit));

        // adding the operation to the end of the list of gates
        add_gate( operation );
}

void Gates_block::add_t_to_front(int target_qbit ) {

        // create the operation
        Gate* gate = static_cast<Gate*>(new T( qbit_num, target_qbit ));

        // adding the operation to the front of the list of gates
        add_gate_to_front( gate );

}




void Gates_block::add_tdg(int target_qbit) {

        // create the operation
        Gate* operation = static_cast<Gate*>(new Tdg( qbit_num, target_qbit));

        // adding the operation to the end of the list of gates
        add_gate( operation );
}

void Gates_block::add_tdg_to_front(int target_qbit ) {

        // create the operation
        Gate* gate = static_cast<Gate*>(new Tdg( qbit_num, target_qbit ));

        // adding the operation to the front of the list of gates
        add_gate_to_front( gate );

}






void Gates_block::add_sx(int target_qbit) {

        // create the operation
        Gate* operation = static_cast<Gate*>(new SX( qbit_num, target_qbit));

        // adding the operation to the end of the list of gates
        add_gate( operation );
}

void Gates_block::add_sx_to_front(int target_qbit ) {

        // create the operation
        Gate* gate = static_cast<Gate*>(new SX( qbit_num, target_qbit ));

        // adding the operation to the front of the list of gates
        add_gate_to_front( gate );

}


void Gates_block::add_sxdg(int target_qbit) {

        // create the operation
        Gate* operation = static_cast<Gate*>(new SXdg( qbit_num, target_qbit));

        // adding the operation to the end of the list of gates
        add_gate( operation );
}

void Gates_block::add_sxdg_to_front(int target_qbit ) {

        // create the operation
        Gate* gate = static_cast<Gate*>(new SXdg( qbit_num, target_qbit ));

        // adding the operation to the front of the list of gates
        add_gate_to_front( gate );

}


void Gates_block::add_syc(  int target_qbit, int control_qbit) {

        // new cz operation
        Gate* gate = static_cast<Gate*>(new SYC(qbit_num, target_qbit, control_qbit ));

        // append the operation to the list
        add_gate(gate);

}



void Gates_block::add_syc_to_front(  int target_qbit, int control_qbit) {

        // new cz operation
        Gate* gate = static_cast<Gate*>(new SYC(qbit_num, target_qbit, control_qbit ));

        // put the operation to tghe front of the list
        add_gate_to_front(gate);

}




void Gates_block::add_ch(  int target_qbit, int control_qbit) {

        // new cz operation
        Gate* gate = static_cast<Gate*>(new CH(qbit_num, target_qbit, control_qbit ));

        // append the operation to the list
        add_gate(gate);

}



void Gates_block::add_ch_to_front(  int target_qbit, int control_qbit) {

        // new cz operation
        Gate* gate = static_cast<Gate*>(new CH(qbit_num, target_qbit, control_qbit ));

        // put the operation to tghe front of the list
        add_gate_to_front(gate);

}

void Gates_block::add_cu(  int target_qbit, int control_qbit) {

        // new cz operation
        Gate* gate = static_cast<Gate*>(new CU(qbit_num, target_qbit, control_qbit ));

        // append the operation to the list
        add_gate(gate);

}



void Gates_block::add_cu_to_front(  int target_qbit, int control_qbit) {

        // new cz operation
        Gate* gate = static_cast<Gate*>(new CU(qbit_num, target_qbit, control_qbit ));

        // put the operation to tghe front of the list
        add_gate_to_front(gate);

}

void Gates_block::add_ccx( int target_qbit, const std::vector<int>& control_qbits) {
        // new ccx operation
        Gate* gate = static_cast<Gate*>(new CCX(qbit_num, target_qbit, control_qbits));
        // append the operation to the list
        add_gate(gate);
}

void Gates_block::add_ccx_to_front( int target_qbit, const std::vector<int>& control_qbits) {
        // new ccx operation
        Gate* gate = static_cast<Gate*>(new CCX(qbit_num, target_qbit, control_qbits));
        // put the operation to the front of the list
        add_gate_to_front(gate);
}

void Gates_block::add_swap( const std::vector<int>& target_qbits) {
        // new swap operation
        Gate* gate = static_cast<Gate*>(new SWAP(qbit_num, target_qbits));
        // append the operation to the list
        add_gate(gate);
}

void Gates_block::add_swap_to_front( const std::vector<int>& target_qbits) {
        // new swap operation
        Gate* gate = static_cast<Gate*>(new SWAP(qbit_num, target_qbits));
        // put the operation to the front of the list
        add_gate_to_front(gate);
}

void Gates_block::add_cswap( const std::vector<int>& target_qbits, const std::vector<int>& control_qbits) {
        // new cswap operation
        Gate* gate = static_cast<Gate*>(new CSWAP(qbit_num, target_qbits, control_qbits));
        // append the operation to the list
        add_gate(gate);
}

void Gates_block::add_cswap_to_front( const std::vector<int>& target_qbits, const std::vector<int>& control_qbits) {
        // new cswap operation
        Gate* gate = static_cast<Gate*>(new CSWAP(qbit_num, target_qbits, control_qbits));
        // put the operation to the front of the list
        add_gate_to_front(gate);
}

void Gates_block::add_gates( std::vector<Gate*> gates_in) {

        for(std::vector<Gate*>::iterator it = gates_in.begin(); it != gates_in.end(); ++it) {
            add_gate( *it );
        }

}


void Gates_block::add_gates_to_front( std::vector<Gate*>  gates_in) {

        // adding gates in reversed order!!
        for(std::vector<Gate*>::iterator it = gates_in.end(); it != gates_in.begin(); --it) {
            add_gate_to_front( *it );
        }

}



void Gates_block::add_adaptive(int target_qbit, int control_qbit)  {

        // create the operation
        Gate* operation = static_cast<Gate*>(new Adaptive( qbit_num, target_qbit, control_qbit));

        // adding the operation to the end of the list of gates
        add_gate( operation );
}


void Gates_block::add_adaptive_to_front(int target_qbit, int control_qbit)  {

        // create the operation
        Gate* gate = static_cast<Gate*>(new Adaptive( qbit_num, target_qbit, control_qbit ));

        // adding the operation to the front of the list of gates
        add_gate_to_front( gate );

}


void Gates_block::add_general_operation(
    const Matrix& operation_mtx,
    const std::vector<int>& target_qbits,
    const std::vector<int>& control_qbits
) {

    if (operation_mtx.rows != matrix_size || operation_mtx.cols != matrix_size) {
        std::string err("Gates_block::add_general_operation: operation matrix has invalid size.");
        throw err;
    }

    Gate* operation = new Gate(qbit_num);
    operation->set_matrix(operation_mtx);
    if (!target_qbits.empty()) {
        operation->set_target_qbits(target_qbits);
    }
    if (!control_qbits.empty()) {
        operation->set_control_qbits(control_qbits);
    }
    add_gate(operation);
}



void Gates_block::add_gate( Gate* gate ) {

        //set the number of qubit in the gate
        gate->set_qbit_num( qbit_num );

        gate->clear_children();
        gate->clear_parents();
        // determine the parents of the gate
        determine_parents( gate );

        // append the gate to the list
        gates.push_back(gate);
        invalidate_structure_cache();
        
        // set the parameter starting index in the parameters array used to execute the circuit.
        gate->set_parameter_start_idx( parameter_num );

        // increase the number of parameters by the number of parameters
        parameter_num = parameter_num + gate->get_parameter_num();

        // increase the number of layers if necessary
        if (gate->get_type() == BLOCK_OPERATION) {
            layer_num = layer_num + 1;
        }

}



 void Gates_block::add_gate_to_front( Gate* gate) {


        // set the number of qubit in the gate
        gate->set_qbit_num( qbit_num );
        gate->clear_children();
        gate->clear_parents();

        // determine the parents of the gate
        determine_children( gate );

        gates.insert( gates.begin(), gate);
        invalidate_structure_cache();

        // increase the number of U3 gate parameters by the number of parameters
        parameter_num = parameter_num + gate->get_parameter_num();

        // increase the number of layers if necessary
        if (gate->get_type() == BLOCK_OPERATION) {
            layer_num = layer_num + 1;
        }
        
        // TODO: develop a more efficient method for large circuits. Now it is updating the whole circuit
        reset_parameter_start_indices();
        reset_dependency_graph();

}



void 
Gates_block::insert_gate( Gate* gate, int idx ) {


        // set the number of qubit in the gate
        gate->set_qbit_num( qbit_num );

        gates.insert( gates.begin()+idx, gate);
        invalidate_structure_cache();

        // increase the number of U3 gate parameters by the number of parameters
        parameter_num = parameter_num + gate->get_parameter_num();

        // increase the number of layers if necessary
        if (gate->get_type() == BLOCK_OPERATION) {
            layer_num = layer_num + 1;
        }
        
        // TODO: develop a more efficient method for large circuits. Now it is updating the whole circuit
        reset_parameter_start_indices();
        reset_dependency_graph();


}


void Gates_block::add_gate_nums( std::map<std::string, int>& gate_nums ) {

    for(std::vector<Gate*>::iterator it = gates.begin(); it != gates.end(); ++it) {
    
        // get the specific gate in the circuit (might be a gate or another subcircuit)
        Gate* gate = *it;
            
            
        if (gate->get_type() == BLOCK_OPERATION) {
            Gates_block* circuit = static_cast<Gates_block*>(gate);
            circuit->add_gate_nums( gate_nums );
        }
        else {
            std::string gate_name = gate->get_name();
            
            if( gate_nums.find(gate_name) == gate_nums.end() ) {
                gate_nums[ gate_name ] = 1;
            }
            else {
                gate_nums[ gate_name ] = gate_nums[ gate_name ] + 1;
            }
        }
    
    }
    
}


std::map<std::string, int> Gates_block::get_gate_nums() {

    std::map<std::string, int> gate_nums;
    
    add_gate_nums( gate_nums );
    
    return gate_nums;

}


int Gates_block::get_parameter_num() {
    return parameter_num;
}




int Gates_block::get_gate_num() {
    return static_cast<int>(gates.size());
}


void Gates_block::list_gates( const Matrix_real &parameters, int start_index ) {
//TODO: rewrite the function with switch statement

     //The stringstream input to store the output messages.
     std::stringstream sstream;
     sstream << std::endl << "The gates in the list of gates:" << std::endl;
     print(sstream, 1);       
          
    int gate_idx = start_index;
    int parameter_idx = 0;
     double *parameters_data = parameters.get_data();
     //const_cast <Matrix_real&>(parameters);


    for(size_t op_idx = 0; op_idx<gates.size(); op_idx++) {

        Gate* gate = gates[op_idx];

        if (gate->get_type() == CNOT_OPERATION) {
            CNOT* cnot_gate = static_cast<CNOT*>(gate);
            std::stringstream sstream;
          sstream << gate_idx << "th gate: CNOT with control qubit: " << cnot_gate->get_control_qbit() << " and target qubit: " << cnot_gate->get_target_qbit() << std::endl;
              print(sstream, 1);             
            gate_idx = gate_idx + 1;
        }
        else if (gate->get_type() == CROT_OPERATION) {
            CROT* crot_gate = static_cast<CROT*>(gate);
            std::stringstream sstream;
            double theta0,phi0;
            theta0 = std::fmod( 2*parameters_data[parameter_idx], 4*M_PI);
            phi0 = std::fmod( parameters_data[parameter_idx+1], 2*M_PI);
            parameter_idx = parameter_idx + 2;
            sstream << gate_idx << "th gate: CROT with control qubit: " << crot_gate->get_control_qbit() << " and target qubit: " << crot_gate->get_target_qbit()<< " and parameters theta=" << theta0 <<" and phi="<< phi0 << std::endl;
              print(sstream, 1);
            gate_idx = gate_idx + 1;
        }
        else if (gate->get_type() == CZ_OPERATION) {
            CZ* cz_gate = static_cast<CZ*>(gate);
              std::stringstream sstream;
              sstream << gate_idx << "th gate: CZ with control qubit: " << cz_gate->get_control_qbit() << " and target qubit: " << cz_gate->get_target_qbit() << std::endl;
              print(sstream, 1);                     
            gate_idx = gate_idx + 1;
        }
        else if (gate->get_type() == CH_OPERATION) {
            CH* ch_gate = static_cast<CH*>(gate);
              std::stringstream sstream;
              sstream << gate_idx << "th gate: CH with control qubit: " << ch_gate->get_control_qbit() << " and target qubit: " << ch_gate->get_target_qbit() << std::endl;
              print(sstream, 1);        
            gate_idx = gate_idx + 1;
        }
        else if (gate->get_type() == SYC_OPERATION) {
            SYC* syc_gate = static_cast<SYC*>(gate);
              std::stringstream sstream;
              sstream << gate_idx << "th gate: Sycamore gate with control qubit: " << syc_gate->get_control_qbit() << " and target qubit: " << syc_gate->get_target_qbit() << std::endl;
              print(sstream, 1);        
            gate_idx = gate_idx + 1;
        }
        else if (gate->get_type() == U1_OPERATION) {
            U1* u1_gate = static_cast<U1*>(gate);
            double lambda = std::fmod(parameters_data[parameter_idx], 2*M_PI);
            parameter_idx = parameter_idx + 1;
            std::stringstream sstream;
            sstream << gate_idx << "th gate: U1 on target qubit: " << u1_gate->get_target_qbit() << " with parameter lambda = " << lambda << std::endl;
            print(sstream, 1);
            gate_idx = gate_idx + 1;
        }
        else if (gate->get_type() == U2_OPERATION) {
            U2* u2_gate = static_cast<U2*>(gate);
            double phi = std::fmod(parameters_data[parameter_idx], 2*M_PI);
            double lambda = std::fmod(parameters_data[parameter_idx+1], 2*M_PI);
            parameter_idx = parameter_idx + 2;
            std::stringstream sstream;
            sstream << gate_idx << "th gate: U2 on target qubit: " << u2_gate->get_target_qbit() << " with parameters phi = " << phi << " and lambda = " << lambda << std::endl;
            print(sstream, 1);
            gate_idx = gate_idx + 1;
        }
        else if (gate->get_type() == U3_OPERATION) {
            U3* u3_gate = static_cast<U3*>(gate);
            double theta = std::fmod(parameters_data[parameter_idx], 4*M_PI);
            double phi = std::fmod(parameters_data[parameter_idx+1], 2*M_PI);
            double lambda = std::fmod(parameters_data[parameter_idx+2], 2*M_PI);
            parameter_idx = parameter_idx + 3;
            std::stringstream sstream;
            sstream << gate_idx << "th gate: U3 on target qubit: " << u3_gate->get_target_qbit() << " with parameters theta = " << theta << ", phi = " << phi << " and lambda = " << lambda << std::endl;
            print(sstream, 1);
            gate_idx = gate_idx + 1;
        }
        else if (gate->get_type() == R_OPERATION) {
             // definig the rotation parameter
            double vartheta,varphi;
            // get the inverse parameters of the U3 rotation
            R* r_gate = static_cast<R*>(gate);         
            vartheta = std::fmod( 2*parameters_data[parameter_idx], 4*M_PI);
            varphi = std::fmod( parameters_data[parameter_idx+1], 2*M_PI);
            parameter_idx = parameter_idx + 2;

              std::stringstream sstream;
              sstream << gate_idx << "th gate: R on target qubit: " << r_gate->get_target_qbit() << " and with parameters theta = " << vartheta << " and phi:" << varphi << std::endl;
              print(sstream, 1);        
            gate_idx = gate_idx + 1;
        }
        else if (gate->get_type() == RX_OPERATION) {
             // definig the rotation parameter
                double vartheta;
                // get the inverse parameters of the U3 rotation
                RX* rx_gate = static_cast<RX*>(gate);       
                vartheta = std::fmod( 2*parameters_data[parameter_idx], 4*M_PI);
                parameter_idx = parameter_idx + 1;

          std::stringstream sstream;
          sstream << gate_idx << "th gate: RX on target qubit: " << rx_gate->get_target_qbit() << " and with parameters theta = " << vartheta << std::endl;
          print(sstream, 1);       
                gate_idx = gate_idx + 1;
            }
        else if (gate->get_type() == RXX_OPERATION) {
             // definig the rotation parameter
                double vartheta;
                // get the inverse parameters of the U3 rotation
                RXX* rxx_gate = static_cast<RXX*>(gate);
                vartheta = std::fmod( 2*parameters_data[parameter_idx], 4*M_PI);
                parameter_idx = parameter_idx + 1;
                std::vector<int> involved_qbits = rxx_gate->get_involved_qubits(true);
                std::stringstream sstream;
                sstream << gate_idx << "th gate: RXX on target qubits: " << involved_qbits[0] <<" and " << involved_qbits[1] << " and with parameters theta = " << vartheta << std::endl;
                print(sstream, 1);
                gate_idx = gate_idx + 1;
            }
        else if (gate->get_type() == RYY_OPERATION) {
             // definig the rotation parameter
                double vartheta;
                // get the inverse parameters of the U3 rotation
                RYY* ryy_gate = static_cast<RYY*>(gate);
                vartheta = std::fmod( 2*parameters_data[parameter_idx], 4*M_PI);
                parameter_idx = parameter_idx + 1;
                std::vector<int> involved_qbits = ryy_gate->get_involved_qubits(true);
                std::stringstream sstream;
                sstream << gate_idx << "th gate: RYY on target qubits: " << involved_qbits[0] <<" and " << involved_qbits[1] << " and with parameters theta = " << vartheta << std::endl;
                print(sstream, 1);
                gate_idx = gate_idx + 1;
            }
        else if (gate->get_type() == RZZ_OPERATION) {
             // definig the rotation parameter
                double vartheta;
                // get the inverse parameters of the U3 rotation
                RZZ* rzz_gate = static_cast<RZZ*>(gate);
                vartheta = std::fmod( 2*parameters_data[parameter_idx], 4*M_PI);
                parameter_idx = parameter_idx + 1;
                std::vector<int> involved_qbits = rzz_gate->get_involved_qubits(true);
                std::stringstream sstream;
                sstream << gate_idx << "th gate: RZZ on target qubits: " << involved_qbits[0] <<" and " << involved_qbits[1] << " and with parameters theta = " << vartheta << std::endl;
                print(sstream, 1);
                gate_idx = gate_idx + 1;
            }
        else if (gate->get_type() == RY_OPERATION) {
                // definig the rotation parameter
                double vartheta;
                // get the inverse parameters of the U3 rotation
                RY* ry_gate = static_cast<RY*>(gate);
                vartheta = std::fmod( 2*parameters_data[parameter_idx], 4*M_PI);
                parameter_idx = parameter_idx + 1;

          std::stringstream sstream;
          sstream << gate_idx << "th gate: RY on target qubit: " << ry_gate->get_target_qbit() << " and with parameters theta = " << vartheta << std::endl; 
          print(sstream, 1);       
                gate_idx = gate_idx + 1;
            }
            else if (gate->get_type() == CRY_OPERATION) {
                // definig the rotation parameter
                double vartheta;
                // get the inverse parameters of the U3 rotation
                CRY* cry_gate = static_cast<CRY*>(gate);
                vartheta = std::fmod( 2*parameters_data[parameter_idx], 4*M_PI);
                parameter_idx = parameter_idx + 1;

          std::stringstream sstream;
          sstream << gate_idx << "th gate: CRY on target qubit: " << cry_gate->get_target_qbit() << ", control qubit" << cry_gate->get_control_qbit() << " and with parameters theta = " << vartheta << std::endl;
          print(sstream, 1);                                
                gate_idx = gate_idx + 1;
            }
            else if (gate->get_type() == CRX_OPERATION) {
                // definig the rotation parameter
                double vartheta;
                // get the inverse parameters of the U3 rotation
                CRX* crx_gate = static_cast<CRX*>(gate);
                vartheta = std::fmod( 2*parameters_data[parameter_idx], 4*M_PI);
                parameter_idx = parameter_idx + 1;

          std::stringstream sstream;
          sstream << gate_idx << "th gate: CRX on target qubit: " << crx_gate->get_target_qbit() << ", control qubit" << crx_gate->get_control_qbit() << " and with parameters theta = " << vartheta << std::endl;
          print(sstream, 1);                                
                gate_idx = gate_idx + 1;
            }
            else if (gate->get_type() == CRZ_OPERATION) {
                // definig the rotation parameter
                double vartheta;
                // get the inverse parameters of the U3 rotation
                CRZ* crz_gate = static_cast<CRZ*>(gate);
                vartheta = std::fmod( 2*parameters_data[parameter_idx], 4*M_PI);
                parameter_idx = parameter_idx + 1;

          std::stringstream sstream;
          sstream << gate_idx << "th gate: CRZ on target qubit: " << crz_gate->get_target_qbit() << ", control qubit" << crz_gate->get_control_qbit() << " and with parameters theta = " << vartheta << std::endl;
          print(sstream, 1);                                
                gate_idx = gate_idx + 1;
            }
            else if (gate->get_type() == CP_OPERATION) {
                // definig the rotation parameter
                double lambda;
                // get the inverse parameters of the U3 rotation
                CP* cp_gate = static_cast<CP*>(gate);
                lambda = std::fmod( parameters_data[parameter_idx], 2*M_PI);
                parameter_idx = parameter_idx + 1;

          std::stringstream sstream;
          sstream << gate_idx << "th gate: CP on target qubit: " << cp_gate->get_target_qbit() << ", control qubit" << cp_gate->get_control_qbit() << " and with parameters lambda = " << lambda << std::endl;
          print(sstream, 1);                                
                gate_idx = gate_idx + 1;
            }
            else if (gate->get_type() == CR_OPERATION) {
                // definig the rotation parameter
                double vartheta,varphi;
                // get the inverse parameters of the U3 rotation
                CR* cr_gate = static_cast<CR*>(gate);
                vartheta = std::fmod( 2*parameters_data[parameter_idx], 4*M_PI);
                varphi = std::fmod(parameters_data[parameter_idx+1],2*M_PI);
                parameter_idx = parameter_idx + 2;

          std::stringstream sstream;
          sstream << gate_idx << "th gate: CR on target qubit: " << cr_gate->get_target_qbit() << ", control qubit" << cr_gate->get_control_qbit() << " and with parameters theta and phi = " << vartheta << " , " << varphi << std::endl;
          print(sstream, 1);                                
                gate_idx = gate_idx + 1;
            }
            else if (gate->get_type() == RZ_OPERATION) {
                // definig the rotation parameter
                double varphi;
                // get the inverse parameters of the U3 rotation
                RZ* rz_gate = static_cast<RZ*>(gate);
                varphi = std::fmod( 2*parameters_data[parameter_idx], 2*M_PI);
                parameter_idx = parameter_idx + 1;

          std::stringstream sstream;
          sstream << gate_idx << "th gate: RZ on target qubit: " << rz_gate->get_target_qbit() << " and with parameters varphi = " << varphi << std::endl;
          print(sstream, 1);       
                gate_idx = gate_idx + 1;
            }
            else if (gate->get_type() == H_OPERATION) {
                // get the inverse parameters of the U3 rotation
                H* h_gate = static_cast<H*>(gate);
          std::stringstream sstream;
          sstream << gate_idx << "th gate: Hadamard on target qubit: " << h_gate->get_target_qbit() << std::endl;
          print(sstream, 1);       
                gate_idx = gate_idx + 1;
            }
            else if (gate->get_type() == X_OPERATION) {
                // get the inverse parameters of the U3 rotation
                X* x_gate = static_cast<X*>(gate);
          std::stringstream sstream;
          sstream << gate_idx << "th gate: X on target qubit: " << x_gate->get_target_qbit() << std::endl;
          print(sstream, 1);       
                gate_idx = gate_idx + 1;
            }
            else if (gate->get_type() == Y_OPERATION) {
                // get the inverse parameters of the U3 rotation
                Y* y_gate = static_cast<Y*>(gate);
          std::stringstream sstream;
          sstream << gate_idx << "th gate: Y on target qubit: " << y_gate->get_target_qbit() << std::endl;
          print(sstream, 1);       
                gate_idx = gate_idx + 1;
            }
            else if (gate->get_type() == Z_OPERATION) {
                // get the inverse parameters of the U3 rotation
                Z* z_gate = static_cast<Z*>(gate);
          std::stringstream sstream;
          sstream << gate_idx << "th gate: Z on target qubit: " << z_gate->get_target_qbit() << std::endl;
          print(sstream, 1);       
                gate_idx = gate_idx + 1;
            }
            else if (gate->get_type() == S_OPERATION) {
                // get the inverse parameters of the U3 rotation
                S* s_gate = static_cast<S*>(gate);
          std::stringstream sstream;
          sstream << gate_idx << "th gate: S on target qubit: " << s_gate->get_target_qbit() << std::endl;
          print(sstream, 1);       
                gate_idx = gate_idx + 1;
            }
            else if (gate->get_type() == SDG_OPERATION) {
                // get the inverse parameters of the U3 rotation
                SDG* sdg_gate = static_cast<SDG*>(gate);
          std::stringstream sstream;
          sstream << gate_idx << "th gate: Sdg on target qubit: " << sdg_gate->get_target_qbit() << std::endl;
          print(sstream, 1);       
                gate_idx = gate_idx + 1;
            }
            else if (gate->get_type() == SX_OPERATION) {
                // get the inverse parameters of the U3 rotation
                SX* sx_gate = static_cast<SX*>(gate);

          std::stringstream sstream;
          sstream << gate_idx << "th gate: SX on target qubit: " << sx_gate->get_target_qbit() << std::endl;
          print(sstream, 1);       
                gate_idx = gate_idx + 1;
            }
            else if (gate->get_type() == SXDG_OPERATION) {
                // get the inverse parameters of the U3 rotation
                SXdg* sxdg_gate = static_cast<SXdg*>(gate);

          std::stringstream sstream;
          sstream << gate_idx << "th gate: SXdg on target qubit: " << sxdg_gate->get_target_qbit() << std::endl;
          print(sstream, 1);       
                gate_idx = gate_idx + 1;
            }
            else if (gate->get_type() == BLOCK_OPERATION) {
                Gates_block* block_gate = static_cast<Gates_block*>(gate);
                const Matrix_real parameters_layer(parameters.get_data() + parameter_idx, 1, gate->get_parameter_num() );
                block_gate->list_gates( parameters_layer, gate_idx );
                parameter_idx = parameter_idx + block_gate->get_parameter_num();
                gate_idx = gate_idx + block_gate->get_gate_num();
            }
            else if (gate->get_type() == ADAPTIVE_OPERATION) {
                // definig the rotation parameter
                double Theta;
                // get the inverse parameters of the U3 rotation
                Adaptive* ad_gate = static_cast<Adaptive*>(gate);
                Theta = std::fmod( 2*parameters_data[parameter_idx], 4*M_PI);
                parameter_idx = parameter_idx + 1;

          std::stringstream sstream;
          sstream << gate_idx << "th gate: Adaptive gate on target qubit: " << ad_gate->get_target_qbit() << ", control qubit " << ad_gate->get_control_qbit() << " and with parameters Theta = " << Theta << std::endl;
          print(sstream, 1);       
                gate_idx = gate_idx + 1;
            }
            else {
                std::string err("Gates_block::list_gates: unimplemented gate"); 
                throw err;
            }

        }


}


Gates_block* 
Gates_block::create_remapped_circuit( const std::map<int, int>& qbit_map ) {

    return create_remapped_circuit( qbit_map, qbit_num );

}


Gates_block* 
Gates_block::create_remapped_circuit( const std::map<int, int>& qbit_map, const int qbit_num_ ) {


    Gates_block* ret = new Gates_block( qbit_num_ );
    for(std::vector<Gate*>::iterator it = gates.begin(); it != gates.end(); ++it) {

        Gate* op = *it;
        switch (op->get_type()) {
        case CNOT_OPERATION: case CZ_OPERATION:
        case CH_OPERATION: case SYC_OPERATION:
        case U1_OPERATION: case U2_OPERATION:
        case U3_OPERATION: case RY_OPERATION:
        case CRY_OPERATION: case RX_OPERATION:
        case CRX_OPERATION: case CRZ_OPERATION:
        case CP_OPERATION: case SWAP_OPERATION:
        case CSWAP_OPERATION: case CCX_OPERATION:
        case RZ_OPERATION: case X_OPERATION:
        case Y_OPERATION: case Z_OPERATION:
        case SX_OPERATION: case BLOCK_OPERATION:
        case GENERAL_OPERATION:
        case ADAPTIVE_OPERATION: case H_OPERATION:
        case S_OPERATION: case SDG_OPERATION:
        case T_OPERATION: case TDG_OPERATION:
        case CU_OPERATION:
        case RXX_OPERATION: case RYY_OPERATION:
        case RZZ_OPERATION: case SXDG_OPERATION:
        {
            Gate* cloned_op = op->clone();

            if ( qbit_num_ > qbit_num ) {
                // qbit num needs to be set in prior to avoid conflict
                cloned_op->set_qbit_num( qbit_num_ );
            }

            // Handle target qubits (both single and vector-based)
            std::vector<int> target_qbits = cloned_op->get_target_qbits();
            if (!target_qbits.empty()) {
                // Vector-based targets (SWAP, CSWAP, etc.)
                std::vector<int> remapped_targets;
                for (int target_qbit : target_qbits) {
                    if (qbit_map.find(target_qbit) != qbit_map.end()) {
                        remapped_targets.push_back(qbit_map.at(target_qbit));
                    } else {
                        std::string err("Gates_block::create_remapped_circuit: Missing target qubit from the qbit map.");
                        throw err;
                    }
                }
                cloned_op->set_target_qbits(remapped_targets);
            } else {
                // Legacy single target qubit
                int target_qbit = cloned_op->get_target_qbit();
                if (target_qbit != -1) {
                    if (qbit_map.find(target_qbit) != qbit_map.end()) {
                        cloned_op->set_target_qbit(qbit_map.at(target_qbit));
                    } else {
                        std::string err("Gates_block::create_remapped_circuit: Missing target qubit from the qbit map.");
                        throw err;
                    }
                }
            }

            // Handle control qubits (both single and vector-based)
            std::vector<int> control_qbits = cloned_op->get_control_qbits();
            if (!control_qbits.empty()) {
                // Vector-based controls (CCX, CSWAP, etc.)
                std::vector<int> remapped_controls;
                for (int control_qbit : control_qbits) {
                    if (qbit_map.find(control_qbit) != qbit_map.end()) {
                        remapped_controls.push_back(qbit_map.at(control_qbit));
                    } else {
                        std::string err("Gates_block::create_remapped_circuit: Missing control qubit from the qbit map.");
                        throw err;
                    }
                }
                cloned_op->set_control_qbits(remapped_controls);
            } else {
                // Legacy single control qubit
                int control_qbit = cloned_op->get_control_qbit();
                if (control_qbit != -1) {
                    if (qbit_map.find(control_qbit) != qbit_map.end()) {
                        cloned_op->set_control_qbit(qbit_map.at(control_qbit));
                    } else {
                        std::string err("Gates_block::create_remapped_circuit: Missing control qubit from the qbit map.");
                        throw err;
                    }
                }
            }

            if ( qbit_num_ < qbit_num ) {
                // qbit num needs to be set post to avoid conflict
                cloned_op->set_qbit_num( qbit_num_ );
            }

            ret->add_gate( cloned_op );

            break;
        }
        default:
            std::string err("Gates_block::create_remapped_circuit: unimplemented gate"); 
            throw err;
        }

    }

    return ret;
}


void Gates_block::reorder_qubits( std::vector<int>  qbit_list ) {

    for(std::vector<Gate*>::iterator it = gates.begin(); it != gates.end(); ++it) {

        Gate* gate = *it;

        gate->reorder_qubits( qbit_list );


    }

    invalidate_structure_cache();
}

std::vector<int> Gates_block::get_involved_qubits(bool only_target) {

    {
        tbb::spin_mutex::scoped_lock lock(*involved_qubits_cache_mutex);
        if (only_target && involved_target_qubits_cache_valid) {
            return involved_target_qubits_cache;
        }
        if (!only_target && involved_qubits_cache_valid) {
            return involved_qubits_cache;
        }
    }

    std::vector<int> involved_qbits;
    std::vector<char> seen(qbit_num > 0 ? qbit_num : 0, 0);

    for(std::vector<Gate*>::iterator it = gates.begin(); it != gates.end(); ++it) {

        Gate* gate = *it;

        // Use the gate's get_involved_qubits() method which properly handles
        // both single and vector-based qubits (SWAP, CSWAP, CCX, etc.)
        std::vector<int> gate_qbits = gate->get_involved_qubits(only_target);
        for (std::vector<int>::iterator in_it = gate_qbits.begin(); in_it != gate_qbits.end(); ++in_it){
            const int qbit = *in_it;
            if (qbit < 0) {
                continue;
            }
            if (qbit >= static_cast<int>(seen.size())) {
                if (std::find(involved_qbits.begin(), involved_qbits.end(), qbit) == involved_qbits.end()) {
                    involved_qbits.push_back(qbit);
                }
            }
            else if (!seen[qbit]) {
                seen[qbit] = 1;
                involved_qbits.push_back(qbit);
            }
        }

    }

    std::sort(involved_qbits.begin(), involved_qbits.end());

    {
        tbb::spin_mutex::scoped_lock lock(*involved_qubits_cache_mutex);
        if (only_target) {
            involved_target_qubits_cache = involved_qbits;
            involved_target_qubits_cache_valid = true;
        }
        else {
            involved_qubits_cache = involved_qbits;
            involved_qubits_cache_valid = true;
        }
    }

    return involved_qbits;
}


std::vector<Gate*> Gates_block::get_gates() {
    return gates;
}


Gate* Gates_block::get_gate(int idx) {

    if (idx > (int)gates.size() ) {
        return NULL;
    }

    return gates[idx];
}


void Gates_block::combine(Gates_block* op_block) {

    // getting the list of gates
    std::vector<Gate*> gates_in = op_block->get_gates();

    int qbit_num_loc = op_block->get_qbit_num();
    if ( qbit_num_loc != qbit_num ) {
        std::string err("Gates_block::combine: number of qubits in the circuits must be the same"); 
        throw err;
    }

    for(std::vector<Gate*>::iterator it = (gates_in).begin(); it != (gates_in).end(); ++it) {
        Gate* op = *it;
        Gate* op_cloned = op->clone();
        op_cloned->clear_children();
        op_cloned->clear_parents();
        add_gate( op_cloned );
    }

}

void Gates_block::set_min_fusion( int min_fusion ) {
    this->min_fusion = min_fusion;
    fusion_block.update(std::shared_ptr<Gates_block>{});
    for (size_t i = 0; i < gates.size(); i++){
        if ( gates[i]->get_type() == BLOCK_OPERATION ) {
            ((Gates_block*)gates[i])->set_min_fusion( min_fusion );
        }
    }
}

void Gates_block::set_qbit_num( int qbit_num_in ) {

    // setting the number of qubits
    Gate::set_qbit_num(qbit_num_in);

    // setting the number of qubit in the gates
    for(std::vector<Gate*>::iterator it = gates.begin(); it != gates.end(); ++it) {
        Gate* op = *it;
        switch (op->get_type()) {
        case CNOT_OPERATION: case CZ_OPERATION:
        case CH_OPERATION: case SYC_OPERATION:
        case U1_OPERATION: case U2_OPERATION: 
        case U3_OPERATION:
        case CP_OPERATION: case SWAP_OPERATION:
        case CSWAP_OPERATION: case CCX_OPERATION:
        case RY_OPERATION: case CRY_OPERATION: 
        case RX_OPERATION: case CR_OPERATION:
        case RZ_OPERATION: case X_OPERATION:
        case Y_OPERATION: case Z_OPERATION:
        case SX_OPERATION: case BLOCK_OPERATION:
        case GENERAL_OPERATION:
        case ADAPTIVE_OPERATION: case CROT_OPERATION:
        case H_OPERATION: case R_OPERATION:
        case CU_OPERATION:
        case S_OPERATION: case SDG_OPERATION:
        case T_OPERATION: case TDG_OPERATION:
        case RXX_OPERATION: case RYY_OPERATION:
        case RZZ_OPERATION: case SXDG_OPERATION:
            op->set_qbit_num( qbit_num_in );
            break;
        default:
            std::string err("Gates_block::set_qbit_num: unimplemented gate"); 
            throw err;
        }
    }
    invalidate_structure_cache();
}


Gates_block* Gates_block::clone() {

    // creatign new instance of class Gates_block
    Gates_block* ret = new Gates_block( qbit_num );
    
    ret->set_parameter_start_idx( get_parameter_start_idx() );
    ret->set_parents( parents );
    ret->set_children( children );

    // extracting the gates from the current class
    if (extract_gates( ret ) != 0 ) {
        std::string err("Gates_block::clone(): extracting gates was not succesfull"); 
        throw err;
    };

    return ret;

}


int Gates_block::extract_gates( Gates_block* op_block ) {

    op_block->release_gates();

    for ( std::vector<Gate*>::iterator it=gates.begin(); it != gates.end(); ++it ) {
        Gate* op = *it;
        switch (op->get_type()) {
        case CNOT_OPERATION: case CZ_OPERATION:
        case CH_OPERATION: case SYC_OPERATION:
        case U1_OPERATION: case U2_OPERATION: 
        case U3_OPERATION: case CP_OPERATION:
        case RY_OPERATION: case CRY_OPERATION: 
        case CRX_OPERATION: case CRZ_OPERATION:
        case RX_OPERATION: case CR_OPERATION:
        case RZ_OPERATION: case X_OPERATION:
        case Y_OPERATION: case Z_OPERATION:
        case SX_OPERATION: case BLOCK_OPERATION:
        case GENERAL_OPERATION:
        case ADAPTIVE_OPERATION: case CROT_OPERATION:
        case H_OPERATION:  case R_OPERATION:
        case CU_OPERATION:
        case S_OPERATION: case SDG_OPERATION:
        case T_OPERATION: case TDG_OPERATION:
        case SWAP_OPERATION: case RXX_OPERATION:
        case RYY_OPERATION: case RZZ_OPERATION:
        case CSWAP_OPERATION: case CCX_OPERATION:
        case SXDG_OPERATION:
        {
            Gate* op_cloned = op->clone();
            op_block->add_gate( op_cloned );
            break; }
        default:
            std::string err("Gates_block::extract_gates: unimplemented gate"); 
            throw err;
        }

    }
    
    return 0;

}



bool Gates_block::contains_adaptive_gate() {

    for ( std::vector<Gate*>::iterator it=gates.begin(); it != gates.end(); ++it ) {
        Gate* op = *it;

        if (op->get_type() == ADAPTIVE_OPERATION) {
            return true;
        }
        else if (op->get_type() == BLOCK_OPERATION) {
            Gates_block* block_op = static_cast<Gates_block*>( op );
            bool ret = block_op->contains_adaptive_gate();
            if ( ret ) return true;
        }
    }

    return false;

}




bool Gates_block::contains_adaptive_gate(int idx) {

    
    Gate* op = gates[idx];

    if (op->get_type() == ADAPTIVE_OPERATION) {
        return true;
    }
    else if (op->get_type() == BLOCK_OPERATION) {
        Gates_block* block_op = static_cast<Gates_block*>( op );
        return block_op->contains_adaptive_gate();
    }

    return false;

}



Matrix
Gates_block::get_reduced_density_matrix( Matrix_real& parameters_mtx, Matrix& input_state, matrix_base<int>& qbit_list_subset ) {


    if (input_state.cols != 1) {
        std::string error("Gates_block::get_reduced_density_matrix: The number of columns in input state should be 1");
        throw error;
    }


    if (input_state.rows != matrix_size) {
        std::string error("Gates_block::get_reduced_density_matrix: The number of rows in input state should be 2^qbit_num");
        throw error;
    }

    // determine the transformed state
    Matrix transformed_state = input_state.copy();  
    if ( parameters_mtx.size() > 0 ) {
        bool parallel = true;
        apply_to( parameters_mtx, transformed_state, parallel );
    }


    int subset_qbit_num = qbit_list_subset.size();
    int complementary_qbit_num = qbit_num - subset_qbit_num;


    // list of complementary qubits
    matrix_base<int> qbit_list_complementary( complementary_qbit_num, 1);
    int qbit_idx_count = 0;
    for (int qbit_idx=0; qbit_idx<qbit_num; qbit_idx++) {

        bool qbit_idx_in_subset = false;

        for (int subset_qbit_idx=0; subset_qbit_idx<subset_qbit_num; subset_qbit_idx++) {
            if ( qbit_idx == qbit_list_subset[subset_qbit_idx] ) {
                qbit_idx_in_subset = true;
                break;
            }
        }


        if ( qbit_idx_in_subset ) {
            continue;
        }

        qbit_list_complementary[qbit_idx_count] = qbit_idx;
        qbit_idx_count++;
    }


    // 000010000 one-hot encoded numbers indicating the bit position of the qubits in the register
    matrix_base<int> qbit_masks(qbit_num, 1);
    for (int qbit_idx=0; qbit_idx<qbit_num; qbit_idx++) {
        qbit_masks[ qbit_idx ] = 1 << qbit_idx;
    }



    // retrieve the reduced density matrix
    int rho_matrix_size = 1 << subset_qbit_num;

    Matrix rho(rho_matrix_size, rho_matrix_size);
    memset( rho.get_data(), 0, rho.size()*sizeof(QGD_Complex16) );




    int complementary_basis_num = 1 << complementary_qbit_num;
    for ( int row_idx=0; row_idx<rho_matrix_size; row_idx++ ) {

        // index of the amplitude in the state vector
        int idx = 0;
        for (int qbit_idx=0; qbit_idx<subset_qbit_num; qbit_idx++) {
            if ( row_idx & qbit_masks[ qbit_idx ] ) {
                idx = idx | qbit_masks[ qbit_list_subset[qbit_idx] ]; 
            }
        }

        for ( int col_idx=row_idx; col_idx<rho_matrix_size; col_idx++ ) {
        
            // index of the amplitude in the state vector
            int jdx = 0;
            for (int qbit_idx=0; qbit_idx<subset_qbit_num; qbit_idx++) {
                if ( col_idx & qbit_masks[ qbit_idx ] ) {
                    jdx = jdx | qbit_masks[ qbit_list_subset[qbit_idx] ]; 
                }
            }


            // thread local storage for partial permanent
            tbb::combinable<QGD_Complex16> priv_addend {[](){QGD_Complex16 ret; ret.real = 0.0; ret.imag = 0.0; return ret;}};


            tbb::parallel_for( tbb::blocked_range<int>(0, complementary_basis_num, 1024), [&](tbb::blocked_range<int> r) {

                QGD_Complex16& rho_element_priv = priv_addend.local();

                for (int compl_idx=r.begin(); compl_idx<r.end(); compl_idx++) {

 
                    int idx_loc = idx;
                    int jdx_loc = jdx;

                    for (int qbit_idx=0; qbit_idx<complementary_qbit_num; qbit_idx++) {
                        if ( compl_idx & qbit_masks[ qbit_idx ] ) {
                            idx_loc = idx_loc | qbit_masks[ qbit_list_complementary[qbit_idx] ]; 
                            jdx_loc = jdx_loc | qbit_masks[ qbit_list_complementary[qbit_idx] ]; 
                        }
                    }

                    QGD_Complex16  element_idx = transformed_state[ idx_loc ];
                    QGD_Complex16& element_jdx = transformed_state[ jdx_loc ];

                    // conjugate because of the bra-vector
                    element_idx.imag = -element_idx.imag;

                    QGD_Complex16 addend = mult( element_idx, element_jdx );

                    rho_element_priv.real = rho_element_priv.real + addend.real;
                    rho_element_priv.imag = rho_element_priv.imag + addend.imag;


                }
            });

        

            QGD_Complex16 rho_element;
            rho_element.real = 0.0;
            rho_element.imag = 0.0;

            priv_addend.combine_each([&](QGD_Complex16 &a) {
                rho_element.real = rho_element.real + a.real;
                rho_element.imag = rho_element.imag + a.imag;
            });

            rho[ row_idx * rho.stride + col_idx ].real += rho_element.real;
            rho[ row_idx * rho.stride + col_idx ].imag += rho_element.imag;

            if ( row_idx == col_idx ) {
                continue;
            }

            rho[ col_idx * rho.stride + row_idx ].real += rho_element.real;
            rho[ col_idx * rho.stride + row_idx ].imag -= rho_element.imag;

        }

        
    }


    // test the trace of the reduced density matrix
    double trace = 0.0;
    for( int idx=0; idx<rho_matrix_size; idx++) {
        trace = trace + rho[idx*rho.stride+idx].real;
    }

    if ( abs( trace-1.0 ) > 1e-6 ) {
        std::stringstream sstream;
          sstream << "Warning: Gates_block::get_reduced_density_matrix: The error of the trace of the reduced density matrix is " << abs( trace-1.0 ) << std::endl;
         print(sstream, 1);    
    }

    return rho;


}


double Gates_block::get_second_Renyi_entropy( Matrix_real& parameters_mtx, Matrix& input_state, matrix_base<int>& qbit_list_subset ) {

    // determine the reduced density matrix
    Matrix rho = get_reduced_density_matrix( parameters_mtx, input_state, qbit_list_subset );


    // calculate the second Rényi entropy 
    // Tr( rho @ rho )
    double trace_rho_square = 0.0;

    for (int idx=0; idx<rho.rows; idx++) {

        double trace_tmp = 0.0;

        for (int jdx=0; jdx<rho.rows; jdx++) {
            QGD_Complex16& element = rho[ idx*rho.stride + jdx ];
            double tmp = element.real * element.real + element.imag * element.imag;

            trace_tmp = trace_tmp + tmp;
        }


        trace_rho_square = trace_rho_square + trace_tmp;


    }




    double entropy = -log(trace_rho_square);


    return entropy;

}


std::vector<int> remove_list_intersection( std::vector<int>& list1, std::vector<int>& list2 ) {

    std::vector<int> ret = list1;

    for( std::vector<int>::iterator it2 = list2.begin(); it2 != list2.end(); it2++ ) {
    
        std::vector<int>::iterator element_found = std::find(ret.begin(), ret.end(), *it2);
        
        if( element_found != ret.end() ) {
            ret.erase( element_found );
        }    
        
    }
    
    return ret;



}


void 
Gates_block::determine_parents( Gate* gate ) {


    std::vector<int>&& involved_qubits = gate->get_involved_qubits();
/*
    std::cout << "involved qubits in the current gate: " << std::endl;
    for( int idx=0; idx<involved_qubits.size(); idx++ ) {
    std::cout << involved_qubits[idx] << ", ";
    }
    std::cout << std::endl;
  */  
    // iterate over gates in the circuit
    for( int idx=static_cast<int>(gates.size()-1); idx>=0; idx-- ) {
        Gate* gate_loc = gates[idx];
        std::vector<int>&& involved_qubits_loc = gate_loc->get_involved_qubits();
        
        std::vector<int>&& reduced_qbit_list = remove_list_intersection( involved_qubits, involved_qubits_loc );
        
        if( reduced_qbit_list.size() < involved_qubits.size() ) {
            // parent gate found, setting parent-child relation
            
            gate->add_parent( gate_loc );
            gate_loc->add_child( gate );
            
            involved_qubits = std::move(reduced_qbit_list);
            
            
        }
        
        
        // return if no further involved qubits left
        if( involved_qubits.size() == 0 ) {
            break;
        }
    }
    
    
}
    
    
    
void 
Gates_block::determine_children( Gate* gate ) {


    std::vector<int>&& involved_qubits = gate->get_involved_qubits();
    
    // iterate over gates in the circuit
    for( size_t idx=0; idx<gates.size(); idx++ ) {
        Gate* gate_loc = gates[idx];
        std::vector<int>&& involved_qubits_loc = gate_loc->get_involved_qubits();
        
        std::vector<int>&& reduced_qbit_list = remove_list_intersection( involved_qubits, involved_qubits_loc );
        
        if( reduced_qbit_list.size() < involved_qubits.size() ) {
            // child gate found, setting parent-child relation
            
            gate->add_child( gate_loc );
            gate_loc->add_parent( gate );
            
            involved_qubits = std::move(reduced_qbit_list);
            
            
        }
        
        
        // return if no further involved qubits left
        if( involved_qubits.size() == 0 ) {
            break;
        }
    }    

}


void 
Gates_block::reset_parameter_start_indices() {

    int parameter_idx = 0;
    
    for( std::vector<Gate*>::iterator gate_it = gates.begin(); gate_it != gates.end(); gate_it++ ) {
    
        Gate* gate = *gate_it;
        
        gate->set_parameter_start_idx( parameter_idx );
        
        parameter_idx = parameter_idx + gate->get_parameter_num();
    
    }

}


void  
Gates_block::reset_dependency_graph() {

    // first clear parent/children data from the gates
    for( std::vector<Gate*>::iterator gate_it = gates.begin(); gate_it != gates.end(); gate_it++ ) {    
        Gate* gate = *gate_it;
        gate->clear_children();
        gate->clear_parents();    
    }


    // first clear parent/children data from the gates
    for( std::vector<Gate*>::iterator gate_it = gates.begin(); gate_it != gates.end(); gate_it++ ) {    
        Gate* gate = *gate_it;

        // determine the parents of the gate
        determine_parents( gate );

    }




}


Gates_block* 
Gates_block::get_flat_circuit() {

    Gates_block* flat_circuit = new Gates_block( qbit_num );

    for( std::vector<Gate*>::iterator gate_it=gates.begin(); gate_it != gates.end(); gate_it++ ) {

        Gate* gate = *gate_it;

        if( gate->get_type() == BLOCK_OPERATION ) {

            Gates_block* circuit_inner      = static_cast<Gates_block*>( gate );
            Gates_block* flat_circuit_inner = circuit_inner->get_flat_circuit();

            flat_circuit->combine( flat_circuit_inner );

            delete( flat_circuit_inner );
        }
        else {
            flat_circuit->add_gate( gate->clone() );
        }

    }



    return flat_circuit;


}



Matrix_real 
Gates_block::extract_parameters( Matrix_real& parameters ) {

    if ( get_parameter_start_idx() + get_parameter_num() < parameters.size()  ) {
        std::string err("Gates_block::extract_parameters: Cant extract parameters, since th einput arary has not enough elements.");
        throw err;     
    }

    Matrix_real extracted_parameters(1, get_parameter_num());

    memcpy( extracted_parameters.get_data(), parameters.get_data() + get_parameter_start_idx(), get_parameter_num()*sizeof(double) );

    return extracted_parameters;

}


#ifdef __DFE__


DFEgate_kernel_type* Gates_block::convert_to_DFE_gates_with_derivates( Matrix_real& parameters_mtx, int& gatesNum, int& gateSetNum, int& redundantGateSets, bool only_derivates ) {

    int parameter_num = get_parameter_num();
    if ( parameter_num != parameters_mtx.size() ) {
        std::string error("Gates_block::convert_to_DFE_gates: wrong number of parameters");
        throw error;
    }

    std::map<std::string, int> gate_nums  = get_gate_nums();
    int gates_total_num = 0;
    for( auto it=gate_nums.begin(); it != gate_nums.end(); it++ ) {
        gates_total_num = gates_total_num + it->second;
    }

    int chained_gates_num = get_chained_gates_num();
    int gate_padding      = gates_total_num % chained_gates_num == 0 ? 0 : chained_gates_num - (gates_total_num % chained_gates_num);
    gatesNum              = gates_total_num+gate_padding;
/*
std::cout << "chained gates num: " << chained_gates_num << std::endl;
std::cout << "number of gates: " << gatesNum << std::endl;
*/


    gateSetNum = only_derivates ? parameter_num : parameter_num+1;

#ifdef __MPI__
    int rem = gateSetNum % (4 * world_size );
    if ( rem == 0 ) {
        redundantGateSets = 0;
    }
    else {
        redundantGateSets = (4 * world_size ) - (gateSetNum % (4 * world_size ));
        gateSetNum = gateSetNum + redundantGateSets;
    }
#else
    int rem = gateSetNum % 4;
    if ( rem == 0 ) {
        redundantGateSets = 0;
    }
    else {
        redundantGateSets = 4 - (gateSetNum % 4);
        gateSetNum = gateSetNum + redundantGateSets;
    }
#endif


    DFEgate_kernel_type* DFEgates = new DFEgate_kernel_type[gatesNum*gateSetNum];
    

    int gate_idx = 0;
    convert_to_DFE_gates( parameters_mtx, DFEgates, gate_idx );


    // padding with identity gates
    for (int idx=gate_idx; idx<gatesNum; idx++ ){

        DFEgate_kernel_type& DFEGate = DFEgates[idx];

        
        DFEGate.target_qbit = 0;
        DFEGate.control_qbit = -1;
        DFEGate.gate_type = U3_OPERATION;
        DFEGate.ThetaOver2 = (int32_t)(0);
        DFEGate.Phi = (int32_t)(0);
        DFEGate.Lambda = (int32_t)(0); 
        DFEGate.metadata = 0;

    }
/*
    for ( int idx=0; idx<gatesNum; idx++ ) {

        std::cout << "target qubit: " << (int)DFEgates[idx].target_qbit << " control qubit: " << (int)DFEgates[idx].control_qbit << " gate type: " << (int)DFEgates[idx].gate_type << std::endl; 
    }
*/

    // adjust parameters for derivation
    if (only_derivates ) {
        for (int idx=1; idx<(gateSetNum-1); idx++) {
           memcpy(DFEgates+idx*gatesNum, DFEgates, gatesNum*sizeof(DFEgate_kernel_type));
        }
    }
    else {
        for (int idx=0; idx<(gateSetNum-1); idx++) {
           memcpy(DFEgates+(idx+1)*gatesNum, DFEgates, gatesNum*sizeof(DFEgate_kernel_type));
        }
    }

    gate_idx = 0;
    int gate_set_index = parameter_num-1;
    if (only_derivates) {
        adjust_parameters_for_derivation( DFEgates, gatesNum, gate_idx, gate_set_index );
    } 
    else {
        adjust_parameters_for_derivation( DFEgates+gatesNum, gatesNum, gate_idx, gate_set_index );
    }

/*
    for ( int idx=0; idx<gatesNum*(parameter_num+1); idx++ ) {

        std::cout << "target qubit: " << (int)DFEgates[idx].target_qbit << " control qubit: " << (int)DFEgates[idx].control_qbit << " theta: " << (int)DFEgates[idx].ThetaOver2 << " lambda: " << (int)DFEgates[idx].Lambda << std::endl; 
    }
*/
    return DFEgates;

}


void Gates_block::adjust_parameters_for_derivation( DFEgate_kernel_type* DFEgates, const int gatesNum, int& gate_idx, int& gate_set_index) {

        int parameter_idx = parameter_num;
        //int gate_set_index = parameter_num-1;

        int32_t parameter_shift = (int32_t)(M_PI/2*(1<<25));

        for(int op_idx = gates.size()-1; op_idx>=0; op_idx--) { 

            Gate* gate = gates[op_idx];
//std::cout <<   gate_idx << " " <<   gate_set_index << " " << gate->get_type() << std::endl;        

            if (gate->get_type() == CNOT_OPERATION) {
                gate_idx = gate_idx + 1;
            }
            else if (gate->get_type() == CZ_OPERATION) {            
                gate_idx = gate_idx + 1;
            }
            else if (gate->get_type() == CH_OPERATION) {         
                gate_idx = gate_idx + 1;
            }
            else if (gate->get_type() == SYC_OPERATION) {
                std::string error("Gates_block::convert_to_DFE_gates: SYC_gate not implemented");
                throw error;
            }
            else if (gate->get_type() == U1_OPERATION) {
                // Lambda parameter derivative
                DFEgate_kernel_type& DFEGate = DFEgates[gate_set_index*gatesNum + gate_idx];
                DFEGate.Lambda = DFEGate.Lambda + parameter_shift;
                DFEGate.metadata = 5 + (1<<7); // Lambda parameter derivative// The 0st and 3nd element in kernel matrix should be zero for derivates and 5 = 0101, plus the leading bit indicates that derivate is processed
                gate_set_index = gate_set_index - 1;
                
                parameter_idx = parameter_idx - 1;
                gate_idx = gate_idx + 1;
            }
            else if (gate->get_type() == U2_OPERATION) {
                // Lambda parameter derivative
                DFEgate_kernel_type& DFEGate2 = DFEgates[gate_set_index*gatesNum + gate_idx];
                DFEGate2.Lambda = DFEGate2.Lambda + parameter_shift;
                DFEGate2.metadata = 5 + (1<<7); // The 0st and 3nd element in kernel matrix should be zero for derivates and 5 = 0101, plus the leading bit indicates that derivate is processed
                gate_set_index = gate_set_index - 1;

                // Phi parameter derivative
                DFEgate_kernel_type& DFEGate = DFEgates[gate_set_index*gatesNum + gate_idx];
                DFEGate.Phi = DFEGate.Phi + parameter_shift;
                DFEGate.metadata = 3 + (1<<7); // The 0st and 3nd element in kernel matrix should be zero for derivates and 5 = 0101, plus the leading bit indicates that derivate is processed
                gate_set_index = gate_set_index - 1;

                parameter_idx = parameter_idx - 2;
                gate_idx = gate_idx + 1;
            }
            else if (gate->get_type() == U3_OPERATION) {
                // Lambda parameter derivative
                DFEgate_kernel_type& DFEGate3 = DFEgates[gate_set_index*gatesNum + gate_idx];
                DFEGate3.Lambda = DFEGate3.Lambda + parameter_shift;
                DFEGate3.metadata = 5 + (1<<7); // The 0st and 3nd element in kernel matrix should be zero for derivates and 5 = 0101, plus the leading bit indicates that derivate is processed
                gate_set_index = gate_set_index - 1;

                // Phi parameter derivative
                DFEgate_kernel_type& DFEGate2 = DFEgates[gate_set_index*gatesNum + gate_idx];
                DFEGate2.Phi = DFEGate2.Phi + parameter_shift;
                DFEGate2.metadata = 3 + (1<<7); // The 0st and 3nd element in kernel matrix should be zero for derivates and 5 = 0101, plus the leading bit indicates that derivate is processed
                gate_set_index = gate_set_index - 1;

                // Theta parameter derivative
                DFEgate_kernel_type& DFEGate = DFEgates[gate_set_index*gatesNum + gate_idx];
                DFEGate.ThetaOver2 = DFEGate.ThetaOver2 + parameter_shift;
                DFEGate.metadata = (1<<7); // the leading bit indicates that derivate is processed
                gate_set_index = gate_set_index - 1;

                parameter_idx = parameter_idx - 3;
                gate_idx = gate_idx + 1;
            }
            else if (gate->get_type() == RX_OPERATION) { // Did not cehcked

                DFEgate_kernel_type& DFEGate = DFEgates[gate_set_index*gatesNum + gate_idx];
                DFEGate.ThetaOver2           = DFEGate.ThetaOver2 + parameter_shift;
                DFEGate.metadata             = (1<<7); // the leading bit indicates that derivate is processed
                gate_set_index               = gate_set_index - 1;
     
                parameter_idx = parameter_idx - 1;    
           
                gate_idx = gate_idx + 1;
            }
            else if (gate->get_type() == RY_OPERATION) { // Did not cehcked

                DFEgate_kernel_type& DFEGate = DFEgates[gate_set_index*gatesNum + gate_idx];
                DFEGate.ThetaOver2           = DFEGate.ThetaOver2 + parameter_shift;
                DFEGate.metadata             = (1<<7); // the leading bit indicates that derivate is processed
                gate_set_index               = gate_set_index - 1;

                parameter_idx = parameter_idx - 1;

     
                gate_idx = gate_idx + 1;
            }
            else if (gate->get_type() == CRY_OPERATION) {

                DFEgate_kernel_type& DFEGate = DFEgates[gate_set_index*gatesNum + gate_idx];
                DFEGate.ThetaOver2           = DFEGate.ThetaOver2 + parameter_shift;
                DFEGate.metadata             = (1<<7); // the leading bit indicates that derivate is processed
                gate_set_index               = gate_set_index - 1;

                parameter_idx = parameter_idx - 1;
                                   
                gate_idx = gate_idx + 1;
            }
            else if (gate->get_type() == RZ_OPERATION) { // Did not cehcked

                std::string error("Gates_block::adjust_parameters_for_derivation: RZ gate not implemented for DFE");
                throw error;       
            }
            else if (gate->get_type() == CU_OPERATION) {
                std::string error("Gates_block::convert_to_DFE_gates: CU_gate not implemented");
                throw error;   

            }
            else if (gate->get_type() == S_OPERATION) {
                std::string error("Gates_block::convert_to_DFE_gates: S_gate not implemented");
                throw error;   

            }
            else if (gate->get_type() == SDG_OPERATION) {
                std::string error("Gates_block::convert_to_DFE_gates: Sdg_gate not implemented");
                throw error;   

            }
            else if (gate->get_type() == H_OPERATION) {
                std::string error("Gates_block::convert_to_DFE_gates: H_gate not implemented");
                throw error;   

            }
            else if (gate->get_type() == X_OPERATION) {
                std::string error("Gates_block::convert_to_DFE_gates: X_gate not implemented");
                throw error;   

            }
            else if (gate->get_type() == Y_OPERATION) {
                std::string error("Gates_block::convert_to_DFE_gates: Y_gate not implemented");
                throw error;   

            }
            else if (gate->get_type() == Z_OPERATION) {
                std::string error("Gates_block::convert_to_DFE_gates: Z_gate not implemented");
                throw error;   

            }
            else if (gate->get_type() == T_OPERATION) {
                std::string error("Gates_block::convert_to_DFE_gates: T_gate not implemented");
                throw error;   

            }
            else if (gate->get_type() == TDG_OPERATION) {
                std::string error("Gates_block::convert_to_DFE_gates: Tdg_gate not implemented");
                throw error;   
            }
            else if (gate->get_type() == SX_OPERATION) {
                std::string error("Gates_block::convert_to_DFE_gates: SX_gate not implemented");
                throw error;       
            }
            else if (gate->get_type() == SXDG_OPERATION) {
                std::string error("Gates_block::convert_to_DFE_gates: SXdg_gate not implemented");
                throw error;       
            }
            else if (gate->get_type() == BLOCK_OPERATION) {

                Gates_block* block_gate = static_cast<Gates_block*>(gate);
                block_gate->adjust_parameters_for_derivation( DFEgates, gatesNum, gate_idx, gate_set_index);         
                //gate_set_index               = gate_set_index - block_gate->get_parameter_num();

                parameter_idx = parameter_idx - block_gate->get_parameter_num();
            }
            else if (gate->get_type() == GENERAL_OPERATION) {
                std::string error("Gates_block::convert_to_DFE_gates: general_gate not implemented");
                throw error;
            }
            else if (gate->get_type() == ADAPTIVE_OPERATION) {

                DFEgate_kernel_type& DFEGate = DFEgates[gate_set_index*gatesNum + gate_idx];
                DFEGate.metadata              = (1<<7); // the leading bit indicates that derivate is processed
                DFEGate.ThetaOver2           = DFEGate.ThetaOver2 + parameter_shift;
                gate_set_index               = gate_set_index - 1;

                parameter_idx = parameter_idx - 1;
                                   
                gate_idx = gate_idx + 1;
            }
            else {
                std::string err("Gates_block::adjust_parameters_for_derivation: unimplemented gate"); 
                throw err;
            }

        }


}


DFEgate_kernel_type* 
Gates_block::convert_to_batched_DFE_gates( std::vector<Matrix_real>& parameters_mtx_vec, int& gatesNum, int& gateSetNum, int& redundantGateSets ) {


    std::map<std::string, int> gate_nums  = get_gate_nums();
    int gates_total_num = 0;
    for( auto it=gate_nums.begin(); it != gate_nums.end(); it++ ) {
        gates_total_num = gates_total_num + it->second;
    }

    int chained_gates_num = get_chained_gates_num();
    int gate_padding      = gates_total_num % chained_gates_num == 0 ? 0 : chained_gates_num - (gates_total_num % chained_gates_num);
    gatesNum              = gates_total_num+gate_padding;
/*
std::cout << "chained gates num: " << chained_gates_num << std::endl;
std::cout << "number of gates: " << gatesNum << std::endl;
*/


    gateSetNum = parameters_mtx_vec.size();

#ifdef __MPI__
    int rem = gateSetNum % (4 * world_size );
    if ( rem == 0 ) {
        redundantGateSets = 0;
    }
    else {
        redundantGateSets = (4 * world_size ) - (gateSetNum % (4 * world_size ));
        gateSetNum = gateSetNum + redundantGateSets;
    }
#else
    int rem = gateSetNum % 4;
    if ( rem == 0 ) {
        redundantGateSets = 0;
    }
    else {
        redundantGateSets = 4 - (gateSetNum % 4);
        gateSetNum = gateSetNum + redundantGateSets;
    }
#endif

    DFEgate_kernel_type* DFEgates = new DFEgate_kernel_type[gatesNum*gateSetNum];


    tbb::parallel_for( 0, gateSetNum, 1, [&](int gateset_idx) {
    
        int gate_idx = gateset_idx * gatesNum;

        if ( gateset_idx < parameters_mtx_vec.size() ) {
            Matrix_real& parameters_mtx = parameters_mtx_vec[gateset_idx];
            convert_to_DFE_gates( parameters_mtx, DFEgates, gate_idx );
        }

        // padding with identity gates
        for (int idx=gate_idx; idx<(gateset_idx+1)*gatesNum; idx++ ){

            DFEgate_kernel_type& DFEGate = DFEgates[idx];

        
            DFEGate.target_qbit = 0;
            DFEGate.control_qbit = -1;
            DFEGate.gate_type = U3_OPERATION;
            DFEGate.ThetaOver2 = (int32_t)(0);
            DFEGate.Phi = (int32_t)(0);
            DFEGate.Lambda = (int32_t)(0); 
            DFEGate.metadata = 0;

            gate_idx++;

        }

    });
    

    return DFEgates;

}

DFEgate_kernel_type* Gates_block::convert_to_DFE_gates( Matrix_real& parameters_mtx, int& gatesNum ) {

    int parameter_num = get_parameter_num();
    if ( parameter_num != parameters_mtx.size() ) {
        std::string error("Gates_block::convert_to_DFE_gates: wrong number of parameters");
        throw error;
    }


    std::map<std::string, int> gate_nums  = get_gate_nums();
    int gates_total_num = 0;
    for( auto it=gate_nums.begin(); it != gate_nums.end(); it++ ) {
        gates_total_num = gates_total_num + it->second;
    }

    int chained_gates_num = get_chained_gates_num();
    int gate_padding      = chained_gates_num - (gates_total_num % chained_gates_num);
    gatesNum              = gates_total_num+gate_padding;



    DFEgate_kernel_type* DFEgates = new DFEgate_kernel_type[gates_total_num+gate_padding];
    
    int gate_idx = 0;
    convert_to_DFE_gates( parameters_mtx, DFEgates, gate_idx );


    // padding with identity gates
    for (int idx=gate_idx; idx<gatesNum; idx++ ){

        DFEgate_kernel_type& DFEGate = DFEgates[idx];

        
        DFEGate.target_qbit = 0;
        DFEGate.control_qbit = -1;
        DFEGate.gate_type = U3_OPERATION;
        DFEGate.ThetaOver2 = (int32_t)(0);
        DFEGate.Phi = (int32_t)(0);
        DFEGate.Lambda = (int32_t)(0); 
        DFEGate.metadata = 0;

    }
/*
    for ( int idx=0; idx<gatesNum; idx++ ) {

        std::cout << "target qubit: " << (int)DFEgates[idx].target_qbit << " control qubit: " << (int)DFEgates[idx].control_qbit << " gate type: " << (int)DFEgates[idx].gate_type << std::endl; 
    }
*/


    return DFEgates;

}



void Gates_block::convert_to_DFE_gates( const Matrix_real& parameters_mtx, DFEgate_kernel_type* DFEgates, int& start_index ) {

     
        int& gate_idx = start_index;
        int parameter_idx = parameter_num;
     double *parameters_data = parameters_mtx.get_data();
     //const_cast <Matrix_real&>(parameters);


        for(int op_idx = gates.size()-1; op_idx>=0; op_idx--) { 

            Gate* gate = gates[op_idx];
            DFEgate_kernel_type& DFEGate = DFEgates[gate_idx];

            if (gate->get_type() == CNOT_OPERATION) {
                CNOT* cnot_gate = static_cast<CNOT*>(gate);
                DFEGate.target_qbit = cnot_gate->get_target_qbit();
                DFEGate.control_qbit = cnot_gate->get_control_qbit();
                DFEGate.gate_type = CNOT_OPERATION;
                DFEGate.ThetaOver2 = (int32_t)(M_PI/2*(1<<25));
                DFEGate.Phi = (int32_t)(0);
                DFEGate.Lambda = (int32_t)(M_PI*(1<<25)); 
                DFEGate.metadata = 0;
                gate_idx = gate_idx + 1;
            }
            else if (gate->get_type() == CZ_OPERATION) {
                CZ* cz_gate = static_cast<CZ*>(gate);    
                DFEGate.target_qbit = cz_gate->get_target_qbit();
                DFEGate.control_qbit = cz_gate->get_control_qbit();
                DFEGate.gate_type = CZ_OPERATION;
                DFEGate.ThetaOver2 = (int32_t)(0); // TODO: check !!!!!!!!!!!!!!!!!!!!!!!!!!!
                DFEGate.Phi = (int32_t)(0); // TODO: check !!!!!!!!!!!!!!!!!!!!!!!!!!!
                DFEGate.Lambda = (int32_t)(M_PI*(1<<25));  // TODO: check !!!!!!!!!!!!!!!!!!!!!!!!!!!    
                DFEGate.metadata = 0;             
                gate_idx = gate_idx + 1;
            }
            else if (gate->get_type() == CH_OPERATION) {
                CH* ch_gate = static_cast<CH*>(gate);   
                DFEGate.target_qbit = ch_gate->get_target_qbit();
                DFEGate.control_qbit = ch_gate->get_control_qbit();
                DFEGate.gate_type = CH_OPERATION;
                DFEGate.ThetaOver2 = (int32_t)(M_PI/4*(1<<25)); // TODO: check !!!!!!!!!!!!!!!!!!!!!!!!!!!
                DFEGate.Phi = (int32_t)(0); // TODO: check !!!!!!!!!!!!!!!!!!!!!!!!!!!
                DFEGate.Lambda = (int32_t)(M_PI*(1<<25));  // TODO: check !!!!!!!!!!!!!!!!!!!!!!!!!!!     
                DFEGate.metadata = 0;             
                gate_idx = gate_idx + 1;
            }
            else if (gate->get_type() == SYC_OPERATION) {
                std::string error("Gates_block::convert_to_DFE_gates: SYC_gate not implemented");
                throw error;
            }
            else if (gate->get_type() == U1_OPERATION) {
                double lambda = std::fmod(parameters_data[parameter_idx-1], 2*M_PI);
                parameter_idx = parameter_idx - 1;

                DFEGate.target_qbit = gate->get_target_qbit();
                DFEGate.control_qbit = -1;
                DFEGate.gate_type = U1_OPERATION;
                DFEGate.ThetaOver2 = (int32_t)(0);
                DFEGate.Phi = (int32_t)(0);
                DFEGate.Lambda = (int32_t)(lambda*(1<<25));
                DFEGate.metadata = 0;
                gate_idx = gate_idx + 1;
            }
            else if (gate->get_type() == U2_OPERATION) {
                double phi = std::fmod(parameters_data[parameter_idx-2], 2*M_PI);
                double lambda = std::fmod(parameters_data[parameter_idx-1], 2*M_PI);
                parameter_idx = parameter_idx - 2;

                DFEGate.target_qbit = gate->get_target_qbit();
                DFEGate.control_qbit = -1;
                DFEGate.gate_type = U2_OPERATION;
                DFEGate.ThetaOver2 = (int32_t)(M_PI/4*(1<<25)); // theta/2 = pi/4 for U2
                DFEGate.Phi = (int32_t)(phi*(1<<25));
                DFEGate.Lambda = (int32_t)(lambda*(1<<25));
                DFEGate.metadata = 0;
                gate_idx = gate_idx + 1;
            }
            else if (gate->get_type() == U3_OPERATION) {
                double theta = std::fmod(parameters_data[parameter_idx-3], 4*M_PI);
                double phi = std::fmod(parameters_data[parameter_idx-2], 2*M_PI);
                double lambda = std::fmod(parameters_data[parameter_idx-1], 2*M_PI);
                parameter_idx = parameter_idx - 3;

                DFEGate.target_qbit = gate->get_target_qbit();
                DFEGate.control_qbit = -1;
                DFEGate.gate_type = U3_OPERATION;
                DFEGate.ThetaOver2 = (int32_t)(theta/2.0*(1<<25)); // theta/2
                DFEGate.Phi = (int32_t)(phi*(1<<25));
                DFEGate.Lambda = (int32_t)(lambda*(1<<25));
                DFEGate.metadata = 0;
                gate_idx = gate_idx + 1;
            }
            else if (gate->get_type() == RX_OPERATION) {
             // definig the rotation parameter
                double varthetaOver2;
                // get the inverse parameters of the U3 rotation
                RX* rx_gate = static_cast<RX*>(gate);       
                varthetaOver2 = std::fmod( parameters_data[parameter_idx-1], 2*M_PI);
                parameter_idx = parameter_idx - 1;

                DFEGate.target_qbit = rx_gate->get_target_qbit();
                DFEGate.control_qbit = -1;
                DFEGate.gate_type = RX_OPERATION;
                DFEGate.ThetaOver2 = (int32_t)(varthetaOver2*(1<<25)); // TODO: check !!!!!!!!!!!!!!!!!!!!!!!!!!!
                DFEGate.Phi = (int32_t)(-M_PI/2*(1<<25));  // TODO: check !!!!!!!!!!!!!!!!!!!!!!!!!!!
                DFEGate.Lambda = (int32_t)(M_PI/2*(1<<25)); // TODO: check !!!!!!!!!!!!!!!!!!!!!!!!!!!   
                DFEGate.metadata = 0;    
           
                gate_idx = gate_idx + 1;
            }
            else if (gate->get_type() == RY_OPERATION) {
                // definig the rotation parameter
                double varthetaOver2;
                // get the inverse parameters of the U3 rotation
                RY* ry_gate = static_cast<RY*>(gate);
                varthetaOver2 = std::fmod( parameters_data[parameter_idx-1], 2*M_PI);
                parameter_idx = parameter_idx - 1;

                DFEGate.target_qbit = ry_gate->get_target_qbit();
                DFEGate.control_qbit = -1;
                DFEGate.gate_type = RY_OPERATION;
                DFEGate.ThetaOver2 = (int32_t)(varthetaOver2*(1<<25)); // TODO: check !!!!!!!!!!!!!!!!!!!!!!!!!!!
                DFEGate.Phi = (int32_t)(0);  // TODO: check !!!!!!!!!!!!!!!!!!!!!!!!!!!
                DFEGate.Lambda = (int32_t)(0); // TODO: check !!!!!!!!!!!!!!!!!!!!!!!!!!! 
                DFEGate.metadata = 0;
     
                gate_idx = gate_idx + 1;
            }
            else if (gate->get_type() == CRY_OPERATION) {
                // definig the rotation parameter
                double Phi;
                // get the inverse parameters of the U3 rotation
                CRY* cry_gate = static_cast<CRY*>(gate);
                double varthetaOver2 = std::fmod( parameters_data[parameter_idx-1], 2*M_PI);
                parameter_idx = parameter_idx - 1;
                DFEGate.target_qbit = cry_gate->get_target_qbit();
                DFEGate.control_qbit = cry_gate->get_control_qbit();
                DFEGate.gate_type = CRY_OPERATION;
                DFEGate.ThetaOver2 = (int32_t)(varthetaOver2*(1<<25)); // TODO: check !!!!!!!!!!!!!!!!!!!!!!!!!!!
                DFEGate.Phi = (int32_t)(0);  // TODO: check !!!!!!!!!!!!!!!!!!!!!!!!!!!
                DFEGate.Lambda = (int32_t)(0); // TODO: check !!!!!!!!!!!!!!!!!!!!!!!!!!!
                DFEGate.metadata = 0;
                                   
                gate_idx = gate_idx + 1;
            }
            else if (gate->get_type() == RZ_OPERATION) { // Did not cehcked

                std::string error("Gates_block::convert_to_DFE_gates: RZ gate not implemented for DFE");
                throw error;       
            }
            else if (gate->get_type() == H_OPERATION) {
                // get the inverse parameters of the Hadamard rotation
                H* h_gate = static_cast<H*>(gate);

                std::string error("Gates_block::convert_to_DFE_gates: Hadamard gate not implemented");
                throw error;
          
                gate_idx = gate_idx + 1;
            }
            else if (gate->get_type() == X_OPERATION) {
                // get the inverse parameters of the U3 rotation
                X* x_gate = static_cast<X*>(gate);

                DFEGate.target_qbit = x_gate->get_target_qbit();
                DFEGate.control_qbit = -1;
                DFEGate.gate_type = X_OPERATION;
                DFEGate.ThetaOver2 = (int32_t)(M_PI/2*(1<<25));
                DFEGate.Phi = (int32_t)(0);
                DFEGate.Lambda = (int32_t)(M_PI*(1<<25)); 
                DFEGate.metadata = 0;
          
                gate_idx = gate_idx + 1;
            }
            else if (gate->get_type() == Y_OPERATION) {
                // get the inverse parameters of the U3 rotation
                Y* y_gate = static_cast<Y*>(gate);

                DFEGate.target_qbit = y_gate->get_target_qbit();
                DFEGate.control_qbit = -1;
                DFEGate.gate_type = Y_OPERATION;
                DFEGate.ThetaOver2 = (int32_t)(M_PI/2*(1<<25));
                DFEGate.Phi = (int32_t)(M_PI/2*(1<<25));
                DFEGate.Lambda = (int32_t)(M_PI/2*(1<<25)); 
                DFEGate.metadata = 0;
          
                gate_idx = gate_idx + 1;
            }
            else if (gate->get_type() == Z_OPERATION) {
                // get the inverse parameters of the U3 rotation
                Z* z_gate = static_cast<Z*>(gate);

                DFEGate.target_qbit = z_gate->get_target_qbit();
                DFEGate.control_qbit = -1;
                DFEGate.gate_type = Z_OPERATION;
                DFEGate.ThetaOver2 = (int32_t)(0);
                DFEGate.Phi = (int32_t)(0);
                DFEGate.Lambda = (int32_t)(M_PI*(1<<25)); 
                DFEGate.metadata = 0;
          
                gate_idx = gate_idx + 1;
            }
            else if (gate->get_type() == CU_OPERATION) {
                // get the inverse parameters of the U3 rotation
                CU* s_gate = static_cast<CU*>(gate);
                std::string error("Gates_block::convert_to_DFE_gates: CU_gate not implemented");
                throw error;       
            }
            else if (gate->get_type() == S_OPERATION) {
                // get the inverse parameters of the U3 rotation
                S* s_gate = static_cast<S*>(gate);
                std::string error("Gates_block::convert_to_DFE_gates: S_gate not implemented");
                throw error;       
            }
            else if (gate->get_type() == SDG_OPERATION) {
                // get the inverse parameters of the U3 rotation
                SDG* s_gate = static_cast<SDG*>(gate);
                std::string error("Gates_block::convert_to_DFE_gates: Sdg_gate not implemented");
                throw error;       
            }
            else if (gate->get_type() == SX_OPERATION) {
                // get the inverse parameters of the U3 rotation
                SX* sx_gate = static_cast<SX*>(gate);
                std::string error("Gates_block::convert_to_DFE_gates: SX_gate not implemented");
                throw error;       
            }
            else if (gate->get_type() == BLOCK_OPERATION) {
                Gates_block* block_gate = static_cast<Gates_block*>(gate);
                const Matrix_real parameters_layer_mtx(parameters_mtx.get_data() + parameter_idx - gate->get_parameter_num(), 1, gate->get_parameter_num() );
                block_gate->convert_to_DFE_gates( parameters_layer_mtx, DFEgates, gate_idx );
                parameter_idx = parameter_idx - block_gate->get_parameter_num();
            }
            else if (gate->get_type() == GENERAL_OPERATION) {
                std::string error("Gates_block::convert_to_DFE_gates: general_gate not implemented");
                throw error;
            }
            else if (gate->get_type() == ADAPTIVE_OPERATION) {
                // definig the rotation parameter
                double Phi;
                // get the inverse parameters of the U3 rotation
                Adaptive* ad_gate = static_cast<Adaptive*>(gate);
                double varthetaOver2 = std::fmod( activation_function(parameters_data[parameter_idx-1], ad_gate->get_limit()), 2*M_PI);
                parameter_idx = parameter_idx - 1;
                DFEGate.target_qbit = ad_gate->get_target_qbit();
                DFEGate.control_qbit = ad_gate->get_control_qbit();
                DFEGate.gate_type = ADAPTIVE_OPERATION;
                DFEGate.ThetaOver2 = (int32_t)(varthetaOver2*(1<<25)); 
                DFEGate.Phi = (int32_t)(0);  
                DFEGate.Lambda = (int32_t)(0); 
                DFEGate.metadata = 0;
                                   
                gate_idx = gate_idx + 1;
            }
            else {
                std::string err("Gates_block::convert_to_DFE_gates: unimplemented gate"); 
                throw err;
            }

        }


    return;

}


#endif


#ifdef __GROQ__

void 
Gates_block::extract_gate_kernels_target_and_control_qubits(std::vector<Matrix> &u3_qbit, std::vector<int> &target_qbit, std::vector<int> &control_qbit, Matrix_real& parameters_mtx)
{   

    if ( u3_qbit.size() == 0 ) {
        std::map<std::string, int> gate_nums  = get_gate_nums();
        int gates_total_num = 0;
        for( auto it=gate_nums.begin(); it != gate_nums.end(); it++ ) {
            gates_total_num = gates_total_num + it->second;
        }

        u3_qbit.reserve( gates_total_num );
        target_qbit.reserve( gates_total_num );
        control_qbit.reserve( gates_total_num );
    }
/*
    u3_qbit.reserve(u3_qbit.capacity() + gates.size());
    target_qbit.reserve(target_qbit.capacity() + gates.size());
    control_qbit.reserve(control_qbit.capacity() + gates.size());
*/
    double* parameters = parameters_mtx.get_data();


    for( int idx=0; idx<gates.size(); idx++) {
        Gate* operation = gates[idx];
        //parameters = parameters + operation->get_parameter_num();

        Matrix_real params_mtx(parameters + operation->get_parameter_start_idx(), 1, operation->get_parameter_num());

        switch (operation->get_type()) {
        case CNOT_OPERATION: case CZ_OPERATION:
        case CH_OPERATION: {
            Matrix_real empty_params(0, 0);
            u3_qbit.push_back(operation->gate_kernel(empty_params));
            break;    
        }
        case H_OPERATION: {
            Matrix_real empty_params(0, 0);
            u3_qbit.push_back(operation->gate_kernel(empty_params));
            break;
        }
        case X_OPERATION: {
            Matrix_real empty_params(0, 0);
            u3_qbit.push_back(operation->gate_kernel(empty_params));
            break;
        }
        case Y_OPERATION: {
            Matrix_real empty_params(0, 0);
            u3_qbit.push_back(operation->gate_kernel(empty_params));
            break;
        }
        case Z_OPERATION: {
            Matrix_real empty_params(0, 0);
            u3_qbit.push_back(operation->gate_kernel(empty_params));
            break;
        }
        case T_OPERATION: {
            Matrix_real empty_params(0, 0);
            u3_qbit.push_back(operation->gate_kernel(empty_params));
            break;
        }
        case S_OPERATION: {
            Matrix_real empty_params(0, 0);
            u3_qbit.push_back(operation->gate_kernel(empty_params));
            break;
        }
        case SDG_OPERATION: {
            Matrix_real empty_params(0, 0);
            u3_qbit.push_back(operation->gate_kernel(empty_params));
            break;
        }
        case TDG_OPERATION: {
            Matrix_real empty_params(0, 0);
            u3_qbit.push_back(operation->gate_kernel(empty_params));
            break;
        }
        case SX_OPERATION: {
            Matrix_real empty_params(0, 0);
            u3_qbit.push_back(operation->gate_kernel(empty_params));
            break;
        }
        case SXDG_OPERATION: {
            Matrix_real empty_params(0, 0);
            u3_qbit.push_back(operation->gate_kernel(empty_params));
            break;
        }
        case U1_OPERATION: {
            u3_qbit.push_back(Gate::calc_one_qubit_u3(0.0, 0.0, params_mtx[0]));
            break;
        }
        case U2_OPERATION: {
            u3_qbit.push_back(Gate::calc_one_qubit_u3(M_PI / 4.0, params_mtx[0], params_mtx[1]));
            break;
        }
        case U3_OPERATION: {
            u3_qbit.push_back(Gate::calc_one_qubit_u3(params_mtx[0], params_mtx[1], params_mtx[2]));
            break;
        }
        case RX_OPERATION: {
            // set static values for the angles
            double ThetaOver2 = params_mtx[0];
            double Phi    = -M_PI/2;
            double Lambda = M_PI/2;
            u3_qbit.push_back(Gate::calc_one_qubit_u3(ThetaOver2, Phi, Lambda));
            break;
        }
        case RY_OPERATION: {
            double ThetaOver2 = params_mtx[0];
            double Phi    = 0.0;
            double Lambda = 0.0;
            u3_qbit.push_back(Gate::calc_one_qubit_u3(ThetaOver2, Phi, Lambda));
            break;
        }
        case CRY_OPERATION: {
            double ThetaOver2 = params_mtx[0];
            double Phi = 0.0;
            double Lambda = 0.0;
            u3_qbit.push_back(Gate::calc_one_qubit_u3(ThetaOver2, Phi, Lambda));
            break;
        }
        case RZ_OPERATION: {
            u3_qbit.push_back(Gate::calc_one_qubit_u3(0.0, 0.0, params_mtx[0]));
            break;
        }
        case BLOCK_OPERATION: {
        
            Gates_block* block_operation = static_cast<Gates_block*>(operation);
            //std::vector<Gate*> involved_gates = block_operation->get_gates();

            std::vector<int>&& involved_qubits = block_operation->get_involved_qubits();
            
            if  ( involved_qubits.size() == 1 && block_operation->gates.size() > 1 && block_operation->get_qbit_num() > 1 ) {           
                // possibly merge successive single qubit gates        
    
                Gates_block gates_block_mini = Gates_block(1);
  
                for (int idx=0; idx<block_operation->gates.size(); idx++){
                    Gate* gate = block_operation->gates[idx]->clone();
                    gate->set_target_qbit(0);
                    gate->set_qbit_num(1);

                    gates_block_mini.add_gate(gate);
                }
                
                Matrix merged_kernel = create_identity(2);
             
                gates_block_mini.apply_to(params_mtx, merged_kernel);
                
                u3_qbit.push_back( merged_kernel );
                target_qbit.push_back( involved_qubits[0] );
                control_qbit.push_back( -1 );

            }
            else {                
                block_operation->extract_gate_kernels_target_and_control_qubits(u3_qbit, target_qbit, control_qbit, params_mtx);
            }
            continue;
        }
        //case ADAPTIVE_OPERATION:
        //case SYC_OPERATION:
        //case GENERAL_OPERATION:
        default:
            std::string err("Optimization_Interface::apply_to: unimplemented gate (" + std::to_string(operation->get_type()) + ")"); 
            throw err;
        }

        target_qbit.push_back(operation->get_target_qbit());
        control_qbit.push_back(operation->get_control_qbit());
    }    



}
#endif



void 
export_gate_list_to_binary(Matrix_real& parameters, Gates_block* gates_block, const std::string& filename, int verbosity) {

    std::stringstream sstream;
    sstream << "Exporting circuit into binary format. Filename: " << filename << std::endl;
    logging log;
    log.verbose = verbosity;
    log.print(sstream, 3);    

    FILE* pFile;
    const char* c_filename = filename.c_str();
    
#ifdef _WIN32
    errno_t err = fopen_s(&pFile, c_filename, "wb");
    if (err != 0) {
        pFile = NULL;
    }
#else
    pFile = fopen(c_filename, "wb");
#endif
    if (pFile==NULL) {fputs ("File error",stderr); exit (1);}

    export_gate_list_to_binary( parameters, gates_block, pFile, verbosity );

    fclose(pFile);
    return;

}



void 
export_gate_list_to_binary(Matrix_real& parameters, Gates_block* gates_block, FILE* pFile, int verbosity) {

    int qbit_num = gates_block->get_qbit_num();   
    fwrite(&qbit_num, sizeof(int), 1, pFile);

    int parameter_num = gates_block->get_parameter_num();
    fwrite(&parameter_num, sizeof(int), 1, pFile);


    int gates_num = gates_block->get_gate_num();
    fwrite(&gates_num, sizeof(int), 1, pFile);



    std::vector<Gate*> gates = gates_block->get_gates();
    double* parameters_data = parameters.get_data();

    for ( std::vector<Gate*>::iterator it=gates.begin(); it != gates.end(); ++it ) {
        Gate* op = *it;

        gate_type gt_type = op->get_type();

        fwrite(&gt_type, sizeof(gate_type), 1, pFile);

        int parameter_num = op->get_parameter_num();

        if (gt_type == CNOT_OPERATION || gt_type == CZ_OPERATION || gt_type == CH_OPERATION || gt_type == SYC_OPERATION) {
            int target_qbit = op->get_target_qbit();
            int control_qbit = op->get_control_qbit();
            fwrite(&target_qbit, sizeof(int), 1, pFile);
            fwrite(&control_qbit, sizeof(int), 1, pFile);
        }
        else if (gt_type == U1_OPERATION || gt_type == U2_OPERATION || gt_type == U3_OPERATION ) {
            int target_qbit = op->get_target_qbit();
            fwrite(&target_qbit, sizeof(int), 1, pFile);
            fwrite(parameters_data, sizeof(double), parameter_num, pFile);
        }
        else if (gt_type == RX_OPERATION || gt_type == RY_OPERATION || gt_type == RZ_OPERATION ) {
            int target_qbit = op->get_target_qbit();
            fwrite(&target_qbit, sizeof(int), 1, pFile);

            fwrite(parameters_data, sizeof(double), parameter_num, pFile);
        }
        else if (gt_type == CRY_OPERATION) {
            int target_qbit = op->get_target_qbit();
            int control_qbit = op->get_control_qbit();
            fwrite(&target_qbit, sizeof(int), 1, pFile);
            fwrite(&control_qbit, sizeof(int), 1, pFile);

            fwrite(parameters_data, sizeof(double), parameter_num, pFile);
        }
        
        else if (gt_type == X_OPERATION || gt_type == Y_OPERATION || gt_type == Z_OPERATION || gt_type == S_OPERATION || gt_type == SDG_OPERATION || gt_type == SX_OPERATION || gt_type == H_OPERATION) {
            int target_qbit = op->get_target_qbit();
            fwrite(&target_qbit, sizeof(int), 1, pFile);
        }
        else if (gt_type == BLOCK_OPERATION) {
            Gates_block* block_op = static_cast<Gates_block*>( op );
            Matrix_real parameters_loc(parameters_data, 1, parameter_num);
            export_gate_list_to_binary( parameters_loc, block_op, pFile );

        }
        else if (gt_type == ADAPTIVE_OPERATION) {
            int target_qbit = op->get_target_qbit();
            int control_qbit = op->get_control_qbit();
            fwrite(&target_qbit, sizeof(int), 1, pFile);
            fwrite(&control_qbit, sizeof(int), 1, pFile);

            fwrite(parameters_data, sizeof(double), parameter_num, pFile);
        }
        else {
            std::string err("export_gate_list_to_binary: unimplemented gate"); 
            throw err;
        }


        parameters_data = parameters_data + parameter_num;

    }

}


Gates_block* import_gate_list_from_binary(Matrix_real& parameters, const std::string& filename, int verbosity) {

    std::stringstream sstream;
    sstream << "Importing quantum circuit from binary file " << filename << std::endl;    
    logging log;
    log.verbose = verbosity;
    log.print(sstream, 2);    

    FILE* pFile;
    const char* c_filename = filename.c_str();
    
#ifdef _WIN32
    errno_t err = fopen_s(&pFile, c_filename, "rb");
    if (err != 0) {
        pFile = NULL;
    }
#else
    pFile = fopen(c_filename, "rb");
#endif
    if (pFile==NULL) {fputs ("File error",stderr); exit (1);}

    Gates_block* ret = import_gate_list_from_binary(parameters, pFile, verbosity);

    fclose(pFile);
    return ret;
}




Gates_block* import_gate_list_from_binary(Matrix_real& parameters, FILE* pFile, int verbosity) {

    std::stringstream sstream;

    int qbit_num;

    fread_wrapper(&qbit_num, sizeof(int), 1, pFile);
    sstream << "qbit_num: " << qbit_num << std::endl;
    Gates_block* gate_block = new Gates_block(qbit_num);

    int parameter_num;
    fread_wrapper(&parameter_num, sizeof(int), 1, pFile);
    sstream << "parameter_num: " << parameter_num << std::endl;
    parameters = Matrix_real(1, parameter_num);
    double* parameters_data = parameters.get_data();

    int gates_num;
    fread_wrapper(&gates_num, sizeof(int), 1, pFile);
    sstream << "gates_num: " << gates_num << std::endl;

    std::vector<int> gate_block_level_gates_num;
    std::vector<Gates_block*> gate_block_levels;
    gate_block_level_gates_num.push_back( gates_num );
    gate_block_levels.push_back(gate_block);
    int current_level = 0;

    

    int iter_max = 100000;
    int iter = 0;
    while ( gate_block_level_gates_num[0] > 0 && iter < iter_max) {

        gate_type gt_type;
        fread_wrapper(&gt_type, sizeof(gate_type), 1, pFile);

        //std::cout << "gate type: " << gt_type << std::endl;

        if (gt_type == CNOT_OPERATION) {
            sstream << "importing CNOT gate" << std::endl;

            int target_qbit;
            fread_wrapper(&target_qbit, sizeof(int), 1, pFile);
            sstream << "target_qbit: " << target_qbit << std::endl;

            int control_qbit;
            fread_wrapper(&control_qbit, sizeof(int), 1, pFile);
            sstream << "control_qbit: " << control_qbit << std::endl;

            gate_block_levels[current_level]->add_cnot(target_qbit, control_qbit);
            gate_block_level_gates_num[current_level]--;
        }
        else if (gt_type == CZ_OPERATION) {
            sstream << "importing CZ gate" << std::endl;

            int target_qbit;
            fread_wrapper(&target_qbit, sizeof(int), 1, pFile);
            sstream << "target_qbit: " << target_qbit << std::endl;

            int control_qbit;
            fread_wrapper(&control_qbit, sizeof(int), 1, pFile);
            sstream << "control_qbit: " << control_qbit << std::endl;

            gate_block_levels[current_level]->add_cz(target_qbit, control_qbit);
            gate_block_level_gates_num[current_level]--;
        }
        else if (gt_type == CH_OPERATION) {
            sstream << "importing CH gate" << std::endl;

            int target_qbit;
            fread_wrapper(&target_qbit, sizeof(int), 1, pFile);
            sstream << "target_qbit: " << target_qbit << std::endl;

            int control_qbit;
            fread_wrapper(&control_qbit, sizeof(int), 1, pFile);
            sstream << "control_qbit: " << control_qbit << std::endl;

            gate_block_levels[current_level]->add_ch(target_qbit, control_qbit);
            gate_block_level_gates_num[current_level]--;
        }
        else if (gt_type == SYC_OPERATION) {
            sstream << "importing SYCAMORE gate" << std::endl;

            int target_qbit;
            fread_wrapper(&target_qbit, sizeof(int), 1, pFile);
            sstream << "target_qbit: " << target_qbit << std::endl;

            int control_qbit;
            fread_wrapper(&control_qbit, sizeof(int), 1, pFile);
            sstream << "control_qbit: " << control_qbit << std::endl;

            gate_block_levels[current_level]->add_syc(target_qbit, control_qbit);
            gate_block_level_gates_num[current_level]--;
        }
        else if (gt_type == U1_OPERATION) {
            sstream << "importing U1 gate" << std::endl;
            
            int target_qbit;
            fread_wrapper(&target_qbit, sizeof(int), 1, pFile);
            sstream << "target_qbit: " << target_qbit << std::endl;

            fread_wrapper(parameters_data, sizeof(double), 1, pFile);
            
            parameters_data++;
            gate_block_levels[current_level]->add_u1(target_qbit);
            gate_block_level_gates_num[current_level]--;
        }
        else if (gt_type == U2_OPERATION) {
        sstream << "importing U2 gate" << std::endl;
            
            int target_qbit;
            fread_wrapper(&target_qbit, sizeof(int), 1, pFile);
            sstream << "target_qbit: " << target_qbit << std::endl;

            fread_wrapper(parameters_data, sizeof(double), 2, pFile);
            
            parameters_data+=2;
            gate_block_levels[current_level]->add_u2(target_qbit);
            gate_block_level_gates_num[current_level]--;
        }
        else if (gt_type == U3_OPERATION) {
            sstream << "importing U3 gate" << std::endl;

            int target_qbit;
            fread_wrapper(&target_qbit, sizeof(int), 1, pFile);
            sstream << "target_qbit: " << target_qbit << std::endl;

            fread_wrapper(parameters_data, sizeof(double), 3, pFile);
            
            parameters_data += 3;
            gate_block_levels[current_level]->add_u3(target_qbit);
            gate_block_level_gates_num[current_level]--;
        }
        else if (gt_type == RX_OPERATION) {

            sstream << "importing RX gate" << std::endl;

            int target_qbit;
            fread_wrapper(&target_qbit, sizeof(int), 1, pFile);
            sstream << "target_qbit: " << target_qbit << std::endl;

            fread_wrapper(parameters_data, sizeof(double), 1, pFile);
            parameters_data++;

            gate_block_levels[current_level]->add_rx(target_qbit);
            gate_block_level_gates_num[current_level]--;

        }
        else if (gt_type == RY_OPERATION) {

            sstream << "importing RY gate" << std::endl;

            int target_qbit;
            fread_wrapper(&target_qbit, sizeof(int), 1, pFile);
            sstream << "target_qbit: " << target_qbit << std::endl;

            fread_wrapper(parameters_data, sizeof(double), 1, pFile);
            parameters_data++;

            gate_block_levels[current_level]->add_ry(target_qbit);
            gate_block_level_gates_num[current_level]--;

        }
        else if (gt_type == CRY_OPERATION) {

            sstream << "importing CRY gate" << std::endl;

            int target_qbit;
            fread_wrapper(&target_qbit, sizeof(int), 1, pFile);
            sstream << "target_qbit: " << target_qbit << std::endl;

            int control_qbit;
            fread_wrapper(&control_qbit, sizeof(int), 1, pFile);
            sstream << "control_qbit: " << control_qbit << std::endl;

            fread_wrapper(parameters_data, sizeof(double), 1, pFile);
            parameters_data++;

            gate_block_levels[current_level]->add_cry(target_qbit, control_qbit);
            gate_block_level_gates_num[current_level]--;

        }
        else if (gt_type == RZ_OPERATION) {

            sstream << "importing RZ gate" << std::endl;

            int target_qbit;
            fread_wrapper(&target_qbit, sizeof(int), 1, pFile);
            sstream << "target_qbit: " << target_qbit << std::endl;

            fread_wrapper(parameters_data, sizeof(double), 1, pFile);
            parameters_data++;

            gate_block_levels[current_level]->add_rz(target_qbit);
            gate_block_level_gates_num[current_level]--;

        }
        else if (gt_type == H_OPERATION) {

            sstream << "importing Hadamard gate" << std::endl;

            int target_qbit;
            fread_wrapper(&target_qbit, sizeof(int), 1, pFile);
            sstream << "target_qbit: " << target_qbit << std::endl;

            gate_block_levels[current_level]->add_h(target_qbit);
            gate_block_level_gates_num[current_level]--;

        }
        else if (gt_type == X_OPERATION) {

            sstream << "importing X gate" << std::endl;

            int target_qbit;
            fread_wrapper(&target_qbit, sizeof(int), 1, pFile);
            sstream << "target_qbit: " << target_qbit << std::endl;

            gate_block_levels[current_level]->add_x(target_qbit);
            gate_block_level_gates_num[current_level]--;

        }
        else if (gt_type == Y_OPERATION) {

            sstream << "importing Y gate" << std::endl;

            int target_qbit;
            fread_wrapper(&target_qbit, sizeof(int), 1, pFile);
            sstream << "target_qbit: " << target_qbit << std::endl;

            gate_block_levels[current_level]->add_y(target_qbit);
            gate_block_level_gates_num[current_level]--;

        }
        else if (gt_type == Z_OPERATION) {

            sstream << "importing Z gate" << std::endl;

            int target_qbit;
            fread_wrapper(&target_qbit, sizeof(int), 1, pFile);
            sstream << "target_qbit: " << target_qbit << std::endl;

            gate_block_levels[current_level]->add_z(target_qbit);
            gate_block_level_gates_num[current_level]--;

        }
        else if (gt_type == S_OPERATION) {

            sstream << "importing S gate" << std::endl;

            int target_qbit;
            fread_wrapper(&target_qbit, sizeof(int), 1, pFile);
            sstream << "target_qbit: " << target_qbit << std::endl;

            gate_block_levels[current_level]->add_s(target_qbit);
            gate_block_level_gates_num[current_level]--;

        }
        else if (gt_type == SDG_OPERATION) {

            sstream << "importing Sdg gate" << std::endl;

            int target_qbit;
            fread_wrapper(&target_qbit, sizeof(int), 1, pFile);
            sstream << "target_qbit: " << target_qbit << std::endl;

            gate_block_levels[current_level]->add_sdg(target_qbit);
            gate_block_level_gates_num[current_level]--;

        }
        else if (gt_type == T_OPERATION) {

            sstream << "importing T gate" << std::endl;

            int target_qbit;
            fread_wrapper(&target_qbit, sizeof(int), 1, pFile);
            sstream << "target_qbit: " << target_qbit << std::endl;

            gate_block_levels[current_level]->add_t(target_qbit);
            gate_block_level_gates_num[current_level]--;

        }
        else if (gt_type == TDG_OPERATION) {

            sstream << "importing Tdg gate" << std::endl;

            int target_qbit;
            fread_wrapper(&target_qbit, sizeof(int), 1, pFile);
            sstream << "target_qbit: " << target_qbit << std::endl;

            gate_block_levels[current_level]->add_t(target_qbit);
            gate_block_level_gates_num[current_level]--;

        }
        else if (gt_type == SX_OPERATION) {

            sstream << "importing SX gate" << std::endl;

            int target_qbit;
            fread_wrapper(&target_qbit, sizeof(int), 1, pFile);
            sstream << "target_qbit: " << target_qbit << std::endl;

            gate_block_levels[current_level]->add_sx(target_qbit);
            gate_block_level_gates_num[current_level]--;

        }
        else if (gt_type == BLOCK_OPERATION) {


            sstream << "******* importing gates block ********" << std::endl;

            int qbit_num_loc;
            fread_wrapper(&qbit_num_loc, sizeof(int), 1, pFile);
            //std::cout << "qbit_num_loc: " << qbit_num_loc << std::endl;
            Gates_block* gate_block_inner = new Gates_block(qbit_num_loc);

            int parameter_num_loc;
            fread_wrapper(&parameter_num_loc, sizeof(int), 1, pFile);
            //std::cout << "parameter_num_loc: " << parameter_num_loc << std::endl;
        

            int gates_num_loc;
            fread_wrapper(&gates_num_loc, sizeof(int), 1, pFile);
            //std::cout << "gates_num_loc: " << gates_num_loc << std::endl;
                        
            gate_block_levels.push_back( gate_block_inner );
            gate_block_level_gates_num.push_back(gates_num_loc);
            current_level++;
        }
        else if (gt_type == ADAPTIVE_OPERATION) {

            sstream << "importing adaptive gate" << std::endl;

            int target_qbit;
            fread_wrapper(&target_qbit, sizeof(int), 1, pFile);
            sstream << "target_qbit: " << target_qbit << std::endl;

            int control_qbit;
            fread_wrapper(&control_qbit, sizeof(int), 1, pFile);
            sstream << "control_qbit: " << control_qbit << std::endl;

            fread_wrapper(parameters_data, sizeof(double), 1, pFile);
            parameters_data++;

            gate_block_levels[current_level]->add_adaptive(target_qbit, control_qbit);
            gate_block_level_gates_num[current_level]--;

        }
        else {
            std::string err("import_gate_list_from_binary: unimplemented gate"); 
            throw err;
        }


        if ( gate_block_level_gates_num[current_level] == 0 ) {
            gate_block_levels[ current_level-1 ]->add_gate( static_cast<Gate*>(gate_block_levels[ current_level ]) );
            gate_block_levels.pop_back();
            gate_block_level_gates_num.pop_back();
            current_level--;
            gate_block_level_gates_num[current_level]--;
            sstream << "finishing gates block" << std::endl;
        }


        iter++;
    }

    logging log;
    log.verbose = verbosity;
    log.print(sstream, 4);    

  
    if ( iter == iter_max ) {
        std::string error("Corrupted input file, reached end of the file before contructing the whole gate structure");
        throw error;
    }

    return gate_block;

}



Matrix_real reverse_parameters( const Matrix_real& parameters_in, std::vector<Gate*>::iterator gates_it, int num_of_gates ) {

     if ( parameters_in.cols > 1 && parameters_in.rows > 1 ) {
          std::string error("reverse_parameters: Input array should have a single column or a single row.");
          throw error;
     }
     
//return parameters_in.copy();
     
    // determine the number of parameters
    int parameters_num_total = 0;
    for (int idx=0; idx<num_of_gates; idx++) {

        // The current gate
        Gate* gate = *(gates_it++);
        parameters_num_total = parameters_num_total + gate->get_parameter_num();

    }
    
    
     if ( parameters_num_total  == 0) {
          return Matrix_real(0,0);
     }        
     
     //std::cout << "uu" << std::endl;
     //parameters_in.print_matrix();
    

     Matrix_real parameters_ret(1, parameters_num_total);
     int parameter_num_copied = 0;

    // apply the gate operations on the inital matrix
    for (int idx=num_of_gates-0; idx>0; idx--) {

        // The current gate
        Gate* gate = *(--gates_it);
        
        int parameter_num_gate = gate->get_parameter_num();
        
        if ( parameter_num_gate == 0 ) {
          continue;
        }     
           
        else if (gate->get_type() == BLOCK_OPERATION ) {
        
         //std::cout << "block: " << parameter_num_gate << " " << parameters_num_total << std::endl;
              parameters_num_total = parameters_num_total - gate->get_parameter_num();
                     
            Matrix_real parameters_of_block( parameters_in.get_data()+parameters_num_total, 1, parameter_num_gate );
            //parameters_of_block.print_matrix();
          
            Gates_block* block_gate = static_cast<Gates_block*>( gate );
            
            std::vector<Gate*> gates_loc = block_gate->get_gates();
            
            Matrix_real parameters_of_block_reversed = reverse_parameters( parameters_of_block, gates_loc.begin(), static_cast<int>(gates_loc.size()) );
            
            //parameters_of_block_reversed.print_matrix();
            
            memcpy( parameters_ret.get_data()+parameter_num_copied, parameters_of_block_reversed.get_data(), parameters_of_block_reversed.size()*sizeof(double) );
            parameter_num_copied = parameter_num_copied + parameters_of_block_reversed.size();
               
            //parameters_ret.print_matrix();
               
        }
        
        else {
        
             //std::cout << parameter_num_gate << std::endl;
        
         parameters_num_total = parameters_num_total - gate->get_parameter_num();
            memcpy( parameters_ret.get_data()+parameter_num_copied, parameters_in.get_data()+parameters_num_total, gate->get_parameter_num()*sizeof(double) );
            parameter_num_copied = parameter_num_copied + gate->get_parameter_num();
               
       }


    }



     
     return parameters_ret;
     
     


}




Matrix_real inverse_reverse_parameters( const Matrix_real& parameters_in, std::vector<Gate*>::iterator gates_it, int num_of_gates ) {

     if ( parameters_in.cols > 1 && parameters_in.rows > 1 ) {
          std::string error("reverse_parameters: Input array should have a single column or a single row.");
          throw error;
     }
     
     
    // determine the number of parameters
    int parameters_num_total = 0;
    for (int idx=0; idx<num_of_gates; idx++) {

        // The current gate
        Gate* gate = *(gates_it++);
        parameters_num_total = parameters_num_total + gate->get_parameter_num();

    }
    
    
     if ( parameters_num_total  == 0) {
          return Matrix_real(0,0);
     }        
     
     //std::cout << "uu" << std::endl;
     //parameters_in.print_matrix();
    

     Matrix_real parameters_ret(1, parameters_num_total);
     int parameter_num_copied = 0;

    // apply the gate operations on the inital matrix
    for (int idx=num_of_gates-0; idx>0; idx--) {

        // The current gate
        Gate* gate = *(--gates_it);
        
        int parameter_num_gate = gate->get_parameter_num();
        
        if ( parameter_num_gate == 0 ) {
          continue;
        }     
         
        else if (gate->get_type() == BLOCK_OPERATION ) {
        
            //std::cout << "block: " << parameter_num_gate << " " << parameters_num_total << std::endl;
              parameters_num_total = parameters_num_total - gate->get_parameter_num();
                     
            Matrix_real parameters_of_block( parameters_in.get_data()+parameters_num_total, 1, parameter_num_gate );
            //parameters_of_block.print_matrix();
          
            Gates_block* block_gate = static_cast<Gates_block*>( gate );
            
            std::vector<Gate*> gates_loc = block_gate->get_gates();
            
            Matrix_real parameters_of_block_reversed = reverse_parameters( parameters_of_block, gates_loc.begin(), static_cast<int>(gates_loc.size()) );
            
            //parameters_of_block_reversed.print_matrix();
            
            memcpy( parameters_ret.get_data()+parameter_num_copied, parameters_of_block_reversed.get_data(), parameters_of_block_reversed.size()*sizeof(double) );
            parameter_num_copied = parameter_num_copied + parameters_of_block_reversed.size();
               
           //parameters_ret.print_matrix();
               
        }
        
        else {
        
             //std::cout << parameter_num_gate << std::endl;
        
         parameters_num_total = parameters_num_total - gate->get_parameter_num();
            memcpy( parameters_ret.get_data()+parameters_num_total, parameters_in.get_data()+parameter_num_copied, gate->get_parameter_num()*sizeof(double) );
            parameter_num_copied = parameter_num_copied + gate->get_parameter_num();
               
          }


    }



     
     return parameters_ret;
     
     


}