text:
# Automatically generated using Clang.jl

const SECRET_NUM = -1234

@cenum UNUSED_ENUM_TYPE::UInt32 begin
    UNUSED_DEF_VAL = 0
end

@cenum LAYER_TYPE::UInt32 begin
    CONVOLUTIONAL = 0
    DECONVOLUTIONAL = 1
    CONNECTED = 2
    MAXPOOL = 3
    LOCAL_AVGPOOL = 4
    SOFTMAX = 5
    DETECTION = 6
    DROPOUT = 7
    CROP = 8
    ROUTE = 9
    COST = 10
    NORMALIZATION = 11
    AVGPOOL = 12
    LOCAL = 13
    SHORTCUT = 14
    SCALE_CHANNELS = 15
    SAM = 16
    ACTIVE = 17
    RNN = 18
    GRU = 19
    LSTM = 20
    CONV_LSTM = 21
    CRNN = 22
    BATCHNORM = 23
    NETWORK = 24
    XNOR = 25
    REGION = 26
    YOLO = 27
    GAUSSIAN_YOLO = 28
    ISEG = 29
    REORG = 30
    REORG_OLD = 31
    UPSAMPLE = 32
    LOGXENT = 33
    L2NORM = 34
    EMPTY = 35
    BLANK = 36
end

@cenum ACTIVATION::UInt32 begin
    LOGISTIC = 0
    RELU = 1
    RELU6 = 2
    RELIE = 3
    LINEAR = 4
    RAMP = 5
    TANH = 6
    PLSE = 7
    LEAKY = 8
    ELU = 9
    LOGGY = 10
    STAIR = 11
    HARDTAN = 12
    LHTAN = 13
    SELU = 14
    GELU = 15
    SWISH = 16
    MISH = 17
    NORM_CHAN = 18
    NORM_CHAN_SOFTMAX = 19
    NORM_CHAN_SOFTMAX_MAXVAL = 20
end

@cenum COST_TYPE::UInt32 begin
    SSE = 0
    MASKED = 1
    L1 = 2
    SEG = 3
    SMOOTH = 4
    WGAN = 5
end

@cenum WEIGHTS_TYPE_T::UInt32 begin
    NO_WEIGHTS = 0
    PER_FEATURE = 1
    PER_CHANNEL = 2
end

@cenum WEIGHTS_NORMALIZATION_T::UInt32 begin
    NO_NORMALIZATION = 0
    RELU_NORMALIZATION = 1
    SOFTMAX_NORMALIZATION = 2
end

@cenum IOU_LOSS::UInt32 begin
    IOU = 0
    GIOU = 1
    MSE = 2
    DIOU = 3
    CIOU = 4
end

@cenum NMS_KIND::UInt32 begin
    DEFAULT_NMS = 0
    GREEDY_NMS = 1
    DIOU_NMS = 2
    CORNERS_NMS = 3
end

@cenum YOLO_POINT::UInt32 begin
    YOLO_CENTER = 1
    YOLO_LEFT_TOP = 2
    YOLO_RIGHT_BOTTOM = 4
end

struct tree
    leaf::Ptr{Cint}
    n::Cint
    parent::Ptr{Cint}
    child::Ptr{Cint}
    group::Ptr{Cint}
    name::Ptr{Cstring}
    groups::Cint
    group_size::Ptr{Cint}
    group_offset::Ptr{Cint}
end

struct layer
    type::LAYER_TYPE
    activation::ACTIVATION
    lstm_activation::ACTIVATION
    cost_type::COST_TYPE
    forward::Ptr{Cvoid}
    backward::Ptr{Cvoid}
    update::Ptr{Cvoid}
    forward_gpu::Ptr{Cvoid}
    backward_gpu::Ptr{Cvoid}
    update_gpu::Ptr{Cvoid}
    share_layer::Ptr{layer}
    train::Cint
    avgpool::Cint
    batch_normalize::Cint
    shortcut::Cint
    batch::Cint
    dynamic_minibatch::Cint
    forced::Cint
    flipped::Cint
    inputs::Cint
    outputs::Cint
    mean_alpha::Cfloat
    nweights::Cint
    nbiases::Cint
    extra::Cint
    truths::Cint
    h::Cint
    w::Cint
    c::Cint
    out_h::Cint
    out_w::Cint
    out_c::Cint
    n::Cint
    max_boxes::Cint
    groups::Cint
    group_id::Cint
    size::Cint
    side::Cint
    stride::Cint
    stride_x::Cint
    stride_y::Cint
    dilation::Cint
    antialiasing::Cint
    maxpool_depth::Cint
    out_channels::Cint
    reverse::Cint
    flatten::Cint
    spatial::Cint
    pad::Cint
    sqrt::Cint
    flip::Cint
    index::Cint
    scale_wh::Cint
    binary::Cint
    xnor::Cint
    peephole::Cint
    use_bin_output::Cint
    keep_delta_gpu::Cint
    optimized_memory::Cint
    steps::Cint
    bottleneck::Cint
    time_normalizer::Cfloat
    state_constrain::Cint
    hidden::Cint
    truth::Cint
    smooth::Cfloat
    dot::Cfloat
    deform::Cint
    grad_centr::Cint
    sway::Cint
    rotate::Cint
    stretch::Cint
    stretch_sway::Cint
    angle::Cfloat
    jitter::Cfloat
    resize::Cfloat
    saturation::Cfloat
    exposure::Cfloat
    shift::Cfloat
    ratio::Cfloat
    learning_rate_scale::Cfloat
    clip::Cfloat
    focal_loss::Cint
    classes_multipliers::Ptr{Cfloat}
    label_smooth_eps::Cfloat
    noloss::Cint
    softmax::Cint
    classes::Cint
    coords::Cint
    background::Cint
    rescore::Cint
    objectness::Cint
    does_cost::Cint
    joint::Cint
    noadjust::Cint
    reorg::Cint
    log::Cint
    tanh::Cint
    mask::Ptr{Cint}
    total::Cint
    bflops::Cfloat
    adam::Cint
    B1::Cfloat
    B2::Cfloat
    eps::Cfloat
    t::Cint
    alpha::Cfloat
    beta::Cfloat
    kappa::Cfloat
    coord_scale::Cfloat
    object_scale::Cfloat
    noobject_scale::Cfloat
    mask_scale::Cfloat
    class_scale::Cfloat
    bias_match::Cint
    random::Cfloat
    ignore_thresh::Cfloat
    truth_thresh::Cfloat
    iou_thresh::Cfloat
    thresh::Cfloat
    focus::Cfloat
    classfix::Cint
    absolute::Cint
    assisted_excitation::Cint
    onlyforward::Cint
    stopbackward::Cint
    train_only_bn::Cint
    dont_update::Cint
    burnin_update::Cint
    dontload::Cint
    dontsave::Cint
    dontloadscales::Cint
    numload::Cint
    temperature::Cfloat
    probability::Cfloat
    dropblock_size_rel::Cfloat
    dropblock_size_abs::Cint
    dropblock::Cint
    scale::Cfloat
    receptive_w::Cint
    receptive_h::Cint
    receptive_w_scale::Cint
    receptive_h_scale::Cint
    cweights::Cstring
    indexes::Ptr{Cint}
    input_layers::Ptr{Cint}
    input_sizes::Ptr{Cint}
    layers_output::Ptr{Ptr{Cfloat}}
    layers_delta::Ptr{Ptr{Cfloat}}
    weights_type::WEIGHTS_TYPE_T
    weights_normalization::WEIGHTS_NORMALIZATION_T
    map::Ptr{Cint}
    counts::Ptr{Cint}
    sums::Ptr{Ptr{Cfloat}}
    rand::Ptr{Cfloat}
    cost::Ptr{Cfloat}
    state::Ptr{Cfloat}
    prev_state::Ptr{Cfloat}
    forgot_state::Ptr{Cfloat}
    forgot_delta::Ptr{Cfloat}
    state_delta::Ptr{Cfloat}
    combine_cpu::Ptr{Cfloat}
    combine_delta_cpu::Ptr{Cfloat}
    concat::Ptr{Cfloat}
    concat_delta::Ptr{Cfloat}
    binary_weights::Ptr{Cfloat}
    biases::Ptr{Cfloat}
    bias_updates::Ptr{Cfloat}
    scales::Ptr{Cfloat}
    scale_updates::Ptr{Cfloat}
    weights::Ptr{Cfloat}
    weight_updates::Ptr{Cfloat}
    scale_x_y::Cfloat
    objectness_smooth::Cint
    max_delta::Cfloat
    uc_normalizer::Cfloat
    iou_normalizer::Cfloat
    cls_normalizer::Cfloat
    iou_loss::IOU_LOSS
    iou_thresh_kind::IOU_LOSS
    nms_kind::NMS_KIND
    beta_nms::Cfloat
    yolo_point::YOLO_POINT
    align_bit_weights_gpu::Cstring
    mean_arr_gpu::Ptr{Cfloat}
    align_workspace_gpu::Ptr{Cfloat}
    transposed_align_workspace_gpu::Ptr{Cfloat}
    align_workspace_size::Cint
    align_bit_weights::Cstring
    mean_arr::Ptr{Cfloat}
    align_bit_weights_size::Cint
    lda_align::Cint
    new_lda::Cint
    bit_align::Cint
    col_image::Ptr{Cfloat}
    delta::Ptr{Cfloat}
    output::Ptr{Cfloat}
    activation_input::Ptr{Cfloat}
    delta_pinned::Cint
    output_pinned::Cint
    loss::Ptr{Cfloat}
    squared::Ptr{Cfloat}
    norms::Ptr{Cfloat}
    spatial_mean::Ptr{Cfloat}
    mean::Ptr{Cfloat}
    variance::Ptr{Cfloat}
    mean_delta::Ptr{Cfloat}
    variance_delta::Ptr{Cfloat}
    rolling_mean::Ptr{Cfloat}
    rolling_variance::Ptr{Cfloat}
    x::Ptr{Cfloat}
    x_norm::Ptr{Cfloat}
    m::Ptr{Cfloat}
    v::Ptr{Cfloat}
    bias_m::Ptr{Cfloat}
    bias_v::Ptr{Cfloat}
    scale_m::Ptr{Cfloat}
    scale_v::Ptr{Cfloat}
    z_cpu::Ptr{Cfloat}
    r_cpu::Ptr{Cfloat}
    h_cpu::Ptr{Cfloat}
    stored_h_cpu::Ptr{Cfloat}
    prev_state_cpu::Ptr{Cfloat}
    temp_cpu::Ptr{Cfloat}
    temp2_cpu::Ptr{Cfloat}
    temp3_cpu::Ptr{Cfloat}
    dh_cpu::Ptr{Cfloat}
    hh_cpu::Ptr{Cfloat}
    prev_cell_cpu::Ptr{Cfloat}
    cell_cpu::Ptr{Cfloat}
    f_cpu::Ptr{Cfloat}
    i_cpu::Ptr{Cfloat}
    g_cpu::Ptr{Cfloat}
    o_cpu::Ptr{Cfloat}
    c_cpu::Ptr{Cfloat}
    stored_c_cpu::Ptr{Cfloat}
    dc_cpu::Ptr{Cfloat}
    binary_input::Ptr{Cfloat}
    bin_re_packed_input::Ptr{UInt32}
    t_bit_input::Cstring
    input_layer::Ptr{layer}
    self_layer::Ptr{layer}
    output_layer::Ptr{layer}
    reset_layer::Ptr{layer}
    update_layer::Ptr{layer}
    state_layer::Ptr{layer}
    input_gate_layer::Ptr{layer}
    state_gate_layer::Ptr{layer}
    input_save_layer::Ptr{layer}
    state_save_layer::Ptr{layer}
    input_state_layer::Ptr{layer}
    state_state_layer::Ptr{layer}
    input_z_layer::Ptr{layer}
    state_z_layer::Ptr{layer}
    input_r_layer::Ptr{layer}
    state_r_layer::Ptr{layer}
    input_h_layer::Ptr{layer}
    state_h_layer::Ptr{layer}
    wz::Ptr{layer}
    uz::Ptr{layer}
    wr::Ptr{layer}
    ur::Ptr{layer}
    wh::Ptr{layer}
    uh::Ptr{layer}
    uo::Ptr{layer}
    wo::Ptr{layer}
    vo::Ptr{layer}
    uf::Ptr{layer}
    wf::Ptr{layer}
    vf::Ptr{layer}
    ui::Ptr{layer}
    wi::Ptr{layer}
    vi::Ptr{layer}
    ug::Ptr{layer}
    wg::Ptr{layer}
    softmax_tree::Ptr{tree}
    workspace_size::Cint
    indexes_gpu::Ptr{Cint}
    z_gpu::Ptr{Cfloat}
    r_gpu::Ptr{Cfloat}
    h_gpu::Ptr{Cfloat}
    stored_h_gpu::Ptr{Cfloat}
    bottelneck_hi_gpu::Ptr{Cfloat}
    bottelneck_delta_gpu::Ptr{Cfloat}
    temp_gpu::Ptr{Cfloat}
    temp2_gpu::Ptr{Cfloat}
    temp3_gpu::Ptr{Cfloat}
    dh_gpu::Ptr{Cfloat}
    hh_gpu::Ptr{Cfloat}
    prev_cell_gpu::Ptr{Cfloat}
    prev_state_gpu::Ptr{Cfloat}
    last_prev_state_gpu::Ptr{Cfloat}
    last_prev_cell_gpu::Ptr{Cfloat}
    cell_gpu::Ptr{Cfloat}
    f_gpu::Ptr{Cfloat}
    i_gpu::Ptr{Cfloat}
    g_gpu::Ptr{Cfloat}
    o_gpu::Ptr{Cfloat}
    c_gpu::Ptr{Cfloat}
    stored_c_gpu::Ptr{Cfloat}
    dc_gpu::Ptr{Cfloat}
    m_gpu::Ptr{Cfloat}
    v_gpu::Ptr{Cfloat}
    bias_m_gpu::Ptr{Cfloat}
    scale_m_gpu::Ptr{Cfloat}
    bias_v_gpu::Ptr{Cfloat}
    scale_v_gpu::Ptr{Cfloat}
    combine_gpu::Ptr{Cfloat}
    combine_delta_gpu::Ptr{Cfloat}
    forgot_state_gpu::Ptr{Cfloat}
    forgot_delta_gpu::Ptr{Cfloat}
    state_gpu::Ptr{Cfloat}
    state_delta_gpu::Ptr{Cfloat}
    gate_gpu::Ptr{Cfloat}
    gate_delta_gpu::Ptr{Cfloat}
    save_gpu::Ptr{Cfloat}
    save_delta_gpu::Ptr{Cfloat}
    concat_gpu::Ptr{Cfloat}
    concat_delta_gpu::Ptr{Cfloat}
    binary_input_gpu::Ptr{Cfloat}
    binary_weights_gpu::Ptr{Cfloat}
    bin_conv_shortcut_in_gpu::Ptr{Cfloat}
    bin_conv_shortcut_out_gpu::Ptr{Cfloat}
    mean_gpu::Ptr{Cfloat}
    variance_gpu::Ptr{Cfloat}
    m_cbn_avg_gpu::Ptr{Cfloat}
    v_cbn_avg_gpu::Ptr{Cfloat}
    rolling_mean_gpu::Ptr{Cfloat}
    rolling_variance_gpu::Ptr{Cfloat}
    variance_delta_gpu::Ptr{Cfloat}
    mean_delta_gpu::Ptr{Cfloat}
    col_image_gpu::Ptr{Cfloat}
    x_gpu::Ptr{Cfloat}
    x_norm_gpu::Ptr{Cfloat}
    weights_gpu::Ptr{Cfloat}
    weight_updates_gpu::Ptr{Cfloat}
    weight_deform_gpu::Ptr{Cfloat}
    weight_change_gpu::Ptr{Cfloat}
    weights_gpu16::Ptr{Cfloat}
    weight_updates_gpu16::Ptr{Cfloat}
    biases_gpu::Ptr{Cfloat}
    bias_updates_gpu::Ptr{Cfloat}
    bias_change_gpu::Ptr{Cfloat}
    scales_gpu::Ptr{Cfloat}
    scale_updates_gpu::Ptr{Cfloat}
    scale_change_gpu::Ptr{Cfloat}
    input_antialiasing_gpu::Ptr{Cfloat}
    output_gpu::Ptr{Cfloat}
    output_avg_gpu::Ptr{Cfloat}
    activation_input_gpu::Ptr{Cfloat}
    loss_gpu::Ptr{Cfloat}
    delta_gpu::Ptr{Cfloat}
    rand_gpu::Ptr{Cfloat}
    drop_blocks_scale::Ptr{Cfloat}
    drop_blocks_scale_gpu::Ptr{Cfloat}
    squared_gpu::Ptr{Cfloat}
    norms_gpu::Ptr{Cfloat}
    gt_gpu::Ptr{Cfloat}
    a_avg_gpu::Ptr{Cfloat}
    input_sizes_gpu::Ptr{Cint}
    layers_output_gpu::Ptr{Ptr{Cfloat}}
    layers_delta_gpu::Ptr{Ptr{Cfloat}}
    srcTensorDesc::Ptr{Cvoid}
    dstTensorDesc::Ptr{Cvoid}
    srcTensorDesc16::Ptr{Cvoid}
    dstTensorDesc16::Ptr{Cvoid}
    dsrcTensorDesc::Ptr{Cvoid}
    ddstTensorDesc::Ptr{Cvoid}
    dsrcTensorDesc16::Ptr{Cvoid}
    ddstTensorDesc16::Ptr{Cvoid}
    normTensorDesc::Ptr{Cvoid}
    normDstTensorDesc::Ptr{Cvoid}
    normDstTensorDescF16::Ptr{Cvoid}
    weightDesc::Ptr{Cvoid}
    weightDesc16::Ptr{Cvoid}
    dweightDesc::Ptr{Cvoid}
    dweightDesc16::Ptr{Cvoid}
    convDesc::Ptr{Cvoid}
    fw_algo::UNUSED_ENUM_TYPE
    fw_algo16::UNUSED_ENUM_TYPE
    bd_algo::UNUSED_ENUM_TYPE
    bd_algo16::UNUSED_ENUM_TYPE
    bf_algo::UNUSED_ENUM_TYPE
    bf_algo16::UNUSED_ENUM_TYPE
    poolingDesc::Ptr{Cvoid}
end

@cenum learning_rate_policy::UInt32 begin
    CONSTANT = 0
    STEP = 1
    EXP = 2
    POLY = 3
    STEPS = 4
    SIG = 5
    RANDOM = 6
    SGDR = 7
end

struct network
    n::Cint
    batch::Cint
    seen::Ptr{UInt64}
    cur_iteration::Ptr{Cint}
    loss_scale::Cfloat
    t::Ptr{Cint}
    epoch::Cfloat
    subdivisions::Cint
    layers::Ptr{layer}
    output::Ptr{Cfloat}
    policy::learning_rate_policy
    benchmark_layers::Cint
    learning_rate::Cfloat
    learning_rate_min::Cfloat
    learning_rate_max::Cfloat
    batches_per_cycle::Cint
    batches_cycle_mult::Cint
    momentum::Cfloat
    decay::Cfloat
    gamma::Cfloat
    scale::Cfloat
    power::Cfloat
    time_steps::Cint
    step::Cint
    max_batches::Cint
    num_boxes::Cint
    train_images_num::Cint
    seq_scales::Ptr{Cfloat}
    scales::Ptr{Cfloat}
    steps::Ptr{Cint}
    num_steps::Cint
    burn_in::Cint
    cudnn_half::Cint
    adam::Cint
    B1::Cfloat
    B2::Cfloat
    eps::Cfloat
    inputs::Cint
    outputs::Cint
    truths::Cint
    notruth::Cint
    h::Cint
    w::Cint
    c::Cint
    max_crop::Cint
    min_crop::Cint
    max_ratio::Cfloat
    min_ratio::Cfloat
    center::Cint
    flip::Cint
    gaussian_noise::Cint
    blur::Cint
    mixup::Cint
    label_smooth_eps::Cfloat
    resize_step::Cint
    attention::Cint
    adversarial::Cint
    adversarial_lr::Cfloat
    max_chart_loss::Cfloat
    letter_box::Cint
    angle::Cfloat
    aspect::Cfloat
    exposure::Cfloat
    saturation::Cfloat
    hue::Cfloat
    random::Cint
    track::Cint
    augment_speed::Cint
    sequential_subdivisions::Cint
    init_sequential_subdivisions::Cint
    current_subdivision::Cint
    try_fix_nan::Cint
    gpu_index::Cint
    hierarchy::Ptr{tree}
    input::Ptr{Cfloat}
    truth::Ptr{Cfloat}
    delta::Ptr{Cfloat}
    workspace::Ptr{Cfloat}
    train::Cint
    index::Cint
    cost::Ptr{Cfloat}
    clip::Cfloat
    delta_gpu::Ptr{Cfloat}
    output_gpu::Ptr{Cfloat}
    input_state_gpu::Ptr{Cfloat}
    input_pinned_cpu::Ptr{Cfloat}
    input_pinned_cpu_flag::Cint
    input_gpu::Ptr{Ptr{Cfloat}}
    truth_gpu::Ptr{Ptr{Cfloat}}
    input16_gpu::Ptr{Ptr{Cfloat}}
    output16_gpu::Ptr{Ptr{Cfloat}}
    max_input16_size::Ptr{Cint}
    max_output16_size::Ptr{Cint}
    wait_stream::Cint
    global_delta_gpu::Ptr{Cfloat}
    state_delta_gpu::Ptr{Cfloat}
    max_delta_gpu_size::Cint
    optimized_memory::Cint
    dynamic_minibatch::Cint
    workspace_size_limit::Cint
end

struct network_state
    truth::Ptr{Cfloat}
    input::Ptr{Cfloat}
    delta::Ptr{Cfloat}
    workspace::Ptr{Cfloat}
    train::Cint
    index::Cint
    net::network
end

struct image
    w::Cint
    h::Cint
    c::Cint
    data::Ptr{Cfloat}
end

struct box
    x::Cfloat
    y::Cfloat
    w::Cfloat
    h::Cfloat
end

struct detection
    bbox::box
    classes::Cint
    prob::Ptr{Cfloat}
    mask::Ptr{Cfloat}
    objectness::Cfloat
    sort_class::Cint
    uc::Ptr{Cfloat}
    points::Cint
end

struct matrix
    rows::Cint
    cols::Cint
    vals::Ptr{Ptr{Cfloat}}
end

struct data
    w::Cint
    h::Cint
    X::matrix
    y::matrix
    shallow::Cint
    num_boxes::Ptr{Cint}
    boxes::Ptr{Ptr{box}}
end

@cenum data_type::UInt32 begin
    CLASSIFICATION_DATA = 0
    DETECTION_DATA = 1
    CAPTCHA_DATA = 2
    REGION_DATA = 3
    IMAGE_DATA = 4
    COMPARE_DATA = 5
    WRITING_DATA = 6
    SWAG_DATA = 7
    TAG_DATA = 8
    OLD_CLASSIFICATION_DATA = 9
    STUDY_DATA = 10
    DET_DATA = 11
    SUPER_DATA = 12
    LETTERBOX_DATA = 13
    REGRESSION_DATA = 14
    SEGMENTATION_DATA = 15
    INSTANCE_DATA = 16
    ISEG_DATA = 17
end

struct load_args
    threads::Cint
    paths::Ptr{Cstring}
    path::Cstring
    n::Cint
    m::Cint
    labels::Ptr{Cstring}
    h::Cint
    w::Cint
    c::Cint
    out_w::Cint
    out_h::Cint
    nh::Cint
    nw::Cint
    num_boxes::Cint
    min::Cint
    max::Cint
    size::Cint
    classes::Cint
    background::Cint
    scale::Cint
    center::Cint
    coords::Cint
    mini_batch::Cint
    track::Cint
    augment_speed::Cint
    letter_box::Cint
    show_imgs::Cint
    dontuse_opencv::Cint
    jitter::Cfloat
    resize::Cfloat
    flip::Cint
    gaussian_noise::Cint
    blur::Cint
    mixup::Cint
    label_smooth_eps::Cfloat
    angle::Cfloat
    aspect::Cfloat
    saturation::Cfloat
    exposure::Cfloat
    hue::Cfloat
    d::Ptr{data}
    im::Ptr{image}
    resized::Ptr{image}
    type::data_type
    hierarchy::Ptr{tree}
end

struct metadata
    classes::Cint
    names::Ptr{Cstring}
end

@cenum IMTYPE::UInt32 begin
    PNG = 0
    BMP = 1
    TGA = 2
    JPG = 3
end

@cenum BINARY_ACTIVATION::UInt32 begin
    MULT = 0
    ADD = 1
    SUB = 2
    DIV = 3
end

struct update_args
    batch::Cint
    learning_rate::Cfloat
    momentum::Cfloat
    decay::Cfloat
    adam::Cint
    B1::Cfloat
    B2::Cfloat
    eps::Cfloat
    t::Cint
end

struct boxabs
    left::Cfloat
    right::Cfloat
    top::Cfloat
    bot::Cfloat
end

struct dxrep
    dt::Cfloat
    db::Cfloat
    dl::Cfloat
    dr::Cfloat
end

struct ious
    iou::Cfloat
    giou::Cfloat
    diou::Cfloat
    ciou::Cfloat
    dx_iou::dxrep
    dx_giou::dxrep
end

struct det_num_pair
    num::Cint
    dets::Ptr{detection}
end

const pdet_num_pair = Ptr{det_num_pair}

struct box_label
    id::Cint
    x::Cfloat
    y::Cfloat
    w::Cfloat
    h::Cfloat
    left::Cfloat
    right::Cfloat
    top::Cfloat
    bottom::Cfloat
end


text:
using PkgPage
using Pkg
using Test
import Franklin
using FranklinUtils
using NodeJS
const P = PkgPage
const F = Franklin

m = F.newmodule(F.utils_name())
m.eval(Meta.parse("using PkgPage"))

include("latex/basic.jl")
include("latex/environments.jl")

include("instantiate.jl")

ext: jl


text:
using RecursiveArrayTools: ArrayPartition
using ComponentArrays: ComponentVector
using ..RingPolymers: RingPolymers

get_positions(u::ComponentVector) = u.r
get_velocities(u::ComponentVector) = u.v
get_flat_positions(u::ComponentVector) = @view get_positions(u)[:]
get_flat_velocities(u::ComponentVector) = @view get_velocities(u)[:]

# These must be kept since the symplectic integrators (VelocityVerlet) in DiffEq use ArrayPartitions.
get_positions(u::ArrayPartition) = u.x[2]
get_velocities(u::ArrayPartition) = u.x[1]

function get_quantum_subsystem end

get_mapping_positions(u::ComponentVector) = u.qmap
get_mapping_momenta(u::ComponentVector) = u.pmap
get_mapping_positions(u::ComponentVector, i) = @view get_mapping_positions(u)[:,i]
get_mapping_momenta(u::ComponentVector, i) = @view get_mapping_momenta(u)[:,i]


text:
# DL model for the Prima Indian data

# load packages
using Pkg
Pkg.status()
Pkg.build("Knet")
using Knet
using PyPlot
using LinearAlgebra
using Random
using Printf
import StatsBase.predict

# load help function
include(pwd()*"/utilities.jl")

Random.seed!(1234) # set random numbers

x_train,y_train,x_val,y_val,x_test,y_test = load_data([1,2],0.2)

# define network
n_input = 8
n_hidden_1 = 50
n_hidden_2 = 20
n_out = 2

# predict function

function predict(w,x,p)
    output = zeros(2,size(x,2))
    for i=1:2:length(w)-2
        if i > 1
            x = dropout(x,p)
        end
        x = relu.(w[i]*x .+ w[i+1])
    end
    return softmax(w[end-1]*x .+ w[end];dims=1)
end

loss(w,x,ygold,p, lambda) = Knet.nll(predict(w,x,p), ygold; average=true) + lambda[1]*sum(w[1].^2)  + lambda[2]*sum(w[3].^2)  + lambda[3]*sum(w[5].^2)

lossgradient = grad(loss)

function train(epoch::Int,dropout_percentage::Real,lambda::Vector, w, optim)

    loss_train = zeros(epoch)
    loss_val = zeros(epoch)

    for i in 1:epoch

        idx_train = rand(1:size(x_train,2),size(x_train,2))
        idx_val = rand(1:size(x_val,2),size(x_val,2))

        grads = lossgradient(w,x_train[:,idx_train[:]],y_train[idx_train],dropout_percentage,lambda)
        update!(w, grads, optim)

        loss_train[i] = loss(w,x_train[:,idx_train[:]], y_train[idx_train], 0,lambda)
        loss_val[i] = loss(w,x_val[:,idx_val], y_val[idx_val], 0,lambda)

    end

    return loss_train, loss_val

end

max_stars_repo_name: SamuelWiqvist/dl-compute-project

max_issues_repo_name: SamuelWiqvist/dl-compute-project

max_forks_repo_name: SamuelWiqvist/dl-compute-project


text:
mutable struct NonlinearSolverStatus{XT,YT,AXT,AYT}
    i::Int        # iteration number

    rxₐ::XT       # residual (absolute)
    rxᵣ::XT       # residual (relative)
    rfₐ::YT       # residual (absolute)
    rfᵣ::YT       # residual (relative)

    x::AXT        # initial solution
    x̄::AXT        # previous solution
    δ::AXT        #

    f::AYT        # initial function
    f̄::AYT        # previous function
    γ::AYT        # initial function

    x_converged::Bool
    f_converged::Bool
    g_converged::Bool
    f_increased::Bool

    NonlinearSolverStatus{T}(n) where {T} = new{T,T,Vector{T},Vector{T}}(
        0, 0, 0, 0, 0,
        zeros(T,n), zeros(T,n), zeros(T,n),
        zeros(T,n), zeros(T,n), zeros(T,n),
        false, false, false, false)
end

solution(status::NonlinearSolverStatus) = status.x

function clear!(status::NonlinearSolverStatus{XT,YT}) where {XT,YT}
    status.i = 0

    status.rxₐ = XT(NaN)
    status.rxᵣ = XT(NaN)
    status.rfₐ = YT(NaN)
    status.rfᵣ = YT(NaN)

    status.x̄ .= XT(NaN)
    status.x .= XT(NaN)
    status.δ .= XT(NaN)

    status.f̄ .= YT(NaN)
    status.f .= YT(NaN)
    status.γ .= YT(NaN)

    status.x_converged = false
    status.f_converged = false
    status.f_increased = false
end

Base.show(io::IO, status::NonlinearSolverStatus) = print(io,
                        (@sprintf "    i=%4i" status.i),  ",   ",
                        (@sprintf "rxₐ=%14.8e" status.rxₐ), ",   ",
                        (@sprintf "rxᵣ=%14.8e" status.rxᵣ), ",   ",
                        (@sprintf "rfₐ=%14.8e" status.rfₐ), ",   ",
                        (@sprintf "rfᵣ=%14.8e" status.rfᵣ))

function print_status(status::NonlinearSolverStatus, config::Options)
    if (config.verbosity ≥ 1 && !(assess_convergence!(status, config) && status.i ≤ config.max_iterations)) ||
        config.verbosity > 1
        println(status)
    end
end

increase_iteration_number!(status::NonlinearSolverStatus) = status.i += 1

isconverged(status::NonlinearSolverStatus) = status.x_converged || status.f_converged

function assess_convergence!(status::NonlinearSolverStatus, config::Options)
    x_converged = status.rxₐ ≤ config.x_abstol ||
                  status.rxᵣ ≤ config.x_reltol

    f_converged = status.rfₐ ≤ config.f_abstol ||
                  status.rfᵣ ≤ config.f_reltol

    status.x_converged = x_converged
    status.f_converged = f_converged

    status.f_increased = norm(status.f) > norm(status.f̄)

    return isconverged(status)
end

function meets_stopping_criteria(status::NonlinearSolverStatus, config::Options)
    assess_convergence!(status, config)

    ( isconverged(status) && status.i ≥ config.min_iterations ) ||
    ( status.f_increased && !config.allow_f_increases ) ||
      status.i ≥ config.max_iterations ||
      status.rxₐ > config.x_abstol_break ||
      status.rxᵣ > config.x_reltol_break ||
      status.rfₐ > config.f_abstol_break ||
      status.rfᵣ > config.f_reltol_break ||
      any(isnan, status.x) ||
      any(isnan, status.f)
end

# function check_solver_status(status::NonlinearSolverStatus, config::Options)
#     if any(x -> isnan(x), status.xₚ)
#         throw(NonlinearSolverException("Detected NaN"))
#     end

#     if status.rₐ > config.f_abstol_break
#         throw(NonlinearSolverException("Absolute error ($(status.rₐ)) larger than allowed ($(config.f_abstol_break))"))
#     end

#     if status.rᵣ > config.f_reltol_break
#         throw(NonlinearSolverException("Relative error ($(status.rᵣ)) larger than allowed ($(config.f_reltol_break))"))
#     end
# end

# function get_solver_status!(status::NonlinearSolverStatus, config::Options, status_dict::Dict)
#     status_dict[:nls_niter] = status.i
#     status_dict[:nls_atol] = status.rₐ
#     status_dict[:nls_rtol] = status.rᵣ
#     status_dict[:nls_stol] = status.rₛ
#     status_dict[:nls_converged] = assess_convergence(status, config)
#     return status_dict
# end

function warn_iteration_number(status::NonlinearSolverStatus, config::Options)
    if config.warn_iterations > 0 && status.i ≥ config.warn_iterations
        println("WARNING: Solver took ", status.i, " iterations.")
    end
end

function residual_absolute!(status::NonlinearSolverStatus{T}, y::Vector{T}) where {T}
    @assert length(y) == length(status.y₀)
    status.rₐ = maxnorm(y)
end

function residual_relative!(status::NonlinearSolverStatus{T}, y::Vector{T}) where {T}
    @assert length(y) == length(status.yₚ)
    status.rᵣ = norm(y .- status.yₚ) / norm(status.yₚ)
end

function residual_successive!(status::NonlinearSolverStatus{T}, δx::Vector{T}, x::Vector{T}) where {T}
    @assert length(δx) == length(x)
    status.rₛ = norm(δx) / norm(x)
end

function residual!(status::NonlinearSolverStatus)
    status.rxₐ = norm(status.δ)
    status.rxᵣ = status.rxₐ / norm(status.x)
    status.rfₐ = norm(status.γ)
    status.rfᵣ = status.rfₐ / norm(status.f)
end

function residual!(status::NonlinearSolverStatus, x, f)
    copyto!(status.x, x)
    copyto!(status.f, f)

    status.δ .= status.x - status.x̄
    status.γ .= status.f - status.f̄

    residual!(status)
end

function initialize!(status::NonlinearSolverStatus, x, f)
    clear!(status)
    status.x̄ .= status.x .= x
    status.f̄ .= status.f .= f
    status.δ .= 0
    status.γ .= 0
end

function next_iteration!(status::NonlinearSolverStatus)
    increase_iteration_number!(status)
    status.x̄ .= status.x
    status.f̄ .= status.f
end


text:
push!(LOAD_PATH, "src")

tests = ARGS
if length(ARGS) == 0
    files = readdir("src")

    tests = Array{Integer,1}()
    for file in files
        m = match(r"(\d+)\.jl", file)
        if m != nothing
            append!(tests, parse(Int32, m.captures[1]))
        end
    end
end

for test in tests
    println("Problem ", test, ":")
    include("../src/" * lpad(test, 3, '0') * ".jl")
    println()
end

max_stars_repo_name: gnujosh/julia-euler-project

max_issues_repo_name: gnujosh/julia-euler-project

max_forks_repo_name: gnujosh/julia-euler-project


text:
using LinearAlgebra

include("../datastructs.jl")
include("../kernels/kernel.jl")
include("../bayesopt/incremental.jl")
include("../quadrature/quadrature.jl")
include("../transforms/transforms.jl")

using TimerOutputs
#module buffers0
#export train_buffer, test_buffer, init_train_buffer_dict, update!
if !@isdefined(to)
    const to = TimerOutput()
end

######s
###### Assortment of buffers for storing values in regular BTG workflow
###### Includes initialization and unpacking functions for each buffer type
######

"""
θ-dependent buffer
"""
mutable struct train_buffer
    Σθ::Array{Float64, 2}  #only top n by n block is filled
    Σθ_inv_X::Array{Float64, 2}
    qr_Σθ_inv_X::Union{LinearAlgebra.QRCompactWY{Float64,Array{Float64,2}}, Nothing}
    choleskyΣθ::IncrementalCholesky
    choleskyXΣX::Cholesky{Float64,Array{Float64, 2}}
    logdetΣθ::Float64
    logdetXΣX::Float64
    capacity::Int64 #size of Σθ, maximum value of n
    n::Int64 #number data points incorporated
    k::AbstractCorrelation
    θ::Union{Array{T, 1}, T} where T<:Real
    L_inv_X::Union{Array{Float64, 2}, Nothing}
    function train_buffer(θ::Union{Array{T, 1}, T} where T<:Real, train::AbstractTrainingData, corr::AbstractCorrelation = Gaussian())
        #x = train.x
        #Fx = train.Fx
        #n = train.n
        x = getPosition(train)
        Fx = getCovariates(train)
        n = getNumPts(train)
        capacity = typeof(train)<: extensible_trainingData ? train.capacity : n #if not extensible training type, then set size of buffer to be number of data points
        @timeit to "Σθ" begin
            Σθ = Array{Float64}(undef, capacity, capacity)
            if length(θ)>1
                Σθ[1:n, 1:n] = correlation(corr, θ, x[1:n, :]; jitter = 1e-6) #tell correlation there is single length scale
            else
                #@info "buffers0 51: x[1:n, :]", x[1:n, :]
                Σθ[1:n, 1:n] = correlation(corr, θ[1], x[1:n, :]; jitter = 1e-6) #tell correlation there is single length scale
            end
            choleskyΣθ = incremental_cholesky!(Σθ, n)
        end
        #L = get_chol(choleskyΣθ).L
        #U = get_chol(choleskyΣθ).U
        #L_inv_X = L\Fx
        L_inv_X = nothing #not used for now
        qr_Σθ_inv_X = nothing #we don't use the qr factorization right now
        Σθ_inv_X = choleskyΣθ\Fx
        XΣ_inv_X = Hermitian(Fx'*(Σθ_inv_X))
        @assert issymmetric(XΣ_inv_X)
        @assert isposdef(XΣ_inv_X)
        choleskyXΣX = cholesky(XΣ_inv_X) #regular cholesky because we don't need to extend this factorization
        #catch e
        #    @info "Hermitian(Fx'*(Σθ_inv_X))", Hermitian(Fx'*(Σθ_inv_X))
        #end
        new(Σθ, Σθ_inv_X, qr_Σθ_inv_X, choleskyΣθ, choleskyXΣX, logdet(choleskyΣθ), logdet(choleskyXΣX), capacity, n, corr, θ, L_inv_X)
    end
end

function print(b::train_buffer)
    println("\n ############## TRAIN_BUFFER: ", b.θ)
    p = size(b.Σθ_inv_X, 2)
    println("=============== Size of kernel system ===============")
    display(b.n)
    if n>100
        @warn "Attempting to display large kernel system of size: ", n
    end
    println("=============== Kernel Matrix ===============")
    display(inv(b.choleskyΣθ\Matrix(I, n, n)))
    println("=============== Correlation type ===============")
    display( b.k)
    println("============== Σθ_inv_X =============== ")
    display(b.Σθ_inv_X')
    println(" =============== choleskyXΣX =============== ")
    display(inv(b.choleskyXΣX\Matrix(I, p, p)))
    println("=============== logdetΣθ =============== ")
    display(b.logdetΣθ)
    println("=============== logdetXΣX =============== ")
    display(b.logdetXΣX)
end

"""
(θ, testingData)-dependent quantities, also depends on training data, specifically Σθ
"""
mutable struct test_buffer
    Eθ::Union{Array{Float64, 2}, Nothing}
    Bθ::Union{Array{Float64, 2}, Nothing}
    ΣθinvBθ::Union{Array{Float64, 2}, Nothing}
    Dθ::Union{Array{Float64, 2}, Nothing}
    Hθ::Union{Array{Float64, 2}, Nothing}
    Cθ::Union{Array{Float64, 2}, Nothing}
    θ::Union{Union{Array{T, 1}, T}, Nothing} where T<:Real
    update_bit::Union{Bool, Nothing} #can we update the test buffer without recomputing stuff?
    init_bit::Bool #has the test_buffer been initialized yet?
    test_buffer() = new(nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, false)
end

mutable struct jac_buffer
    jacB::Union{Union{Array{T, 1}, Array{T, 2}}, Nothing} where T<:Float64
    jacD::Union{Union{Array{T, 1}, Array{T, 2}}, Nothing} where T<:Float64
    jacH::Union{Union{Array{T, 1}, Array{T, 2}}, Nothing} where T<:Float64
    jacC::Union{Union{Array{T, 1}, Array{T, 2}}, Nothing} where T<:Float64
    update_bit::Union{Bool, Nothing} #can we update the test buffer without recomputing stuff?
    init_bit::Bool #has the test_buffer been initialized yet?
    jac_buffer() = new(nothing, nothing, nothing, nothing, nothing, false)
end

"""
Store results like Σθ_inv_y, so we have a fast way to do hat(Σθ)_inv_y, where hat(Σθ) is a submatrix
"""
mutable struct θλbuffer
    θ::Union{Array{T, 1}, T} where T<:Real
    λ::T where T<:Real
    βhat::Array{T} where T<:Real
    qtilde::Real
    L_inv_y::Union{Array{T} where T<:Real, Nothing}
    Σθ_inv_y::Array{T} where T<:Real
    remainder::Union{Array{T} where T<:Real, Nothing}
    Σθ_inv_remainder::Array{T} where T<:Real
    θλbuffer() = new()
    function θλbuffer(θ, λ, βhat, qtilde, L_inv_y, Σθ_inv_y, remainder, Σθ_inv_remainder)
        @assert typeof(λ)<:Real
        @assert typeof(θ)<:Union{Array{T}, T} where T<: Real
        if typeof(θ)<:Array{T} where T
            @assert length(θ) > 1
        end
        return new(θ, λ, βhat, qtilde, L_inv_y, Σθ_inv_y, remainder, Σθ_inv_remainder) #Σθ_inv_y - Σθ_inv_X*βhat)
    end
end

function print(b::θλbuffer)
    println("\n ##############  θλ BUFFER: ", (b.θ, b.λ))
    println("=============== βhat ===============")
    display(b.βhat')
    println("=============== qtilde ===============")
    display(b.qtilde)
    println("=============== Σθ_inv_y=============== ")
    display(b.Σθ_inv_y')
end

"""
Stores transformed data
"""
mutable struct λbuffer
    λ::Real
    gλz::Array{T, 1} where T<:Float64
    logjacval::Real where T<:Float64
    dgλz::Array{T, 1} where T<:Float64
    function λbuffer(λ, gλz, logjacval, dgλz)
        return new(λ, gλz, logjacval, dgλz)
    end
end
function print(b::λbuffer)
    println("\n ############## λ BUFFER: ", b.λ)
    println("=============== gλz ===============")
    display(b.gλz')
    println("=============== logjacval ===============")
    display(b.logjacval)
    println("=============== dgλz=============== ")
    display(b.dgλz')
end

"""
update λbuffer with single new y-observation
"""
function update!(λbuffer::λbuffer, y_new::Union{A, B} where A<:Real where B<:Array{H, 1} where H<:Real, g::NonlinearTransform)
    if typeof(y_new)<:Array{T, 1} where T<:Real && length(y_new)==1
        y_new = y_new[1]
    end
    λ = λbuffer.λ;
    #@info "partialx", typeof(partialx)
    dg = x-> partialx(g, x, λ)   #fix lambda
    push!(λbuffer.gλz, g(y_new, λ))
    dgy_new = dg(y_new)
    λbuffer.logjacval = λbuffer.logjacval + log(abs(dgy_new))
    push!(λbuffer.dgλz, dgy_new)
end

"""
Updates θλbuffer
Must have already updated train_buffer, lambda_buffer, and train_data before calling this function
"""
function update!(θλbuffer::θλbuffer, cur_λbuffer::λbuffer, train_buffer::train_buffer, train_data::AbstractTrainingData)
    #(_, Σθ_inv_X, qr_Σθ_inv_X, choleskyΣθ, choleskyXΣX) = unpack(train_buffer)
    #L_inv_X = cur_train_buf.L_inv_X
    (λ, gλz, logjacval) = unpack(cur_λbuffer)
    Σθ_inv_X = train_buffer.Σθ_inv_X
    choleskyXΣX = train_buffer.choleskyXΣX
    choleskyΣθ = train_buffer.choleskyΣθ
    Fx = getCovariates(train_data)
    @timeit to "Σθ_inv_y" Σθ_inv_y = choleskyΣθ\gλz
    #Σθ_inv_y = (choleskyΣθ \ gλz) #O(n^2)
    @timeit to "βhat" βhat = choleskyXΣX\(Fx'*Σθ_inv_y)  #O
    #qtilde = (expr = gλz-Fx*βhat; expr'*(choleskyΣθ\expr))
    @timeit to "qtilde" qtilde =  gλz'*Σθ_inv_y  - 2*gλz'*Σθ_inv_X*βhat + βhat'*Fx'*Σθ_inv_X*βhat #O(np) checks out b/c qtilde = norm(remainder)^2
    #remainder = L_inv_y - L_inv_X*βhat
    θλbuffer.qtilde = qtilde
    θλbuffer.βhat = βhat
    θλbuffer.Σθ_inv_y = Σθ_inv_y
end

"""
Stores Σθ_inv_X_minus_i
"""
mutable struct validation_train_buffer
    θ::Union{Array{T, 1}, T} where T<:Real
    i::Int64
    Σθ_inv_X_minus_i::Array{T} where T<:Real
    logdetΣθ_minus_i::Float64
    logdetXΣX_minus_i::Float64
    choleskyXΣX_X_minus_i::Array{Float64, 2}
    function validation_train_buffer(θ::Union{Array{T, 1}, T} where T<:Float64, i::Int64, train_buf::train_buffer, trainingData::AbstractTrainingData)
        (_, Σθ_inv_X, qr_Σθ_inv_X, choleskyΣθ, _, logdetΣθ, _)= unpack(train_buf) #will need QR factorization and other vals to perform fast LOOCV for least squares
        Σθ_inv_X_minus_i = lin_sys_loocv_IC(Σθ_inv_X, choleskyΣθ, i) #IC stands for incremental_cholesky
        Fx_minus_i = getCovariates(trainingData)[[1:i-1;i+1:end], :]
        sublogdetΣθ = logdet_loocv_IC(choleskyΣθ, logdetΣθ, i)
        sublogdetXΣX = logdet_XΣX_loocv_IC(getCovariates(trainingData), choleskyΣθ, logdetΣθ, i)
        choleskyXΣX_X_minus_i = cholesky(Hermitian(Fx_minus_i'*(Σθ_inv_X_minus_i))) \ (Fx_minus_i')  # for computing beta(-i)
        return new(θ, i, Σθ_inv_X_minus_i, sublogdetΣθ, sublogdetXΣX, choleskyXΣX_X_minus_i)
    end
end

"""
Stores ΣθinvBθ
"""
mutable struct validation_test_buffer
    θ::Union{Array{T, 1}, T} where T<:Real
    ΣθinvBθ::Array{T} where T<:Real
    validation_test_buffer() = new()
    #function validation_test_buffer(θ::Union{Array{T, 1}, T} where T<:Real, train_buf::train_buffer, test_buf::test_buffer)
    #    #ΣθinvBθ = test_buf.ΣθinvBθ
    #    #choleskyΣθ = train_buf.choleskyΣθ
    #    # now these are giving me issues vvv
    #    (_, _, ΣθinvBθ, _, _, _) = unpack(test_buf) #a critical assumption is that the covariates Fx0 remain constant throughout cross-validation
    #    (_, _, _, choleskyΣθ, _, _, _) = unpack(train_buf)
    #    ΣθinvBθ = lin_sys_loocv_IC(ΣθinvBθ, choleskyΣθ, validate) #new ΣθinvBθ of dimension n-1 x 1
    #    return new(θ, ΣθinvBθ)
    #end
end

"""
Stores qtilde_minus_i, βhat_minus_i
"""
mutable struct validation_θλ_buffer
    θ::Union{Array{T, 1}, T} where T<:Real
    λ::Real
    i::Int64
    βhat_minus_i::Array{T} where T<:Real  #depends on theta and lambda
    qtilde_minus_i::Real #depends on theta and lambda
    Σθ_inv_y_minus_i::Array{T} where T<:Real
    # remainder_minus_i
    function validation_θλ_buffer(θ::Union{Array{T, 1}, T} where T<:Float64, λ::Float64, i::Int64, train_buf::train_buffer, λbuf::λbuffer, θλbuf::θλbuffer, val_train_buf::validation_train_buffer, trainingData::AbstractTrainingData)
        # why can't this line below use/find unpack?
        #println(" type of unpacked train_buf in buffers.jl 157:", typeof(unpack(train_buf)))
        #(_, Σθ_inv_X, qr_Σθ_inv_X, choleskyΣθ, _) = unpack(train_buf) #will need QR factorization to perform fast LOOCV for least squares
        Σθ_inv_X =  train_buf.Σθ_inv_X
        # qr_Σθ_inv_X = train_buf.qr_Σθ_inv_X
        choleskyΣθ = train_buf.choleskyΣθ
        # L_inv_X = train_buf.L_inv_X
        #Σθ_inv_X = train_buf.Σθ_inv_X; qr_Σθ_inv_X = train_buf.qr_Σθ_inv_X; choleskyΣθ = train_buf.choleskyΣθ
        #println(" type of unpacked theta_lambda buf in buffers.jl 157:", typeof(unpack(θλbuf)))
        #(_, _, βhat, _, Σθ_inv_y, remainder, Σθ_inv_remainder) = unpack(θλbuf)  #the key here is that βhat, qtilde, etc will already have been computed if we are now doing LOOCV on model
        # βhat = θλbuf.βhat
        Σθ_inv_y = θλbuf.Σθ_inv_y
        # remainder = θλbuf.remainder
        # Σθ_inv_remainder = θλbuf.Σθ_inv_remainder
        Σθ_inv_y_minus_i = lin_sys_loocv_IC(Σθ_inv_y, choleskyΣθ, i)
        Σθ_inv_X_minus_i = lin_sys_loocv_IC(Σθ_inv_X, choleskyΣθ, i)
        gλz_minus_i = λbuf.gλz[[1:i-1;i+1:end]]
        Fx_minus_i = trainingData.Fx[[1:i-1;i+1:end], :]
        # remainder_minus_i, βhat_minus_i = lsq_loocv(L_inv_X, qr_Σθ_inv_X, remainder, βhat, i)
        # qtilde_minus_i = norm(remainder_minus_i[[1:i-1;i+1:end]])^2
        choleskyXΣX_X_minus_i = val_train_buf.choleskyXΣX_X_minus_i
        βhat_minus_i = choleskyXΣX_X_minus_i * Σθ_inv_y_minus_i
        temp = Σθ_inv_y_minus_i - Σθ_inv_X_minus_i*βhat_minus_i
        qtilde_minus_i = ((gλz_minus_i - Fx_minus_i*βhat_minus_i)' * temp)[1]
        #qtilde_minus_i = norm(remainder_minus_i)^2 very blatantly wrong
        #Σθ_inv_remainder_minus_i = lin_sys_loocv_IC(Σθ_inv_remainder, choleskyΣθ, i)
        #qtilde_minus_i = remainder_minus_i[[1:i-1;i+1:end]]'*Σθ_inv_remainder_minus_i
        #@info "Σθ_inv_remainder_minus_i" Σθ_inv_remainder_minus_i
        #@info "remainder_minus_i[[1:i-1;i+1:end]]" remainder_minus_i[[1:i-1;i+1:end]]
        #@info "alternate comp of Σθ_inv_remainder_minus_i" choleskyΣθ\remainder_minus_i
        # @info "qtilde_minus_i" qtilde_minus_i
        @assert typeof(qtilde_minus_i)<:Real
        return new(θ, λ, i, βhat_minus_i, qtilde_minus_i, Σθ_inv_y_minus_i)
    end
end

"""
Stores log jacobian of gλz with ith point deleted. i is not stored because the i is specific
to the call to solve. Rather it's used to construct a validation_λ_buffer.
"""
mutable struct validation_λ_buffer
    logjacval::Real
    function validation_λ_buffer(λbuffer::λbuffer, i::Int64)
        (_, _, logjacval, dgλz) = unpack(λbuffer)
        logjacval = logjacval - log(abs(dgλz[i]))
        return new(logjacval)
    end
end

function print_test_buffer(b::test_buffer)
    println("Eθ: ", b.Eθ)
    println("Bθ: ", b.Bθ)
    println("ΣθinvBθ: ", b.ΣθinvBθ)
    println("Dθ: ", b.Dθ)
    println("Hθ: ", b.Hθ)
    println("Cθ: ", b.Cθ)
end

"""
Checks whether a θλbuffer buffer corresponds to a validation_θλ_buffer
"""
function equals(buf::θλbuffer, val_buf::validation_θλ_buffer)
    try
        @assert buf.θ == val_buf.θ
        @assert buf.λ == val_buf.λ
        @assert isapprox(buf.βhat, val_buf.βhat_minus_i)
        @assert isapprox(buf.qtilde, val_buf.qtilde_minus_i)
    catch e
        return false
    end
    return true
end
"""
Checks if a θλbuffer_dict is equal to a validation_θ_λ_buffer_dict
"""
function equals(θλbuffer_dict::Union{Dict{Tuple{Array{T, 1}, T} where T<:Real , θλbuffer}, Dict{Tuple{T, T} where T<:Real , θλbuffer}}, val_θλbuffer_dict::Union{Dict{Tuple{Array{T, 1}, T} where T<:Real, validation_θλ_buffer}, Dict{Tuple{T, T} where T<:Real, validation_θλ_buffer}})
        @assert keys(θλbuffer_dict) == keys(val_θλbuffer_dict)
        for key in keys(θλbuffer_dict)
            try
                @assert equals(θλbuffer_dict[key], val_θλbuffer_dict[key])
            catch e
                @warn "bad key", key
                return false
            end
        end
    return true
end

###########   what's going on here???? ##############################
function anotherone(b::test_buffer)
    #println("in anotherone")
    return (b.Eθ, b.Bθ, b.ΣθinvBθ, b.Dθ, b.Hθ, b.Cθ)
end
##################### I now deem the unpack functions to be unsafe...it's hard to count underscores...
##################### to get a field, just use dot notation

unpack(b::θλbuffer) = (b.θ, b.λ, b.βhat, b.qtilde, b.L_inv_y, b.Σθ_inv_y, b.remainder, b.Σθ_inv_remainder)
function anotherone(b::θλbuffer)
    return (b.θ, b.λ, b.βhat, b.qtilde, b.L_inv_y, b.Σθ_inv_y, b.remainder, b.Σθ_inv_remainder)
end

function unpack(b::test_buffer)
    #println("in unpack")
    return (b.Eθ, b.Bθ, b.ΣθinvBθ, b.Dθ, b.Hθ, b.Cθ, b.θ)
end
unpack(b::train_buffer) = (b.Σθ, b.Σθ_inv_X, b.qr_Σθ_inv_X, b.choleskyΣθ, b.choleskyXΣX, b.logdetΣθ, b.logdetXΣX, b.θ)
unpack(b::λbuffer) = (b.λ, b.gλz, b.logjacval, b.dgλz)
unpack(b::validation_train_buffer) = (b.θ, b.i, b.Σθ_inv_X_minus_i, b.logdetΣθ_minus_i, b.logdetXΣX_minus_i)
unpack(b::validation_test_buffer) = (b.θ, b.ΣθinvBθ)
unpack(b::validation_θλ_buffer) = (b.θ, b.λ, b.i, b.βhat_minus_i, b.qtilde_minus_i, b.Σθ_inv_y_minus_i) #depends on theta and lambda
unpack(b::validation_λ_buffer) = (b.logjacval)
unpack(b::jac_buffer) = (b.jacC, b.jacB, b.jacH, b.jacD)

#mutable struct validation_buffer_block
#    validation_test_buffer_dict::Dict{Union{Array{T}, T} where T<:Real, validation_test_buffer}
#    validation_train_buffer_dict::Dict{Union{Array{T}, T} where T<:Real, validation_train_buffer}
#    validation_θλ_buffer_dict::Dict{Tuple{Union{Array{T}, T}, T} where T<:Float64, validation_θλ_buffer}
#   function validation_buffer_block()
#         new(Dict{Union{Array{T}, T} where T<:Real, validation_test_buffer}(),
#         {Union{Array{T}, T} where T<:Real, validation_train_buffer}(),
#         Dict{Tuple{Union{Array{T}, T}, T} where T<:Float64, validation_θλ_buffer}())
#    end
#end

"""
Initialize θ to train_buffer dictionary
"""
function init_train_buffer_dict(nw::nodesWeights, trainingData::AbstractTrainingData, corr::AbstractCorrelation = Gaussian(), quadtype::String = "Gaussian")
    if getDimension(nw) == 1 #single theta length scale
        train_buffer_dict = Dict{Float64, train_buffer}()
        for i in 1:size(nw)[2]
            node = nw.nodes[:, i][1]
            push!(train_buffer_dict, node => train_buffer(node, trainingData, corr))
        end
    else
        train_buffer_dict = Dict{Array{Float64, 1}, train_buffer}()
        if quadtype != "Gaussian" # for MC, QMC, SparseGrid
            for i in 1:size(nw)[2]
                node = nw.nodes[:, i]
                push!(train_buffer_dict, node => train_buffer(node, trainingData, corr))
            end
        else
            #println("size nw nodes: ", size(nw.nodes))
            CI = CartesianIndices(Tuple([size(nw)[2] for i = 1:size(nw)[1]]))
            nodeSet = Set(getNodeSequence(nw.nodes, I) for I in CI)
            #println("Iterating over nodeset to build train_buffer_dict...")
            # counter = 1
            for node in nodeSet #this loop is pretty expensive
                #println("Iteration: ", counter); counter += 1
                push!(train_buffer_dict, node => train_buffer(node, trainingData, corr))
            end
        end
    end
    return train_buffer_dict
end

"""
Initialize λbuffer_dict using nodesWeightsλ, trainingData, and nonLinearTransform
"""
function init_λbuffer_dict(nw::nodesWeights, train::AbstractTrainingData, nt:: NonlinearTransform)
    function jac_comp(btg)
        nl2 = nw.num #number of lambda quadrature in each dimension
        Fx = getCovariates(train)
        z = getLabel(train)
        g = (x, λ) -> nt(x, λ); dg = (x, λ) -> partialx(nt, x, λ); lmbda = λ -> g(z, λ); prime =  λ -> dg(z, λ)
        gλvals = Array{Float64, 2}(undef, length(z), nl2) #preallocate space to store gλz arrays
        dgλvals = Array{Float64, 2}(undef, length(z), nl2)
        for i = 1:nl2
            gλvals[:, i] = lmbda(nw.nodes[i])
            dgλvals[:, i] .= prime(nw.nodes[i])
        end
        logjacvals = zeros(1, nl2)  #compute exponents of Jac(z)^(1-p/n)
        for i = 1:nl2
            logjacvals[i] = sum(log.(abs.(map( x-> dg(x, nw.nodes[i]), z))))
        end
        return gλvals, logjacvals, dgλvals
    end
    gλvals, logjacvals, dgλvals = jac_comp(btg)
    # @timeit to "gλvals, logjacvals, dgλvals" gλvals, logjacvals, dgλvals = jac_comp(btg)
    # @timeit to "build λbuffer_dict" begin
    λbuffer_dict = Dict{Real, λbuffer}()
    for i = 1:size(gλvals, 2)
        cur_node = nw.nodes[i]
        @assert typeof(cur_node) <: Real #lambda is always real
        push!(λbuffer_dict, cur_node => λbuffer(cur_node, gλvals[:, i], logjacvals[i], dgλvals[:, i]))
    end
    # end
    return λbuffer_dict
end

"""
Create dictionary of qtilde, βhat, and Σθ_inv_y values.
These are essentially the minimizer and minimum of least squares problem involving the transformed observations, covariance matrix, and covariates
Must be initialized after λbuffer_dict and train_buffer_dict. This function is only called once in the lifetime of the btg object, namely when it is initialized.
"""
function init_θλbuffer_dict(nwθ::nodesWeights, nwλ::nodesWeights, train::AbstractTrainingData, λbuffer_dict::Union{Dict{T, λbuffer}, Dict{Array{T, 1}, λbuffer}} where T<:Real, train_buffer_dict::Union{Dict{Array{T, 1}, train_buffer}, Dict{T, train_buffer}}  where T<: Real, quadtype::Array{String,1})

    #println("type of lambda buffer dict: ", typeof(λbuffer_dict))
    #λbuffer_dict::Dict{Union{Array{T, 1}, T}, λbuffer} where T<:Real
    (R, _) = get_btg_iterator(nwθ, nwλ, quadtype)
    d = getDimension(nwθ)
    if d==1
        θλbuffer_dict = Dict{Tuple{T, T} where T<:Real, θλbuffer}()
    else
        θλbuffer_dict = Dict{Tuple{Array{T, 1}, T} where T<:Real, θλbuffer}()
    end
    for I in R
        (r1, r2, t1, t2) = get_index_slices(nwθ::nodesWeights, nwλ::nodesWeights, quadtype::Array{String,1}, I) #t1, t2 are keys for theta and lambda, respectively
        @assert typeof(t1)<:Union{T, Array{T, 1}} where T<:Real
        @assert typeof(t2)<:Real
        #@info "t1" t1
        if length(t1)==1
            t1 = t1[1]
        end
        cur_train_buf = train_buffer_dict[t1] #get train_buffer
        cur_λ_buf = λbuffer_dict[t2] #get lambda buffer
        θλpair = (t1, t2)::Union{Tuple{Array{T, 1}, T}, Tuple{T, T}} where T<:Real #key used for this buffer
        (_, Σθ_inv_X, qr_Σθ_inv_X, choleskyΣθ, choleskyXΣX) = unpack(cur_train_buf)
        L_inv_X = cur_train_buf.L_inv_X
        (λ, gλz, logjacval) = unpack(cur_λ_buf)
        choleskyXΣX = cur_train_buf.choleskyXΣX
        choleskyΣθ = cur_train_buf.choleskyΣθ
        Fx = getCovariates(train)
        @timeit to "Σθ_inv_y" Σθ_inv_y = choleskyΣθ\gλz
        #Σθ_inv_y = (choleskyΣθ \ gλz) #O(n^2)
        @timeit to "βhat" βhat = choleskyXΣX\(Fx'*Σθ_inv_y)  #O
        #qtilde = (expr = gλz-Fx*βhat; expr'*(choleskyΣθ\expr))
        @timeit to "qtilde" qtilde =  gλz'*Σθ_inv_y  - 2*gλz'*Σθ_inv_X*βhat + βhat'*Fx'*Σθ_inv_X*βhat #O(np) checks out b/c qtilde = norm(remainder)^2
        #remainder = L_inv_y - L_inv_X*βhat
        Σθ_inv_remainder = Σθ_inv_y - Σθ_inv_X*βhat
        cur_θλbuffer = θλbuffer(t1, t2, βhat, qtilde, nothing, Σθ_inv_y, nothing, Σθ_inv_remainder)
        push!(θλbuffer_dict, θλpair => cur_θλbuffer)
    end
        return θλbuffer_dict
end

function init_empty_buffer_dict(nw::nodesWeights, train_buffer_dict::Union{Dict{Array{T, 1}, train_buffer}, Dict{T, train_buffer}} where T<: Real,
    buf_type::Union{Type{jac_buffer}, Type{test_buffer}})
    d = getDimension(nw)
    if d==1 #if single length scale, then key will be Real, else array
        return Dict{Real, buf_type}(arr => buf_type() for arr in keys(train_buffer_dict))
    else
        return Dict{Array{T, 1} where T<:Real, buf_type}(arr => buf_type() for arr in keys(train_buffer_dict))
    end
end

# """
# Initializes test_buffer_dict with empty test buffers, so that the keys match those of train_buffer
# """
# function init_test_buffer_dict(nw::nodesWeights, train_buffer_dict::Union{Dict{Array{T, 1}, train_buffer}, Dict{T, train_buffer}} where T<: Real)
#     d = getDimension(nw)
#     if d==1 #if single length scale, then key will be Real, else array
#         return Dict{Real, test_buffer}(arr => test_buffer() for arr in keys(train_buffer_dict))
#     else
#         return Dict{Array{T, 1} where T<:Real, test_buffer}(arr => test_buffer() for arr in keys(train_buffer_dict))
#     end
# end

function init_validation_train_buffer_dict(nw::nodesWeights)
    d = getDimension(nw)
    if d==1
        return Dict{T where T<: Real, validation_train_buffer}()
    else
        return Dict{Array{T, 1} where T<: Real, validation_train_buffer}()
    end
end

function init_validation_θλ_buffer_dict(nw::nodesWeights)
    d = getDimension(nw)
    if d==1
        return Dict{Tuple{T, T} where T<: Real, validation_θλ_buffer}()
    else
        return Dict{Tuple{Array{T, 1}, T} where T<: Real, validation_θλ_buffer}()
    end
end

"""
Notice that this function is initialized differently from the other validation buffers. It is initialized the same way as
test_buffer. This is because we can't tell what the contents of the buffers will be when the btg object is created, or rather
when we see that we want to do LOOCV. While other buffers are only initialized once, this buffer will be continually updated
with each new data point. Hence we preallocate all the space and then use update! to change the contents of this buffer.
Essentially we are initializing the keys here so we don't have to worry about it in update! If we didn't initialize the keys,
we would have to detect whether they have been initialized or not with each call to update!, or add a field which flags whether
the keys have been initialized.
"""
function init_validation_test_buffer_dict(nw::nodesWeights, test_buffer_dict::Union{Dict{Array{T, 1} where T<: Real, test_buffer}, Dict{T where T<: Real, test_buffer}})
    d = getDimension(nw)
    if d==1
        return Dict{T where T<: Real, validation_test_buffer}(key => validation_test_buffer() for key in keys(test_buffer_dict))
    else
        return Dict{Array{T, 1} where T<: Real, validation_test_buffer}(key => validation_test_buffer() for key in keys(test_buffer_dict))
    end
end

function init_validation_λ_buffer_dict()
    return Dict{T where T<: Real, validation_λ_buffer}()
end

######
###### Update train_buffer when extending kernel system (for Bayesian optimization)
###### Update test_buffer for each new problem, parametrized by triple (x0, Fx0, y0)
######

    #Σθ::Array{Float64, 2}  #only top n by n block is filled
    #Σθ_inv_X::Array{Float64, 2}
    #qr_Σθ_inv_X::LinearAlgebra.QRCompactWY{Float64,Array{Float64,2}}
    #choleskyΣθ::IncrementalCholesky
    #choleskyXΣX::Cholesky{Float64,Array{Float64, 2}}
    #logdetΣθ::Float64
    #logdetXΣX::Float64
    #capacity::Int64 #size of Σθ, maximum value of n
    #n::Int64 #number data points incorporated
    #k::AbstractCorrelation
    #θ::Union{Array{T, 1}, T} where T<:Real
    #L_inv_X::Array{Float64, 2}

"""
Updates training buffer with new training data
-
"""
function update!(train_buffer::train_buffer, trainingData::AbstractTrainingData)  #use incremental cholesky to update training train_buffer

    @assert train_buffer.n < trainingData.n #train_train_buffer must be "older" than trainingData
    k = trainingData.n - train_buffer.n  #typically going to be 1
    @assert k == 1 #we only work with single-point updates for now
    A12, A2 = extend(train_buffer.choleskyΣθ, k) #two view objects

    #peel off "new" values from already updated trainingData
    x_new = getPosition(trainingData)[end:end, :]
    Fx_new = getCovariates(trainingData)[end:end, :]
    y_new = getLabel(trainingData)[end]

    Bθ_new = cross_correlation(Gaussian(), train_buffer.θ, x_new, getPosition(trainingData)[1:end-1, :])
    Eθ_new = [1]

    #A12 = cross_correlation(train_buffer.k(), train_buffer.θ, trainingData.x[1:end-k], trainingData.x[end-k+1:end])
    #A2 = correlation(train_buffer.k(), train_buffer.θ, trainingData.x[end-k+1:end]) #Σθ should get updated automatically, but only upper triangular portion
    A12 .= Bθ_new' #mutates choleskyΣθ
    A2 .= Eθ_new #mutates choleskyΣθ
    update!(train_buffer.choleskyΣθ, k) #extends Cholesky decomposition -- this is where the action happens

    train_Fx = getCovariates(trainingData)
    train_buffer.Σθ_inv_X = train_buffer.choleskyΣθ\train_Fx #potential repeated computations
    train_buffer.choleskyXΣX = cholesky(Hermitian(train_Fx' * train_buffer.Σθ_inv_X))
    train_buffer.logdetΣθ = logdet(train_buffer.choleskyΣθ)
    train_buffer.logdetXΣX = logdet(train_buffer.choleskyXΣX)
    train_buffer.n += k #number of incorporated points
    @assert train_buffer.n == trainingData.n #invariant
    return nothing
end

# function update!(train_buffer::train_buffer, test_buffer::test_buffer, trainingData::AbstractTrainingData)  #use incremental cholesky to update training train_buffer
#     @assert typeof(x0)<:Array{T, 2} where T<:Real
#     @assert typeof(Fx0)<:Array{T, 2} where T<:Real
#     @assert typeof(y0)<:Array{T, 1} where T<:Real
#     @assert train_buffer.n < trainingData.n #train_train_buffer must be "older" than trainingData
#     k = trainingData.n - train_buffer.n  #typically going to be 1
#     @assert k == 1 #we only work with single-point updates for now
#     A12, A2 = extend(train_buffer.choleskyΣθ, k) #two view objects
#     #A12 = cross_correlation(train_buffer.k(), train_buffer.θ, trainingData.x[1:end-k], trainingData.x[end-k+1:end])
#     #A2 = correlation(train_buffer.k(), train_buffer.θ, trainingData.x[end-k+1:end]) #Σθ should get updated automatically, but only upper triangular portion
#     A12 .= test_buffer.Bθ' #mutates choleskyΣθ
#     A2 .= test_buffer.Eθ #mutates choleskyΣθ
#     update!(train_buffer.choleskyΣθ, k) #extends Cholesky decomposition -- this is where the action happens
#     train_buffer.Σθ_inv_X = train_buffer.choleskyΣθ\trainingData.Fx
#     train_buffer.choleskyXΣX = cholesky(trainingData.Fx' * train_buffer.Σθ_inv_X)
#     train_buffer.logdetΣθ = logdet(train_buffer.choleskyΣθ)
#     train_buffer.logdetXΣX = logdet(train_buffer.choleskyXΣX)
#     train_buffer.n .+= k #number of incorporated points
#     return nothing
# end

function refresh_buffer(buffer::Union{test_buffer, jac_buffer})
    buffer.update_bit = true
    return nothing
end

function refresh_buffer_dict(buffer_dict::Union{Union{Dict{Array{T, 1} where T<: Real, test_buffer}, Dict{T where T<: Real, test_buffer}},
    Union{Dict{Array{T, 1} where T<: Real, jac_buffer}, Dict{T where T<: Real, jac_buffer}}})
    for key in keys(buffer_dict)
        refresh_buffer(buffer_dict[key])
    end
    return nothing
end

"""
Update test_buffer, which depends on testing data, training data, and train_buffer.
Note that a test_buffer, attached to a unique theta-value, is only capable of being updated once before being having to be
refreshed again. The logic in this function, which checks init_bit and update_bit, ensure that this invariant holds.
"""
function update!(train_buffer::train_buffer, test_buffer::test_buffer, trainingData::AbstractTrainingData, testingData::AbstractTestingData)
    try
        @assert checkCompatibility(trainingData, testingData) #make sure testingData is compatible with trainingData
    catch e
        @info "train", trainingData
        @info "test", testingData
    end
    function update_me()
        @timeit to "Eθ" test_buffer.Eθ = correlation(train_buffer.k, train_buffer.θ, testingData.x0)
        @timeit to "Bθ" Bθ = cross_correlation(train_buffer.k, train_buffer.θ, testingData.x0, getPosition(trainingData))
        @assert size(Bθ, 2)>= size(Bθ, 1) #row vector, as in the paper
        test_buffer.Bθ = Bθ
        #@info "test_buffer.Bθ'", test_buffer.Bθ'
        @timeit to "ΣθinvBθ" test_buffer.ΣθinvBθ = train_buffer.choleskyΣθ\test_buffer.Bθ'
        @timeit to "Dθ" test_buffer.Dθ = test_buffer.Eθ - test_buffer.Bθ*test_buffer.ΣθinvBθ
        #println("shape of Dtheta: ", size(test_buffer.Dθ ))
        #println("shape of Btheta: ", size(test_buffer.Bθ ))
        #println("shape of Fx0: " , size(testingData.Fx0))
        #println("shape of train_buffer.Σθ_inv_X", size(train_buffer.Σθ_inv_X))
        #@info "buffers0 646: Fx0" testingData.Fx0
        #@info "buffers0 646: Bθ" test_buffer.Bθ
        #@info "buffers0 646: Σθ_inv_X" (train_buffer.Σθ_inv_X)
        @timeit to "Hθ" test_buffer.Hθ = testingData.Fx0 - test_buffer.Bθ*(train_buffer.Σθ_inv_X)
        @timeit to "Cθ" test_buffer.Cθ = test_buffer.Dθ + test_buffer.Hθ*(train_buffer.choleskyXΣX\test_buffer.Hθ')
        test_buffer.θ = train_buffer.θ
        test_buffer.update_bit = false
        test_buffer.init_bit = true
        #@info "D" test_buffer.Dθ
        #@info "H" test_buffer.Hθ
    end
    if test_buffer.init_bit
        if test_buffer.update_bit
            update_me()
        end
    else
        update_me()
    end
    return nothing
end

"""
Update jac_buffer
"""
function update!(jac_buffer::jac_buffer, test_data::testingData, train_data::AbstractTrainingData, test_buffer::test_buffer, train_buffer::train_buffer)
    function update_me()
        x0 = test_data.x0
        θ = test_buffer.θ
        Bθ = test_buffer.Bθ
        Hθ = test_buffer.Hθ
        Σθ_inv_X = train_buffer.Σθ_inv_X
        choleskyXΣX = train_buffer.choleskyXΣX
        ΣθinvBθ = test_buffer.ΣθinvBθ

        x = getPosition(train_data)
        n = getNumPts(train_data)

        @timeit to "derivative intermediates" begin
            #@assert typeof(θ) <: Array{T, 1} where T
            #@assert length(θ) ==d

            S = repeat(x0, n, 1) .- x # n x d
            d = size(x0, 2)
            #println("size of S: ", size(S))
            flattenedθ = typeof(θ)<:Real ? [θ for i=1:d] : θ[:]
            @timeit to "jacB" jacB =   diagm(Bθ[:]) * S * diagm(- flattenedθ) #n x d
            #println("size of jacB: ", size(jacB))
            #println("size of Bθ: ", size(Bθ))
            #println("size of choleskyΣθ: ", size(choleskyΣθ))
            #println("choleskyΣθ inv Bθ': ", size(choleskyΣθ \ Bθ'))
            #jacD = -2* jacB' * (choleskyΣθ \ Bθ') #d x 1
            @timeit to "jacD" jacD = -2* jacB' * ΣθinvBθ
            #println("size of jacD: ", size(jacD))
            #assuming linear polynomial basis
            @timeit to "jacFx0" jacFx0 = vcat(zeros(1, d), diagm(ones(d))) #p x d
            #println("size of Fx0: ", size(Fx0))
            @timeit to "jacH" jacH = jacFx0' - jacB' * Σθ_inv_X #d x p
            #println("size of jacH: ", size(jacH))
            @timeit to "jacC" jacC = jacD + 2 * jacH * (choleskyXΣX \ Hθ') #d x 1
            jac_buffer.jacB = jacB
            jac_buffer.jacD = jacD
            jac_buffer.jacH = jacH
            jac_buffer.jacC = jacC
        end
        jac_buffer.update_bit = false
        jac_buffer.init_bit = true
    end
    if jac_buffer.init_bit
        if jac_buffer.update_bit
            update_me()
        end
    else
        update_me()
    end
    return nothing
end

# """
# update validation_test_buffer
# """
# function update!(validation_test_buffer::validation_test_buffer, train_buf::train_buffer, test_buf::test_buffer, validate)
#     ΣθinvBθ = test_buf.ΣθinvBθ
#     choleskyΣθ = train_buf.choleskyΣθ
#     θ = train_buf.θ
#        # now these are giving me issues vvv
#     (_, _, ΣθinvBθ, _, _, _) = unpack(test_buf) #a critical assumption is that the covariates Fx0 remain constant throughout cross-validation
#        (_, _, _, choleskyΣθ, _, _, _, θ) = unpack(train_buf)
#     ΣθinvBθ = lin_sys_loocv_IC(ΣθinvBθ, choleskyΣθ, validate) #new ΣθinvBθ of dimension n-1 x 1
#     validation_test_buffer.θ = θ
#     validation_test_buffer.ΣθinvBθ = ΣθinvBθ
#     return nothing
# end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
module TestAdd

include("preamble.jl")
using Base.Broadcast: broadcasted

@testset begin
    @test add!!(1, 2) === 3
    @test add!!((1,), (2,)) === (3,)
    @test add!!([1], [2]) ==ₜ [3]
    @test add!!(SVector(1), SVector(2)) === SVector(3)
    @test add!!(SVector(1), [2]) == [3]  # ok if `SVector(3)`
    @test add!!([1], SVector(2)) ==ₜ [3]
    @test add!!([1], [0.5]) ==ₜ Real[1.5]
end

@testset "mutation" begin
    x = [1]
    @test add!!(x, [2]) === x ==ₜ [3]
end

@testset "broadcasted" begin
    @test @inferred(add!!(broadcasted(*, [1], [2]), [2])) ==ₜ [4]
    @test @inferred(add!!([1], broadcasted(/, [1], [2]))) ==ₜ Real[1.5]

    x = [1]
    @test @inferred(add!!(x, broadcasted(*, [1], [2]))) === x ==ₜ [3]
end

end  # module


text:
function add_dual_equality_constraints(
    dual_model::MOI.ModelLike,
    primal_model::MOI.ModelLike,
    primal_dual_map::PrimalDualMap,
    dual_names::DualNames,
    primal_objective::PrimalObjective{T},
    con_types::Vector{Tuple{DataType,DataType}},
    variable_parameters::Vector{VI},
) where {T}
    sense_change =
        MOI.get(dual_model, MOI.ObjectiveSense()) == MOI.MAX_SENSE ? one(T) :
        -one(T)

    all_variables = MOI.get(primal_model, MOI.ListOfVariableIndices())
    restricted_variables = setdiff(all_variables, variable_parameters)

    # Loop at every constraint to get the scalar affine terms
    scalar_affine_terms = get_scalar_affine_terms(
        primal_model,
        primal_dual_map.primal_con_dual_var,
        all_variables,
        con_types,
        T,
    )

    # get RHS from objective coeficients
    scalar_terms = get_scalar_terms(primal_objective)

    # Add terms from objective:
    # Terms from quadratic part
    add_scalar_affine_terms_from_quad_obj(
        scalar_affine_terms,
        primal_dual_map.primal_var_dual_quad_slack,
        primal_objective,
    )
    # terms from mixing variables and parameters
    add_scalar_affine_terms_from_quad_params(
        scalar_affine_terms,
        primal_dual_map.primal_parameter,
        primal_objective,
    )

    # Constrained variables
    for ci in keys(primal_dual_map.constrained_var_dual)
        _add_constrained_variable_constraint(
            dual_model,
            primal_model,
            primal_dual_map.constrained_var_zero,
            primal_dual_map.constrained_var_dual,
            ci,
            scalar_affine_terms,
            scalar_terms,
            sense_change,
            T,
        )
    end

    # Free variables
    for primal_vi in restricted_variables
        if primal_vi in keys(primal_dual_map.constrained_var_idx)
            continue # constrained variable
        end
        # Add equality constraint
        dual_ci = MOIU.normalize_and_add_constraint(
            dual_model,
            MOI.ScalarAffineFunction(scalar_affine_terms[primal_vi], zero(T)),
            MOI.EqualTo(sense_change * get(scalar_terms, primal_vi, zero(T))),
        )
        #Set constraint name with the name of the associated priaml variable
        if !is_empty(dual_names)
            set_dual_constraint_name(
                dual_model,
                primal_model,
                primal_vi,
                dual_ci,
                dual_names.dual_constraint_name_prefix,
            )
        end
        # Add primal variable to dual contraint to the link dictionary
        push!(primal_dual_map.primal_var_dual_con, primal_vi => dual_ci)
    end
    return scalar_affine_terms
end

# `_add_constrained_variable_constraint` is a function barrier to avoid type
# unstability.

function _add_constrained_variable_constraint(
    dual_model,
    primal_model,
    zero_map,
    ci_map,
    ci::MOI.ConstraintIndex{MOI.VectorOfVariables,MOI.Zeros},
    scalar_affine_terms,
    scalar_terms,
    sense_change,
    ::Type{T},
) where {T}
    # The dual is `Reals`, adding a constraint `func`-in-`Reals` is equivalent
    # to not adding any constraint.
    func_primal = MOI.get(primal_model, MOI.ConstraintFunction(), ci)
    return zero_map[ci] = MOIU.vectorize([
        MOI.ScalarAffineFunction(
            MOIU.operate_terms(-, scalar_affine_terms[primal_vi]),
            sense_change * get(scalar_terms, primal_vi, zero(T)),
        ) for primal_vi in func_primal.variables
    ])
end
function _add_constrained_variable_constraint(
    dual_model,
    primal_model,
    zero_map,
    ci_map,
    ci::MOI.ConstraintIndex{MOI.VectorOfVariables},
    scalar_affine_terms,
    scalar_terms,
    sense_change,
    ::Type{T},
) where {T}
    set_primal = MOI.get(primal_model, MOI.ConstraintSet(), ci)
    set_dual = MOI.dual_set(set_primal)
    func_primal = MOI.get(primal_model, MOI.ConstraintFunction(), ci)
    func_dual = MOIU.vectorize([
        MOI.ScalarAffineFunction(
            MOIU.operate_term.(
                *,
                -inv(set_dot(i, set_primal, T)),
                scalar_affine_terms[primal_vi],
            ),
            sense_change *
            inv(set_dot(i, set_primal, T)) *
            get(scalar_terms, primal_vi, zero(T)),
        ) for (i, primal_vi) in enumerate(func_primal.variables)
    ])
    ci_map[ci] = MOI.add_constraint(dual_model, func_dual, set_dual)
    return
end

function _add_constrained_variable_constraint(
    dual_model,
    primal_model,
    zero_map,
    ci_map,
    ci::MOI.ConstraintIndex{MOI.SingleVariable,<:MOI.EqualTo},
    scalar_affine_terms,
    scalar_terms,
    sense_change,
    ::Type{T},
) where {T}
    # Nothing to add as the set is `EqualTo`.
    func_primal = MOI.get(primal_model, MOI.ConstraintFunction(), ci)
    primal_vi = func_primal.variable
    return zero_map[ci] = MOI.ScalarAffineFunction(
        MOIU.operate_terms(-, scalar_affine_terms[primal_vi]),
        sense_change * get(scalar_terms, primal_vi, zero(T)),
    )
end
function _add_constrained_variable_constraint(
    dual_model,
    primal_model,
    zero_map,
    ci_map,
    ci::MOI.ConstraintIndex{MOI.SingleVariable},
    scalar_affine_terms,
    scalar_terms,
    sense_change,
    ::Type{T},
) where {T}
    func_primal = MOI.get(primal_model, MOI.ConstraintFunction(), ci)
    primal_vi = func_primal.variable
    func_dual = MOI.ScalarAffineFunction(
        MOIU.operate_terms(-, scalar_affine_terms[primal_vi]),
        sense_change * get(scalar_terms, primal_vi, zero(T)),
    )
    set_primal = MOI.get(primal_model, MOI.ConstraintSet(), ci)
    set_dual = _dual_set(set_primal)
    ci_map[ci] =
        MOIU.normalize_and_add_constraint(dual_model, func_dual, set_dual)
    return
end

function add_scalar_affine_terms_from_quad_obj(
    scalar_affine_terms::Dict{VI,Vector{MOI.ScalarAffineTerm{T}}},
    primal_var_dual_quad_slack::Dict{VI,VI},
    primal_objective::PrimalObjective{T},
) where {T}
    for term in primal_objective.obj.quadratic_terms
        if term.variable_index_1 == term.variable_index_2
            dual_vi = primal_var_dual_quad_slack[term.variable_index_1]
            push_to_scalar_affine_terms!(
                scalar_affine_terms[term.variable_index_1],
                -MOI.coefficient(term),
                dual_vi,
            )
        else
            dual_vi_1 = primal_var_dual_quad_slack[term.variable_index_1]
            push_to_scalar_affine_terms!(
                scalar_affine_terms[term.variable_index_2],
                -MOI.coefficient(term),
                dual_vi_1,
            )
            dual_vi_2 = primal_var_dual_quad_slack[term.variable_index_2]
            push_to_scalar_affine_terms!(
                scalar_affine_terms[term.variable_index_1],
                -MOI.coefficient(term),
                dual_vi_2,
            )
        end
    end
end

function add_scalar_affine_terms_from_quad_params(
    scalar_affine_terms::Dict{VI,Vector{MOI.ScalarAffineTerm{T}}},
    primal_parameter::Dict{VI,VI},
    primal_objective::PrimalObjective{T},
) where {T}
    for (key, val) in primal_objective.quad_cross_parameters
        for term in val
            dual_vi = primal_parameter[term.variable_index]
            push_to_scalar_affine_terms!(
                scalar_affine_terms[key],
                -MOI.coefficient(term),
                dual_vi,
            )
        end
    end
end

function set_dual_constraint_name(
    dual_model::MOI.ModelLike,
    primal_model::MOI.ModelLike,
    primal_vi::VI,
    dual_ci::CI,
    prefix::String,
)
    MOI.set(
        dual_model,
        MOI.ConstraintName(),
        dual_ci,
        prefix * MOI.get(primal_model, MOI.VariableName(), primal_vi),
    )
    return
end

function get_scalar_terms(primal_objective::PrimalObjective{T}) where {T}
    scalar_terms = Dict{VI,T}()
    for term in get_affine_terms(primal_objective)
        if haskey(scalar_terms, term.variable_index)
            scalar_terms[term.variable_index] += MOI.coefficient(term)
        else
            scalar_terms[term.variable_index] = MOI.coefficient(term)
        end
    end
    return scalar_terms
end

function fill_scalar_affine_terms!(
    scalar_affine_terms::Dict{VI,Vector{MOI.ScalarAffineTerm{T}}},
    primal_con_dual_var::Dict{CI,Vector{VI}},
    primal_model::MOI.ModelLike,
    ::Type{F},
    ::Type{S},
) where {T,F,S}
    for ci in MOI.get(primal_model, MOI.ListOfConstraintIndices{F,S}())
        fill_scalar_affine_terms!(
            scalar_affine_terms,
            primal_con_dual_var,
            primal_model,
            ci,
        )
    end
    return
end

function get_scalar_affine_terms(
    primal_model::MOI.ModelLike,
    primal_con_dual_var::Dict{CI,Vector{VI}},
    variables::Vector{VI},
    con_types::Vector{Tuple{DataType,DataType}},
    ::Type{T},
) where {T}
    scalar_affine_terms = Dict{VI,Vector{MOI.ScalarAffineTerm{T}}}(
        vi => MOI.ScalarAffineTerm{T}[] for vi in variables
    )
    for (F, S) in con_types
        fill_scalar_affine_terms!(
            scalar_affine_terms,
            primal_con_dual_var,
            primal_model,
            F,
            S,
        )
    end
    return scalar_affine_terms
end

function push_to_scalar_affine_terms!(
    scalar_affine_terms::Vector{MOI.ScalarAffineTerm{T}},
    affine_term::T,
    vi::VI,
) where {T}
    if !iszero(affine_term) # if term is different than 0 add to the scalar affine terms vector
        push!(scalar_affine_terms, MOI.ScalarAffineTerm(affine_term, vi))
    end
    return
end

function fill_scalar_affine_terms!(
    scalar_affine_terms::Dict{VI,Vector{MOI.ScalarAffineTerm{T}}},
    primal_con_dual_var::Dict{CI,Vector{VI}},
    primal_model::MOI.ModelLike,
    ci::CI{SAF{T},S},
) where {T,S<:Union{MOI.GreaterThan{T},MOI.LessThan{T},MOI.EqualTo{T}}}
    moi_function = get_function(primal_model, ci)
    for term in moi_function.terms
        dual_vi = primal_con_dual_var[ci][1] # In this case we only have one vi
        push_to_scalar_affine_terms!(
            scalar_affine_terms[term.variable_index],
            MOI.coefficient(term),
            dual_vi,
        )
    end
    return
end

function fill_scalar_affine_terms!(
    scalar_affine_terms::Dict{VI,Vector{MOI.ScalarAffineTerm{T}}},
    primal_con_dual_var::Dict{CI,Vector{VI}},
    primal_model::MOI.ModelLike,
    ci::CI{SVF,S},
) where {T,S<:Union{MOI.GreaterThan{T},MOI.LessThan{T},MOI.EqualTo{T}}}
    dual_var = get(primal_con_dual_var, ci, nothing)
    if dual_var === nothing
        return # No variables created as it's a constrained variable
    end
    moi_function = get_function(primal_model, ci)
    dual_vi = dual_var[1] # In this case we only have one vi
    push_to_scalar_affine_terms!(
        scalar_affine_terms[moi_function.variable],
        one(T),
        dual_vi,
    )
    return
end

function fill_scalar_affine_terms!(
    scalar_affine_terms::Dict{VI,Vector{MOI.ScalarAffineTerm{T}}},
    primal_con_dual_var::Dict{CI,Vector{VI}},
    primal_model::MOI.ModelLike,
    ci::CI{VAF{T},S},
) where {T,S<:MOI.AbstractVectorSet}
    moi_function = get_function(primal_model, ci)
    set = get_set(primal_model, ci)
    for term in moi_function.terms
        dual_vi = primal_con_dual_var[ci][term.output_index]
        # term.output_index is the row of the VAF,
        # it corresponds to the dual variable associated with this constraint
        push_to_scalar_affine_terms!(
            scalar_affine_terms[term.scalar_term.variable_index],
            set_dot(term.output_index, set, T) * MOI.coefficient(term),
            dual_vi,
        )
    end
    return
end

function fill_scalar_affine_terms!(
    scalar_affine_terms::Dict{VI,Vector{MOI.ScalarAffineTerm{T}}},
    primal_con_dual_var::Dict{CI,Vector{VI}},
    primal_model::MOI.ModelLike,
    ci::CI{VVF,S},
) where {T,S<:MOI.AbstractVectorSet}
    dual_vars = get(primal_con_dual_var, ci, nothing)
    if dual_vars === nothing
        return # No variables created as it's part of constrained variables
    end
    moi_function = get_function(primal_model, ci)
    set = get_set(primal_model, ci)
    for (i, variable) in enumerate(moi_function.variables)
        dual_vi = dual_vars[i]
        push_to_scalar_affine_terms!(
            scalar_affine_terms[variable],
            set_dot(i, set, T) * one(T),
            dual_vi,
        )
    end
    return
end

function set_dot(i::Integer, s::MOI.AbstractVectorSet, T::Type)
    vec = zeros(T, MOI.dimension(s))
    vec[i] = one(T)
    return MOIU.set_dot(vec, vec, s)
end
function set_dot(::Integer, ::MOI.AbstractScalarSet, T::Type)
    return one(T)
end


text:
abstract AbstractAffineMap <: Transformation

"""
    Translation(v) <: AbstractAffineMap
    Translation(dx, dy)       (2D)
    Translation(dx, dy, dz)   (3D)

Construct the `Translation` transformation for translating Cartesian points by
an offset `v = (dx, dy, ...)`
"""
immutable Translation{V <: AbstractVector} <: AbstractAffineMap
    v::V
end
Translation(x::Tuple) = Translation(SVector(x))
Translation(x,y) = Translation(SVector(x,y))
Translation(x,y,z) = Translation(SVector(x,y,z))
Base.show(io::IO, trans::Translation) = print(io, "Translation$((trans.v...))")

function (trans::Translation{V}){V}(x)
    x + trans.v
end

Base.inv(trans::Translation) = Translation(-trans.v)

function compose(trans1::Translation, trans2::Translation)
    Translation(trans1.v + trans2.v)
end

transform_deriv(trans::Translation, x) = I
transform_deriv_params(trans::Translation, x) = I

function Base.isapprox(t1::Translation, t2::Translation; kwargs...)
    isapprox(t1.v, t2.v; kwargs...)
end

"""
    LinearMap <: AbstractAffineMap
    LinearMap(M)

A general linear transformation, constructed using `LinearMap(M)`
for any `AbstractMatrix` `M`.
"""
immutable LinearMap{M <: AbstractMatrix} <: AbstractAffineMap
    m::M
end
Base.show(io::IO, trans::LinearMap)   = print(io, "LinearMap($(trans.m))") # TODO make this output more petite

function (trans::LinearMap{M}){M}(x)
    trans.m * x
end

Base.inv(trans::LinearMap) = LinearMap(inv(trans.m))

compose(t1::LinearMap, t2::LinearMap) = LinearMap(t1.m * t2.m)

function Base.isapprox(t1::LinearMap, t2::LinearMap; kwargs...)
    isapprox(t1.m, t2.m; kwargs...)
end

function Base.isapprox(t1::LinearMap, t2::Translation; kwargs...)
    isapprox(vecnorm(t1.m), 0; kwargs...) &&
        isapprox(vecnorm(t2.v),0; kwargs...)
end

function Base.isapprox(t1::Translation, t2::LinearMap; kwargs...)
    isapprox(vecnorm(t1.v), 0; kwargs...) &&
        isapprox(vecnorm(t2.m),0; kwargs...)
end

function Base.:(==)(t1::LinearMap, t2::Translation)
    vecnorm(t1.m) == 0 &&
        0 == vecnorm(t2.v)
end

function Base.:(==)(t1::Translation, t2::LinearMap)
    vecnorm(t1.v) == 0 &&
        vecnorm(t2.m) == 0
end

transform_deriv(trans::LinearMap, x) = trans.m
# TODO transform_deriv_params

"""
    AffineMap <: AbstractAffineMap

A concrete affine transformation.  To construct the mapping `x -> M*x + v`, use

    AffineMap(M, v)

where `M` is a matrix and `v` a vector.  An arbitrary `Transformation` may be
converted into an affine approximation by linearizing about a point `x` using

    AffineMap(trans, [x])

For transformations which are already affine, `x` may be omitted.
"""
immutable AffineMap{M <: AbstractMatrix, V <: AbstractVector} <: AbstractAffineMap
    m::M
    v::V
end

function (trans::AffineMap{M,V}){M,V}(x)
    trans.m * x + trans.v
end

# Note: the expression `Tx - dT*Tx` will have large cancellation error for
# large Tx!  However, changing the order of applying the matrix and
# translation won't fix things, because then we'd have `Tx*(x-x0)` which
# also can incur large cancellation error in `x-x0`.
"""
    AffineMap(trans::Transformation, x0)

Create an Affine transformation corresponding to the differential transformation
of `x0 + dx` according to `trans`, i.e. the Affine transformation that is
locally most accurate in the vicinity of `x0`.
"""
function AffineMap(trans::Transformation, x0)
    dT = transform_deriv(trans, x0)
    Tx = trans(x0)
    AffineMap(dT, Tx - dT*x0)
end

Base.show(io::IO, trans::AffineMap) = print(io, "AffineMap($(trans.m), $(trans.v))") # TODO make this output more petite

function compose(t1::Translation, t2::LinearMap)
    AffineMap(t2.m, t1.v)
end

function compose(t1::LinearMap, t2::Translation)
    AffineMap(t1.m, t1.m * t2.v)
end

function compose(t1::AffineMap, t2::AffineMap)
    AffineMap(t1.m * t2.m, t1.v + t1.m * t2.v)
end

function compose(t1::AffineMap, t2::LinearMap)
    AffineMap(t1.m * t2.m, t1.v)
end

function compose(t1::LinearMap, t2::AffineMap)
    AffineMap(t1.m * t2.m, t1.m * t2.v)
end

function compose(t1::AffineMap, t2::Translation)
    AffineMap(t1.m, t1.v + t1.m * t2.v)
end

function compose(t1::Translation, t2::AffineMap)
    AffineMap(t2.m, t1.v + t2.v)
end

function Base.inv(trans::AffineMap)
    m_inv = inv(trans.m)
    AffineMap(m_inv, m_inv * (-trans.v))
end

function Base.isapprox(t1::AffineMap, t2::AffineMap; kwargs...)
    isapprox(t1.m, t2.m; kwargs...) &&
        isapprox(t1.v, t2.v; kwargs...)
end

function Base.isapprox(t1::AffineMap, t2::Translation; kwargs...)
    isapprox(vecnorm(t1.m), 0; kwargs...) &&
        isapprox(t1.v, t2.v; kwargs...)
end

function Base.isapprox(t1::Translation, t2::AffineMap; kwargs...)
    isapprox(vecnorm(t2.m), 0; kwargs...) &&
        isapprox(t1.v, t2.v; kwargs...)
end

function Base.isapprox(t1::AffineMap, t2::LinearMap; kwargs...)
    isapprox(t1.m, t2.m; kwargs...) &&
        isapprox(vecnorm(t1.v), 0; kwargs...)
end

function Base.isapprox(t1::LinearMap, t2::AffineMap; kwargs...)
    isapprox(t1.m, t2.m; kwargs...) &&
        isapprox(0, vecnorm(t2.v); kwargs...)
end

function Base.:(==)(t1::AffineMap, t2::Translation)
    vecnorm(t1.m) == 0 &&
        t1.v == t2.v
end

function Base.:(==)(t1::Translation, t2::AffineMap)
    vecnorm(t2.m) == 0 &&
        t1.v == t2.v
end

function Base.:(==)(t1::AffineMap, t2::LinearMap)
    t1.m == t2.m &&
        vecnorm(t1.v) == 0
end

function Base.:(==)(t1::LinearMap, t2::AffineMap)
    t1.m == t2.m &&
        0 == vecnorm(t2.v)
end

transform_deriv(trans::AffineMap, x) = trans.m
# TODO transform_deriv_params


text:
using CUDA
using CUDA.CUSPARSE
using ExaPF
using KernelAbstractions
using LinearAlgebra
using Random
using SparseArrays
using Test
using TimerOutputs
using FiniteDiff

import ExaPF: ParsePSSE, PowerSystem, IndexSet

Random.seed!(2713)

@testset "Problem formulations" begin
    # Test PowerNetwork object and parsers
    include("powersystem.jl")
    # Test that behavior matches Matpower
    include("matpower.jl")
    # Test polar formulation
    include("polar_form.jl")
end

@testset "Iterative solvers" begin
    include("iterative_solvers.jl")
end

@testset "Optimization evaluators" begin
    # Resolution of powerflow equations with NLPEvaluators
    include("powerflow.jl")
    # Reduced gradient
    include("reduced_gradient.jl")
    # ReducedSpaceEvaluator API
    include("reduced_evaluator.jl")
    include("penalty.jl")
    include("auglag.jl")
    include("MOI_wrapper.jl")
end

@testset "Reduced space algorithms" begin
    include("test_rgm.jl")
end


text:
using BayesianOptimization, GaussianProcesses
using Test, Random
include("acquisitionfunctions.jl")
include("branin.jl")
include("warmstart.jl")
include("utils.jl")


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    StructuredGridData(data, X, Y, Z, ...)

Data spatially distributed on a structured grid where points are
georeferenced by coordinates `X`, `Y`, `Z`, ...

The `data` argument is a dictionary mapping variable names to Julia
arrays with the actual data.

## Examples

A very popular structured grid data format is NetCDF. Given 2D arrays
`LAT` and `LON` with coordinates and arrays with climate data `precipitation`,
`temperature`, the following code can be used to construct a structured grid:

```julia
julia> data = OrderedDict(:precipitation => precipitation, :temperature => temperature)
julia> StructuredGridData(data, LAT, LON)
```

See also: [`StructuredGrid`](@ref)
"""
struct StructuredGridData{T,N} <: AbstractData{T,N}
  data::OrderedDict{Symbol,<:AbstractArray}
  domain::StructuredGrid{T,N}

  function StructuredGridData{T,N}(data, domain) where {N,T}
    sizes = [size(array) for array in values(data)]
    @assert length(unique(sizes)) == 1 "data dimensions must be the same for all variables"
    @assert prod(sizes[1]) == npoints(domain) "data and coordinates arrays must have the same number of indices"
    new(data, domain)
  end
end

StructuredGridData(data::OrderedDict{Symbol,<:AbstractArray},
                   coordarrays::Vararg{<:AbstractArray{T},N}) where {N,T} =
  StructuredGridData{T,N}(data, StructuredGrid(coordarrays...))

Base.size(geodata::StructuredGridData) = size(geodata.domain)

function Base.getindex(geodata::StructuredGridData, var::Symbol)
  vals = [getindex(geodata, ind, var) for ind in 1:npoints(geodata)]
  reshape(vals, size(geodata))
end

function Base.getindex(geodata::StructuredGridData,
                       icoords::Vararg{Int,N}) where {N}
  vars = [var for (var,V) in variables(geodata)]
  vals = [geodata.data[var][icoords...] for var in vars]
  NamedTuple{tuple(vars...)}(vals)
end

# ------------
# IO methods
# ------------
function Base.show(io::IO, geodata::StructuredGridData{T,N}) where {N,T}
  dims = join(size(geodata), "×")
  print(io, "$dims StructuredGridData{$T,$N}")
end


text:
using BinaryBuilder

# Collection of sources required to build libcalceph
sources = [
    "https://www.imcce.fr/content/medias/recherche/equipes/asd/calceph/calceph-3.3.0.tar.gz" =>
    "be0ee7e838cbd60bb4a90834a1deb4e657d2d83f8dfe7f34774dc41a7c1cb1b4",
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd calceph-3.3.0/
if [[ ${target} == i686-w64* ]] || [[ ${target} == x86_64-w64* ]]; then
    echo 'LT_INIT([win32-dll])' >> configure.ac
    sed -i '/^libcalceph_la_LDFLAGS/ s/$/ -no-undefined/' src/Makefile.am
fi
autoreconf -fi
./configure --prefix=$prefix --host=$target --enable-fortran=no --enable-python=no --disable-static
make -j${nproc}
make install
"""

# Build on all default platforms
platforms = supported_platforms()

# No dependencies
dependencies = []

products = prefix -> [
    LibraryProduct(prefix, "libcalceph", :libcalceph)
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, "libcalceph", v"3.3.0", sources, script, platforms, products, dependencies)


text:
export Mpo,
       add_ops!,
       number_projector

"""
    Mpo{T<:Number}

Matrix product operator with type T.

# Attributes:
- `W::Vector{Array{T, 4}}`: tensor representation of the Mpo.
- `L::Int`: length of the Mpo.
- `d::Int`: physical bond dimension.
"""
struct Mpo{T<:Number}
    W::Vector{Array{T, 4}}
    L::Int
    d::Int
end

Mpo(L::Int, d::Int) = Mpo(Float64, L, d)

"""
    Mpo(::Type{T}, L::Int, d::Int) where T<:Number

Create an empty Mpo of length `L` with bond dimension `d`.
"""
function Mpo(T::Type, L::Int, d::Int)
    Id = Matrix{T}(I, d, d)
    W = Vector{Array{T, 4}}(undef, L)
    W1 = zeros(T, 1, d, d, 2)
    W1[1, :, :, 1] = Id
    W[1] = W1
    Wi = zeros(T, 2, d, d, 2)
    Wi[1, :, :, 1] = Id
    Wi[2, :, :, 2] = Id
    for i=2:L-1
        W[i] = deepcopy(Wi)
    end
    Wend = zeros(T, 2, d, d, 1)
    Wend[2, :, :, 1] = Id
    W[end] = Wend
    return Mpo{T}(W, L, d)
end

"""
    add_ops!(Op::Mpo{T}, op_i::AbstractMatrix{<:Number},
             weights::Vector{T}) where T<:Number

Add local/on-site operators to an Mpo.

New operator is: Op + ∑ weights[i]*op_i.
"""
function add_ops!(Op::Mpo{T}, op_i::AbstractMatrix{<:Number},
                  weights::Vector{T}) where T<:Number
    # Convert op_i to have the same type as Op.
    c_op_i = convert.(T, op_i)

    for i=1:Op.L-1
        Op.W[i][1, :, :, 2] .+= weights[i]*c_op_i
    end
    Op.W[end][1, :, :, 1] .+= weights[end]*c_op_i
    return Op
end

"""
    add_ops!(Op::Mpo{T}, str_op_i::String,
             weights::Vector{T}) where T<:Number

Add local/on-site operators to an Mpo.

New operator is: Op + ∑ weights[i]*op_i.
"""
function add_ops!(Op::Mpo{T}, str_op_i::String,
                  weights::Vector{T}) where T<:Number
    return add_ops!(Op, str_to_op(str_op_i), weights)
end

"""
    add_ops!(Op::Mpo{T}, op_i::AbstractMatrix{<:Number},
             op_j::AbstractMatrix{<:Number},
             weights::AbstractMatrix{T};
             ferm_op::AbstractMatrix{<:Number}=Matrix{Float64}(I, Op.d, Op.d)
             ) where T<:Number

Add operators to an Mpo acting on two sites.

New operator is: Op + ∑ weights[i, j]*op_i*op_j. `ferm_op` is the fermion parity
operator: 1-2n that goes between the fermionic creation and annihilation
operators. If the operators `op_i`, `op_j` are not fermionic, `ferm_op` defaults
to the identity matrix.
"""
function add_ops!(Op::Mpo{T}, op_i::AbstractMatrix{<:Number},
                  op_j::AbstractMatrix{<:Number},
                  weights::AbstractMatrix{T};
                  ferm_op::AbstractMatrix{<:Number}=Matrix{Float64}(I, Op.d, Op.d)
                  ) where T<:Number
    # Extract on-site operators.
    onsite_weights = diag(weights)
    add_ops!(Op, op_j*op_i, onsite_weights)

    # Convert op_i, op_j to have the same type as Op.
    c_op_i = convert.(T, op_i)
    c_op_j = convert.(T, op_j)

    # Current bond dimensions of the Mpo.
    max_w = maximum(size.(Op.W, 4))
    L = Op.L

    # Allocate space for the new operators. Let the operators at i=1 and i=L
    # have dimensions (1, d, d, max_w) and (max_w, d, d, 1).
    for i=1:L
        Wi = Op.W[i]
        tmp_wi = zeros(eltype(Wi), i > 1 ? max_w + L : 1,
                                   size(Wi, 2),
                                   size(Wi, 3),
                                   i < L ? max_w + L : 1)
        tmp_wi[1:size(Wi, 1), 1:size(Wi, 2), 1:size(Wi, 3), 1:size(Wi, 4)] += Wi
        Op.W[i] = tmp_wi
    end

    # Write the new operators. Each iteration corresponds to writing all
    # operators in weights[i, :]*op_i*op_j.
    for i=1:L
        # Continue if there are no operators to write.
        all(abs.(weights[i, :]) .< 1e-8) && continue

        # Write op_i.
        if i > 1
            Op.W[i][max_w+i, :, :, i < L ? 2 : 1] = op_i
        end
        if i < L
            Op.W[i][1, :, :, max_w+i] = op_i
        end

        # Write J[i, j]*op_j, j != i.
        for j=1:i-1
            Op.W[j][1, :, :, max_w+i] = op_j*weights[i, j]
        end
        for j=i+1:L-1
            Op.W[j][max_w+i, :, :, 2] = op_j*weights[i, j]
        end
        if i < L
            Op.W[L][max_w+i, :, :, 1] += op_j*weights[i, L]
        end

        # Write operators between op_i and op_j.
        if any(abs.(weights[i, 1:i-1]) .> 1e-8)
            for j=findfirst(abs.(weights[i, :]) .> 1e-8)+1:i-1
                Op.W[j][max_w+i, :, :, max_w+i] = ferm_op
            end
        end
        if any(abs.(weights[i, i+1:L]) .> 1e-8)
            for j=i+1:findlast(abs.(weights[i, :]) .> 1e-8)-1
                Op.W[j][max_w+i, :, :, max_w+i] = ferm_op
            end
        end
    end
    compress_mpo!(Op)
    return Op
end

"""
    add_ops!(Op::Mpo{T}, op_i::String, op_j::String,
             weights::AbstractMatrix{T};
             ferm_op::String="Id") where T<:Number

Add operators to an Mpo acting on two sites.

New operator is: Op + ∑ weights[i, j]*op_i*op_j.
"""
function add_ops!(Op::Mpo{T}, op_i::String, op_j::String,
                  weights::AbstractMatrix{T};
                  ferm_op::String="Id") where T<:Number
    return add_ops!(Op, str_to_op(op_i), str_to_op(op_j), weights,
                    ferm_op=str_to_op(ferm_op, Op.d))
end

"""
    compress_mpo!(Op::Mpo{<:Number})

Remove unnecessary bond dimensions in an Mpo.
"""
function compress_mpo!(Op::Mpo{<:Number})
    # Compress the operators of the Mpo to remove unnecessary zeros. We go over
    # each pair of contiguous tensors, e.g. A and B, with A left of B, and
    # remove the column A[:, :, :, i] and the row B[i, :, :, :] if both the
    # row and the column are empty (have zero norm).
    for i=1:Op.L-1
        A = Op.W[i]
        B = Op.W[i+1]
        to_rm = Int[]
        for j=1:size(A, 4)
            if norm(A[:, :, :, j]) < 1e-8 && norm(B[j, :, :, :]) < 1e-8
                push!(to_rm, j)
            end
        end
        # Make a copy with removed indices and columns.
        cA = zeros(
            eltype(A),
            size(A, 1), size(A, 2), size(A, 3), size(A, 4)-length(to_rm)
        )
        cB = zeros(
            eltype(B),
            size(B, 1)-length(to_rm), size(B, 2), size(B, 3), size(B, 4)
        )
        cont = 1
        for j=1:size(A, 4)
            j ∈ to_rm && continue
            cA[:, :, :, cont] = A[:, :, :, j]
            cB[cont, :, :, :] = B[j, :, :, :]
            cont += 1
        end
        # Overwrite the original tensors.
        Op.W[i] = cA
        Op.W[i+1] = cB
    end
    return Op
end

"""
    str_to_op(str_op::String)

Return the matrix that corresponds to an operator input as a string. Example:
"n" -> [[0. 0.];
        [0. 1.]].
"""
function str_to_op(str_op::String, d::Int=0)
    # Operators for 2 physical dimensions: fermions and hard-core bosons.
    if str_op == "n"
        # Number operator.
        return [[0. 0.];
                [0. 1.]]
    elseif str_op == "sz"
        # Spin 1/2 sz operator.
        return [[-0.5 0.];
                [0. 0.5]]
    elseif str_op == "a+" || str_op == "b+" || str_op == "c+" || str_op == "s+"
        # Fermion/hard-core boson creation or spin 1/2 raising operator.
        return [[0. 1.];
                [0. 0.]]
    elseif str_op == "a" || str_op == "b" || str_op == "c" || str_op == "s-"
        # Fermion/hard-core boson annihilation or spin 1/2 lowering operator.
        return [[0. 0.];
                [1. 0.]]
    elseif str_op == "Z"
        # Parity sign operator.
        return [[1. 0.];
                [0. -1.]]
    # Operators for 3 physical dimensions: spin 1 particles. State 1 is |1, -1>,
    # state 2 is |1, 0> and state 3 is |1, +1>, with |s, m>.
    elseif str_op == "Sz"
        return [[-1. 0. 0.];
                [0. 0. 0.];
                [0. 0. 1.]]
    elseif str_op == "S+"
        return [[0. sqrt(2) 0.];
                [0. 0. sqrt(2)];
                [0. 0. 0.]]
    elseif str_op == "S-"
        return [[0. 0. 0.];
                [sqrt(2) 0. 0.];
                [0. sqrt(2) 0.]]
    # General identity matrix with the same physical dimensions as the Mpo.
    elseif str_op == "Id"
        return Matrix{Float64}(I, d, d)
    else
        throw("Operator $str_op is not defined.")
    end

end

"""
    number_projector(::Type{T}, L::Int, N::Int) where T<:Number

Build a projector over the number of particles N for MPS with d=2.

The projector P consists on the sum of every combination
n_{i1} * n_{i2} * ... * n_{iN} * (1-n_{iN+1}) * ... * (1-n_{iL})
This operator clearly satisfies P*P = P and projects any state into the
space of states with N particles and L-N holes.
"""
function number_projector(T::Type, L::Int, N::Int)
    d = 2

    # Particle number operator, equivalent to a^dagger*a = n.
    np = zeros(T, d, d)
    np[2, 2] = one(T)
    # Hole number operator, equivalent to a*a^dagger = 1-n.
    nh = zeros(T, d, d)
    nh[1, 1] = one(T)

    W = Vector{Array{T, 4}}(undef, L)
    for i=1:N
        nstates_left = min(i, L-N+1)
        nstates_right = min(i+1, L-N+1)
        Wi = zeros(T, nstates_left, d, d, nstates_right)
        for j=1:nstates_right-1
            Wi[j, :, :, j] = np
            Wi[j, :, :, j+1] = nh
        end
        W[i] = Wi
    end
    for i=N+1:L
        nstates_left = min(L-i+2, N+1)
        nstates_right = min(L-i+1, N+1)
        Wi = zeros(T, nstates_left, d, d, nstates_right)
        for j=1:nstates_right
            Wi[j, :, :, j] = nh
            if nstates_left > nstates_right
                Wi[j+1, :, :, j] = np
            end
        end
        W[i] = Wi
    end

    return Mpo{T}(W, L, d)
end

import Base.display

function Base.display(O::Mpo{T}) where T<:Number
    println("MPO:")
    println("   Type: $(eltype(O.W[1]))")
    println("   Length: $(O.L)")
    println("   Physical dims: $(O.d)")
    println("   Max bond dim: $(maximum(size.(O.W, 4)))")
    return
end

function show_bond_dims(O::Mpo{T}) where T<:Number
    bdims = vcat(1, size.(O.W[:], 4))
    println(join(bdims, "-"))
    return
end


text:
"""
The `ZChop` module provides the functions `zchop` and `zchop!`, which replace
numbers that are close to zero with zero. These functions act recursivley on
collections.
"""
module ZChop

export zchop, zchop!

"""
    zeps::Float64

The default lower threshold for a number to be replaced by zero.
"""
const zeps = 1e-14

"""
    zchop(x::T, eps::Real = zeps)

Replace `x` by zero if `abs(x) < eps`. `zchop` acts recursively on
mappable objects. Objects that cannot be sensibly compared to a real
number are passed unaltered.
"""
zchop(x::Real, eps::Real = zeps) = abs(x) > eps ? x : zero(typeof(x))
zchop(x::Complex, eps::Real = zeps) = complex(zchop(real(x), eps), zchop(imag(x), eps))
# Following is not type stable
zchop(x::Irrational, eps::Real = zeps) = zchop(float(x), eps)
# Do not iterate over strings
zchop(x::Union{AbstractString,AbstractChar}, eps::Real = zeps) = x

"""
    zchop!(x::T, eps::Real = zeps)

Perform `zchop` in place.
"""
function zchop!(a::AbstractArray, eps::Real = zeps)
    @inbounds for i in firstindex(a):lastindex(a)
        a[i] = zchop(a[i], eps)
    end
    return a
end

zchop(a::AbstractArray, eps::Real = zeps) = zchop!(copy(a), eps)
zchop(a::Base.Generator, eps::Real = zeps) = zchop!(collect(a), eps)
zchop(a::Tuple, eps::Real = zeps) = zchop.(a, eps)
zchop(x::Expr, eps::Real = zeps) = Expr(x.head, zchop(x.args)...)
zchop(x, eps::Real) = Base.isiterable(typeof(x)) ? map((x)->zchop(x, eps), x) : x
zchop(x) = Base.isiterable(typeof(x)) ? map(zchop, x) : x
zchop(x::Number) = x

end # module ZChop


text:
import HypothesisTests: OneSampleTTest

jtest = JointDistanceDistributionTTest()

@test causality(x, y, jtest) isa OneSampleTTest