text:
using LinearAlgebra

@testset "Test QXZoo gates" begin

    @testset "Test fundamental gate structs" begin
        @testset "GateSymbol" begin
            gs1 = QXZoo.GateSymbol(:label)
            gs2 = QXZoo.GateSymbol(:label, 2)
            gs3 = QXZoo.GateSymbol(:label, 1, true)
            gs4 = sqrt(gs1)
            gs5 = adjoint(gs1)
            @test begin
                gs1.label == :label
                gs1.rt_depth == 0
                gs1.is_adj == false
            end
            @test begin
                gs2.label == :label
                gs2.rt_depth == 2
                gs2.is_adj == false
            end
            @test begin
                gs3.label == :label
                gs3.rt_depth == 1
                gs3.is_adj == true
            end
            @test begin
                gs4.label == gs1.label
                gs4.rt_depth == 1
                gs4.is_adj == false
            end
            @test begin
                gs5.label == gs1.label
                gs5.rt_depth == 0
                gs5.is_adj == true
            end
        end
        @testset "GateSymbolP" begin
            gs1 = QXZoo.GateSymbolP(:label)
            gs2 = QXZoo.GateSymbolP(:label, 2)
            gs3 = QXZoo.GateSymbolP(:label, 1, true)
            gs4 = QXZoo.GateSymbolP(:label, 3, true, pi/7)
            @test begin
                gs1.label == :label
                gs1.rt_depth == 0
                gs1.is_adj == false
                gs1.param === 0
            end
            @test begin
                gs2.label == :label
                gs2.rt_depth == 2
                gs2.is_adj == false
                gs2.param === 0
            end
            @test begin
                gs3.label == :label
                gs3.rt_depth == 1
                gs3.is_adj == true
                gs3.param === 0
            end
            @test begin
                gs4.label == :label
                gs4.rt_depth == 1
                gs4.is_adj == true
                gs4.param === pi/7
            end
        end
        @testset "GateCall1" begin
            gs1 = QXZoo.GateSymbol(:label)
            gc1 = QXZoo.GateCall1(gs1, 3)

            @test begin
                gc1.gate_symbol == gs1
                gc1.target == 3
            end
        end
        @testset "GateCall2" begin
            gs1 = QXZoo.GateSymbol(:c_x, 0, false)
            gs2 = QXZoo.GateSymbolP(:asdf, 0, false, pi/5)
            gc1 = QXZoo.GateCall2(gs1, 2, 3, QXZoo.GateSymbol(:x))
            gc2 = QXZoo.GateCall2(gs2, 1, 2 )

            @test begin
                gc1.gate_symbol == gs1
                gc1.ctrl == 2
                gc1.target == 3
                gc1.base_gate == QXZoo.DefaultGates.GateSymbols.x
                gc2.gate_symbol == gs2
                gc2.ctrl == 1
                gc2.target == 2
                gc2.base_gate === nothing
            end
        end
    end

    @testset "Test default gate structs" begin
        @testset "Default GateSymbol structs" begin
            @test begin
                DefaultGates.GateSymbols.p00 === GateOps.GateSymbol(:p00)
                DefaultGates.GateSymbols.p01 === GateOps.GateSymbol(:p01)
                QXZoo.DefaultGates.GateSymbols.p10 === GateOps.GateSymbol(:p10)
                QXZoo.DefaultGates.GateSymbols.p11 === GateOps.GateSymbol(:p11)
                QXZoo.DefaultGates.GateSymbols.x === GateOps.GateSymbol(:x)
                QXZoo.DefaultGates.GateSymbols.y === GateOps.GateSymbol(:y)
                QXZoo.DefaultGates.GateSymbols.z === GateOps.GateSymbol(:z)

                QXZoo.DefaultGates.GateSymbols.s === GateOps.GateSymbol(:z, 1, false)
                QXZoo.DefaultGates.GateSymbols.t === GateOps.GateSymbol(:z, 2, false)
                QXZoo.DefaultGates.GateSymbols.h === GateOps.GateSymbol(:h)
                QXZoo.DefaultGates.GateSymbols.I === GateOps.GateSymbol(:I)
            end
        end
    end

    @testset "Test single-qubit gate values" begin
        @testset "Test static gates" begin
            @test begin # Individual tests
                QXZoo.GateMap.x() == [0 1; 1 0]
                QXZoo.GateMap.y() == [0 -1im; 1im 0]
                QXZoo.GateMap.z() == [1 0; 0 -1]
                QXZoo.GateMap.h() == [1 1; 1 -1].*(1/sqrt(2))
                QXZoo.GateMap.I() == [1 0; 0 1]

                QXZoo.GateMap.s() == sqrt(QXZoo.GateMap.z())
                QXZoo.GateMap.t() == sqrt(QXZoo.GateMap.s())

                QXZoo.GateMap.p00() == [1 0; 0 0]
                QXZoo.GateMap.p01() == [0 1; 0 0]
                QXZoo.GateMap.p10() == [0 0; 1 0]
                QXZoo.GateMap.p11() == [0 0; 0 1]
            end

            @test begin # identity tests
                QXZoo.GateMap.z() ≈ QXZoo.GateMap.h() * QXZoo.GateMap.x() * QXZoo.GateMap.h()
                -QXZoo.GateMap.y() ≈ QXZoo.GateMap.h() * QXZoo.GateMap.y() * QXZoo.GateMap.h()
                QXZoo.GateMap.x() ≈ QXZoo.GateMap.h() * QXZoo.GateMap.z() * QXZoo.GateMap.h()
            end

            @test begin #
                QXZoo.GateMap.r_x(0) == QXZoo.GateMap.I()
                QXZoo.GateMap.r_y(0) == QXZoo.GateMap.I()
                QXZoo.GateMap.r_z(0) == QXZoo.GateMap.I()
                QXZoo.GateMap.r_phase(0) == QXZoo.GateMap.I()
            end
        end

        @testset "Test rotation gates" begin
            for θ in range(0, 2*pi, step=pi/5)
                @test QXZoo.GateMap.r_x(θ) ≈ cos(θ/2)*[1 0; 0 1] - 1im*sin(θ/2)*QXZoo.GateMap.x()
                @test QXZoo.GateMap.r_y(θ) ≈ cos(θ/2)*[1 0; 0 1] - 1im*sin(θ/2)*QXZoo.GateMap.y()
                @test QXZoo.GateMap.r_z(θ) ≈ cos(θ/2)*[1 0; 0 1] - 1im*sin(θ/2)*QXZoo.GateMap.z()
            end

            for θ in range(0, 2*pi, step=pi/5)
                @test QXZoo.GateMap.r_z(θ) ≈ QXZoo.GateMap.r_y(-pi/2) * QXZoo.GateMap.r_x(θ) * QXZoo.GateMap.r_y(pi/2)
            end
        end
    end

    @testset "Test two-qubit gate values" begin

        @test begin # Individual tests
            QXZoo.GateMap.c_x() == [1 0 0 0; 0 1 0 0; 0 0 0 1; 0 0 1 0]
            QXZoo.GateMap.c_y() == [1 0 0 0; 0 1 0 0; 0 0 0 -1im; 0 0 1im 0]
            QXZoo.GateMap.c_z() == LinearAlgebra.Diagonal([1,1,1,-1])
        end

        @testset "Test rotation gates" begin
            for θ in range(0, 2*pi, step=pi/5)
                @test (QXZoo.GateMap.c_r_x(θ) - (kron(QXZoo.GateMap.p00(), QXZoo.GateMap.I()) + kron(QXZoo.GateMap.p11(), QXZoo.GateMap.r_x(θ)))) ≈ zeros(4,4)

                @test  (QXZoo.GateMap.c_r_y(θ) - (kron(QXZoo.GateMap.p00(), QXZoo.GateMap.I()) + kron(QXZoo.GateMap.p11(), QXZoo.GateMap.r_y(θ)))) ≈ zeros(4,4)

                @test  (QXZoo.GateMap.c_r_z(θ) - (kron(QXZoo.GateMap.p00(), QXZoo.GateMap.I()) + kron(QXZoo.GateMap.p11(), QXZoo.GateMap.r_z(θ)))) ≈ zeros(4,4)

                @test (QXZoo.GateMap.c_r_phase(θ) - (kron(QXZoo.GateMap.p00(), QXZoo.GateMap.I()) + kron(QXZoo.GateMap.p11(), QXZoo.GateMap.r_phase(θ)))) ≈ zeros(4,4)
            end
        end
    end

    @testset "Test custom gates" begin

        @test begin # Individual tests
            gate_func_1q = () -> [0 -1; 1 0].+0im
            q1 = GateMap.create_gate_1q("test_1q", gate_func_1q)
            GateMap.gates[q1]() == gate_func_1q()
            GateMap.gates[GateSymbol(:test_1q, 0, false)]() == gate_func_1q()
        end

        @test begin # Individual tests
            gate_func_1q_p = (a) -> [cos(a) -1; 1 sin(a)].+0im
            q1 = GateMap.create_gate_1q("test_1q_p", gate_func_1q_p)
            GateMap.gates[q1](pi/2) == gate_func_1q_p(pi/2)
        end

    end

end


text:
include("common.jl")
include("compression.jl")

ext: jl


text:
"""
    sum_of_bit_difference(array::Array{Int64,1})

For finding the sum of bit differences in all pairs of an array, the array will be
transformed into a series of subarrays  for a given bit-length.

# Arguments
- `array::Array{Int64,1}`: Unsorted array

# Examples
```julia-repl
julia> import ClassicAlgorithmsCollections
julia> arr = [8, 1, 2, 12, 45, 100, 43]
julia> ClassicAlgorithmsCollections.maximum_xor_subarray(arr)
104
```
"""
function sum_of_bit_difference(array::Array{Int64,1}, bit = 32)
    n = length(array)
    result = 0
    for i in 1:bit
        count = 0
        for j in 1:n
            if (array[j] & (1 << i)) != 0
                count += 1
            end
        end
        result += count * (n - count) * 2
    end
    return result
end

max_stars_repo_name: Anselmoo/ClassicAlgorithmsCollections

max_issues_repo_name: Anselmoo/ClassicAlgorithmsCollections

max_forks_repo_name: Anselmoo/ClassicAlgorithmsCollections


text:
using Test

using Distributed
using Random
using SymbolicUtils

addprocs(2)

@testset "Basic test" begin
    @everywhere include("src/SymbolicRegression.jl")
    @everywhere using .SymbolicRegression

    Random.seed!(0)

    X = randn(Float32, 100, 2)
    y = 2 * cos.(X[:, 2]) + X[:, 1] .^ 2 .- 2

    options = SymbolicRegression.Options(
        binary_operators=(+, *),
        unary_operators=(cos,),
        npopulations=3,
        maxsize=15
    )

    niterations = 5
    hallOfFame = RunSR(X, y, niterations, options)
    dominating = calculateParetoFrontier(X, y, hallOfFame, options)

    @test dominating[end].score < 1e-3
    @test length(stringTree(dominating[end].tree, options)) > 10
end

rmprocs()


text:
"""
    Block

Data Structure representing a Block in the Bitcoin blockchain.

Consists of a `block.header::Header` and
`block.transactions::Vector{Transaction}`.

To get the hash of a `Block`
```julia
double_sha256(block)
```
"""
struct Block
    size                :: UInt32
    header              :: Header
    transaction_counter :: UInt64
    transactions        :: Vector{Transaction}
end

function Block(io::IO)

    block_size = read(io, UInt32)

    block_header = Header(io)

    n_trans = read_varint(io)
    @assert n_trans > zero(n_trans)
    transactions = Transaction[
        Transaction(io)
        for i in 1:n_trans
    ]

    return Block(
        block_size,
        block_header,
        n_trans,
        transactions
    )
end

function Block(bcio::BCIterator)

    check_magic_bytes(bcio)
    block = Block(bcio.io)

    if eof(bcio)
        open_next_file(bcio)
    end

    return block
end

function Block(x::Array{UInt8})
    block_size = length(x)

    io = IOBuffer(x)

    block_header = Header(io)

    n_trans = read_varint(io)
    @assert n_trans > zero(n_trans)
    transactions = Transaction[
        Transaction(io)
        for i in 1:n_trans
    ]

    return Block(
        block_size,
        block_header,
        n_trans,
        transactions
    )
end

function showcompact(io::IO, block::Block)
    @printf(io, "Block, %d bytes, %d transactions:\n",
            block.size, block.transaction_counter)
end

function Base.show(io::IO, block::Block)
    showcompact(io, block)
    if !get(io, :compact, false)
        show(io, block.header)
    end
end

# function Base.showall(io::IO, block::Block)
#     show(io, block)
#     println("Transactions:")
#     show(io, block.transactions)
# end

"""
    double_sha256(x::Block)::UInt256
    double_sha256(x::Header)::UInt256
    double_sha256(x::Transaction)::UInt256

Hash a `Block`, `Header`, or `Transaction`

```julia
double_sha256(block)
double_sha256(header)
double_sha256(transaction)
```
"""
double_sha256(x::Block) = double_sha256(x.header)

function dump_block_data(io::IO)
    block_size = read(io, UInt32)
    read!(io, Array{UInt8}(undef, block_size))
end

function dump_block_data(bcio::BCIterator)
    check_magic_bytes(bcio)
    dump_block_data(bcio.io)
end

function dump_block_data(block::Block)

    data = Array{UInt8}(undef, block.size)

    copyto!(data, 1, block.header.data)

    transaction_counter = to_varint(block.transaction_counter)
    copyto!(data, 81, transaction_counter)
    idx = 81 + length(transaction_counter)

    for i in 1:block.transaction_counter
        idx = dump_tx_data!(data, idx, block.transactions[i])
    end

    @assert idx == length(data) + 1

    return data
end

function Block(fp::FilePointer)
    fp.file_number |>
        get_block_chain_file_path |>
        x -> open(x) do fh
            seek(fh, fp.file_position)
            Block(fh)
        end
end


text:
abstract type FarmField <: AgComponent end

@with_kw mutable struct CropField <: FarmField
    name::String
    total_area_ha::Union{Float64, AgParameter}
    irrigation::Irrigation
    crop::Crop  # Initial value can just be the first item in crop_rotation
    # crop_choices::Array{Crop}  # This can be the unique values in crop_rotation
    crop_rotation::Array{Crop}
    soil_TAW::Union{Float64, AgParameter}
    soil_SWD::Union{Float64, AgParameter}  # soil water deficit
    ssm::Union{Float64, AgParameter} = 0.0
    irrigated_area::Union{Nothing, Float64} = 0.0
    sowed::Bool = false
    _irrigated_volume::DefaultDict = DefaultDict{String, Float64}(0.0)
    _num_irrigation_events::Int64 = 0
    _irrigation_cost::Float64 = 0.0

    # next crop will be 2nd item in crop_rotation and this
    # counter will increment and reset as the seasons go by
    _next_crop_idx::Int64 = 2
    _fname::String = replace("$(name)-", " " => "_")
    _seasonal_log::DataFrame = setup_log()

end

"""Getter for Field"""
function Base.getproperty(f::FarmField, v::Symbol)
    if v == :irrigated_volume
        return sum(values(f._irrigated_volume))

    elseif v == :irrigation_from_source
        return f._irrigated_volume

    elseif v == :irrigated_vol_mm
        irrig_area = f.irrigated_area
        if irrig_area == 0.0
            @assert f.irrigated_volume == 0.0 "Irrigation occured but irrigated area is 0!"
            return 0.0
        end

        val::Float64 = (f.irrigated_volume / irrig_area) * ML_to_mm

        return val

    elseif v == :irrigation_cost
        return f._irrigation_cost

    elseif v == :dryland_area
        return f.total_area_ha - f.irrigated_area

    elseif v == :plant_date
        return f.crop.plant_date

    elseif v == :harvest_date
        return f.crop.harvest_date
    elseif v == :num_crops
        return length(f.crop_rotation)
    else
        return getfield(f, v)

    end
end

"""Setter for Field"""
function Base.setproperty!(f::FarmField, v::Symbol, value)::Nothing
    if v == :irrigated_volume
        if isnothing(value) || length(value) == 1
            for (k, _v) in f._irrigated_volume
                f._irrigated_volume[k] = value
            end
        elseif length(value) == 2
            f._irrigated_volume[value[1]] += value[2]
        end
    elseif v == :plant_date
        f.crop.plant_date = value
    elseif v == :harvest_date
        f.crop.harvest_date = value
    else
        setfield!(f, Symbol(v), value)
    end

    return nothing
end

"""Volume used from a water source in ML"""
function volume_used_by_source(f::FarmField, ws_name::String)::Float64
    if haskey(f._irrigated_volume, Symbol(ws_name))
        return f._irrigated_volume[Symbol(ws_name)] / mm_to_ML
    end

    return 0.0
end

"""Log the (total) cost of irrigation."""
function log_irrigation_cost(f::FarmField, costs::Float64)::Nothing
    f._irrigation_cost += costs
    return nothing
end

"""
Log seasonal results.

*
"""
function seasonal_field_log!(f::FarmField, dt::Date, income::Float64,
                             irrig_vol::Float64, irrig_yield::Float64,
                             dry_yield::Float64, seasonal_rainfall::Float64)::Nothing
    push!(f._seasonal_log, [dt, income, irrig_vol, irrig_yield, dry_yield, seasonal_rainfall, f.irrigated_area, f.dryland_area])

    return nothing
end

@doc """
    update_SWD!(f::FarmField, rainfall::Float64, ET::Float64)

Calculate soil water deficit, and update `f.soil_SWD`

Water deficit is represented as positive values.

`rainfall` and `ET` parameters are expected in mm.
"""
function update_SWD!(f::FarmField, rainfall::Float64, ET::Float64)::Nothing
    tmp::Float64 = f.soil_SWD::Float64 - (rainfall - ET)::Float64

    tmp = max(0.0, min(tmp, f.soil_TAW))
    f.soil_SWD = round(tmp, digits=4)

    return nothing
end

@doc """
    nid(f::FarmField, dt::Date)

Calculate net irrigation depth in mm, 0.0 or above.

Equation taken from [Agriculture Victoria](http://agriculture.vic.gov.au/agriculture/horticulture/vegetables/vegetable-growing-and-management/estimating-vegetable-crop-water-use)

See also:
* http://www.fao.org/docrep/x5560e/x5560e03.htm
* https://www.bae.ncsu.edu/programs/extension/evans/ag452-1.html
* http://dpipwe.tas.gov.au/Documents/Soil-water_factsheet_14_12_2011a.pdf
* https://www.agric.wa.gov.au/water-management/calculating-readily-available-water?nopaging=1

``NID`` = Effective root depth (``D_{rz}``) ``*`` Readily Available Water (``RAW``)

where:

* ``D_{rz}`` = ``Crop_{root depth} * Crop_{e_rz}``, where ``Crop_{root depth}``
  is the estimated root depth for current stage of crop (initial, late, etc.)
  and ``Crop_{e rz}`` is the effective root zone coefficient for the crop.
* ``Crop_{e rz}`` is said to be between 1 and 2/3rds of total root depth
* ``RAW = p * TAW``, ``p`` is depletion fraction of crop, ``TAW`` is
  Total Available Water in Soil

As an example, if a crop has a root depth (``RD_{r}``) of 1m, an effective
root zone (``RD_{erz}``) coefficient of 0.55, a depletion fraction (`p`) of 0.4
and the soil has a ``TAW`` of 180mm:

``(RD_{r} * RD_{erz}) * (p * TAW)``

Works out to be:

``(1 * 0.55) * (0.4 * 180)``

Returns
-------
    * float : net irrigation depth
"""
function nid(f::FarmField, dt::Date)::Float64
    crop::Crop = f.crop
    coefs::NamedTuple = get_stage_coefs(crop, dt)

    e_rootzone_m::Float64 = crop.root_depth_m * crop.effective_root_zone
    soil_RAW::Float64 = f.soil_TAW * coefs.depletion_fraction

    return (e_rootzone_m * soil_RAW)
end

@doc """
    calc_required_water(f::FarmField, dt::Date)

Volume of water to maintain moisture above net irrigation depth (`nid`).

Calculates volume of water needed to replenish soil water deficit (`swd`)
when SWD falls below `nid`, considering irrigation efficiency.

Values are given in mm.

- link to [`nid(f::FarmField, dt::Date)`](@ref)
"""
function calc_required_water(f::FarmField, dt::Date)::Float64
    swd::Float64 = f.soil_SWD::Float64
    to_nid::Float64 = swd - nid(f, dt)
    if to_nid < 0.0
        return 0.0
    end

    return (swd / f.irrigation.efficiency)
end

@doc """Possible irrigation area in hectares."""
function possible_irrigation_area(f::FarmField, vol_ML::Float64, req_ML::Float64)::Float64
    if vol_ML == 0.0
        return 0.0
    end

    area::Float64 = 0.0
    if f.irrigation.name != "dryland"
        area = f.irrigated_area == 0.0 ? f.total_area_ha : f.irrigated_area
    end

    if area == 0.0
        return 0.0
    end

    if req_ML == 0.0
        return area
    end

    perc_area::Float64 = vol_ML / (req_ML * area)

    return min(perc_area * area, area)
end

function set_next_crop!(f::FarmField)::Nothing
    crop_id::Int64 = f._next_crop_idx > f.num_crops ? 1 : f._next_crop_idx

    # Update crop and next crop id
    f.crop, f._next_crop_idx = iterate(f.crop_rotation, crop_id)

    return nothing
end

function reset_state!(f::FarmField)::Nothing
    f.sowed = false

    f.irrigated_volume = 0.0  # This clears the underlying log as well.
    f.irrigated_area = 0.0
    f._irrigation_cost = 0.0
    f._num_irrigation_events = 0
    f.soil_SWD = 0.0

    return nothing
end

@doc """
    french_schultz_cy(; ssm_mm::Float64, gsr_mm::Float64, crop::AgComponent)

Potential crop yield calculation based on a modified French-Schultz equation.
The implemented method is the farmer-friendly version as described by
Oliver et al., (2008) (see [1]).

``YP = (SSM + GSR - E) * WUE``

where

* ``YP`` is yield potential in tonnes per hectare
* ``SSM`` is Stored Soil Moisture (at start of season) in mm, assumed to be 30% of summer rainfall
* ``GSR`` is Growing Season Rainfall in mm
* ``E`` is Crop Evaporation coefficient in mm, the amount of rainfall required before the crop will start
    to grow, commonly 110mm, but can range from 30-170mm (see [2]),
* ``WUE`` is Water Use Efficiency coefficient in kg/mm

References
----------
1. [Oliver et al. (2008) (see Equation 1)](http://www.regional.org.au/au/asa/2008/concurrent/assessing-yield-potential/5827_oliverym.htm)
2. [Whitbread and Hancock (2008)](http://www.regional.org.au/au/asa/2008/concurrent/assessing-yield-potential/5803_whitbreadhancock.htm)

Parameters
----------
* ssm_mm : float, Stored Soil Moisture (mm) at start of season.
* gsr_mm : float, Growing Season Rainfall (mm)
* crop : object, Crop component object

Returns
-----------
* Potential yield in tonnes/Ha
"""
function french_schultz_cy(; ssm_mm::Float64, gsr_mm::Float64,
                             crop::AgComponent)::Float64
    evap_coef_mm::Float64 = crop.et_coef  # Crop evapotranspiration coefficient (mm)
    wue_coef_mm::Float64 = crop.wue_coef  # Water Use Efficiency coefficient (kg/mm)

    # maximum rainfall threshold in mm
    # water above this amount does not contribute to crop yield
    max_thres::Float64 = crop.rainfall_threshold

    gsr_mm::Float64 = min(gsr_mm, max_thres)
    return max(0.0, ((ssm_mm + gsr_mm - evap_coef_mm) * wue_coef_mm) / 1000.0)
end

"""
    calc_potential_crop_yield(method::Function; kwargs...)

Calculate potential crop yield using an arbitrary function.
This function enables a consistent interface to be provided.

Function must return a Tuple of results or a Float64
"""
function calc_potential_crop_yield(method::Function; kwargs...)::Union{Tuple, Float64}
    return method(; kwargs...)
end

"""
    total_income(f::FarmField, ssm::Float64, gsr::Float64, irrig::Float64, comps)::Tuple

Calculate net income considering crop yield and costs incurred.

Returns
-------
    Tuple:
        * float : net income
        * float : irrigated crop yield [t/ha]
        * float : dryland crop yield [t/ha]
"""
function total_income(f::FarmField, ssm::Float64, gsr::Float64,
                      irrig::Float64, comps)::Tuple
    inc::Float64, irrigated::Float64, dryland::Float64 = gross_income(f, ssm, gsr, irrig)
    return inc - total_costs(f, comps...), irrigated, dryland
end

"""
    gross_income(f::FarmField, ssm::Float64, gsr::Float64,
        irrig::Float64, func::Function=french_schultz_cy)::Tuple

Calculate gross income, potential irrigated yield, and potential dryland yield.

Returns
-------
Tuple{Float64} : income, irrigated yield [t], dryland yield [t]
"""
function gross_income(f::FarmField, ssm_mm::Float64, gsr_mm::Float64,
                      irrig::Float64, func::Function=french_schultz_cy)::Tuple
    crop::Crop = f.crop

    income::Float64 = 0.0
    total_irrig_yield::Float64 = 0.0
    if f.irrigated_area > 0.0
        irrigated_yield::Float64 = calc_potential_crop_yield(func; ssm_mm, gsr_mm=gsr_mm+irrig, crop)
        total_irrig_yield = irrigated_yield * f.irrigated_area
        income = irrigated_yield * f.irrigated_area * crop.price_per_yield
    end

    dryland_yield::Float64 = calc_potential_crop_yield(func; ssm_mm, gsr_mm, crop)
    total_dry_yield::Float64 = dryland_yield * f.dryland_area

    income += dryland_yield * f.dryland_area * crop.price_per_yield

    return income, total_irrig_yield, total_dry_yield
end

"""
    gross_income(f::FarmField, area::Float64, func::Function=french_schultz_cy; kwargs...)::Tuple

Calculate gross income and crop yield using an arbitrary crop yield function.

Returns
-------
Tuple{Float64} : income [\$/yield], total_yield [t]
"""
function gross_income(f::FarmField, area::Float64, func::Function=french_schultz_cy; kwargs...)::Tuple
    crop::Crop = f.crop
    crop_yield_per_area::Float64 = calc_potential_crop_yield(func; kwargs...)

    total_yield::Float64 = crop_yield_per_area * area
    income::Float64 = total_yield * crop.price_per_yield

    return income, total_yield
end

"""Calculate total costs for a field.

Maintenance costs can be spread out across a number of fields if desired.
"""
function total_costs(f::FarmField, dt::Date, water_sources::Array{WaterSource}, num_fields::Int64=1)::Float64
    year_val::Int64 = year(dt)
    irrig_area::Float64 = f.irrigated_area
    h20_usage_cost::Float64 = 0.0
    maint_cost::Float64 = 0.0
    for w::WaterSource in water_sources
        water_used::Float64 = volume_used_by_source(f, w.name)
        ws_cost::Float64 = subtotal_costs(w, irrig_area, water_used)
        h20_usage_cost = h20_usage_cost + ws_cost

        pump_cost::Float64 = subtotal_costs(w.pump, year_val) / num_fields
        maint_cost = maint_cost + pump_cost
    end

    irrig_app_cost::Float64 = f._irrigation_cost

    if irrig_app_cost > 0.0
        @assert f.irrigated_volume > 0.0 "Irrigation had to occur for costs to be incurred!"
    elseif f.irrigated_volume > 0.0
        @assert irrig_app_cost > 0.0 "If irrigation occured, costs have to be incurred!"
    end

    h20_usage_cost = h20_usage_cost + irrig_app_cost

    maint_cost = maint_cost + subtotal_costs(f.irrigation, year_val)
    crop_costs::Float64 = subtotal_costs(f.crop, year_val)
    total_costs::Float64 = h20_usage_cost + maint_cost + crop_costs

    return total_costs
end

function setup_log()::DataFrame
    return DataFrame([Date, Float64, Float64, Float64, Float64, Float64, Float64, Float64],
                     [:Date, :income, :irrigated_volume, :irrigated_yield, :dryland_yield,
                      :growing_season_rainfall, :irrigated_area, :dryland_area])
end

function reset!(f::FarmField)::Nothing
    f._seasonal_log = setup_log()

    reset_state!(f)

    return
end

max_stars_repo_name: ConnectedSystems/agtor-jl

max_issues_repo_name: ConnectedSystems/agtor-jl

max_forks_repo_name: ConnectedSystems/agtor-jl


text:
using Test
using LogicCircuits
using Random: bitrand, rand
using DataFrames: DataFrame

include("../helper/gpu.jl")

@testset "Binary flows test" begin

    r = fully_factorized_circuit(PlainLogicCircuit, 10)
    input = [1 0 1 0 1 0 1 0 1 0;
            1 1 1 1 1 1 1 1 1 1;
            0 0 0 0 0 0 0 0 0 0;
            0 1 1 0 1 0 0 1 0 1]

    input = DataFrame(BitArray(input), :auto)
    @test r(input) == BitVector([1,1,1,1])

    vtree = PlainVtree(10, :balanced)
    r = fully_factorized_circuit(StructLogicCircuit, vtree)
    @test r(input) == BitVector([1,1,1,1])

    v, f, node2id= satisfies_flows(r, input)
    foreach(literal_nodes(r)) do n
        id = node2id[n].node_id
        @test v[:,id] == f[:,id] # invariant of logically valid circuits
    end

    cpu_gpu_agree(input) do d
        satisfies_flows(r, d)[1] # same value
    end

    cpu_gpu_agree(input) do d
        satisfies_flows(r, d)[2][3:end,:] # same flows (except constants)
    end

    num_vars = 7
    mgr = SddMgr(num_vars, :balanced)
    v = Dict([(i => compile(mgr, Lit(i))) for i=1:num_vars])
    c = (v[1] | !v[2] | v[3]) &
        (v[2] | !v[7] | v[6]) &
        (v[3] | !v[4] | v[5]) &
        (v[1] | !v[4] | v[6])
    input = DataFrame(bitrand(25,num_vars), :auto)

    r = smooth(PlainLogicCircuit(c)) # flows don't make sense unless the circuit is smooth; cannot smooth trimmed SDDs

    v, f, node2id = satisfies_flows(r, input)
    foreach(literal_nodes(r)) do n
        id = node2id[n].node_id
        @test v[:,id] .& v[:,end] == f[:,id] # invariant of all circuits
    end

    cpu_gpu_agree(input) do d
        satisfies_flows(r, d)[1] # same value
    end

    cpu_gpu_agree(input) do d
        satisfies_flows(r, d)[2][3:end,:] # same flows (except constants)
    end

    l1 = LogicCircuit(Lit(1))
    l2 = LogicCircuit(Lit(2))
    l3 = LogicCircuit(Lit(-1))
    l4 = LogicCircuit(Lit(-2))
    r = (l1 & l2) | (l3 & l4)
    input = DataFrame(bitrand(4,2), :auto)

    v, f, node2id = satisfies_flows(r, input)
    foreach(literal_nodes(r)) do n
        id = node2id[n].node_id
        @test v[:,id] .& v[:,end] == f[:,id] # invariant of all circuits
    end

    cpu_gpu_agree(input) do d
        satisfies_flows(r, d)[1] # same value
    end

    cpu_gpu_agree(input) do d
        satisfies_flows(r, d)[2][3:end,:] # same flows (except constants)
    end

end

@testset "Weighted binary flows test" begin

    r = fully_factorized_circuit(PlainLogicCircuit, 10)
    input = [1 0 1 0 1 0 1 0 1 0;
            1 1 1 1 1 1 1 1 1 1;
            0 0 0 0 0 0 0 0 0 0;
            0 1 1 0 1 0 0 1 0 1]

    input = DataFrame(BitArray(input), :auto)
    @test r(input) == BitVector([1,1,1,1])

    weights = [0.6, 0.6, 0.6, 0.6]
    weights = Array{Float32, 1}(weights)

    vtree = PlainVtree(10, :balanced)
    r = fully_factorized_circuit(StructLogicCircuit, vtree)
    @test r(input) == BitVector([1,1,1,1])

    v, f, node2id = satisfies_flows(r, input; weights = weights)
    foreach(literal_nodes(r)) do n
        id = node2id[n].node_id
        @test v[:,id] == f[:,id] # invariant of logically valid circuits
    end

    if CUDA.functional()
        @test all(satisfies_flows(r, input; weights = weights)[1] .≈ to_cpu(satisfies_flows(r, to_gpu(input); weights = to_gpu(weights))[1]))

        @test all(satisfies_flows(r, input; weights)[2][3:end,:] .≈ to_cpu(satisfies_flows(r, to_gpu(input); weights = to_gpu(weights))[2][3:end,:]))
    end

    num_vars = 7
    mgr = SddMgr(num_vars, :balanced)
    v = Dict([(i => compile(mgr, Lit(i))) for i=1:num_vars])
    c = (v[1] | !v[2] | v[3]) &
        (v[2] | !v[7] | v[6]) &
        (v[3] | !v[4] | v[5]) &
        (v[1] | !v[4] | v[6])
    input = DataFrame(bitrand(25,num_vars), :auto)

    weights = Array{Float32, 1}(ones(25) * 0.6)

    r = smooth(PlainLogicCircuit(c)) # flows don't make sense unless the circuit is smooth; cannot smooth trimmed SDDs

    v, f, node2id = satisfies_flows(r, input; weights = weights)
    foreach(literal_nodes(r)) do n
        id = node2id[n].node_id
        @test v[:,id] .& v[:,end] == f[:,id] # invariant of all circuits
    end

    if CUDA.functional()
        @test all(satisfies_flows(r, input; weights = weights)[1] .≈ to_cpu(satisfies_flows(r, to_gpu(input); weights = to_gpu(weights))[1]))

        @test all(satisfies_flows(r, input; weights)[2][3:end,:] .≈ to_cpu(satisfies_flows(r, to_gpu(input); weights = to_gpu(weights))[2][3:end,:]))
    end

    l1 = LogicCircuit(Lit(1))
    l2 = LogicCircuit(Lit(2))
    l3 = LogicCircuit(Lit(-1))
    l4 = LogicCircuit(Lit(-2))
    r = (l1 & l2) | (l3 & l4)
    input = DataFrame(bitrand(4,2), :auto)

    v, f, node2id = satisfies_flows(r, input; weights = weights)
    foreach(literal_nodes(r)) do n
        id = node2id[n].node_id
        @test v[:,id] .& v[:,end] == f[:,id] # invariant of all circuits
    end

    if CUDA.functional()
        @test all(satisfies_flows(r, input; weights = weights)[1] .≈ to_cpu(satisfies_flows(r, to_gpu(input); weights = to_gpu(weights))[1]))

        @test all(satisfies_flows(r, input; weights = weights)[2][3:end,:] .≈ to_cpu(satisfies_flows(r, to_gpu(input); weights = to_gpu(weights))[2][3:end,:]))
    end

end

@testset "Soft Flows test" begin
    l_a = LogicCircuit(Lit(1))
    l_na = LogicCircuit(Lit(-1))
    l_b = LogicCircuit(Lit(2))
    l_nb = LogicCircuit(Lit(-2))

    o_c = (l_a & l_b) | (l_na & l_b)
    o_d = (l_a & l_b) | (l_a & l_nb)
    r = o_c | o_d

    weights = [0.6, 0.6, 0.6]
    weights = Array{Float32, 1}(weights)

    df = DataFrame(BitMatrix([true true; false true; false false]), :auto)
    sdf = soften(df, 0.001; scale_by_marginal = false)

    v, f, node2id = satisfies_flows(r, sdf)
    @test all(v[:,node2id[o_c].node_id] .≈ sdf[:, 2])
    @test all(v[:,node2id[o_c].node_id] .≈ f[:,node2id[o_c].node_id])
    @test all(v[:,node2id[o_d].node_id] .≈ sdf[:, 1])
    @test all(v[:,node2id[o_d].node_id] .≈ f[:,node2id[o_d].node_id])

    v, f, node2id = satisfies_flows(r, sdf; weights = weights)
    @test all(v[:,node2id[o_c].node_id] .≈ sdf[:, 2])
    @test all(v[:,node2id[o_c].node_id] .≈ f[:,node2id[o_c].node_id])
    @test all(v[:,node2id[o_d].node_id] .≈ sdf[:, 1])
    @test all(v[:,node2id[o_d].node_id] .≈ f[:,node2id[o_d].node_id])

    if CUDA.functional()
        @test all(satisfies_flows(r, sdf; weights = weights)[1] .≈ to_cpu(satisfies_flows(r, to_gpu(sdf); weights = to_gpu(weights))[1]))
    end
end

@testset "Probabilistic Flows test" begin

    r = fully_factorized_circuit(PlainLogicCircuit, 10)

    input = DataFrame(Float64[1 0 1 0 1 0 1 0 1 0;
                    1 1 1 1 1 1 1 1 1 1;
                    0 0 0 0 0 0 0 0 0 0;
                    0 1 1 0 1 0 0 1 0 1], :auto)

    @test r(input) ≈ [1.0, 1.0, 1.0, 1.0]

    v, f, node2id = satisfies_flows(r, input)
    foreach(literal_nodes(r)) do n
        id = node2id[n].node_id
        @test v[:,id] ≈ f[:,id] # invariant of logically valid circuits
    end

    cpu_gpu_agree_approx(input) do d
        satisfies_flows(r, d)[1] # same value
    end

    cpu_gpu_agree_approx(input) do d
        satisfies_flows(r, d)[2][3:end,:] # same flows (except constants)
    end

    l_a = LogicCircuit(Lit(1))
    l_na = LogicCircuit(Lit(-1))
    l_b = LogicCircuit(Lit(2))
    l_nb = LogicCircuit(Lit(-2))
    l_c = LogicCircuit(Lit(3))
    l_nc = LogicCircuit(Lit(-3))

    o_c = (l_a & l_b & l_c) | (l_a & l_nb & l_c) | (l_na & l_b & l_c) | (l_na & l_nb & l_c)
    o_nc = (l_a & l_b & l_nc) | (l_a & l_nb & l_nc) | (l_na & l_b & l_nc) | (l_na & l_nb & l_nc)
    r = (o_c | o_nc)

    input = DataFrame(rand(Float64, (10,3)), :auto)

    @test all(r(input) .≈ 1.0)

    v, f, node2id = satisfies_flows(r, input)
    @test all(v[:,end] .≈ 1.0)
    @test all(f[:,end].≈ 1.0)
    @test all(v[:,node2id[o_c].node_id] .≈ input[:, 3])

    @test all(f[:,node2id[o_c].node_id] .≈ input[:, 3])
    @test all(v[:,node2id[o_nc].node_id] .+ input[:, 3] .≈ 1.0)
    @test all(f[:,node2id[o_nc].node_id] .+ input[:, 3] .≈ 1.0)

    foreach(literal_nodes(r)) do n
        id = node2id[n].node_id
        @test v[:,id] ≈ f[:,id] # invariant of logically valid circuits
    end

    cpu_gpu_agree_approx(input) do d
        satisfies_flows(r, d)[1] # same value
    end

    cpu_gpu_agree_approx(input; atol=1e-6) do d
        satisfies_flows(r, d)[2][3:end,:] # same flows (except constants)
    end

    o_c = (l_a & l_b & l_c) | (l_na & l_b & l_c)
    o_nc = (l_a & l_nb & l_nc) | (l_na & l_nb & l_nc)
    r = (o_c | o_nc)

    input = DataFrame(rand(Float64, (10,3)), :auto)
    input[1:5, 3] .= 0.0

    v, f, node2id = satisfies_flows(r, input)
    @test all(v[:,node2id[o_c].node_id] .≈ input[:, 2] .* input[:, 3])
    @test all(f[:,node2id[o_c].node_id] .≈ v[:,node2id[o_c].node_id])
    @test all(v[:,node2id[o_nc].node_id] .≈ (1.0 .- input[:, 2]) .* (1.0 .- input[:, 3]))
    @test all(f[:,node2id[o_nc].node_id] .≈ v[:,node2id[o_nc].node_id])
    @test all(f[:,node2id[l_a].node_id] .≈ v[:,node2id[l_a].node_id] .* f[:,node2id[r].node_id])
    @test all(f[:,node2id[l_na].node_id] .≈ v[:,node2id[l_na].node_id] .* f[:,node2id[r].node_id])
    @test all(f[:,node2id[l_b].node_id] .≈ input[:, 2] .* input[:, 3])
    @test all(f[:,node2id[l_b].node_id] .≈ f[:,node2id[l_c].node_id])
    @test all(f[:,node2id[l_nb].node_id] .≈ (1.0 .- input[:, 2]) .* (1.0 .- input[:, 3]))
    @test all(f[:,node2id[l_nb].node_id] .≈ f[:,node2id[l_nc].node_id])

    cpu_gpu_agree_approx(input) do d
        satisfies_flows(r, d)[1] # same value
    end

    cpu_gpu_agree_approx(input) do d
        satisfies_flows(r, d)[2][3:end,:] # same flows (except constants)
    end

end

@testset "Downflow API test" begin
    r = fully_factorized_circuit(PlainLogicCircuit, 10)

    input_b = DataFrame(BitArray([1 0 1 0 1 0 1 0 1 0;
                    1 1 1 1 1 1 1 1 1 1;
                    0 0 0 0 0 0 0 0 0 0;
                    0 1 1 0 1 0 0 1 0 1]), :auto)

    input_f = DataFrame(Float64.(Matrix(input_b)), :auto)

    f_b, v_b, node2id = satisfies_flows(r, input_b)
    f_f, v_f, node2id = satisfies_flows(r, input_f)

    N = num_examples(input_b)
    foreach(r) do n
        if is⋁gate(n)
            @test all(Float64.(downflow_all(f_b, v_b, N, n, node2id)) .≈ downflow_all(f_f, v_f, N, n, node2id))
            @test Float64.(count_downflow(f_b, v_b, N, n, node2id)) ≈ count_downflow(f_f, v_f, N, n, node2id)
            for c in children(n)
                @test all(Float64.(downflow_all(f_b, v_b, N, n, c, node2id)) .≈ downflow_all(f_f, v_f, N, n, c, node2id))
                @test Float64.(count_downflow(f_b, v_b, N, n, c, node2id)) ≈ count_downflow(f_f, v_f, N, n, c, node2id)
            end
        end
    end

    input1 = DataFrame(BitArray([1 0 1 0 1 0 1 0 1 0;
                                1 1 1 1 1 1 1 1 1 1;
                                0 0 0 0 0 0 0 0 0 0;
                                0 1 1 0 1 0 0 1 0 1]), :auto)

    input2 = DataFrame(BitArray([0 1 1 1 0 1 1 1 0 0;
                                1 1 0 0 1 1 0 1 0 1;
                                1 0 0 0 1 0 1 0 1 0;
                                0 0 1 1 0 1 1 0 0 1]), :auto)

    input3 = DataFrame(BitArray([1 0 1 0 1 0 1 0 1 0;
                                0 1 0 1 0 1 1 1 0 1;
                                0 1 0 1 0 1 1 0 1 0;
                                0 0 1 1 1 0 1 0 0 1]), :auto)

    input4 = DataFrame(0.3 .* Float64.(Matrix(input1)) .+
                        0.1 .* Float64.(Matrix(input2)) .+
                        0.6 .* Float64.(Matrix(input3)), :auto)

    inputs = [input1, input2, input3, input4]
    vfs = [satisfies_flows(r, input) for input in inputs]

    N = num_examples(input_b)
    foreach(r) do n
        if is⋁gate(n)
            df = [count_downflow(v, f, N, n, node2id) for (v, f) in vfs]
            @test 0.3 * df[1] + 0.1 * df[2] + 0.6 * df[3] ≈ df[4]

            df_all = [downflow_all(v, f, N, n, node2id) for (v, f) in vfs]
            @test all(0.3 * Float64.(df_all[1]) + 0.1 * Float64.(df_all[2]) + 0.6 * Float64.(df_all[3]) .≈ df_all[4])
            for c in children(n)
                df = [count_downflow(v, f, N, n, c, node2id) for (v, f) in vfs]
                @test 0.3 * df[1] + 0.1 * df[2] + 0.6 * df[3] ≈ df[4]
                df_all = [downflow_all(v, f, N, n, c, node2id) for (v, f) in vfs]
                @test all(0.3 * Float64.(df_all[1]) + 0.1 * Float64.(df_all[2]) + 0.6 * Float64.(df_all[3]) .≈ df_all[4])
            end
        end
    end
end

max_issues_repo_name: Juice-jl/LogicCircuits

max_forks_repo_name: Juice-jl/LogicCircuits


text:
module Grover

using QXZoo.NCU
using QXZoo.Circuit
using QXZoo.DefaultGates
using QXZoo.GateOps

export run_grover

"""
    calc_iterations(num_qubits::Integer)

Calculate the required number of iterations to maximise the state's amplitude.
"""
function calc_iterations(num_qubits::Integer)
    return pi*sqrt(num_qubits)/4
end

"""
    state_init!(cct::Circuit.Circ, qubit_indices::Vector)

Initialises the state to the required target; defaults to ``H^{\\otimes n}\\vert \\psi \\rangle`` . Override for custom initialisation.
"""
function state_init!(cct::Circuit.Circ, qubit_indices::Vector)
    for i in qubit_indices
        Circuit.add_gatecall!(cct, DefaultGates.h(i))
    end
end

"""
    bitstring_ncu(cct::Circuit.Circ, bitstring::Integer, ctrl_indices::Vector, tgt_idx::Int, U::GateOps.GateSymbol, aux_indices::Vector=Int[])

Takes bitstring as the binary pattern and indices as the qubits to operate upon. Applies the appropriate PauliX gates to the control lines to call the NCU with the given matrix. Uses aux_indices to reduce Circuit depth by expanding width if |aux|+2 >= |ctrl|. If not specified, uses default quadratic gate expansion.
"""
function bitstring_ncu!(cct::Circuit.Circ, bitstring::Integer, ctrl_indices::Vector, tgt_idx::Int, U::GateOps.GateSymbol, aux_indices::Vector=Int[])
    bitmask =  0x1

    # Filter qubit values to mark specified pattern
    for idx in 0:length(ctrl_indices)
        if ~( (bitstring & (bitmask << idx) ) != 0)
            if idx < length(ctrl_indices)
                cct << DefaultGates.x(ctrl_indices[idx+1])
            else
                cct << DefaultGates.x(tgt_idx)
            end
        end
    end
    NCU.apply_ncu!(cct, ctrl_indices, aux_indices, tgt_idx, U)
    for idx in length(ctrl_indices):-1:0
        if ~( (bitstring & (bitmask << idx) ) != 0)
            if idx < length(ctrl_indices)
                cct << DefaultGates.x(ctrl_indices[idx+1])
            else
                cct << DefaultGates.x(tgt_idx)
            end
        end
    end
end

"""
    bitstring_phase_oracle(cct::Circuit.Circ, bitstring::Integer, ctrl_indices::Vector, tgt_idx::Int, aux_indices::Vector=Int[])

Applies PauliX gates to the appropriate lines in the Circuit, then applies a n-controlled PauliZ to mark the state. Uses aux qubits to reduce circuit depth if |aux|+2 >= |ctrl|. If not provided, uses quadratic gate expansion
"""
function bitstring_phase_oracle!(cct::Circuit.Circ, bitstring::Integer, ctrl_indices::Vector, tgt_idx::Int, aux_indices::Vector=Int[])
    return bitstring_ncu!(cct, bitstring, ctrl_indices, tgt_idx, DefaultGates.GateSymbols.c_z, aux_indices)
end

"""
    apply_diffusion(cct::Circuit.Circ, ctrl_indices::Vector, tgt_index::Int, aux_indices::Vector=Int[] )

Application of the Grover diffusion operator to marked register. Uses additionally provided auxiliary qubits to reduce depth.
"""
function apply_diffusion!(cct::Circuit.Circ, ctrl_indices::Vector, tgt_index::Int, aux_indices::Vector=Int[])

    for ctrl in vcat(ctrl_indices, tgt_index)
        Circuit.add_gatecall!(cct, DefaultGates.h(ctrl) )
        Circuit.add_gatecall!(cct, DefaultGates.x(ctrl) )
    end

    NCU.apply_ncu!(cct, ctrl_indices, aux_indices, tgt_index, DefaultGates.GateSymbols.c_z)

    for ctrl in vcat(ctrl_indices, tgt_index)
        Circuit.add_gatecall!(cct, DefaultGates.x(ctrl) )
        Circuit.add_gatecall!(cct, DefaultGates.h(ctrl) )
    end
end

"""
    apply_grover_iteration!(cct::Circuit.Circ, qubit_indices::Vector, qubit_aux_indices::Vector)

Applies a single Grover iteration. To be used following `mark_state!`
"""
function apply_grover_iteration!(cct::Circuit.Circ, qubit_indices::Vector, qubit_aux_indices::Vector=Int[])
    return apply_diffusion!(cct, qubit_indices[1:end-1], qubit_indices[end], qubit_aux_indices)
end

"""
    mark_state!(cct::Circuit.Circ, state::Integer, qubit_indices::Vector, qubit_aux_indices::Vector)

Applies the state marking procedure of the Grover iteration.
"""
function mark_state!(cct::Circuit.Circ, state::Integer, qubit_indices::Vector, qubit_aux_indices::Vector=Int[])
    return bitstring_phase_oracle!(cct, state, qubit_indices[1:end-1], qubit_indices[end], qubit_aux_indices)
end

"""
    run_grover!(cct::Circuit.Circ, qubit_indices::Vector, state::Integer=1, qubit_aux_indices::Vector=Int[])

Generates a Grover search Circuit sample, marking the state defined by `state`
and performing iterations to amplify the desired result upon measurement.
"""
function run_grover!(cct::Circuit.Circ, qubit_indices::Vector, state::Integer=1, qubit_aux_indices::Vector=Int[])
    num_states = 2^(length(qubit_indices))
    state_init!(cct, qubit_indices)
    for i in range(1,stop=ceil(Integer, calc_iterations(num_states)))
        mark_state!(cct, state, qubit_indices, qubit_aux_indices)
        apply_grover_iteration!(cct, qubit_indices, qubit_aux_indices)
    end
end

end


text:
# Helper routines

function reverseeven!(x::AbstractVector)
    n = length(x)
    if iseven(n)
        @inbounds @simd for k=2:2:n÷2
            x[k],x[n+2-k] = x[n+2-k],x[k]
        end
    else
        @inbounds @simd for k=2:2:n÷2
            x[k],x[n+1-k] = x[n+1-k],x[k]
        end
    end
    x
end

function negateeven!(x::AbstractVector)
    @inbounds @simd for k = 2:2:length(x)
        x[k] *= -1
    end
    x
end

### In-place O(n) interlacing

function highestleader(n::Int)
    i = 1
    while 3i < n i *= 3 end
    i
end

function nextindex(i::Int,n::Int)
    i <<= 1
    while i > n
        i -= n + 1
    end
    i
end

function cycle_rotate!(v::AbstractVector, leader::Int, it::Int, twom::Int)
    i = nextindex(leader, twom)
    while i != leader
        idx1, idx2 = it + i - 1, it + leader - 1
        @inbounds v[idx1], v[idx2] = v[idx2], v[idx1]
        i = nextindex(i, twom)
    end
    v
end

function right_cyclic_shift!(v::AbstractVector, it::Int, m::Int, n::Int)
    itpm = it + m
    itpmm1 = itpm - 1
    itpmpnm1 = itpmm1 + n
    reverse!(v, itpm, itpmpnm1)
    reverse!(v, itpm, itpmm1 + m)
    reverse!(v, itpm + m, itpmpnm1)
    v
end

"""
This function implements the algorithm described in:

    P. Jain, "A simple in-place algorithm for in-shuffle," arXiv:0805.1598, 2008.
"""
function interlace!(v::AbstractVector,offset::Int)
    N = length(v)
    if N < 2 + offset
        return v
    end

    it = 1 + offset
    m = 0
    n = 1

    while m < n
        twom = N + 1 - it
        h = highestleader(twom)
        m = h > 1 ? h÷2 : 1
        n = twom÷2

        right_cyclic_shift!(v,it,m,n)

        leader = 1
        while leader < 2m
            cycle_rotate!(v, leader, it, 2m)
            leader *= 3
        end

        it += 2m
    end
    v
end


text:
##

# Convert Julia types to LLVM types
const LLVMTYPE = Dict{DataType,String}(
    Bool => "i8",   # Julia represents Tuple{Bool} as [1 x i8]
    Int8 => "i8",
    Int16 => "i16",
    Int32 => "i32",
    Int64 => "i64",
    Int128 => "i128",
    UInt8 => "i8",
    UInt16 => "i16",
    UInt32 => "i32",
    UInt64 => "i64",
    UInt128 => "i128",
    Float16 => "half",
    Float32 => "float",
    Float64 => "double",
    Nothing => "void"
)
llvmtype(x)::String = LLVMTYPE[x]
const JuliaPointerType = LLVMTYPE[Int]

const LLVMCompatible = Union{Bool,Int8,Int16,Int32,Int64,Int128,UInt8,UInt16,UInt32,UInt64,UInt128,Float16,Float32,Float64}

@generated function vload(ptr::Ptr{T}) where {T <: LLVMCompatible}
    ptyp = JuliaPointerType
    typ = llvmtype(T)
    instrs = String[]
    decl = "!1 = !{!\"noaliasdomain\"}\n!2 = !{!\"noaliasscope\", !1}\n!3 = !{!2}"
    alignment = Base.datatype_alignment(T)
    # push!(flags, "!noalias !0")
    push!(instrs, "%ptr = inttoptr $ptyp %0 to $typ*")
    push!(instrs, "%res = load $typ, $typ* %ptr, align $alignment, !alias.scope !3")
    push!(instrs, "ret $typ %res")
    quote
        $(Expr(:meta, :inline))
        Base.llvmcall(($decl,$(join(instrs, "\n"))), $T, Tuple{Ptr{$T}}, ptr)
    end
end
@generated function vstore!(ptr::Ptr{T}, v::T) where {T <: LLVMCompatible}
    ptyp = JuliaPointerType
    typ = llvmtype(T)
    instrs = String[]
    alignment = Base.datatype_alignment(T)
    # push!(flags, "!noalias !0")
    push!(instrs, "%ptr = inttoptr $ptyp %0 to $typ*")
    push!(instrs, "store $typ %1, $typ* %ptr, align $alignment")
    push!(instrs, "ret void")
    quote
        $(Expr(:meta, :inline))
        Base.llvmcall($(join(instrs, "\n")), Cvoid, Tuple{Ptr{$T}, $T}, ptr, v)
    end
end
@generated function vload(ptr::Ptr{T}, i::I) where {T <: LLVMCompatible, I<:Integer}
    ityp = llvmtype(I)
    ptyp = JuliaPointerType
    typ = llvmtype(T)
    instrs = String[]
    decl = "!1 = !{!\"noaliasdomain\"}\n!2 = !{!\"noaliasscope\", !1}\n!3 = !{!2}"
    alignment = Base.datatype_alignment(T)
    # push!(flags, "!noalias !0")
    push!(instrs, "%typptr = inttoptr $ptyp %0 to $typ*")
    push!(instrs, "%ptr = getelementptr inbounds $typ, $typ* %typptr, $ityp %1")
    push!(instrs, "%res = load $typ, $typ* %ptr, align $alignment, !alias.scope !3")
    push!(instrs, "ret $typ %res")
    quote
        $(Expr(:meta, :inline))
        Base.llvmcall($(decl, join(instrs, "\n")), $T, Tuple{Ptr{$T}, $I}, ptr, i)
    end
end
@generated function vstore!(ptr::Ptr{T}, v::T, i::I) where {T <: LLVMCompatible, I<:Integer}
    ityp = llvmtype(I)
    ptyp = JuliaPointerType
    typ = llvmtype(T)
    instrs = String[]
    alignment = Base.datatype_alignment(T)
    # push!(flags, "!noalias !0")
    push!(instrs, "%typptr = inttoptr $ptyp %0 to $typ*")
    push!(instrs, "%ptr = getelementptr inbounds $typ, $typ* %typptr, $ityp %2")
    push!(instrs, "store $typ %1, $typ* %ptr, align $alignment")
    push!(instrs, "ret void")
    quote
        $(Expr(:meta, :inline))
        Base.llvmcall($(join(instrs, "\n")), Cvoid, Tuple{Ptr{$T}, $T, $I}, ptr, v, i)
    end
end
@generated function vnoaliasstore!(ptr::Ptr{T}, v::T, i::I) where {T <: LLVMCompatible, I<:Integer}
    ityp = llvmtype(I)
    ptyp = JuliaPointerType
    typ = llvmtype(T)
    instrs = String[]
    decl = "!1 = !{!\"noaliasdomain\"}\n!2 = !{!\"noaliasscope\", !1}\n!3 = !{!2}"
    alignment = Base.datatype_alignment(T)
    # push!(flags, "!noalias !0")
    push!(instrs, "%typptr = inttoptr $ptyp %0 to $typ*")
    push!(instrs, "%ptr = getelementptr inbounds $typ, $typ* %typptr, $ityp %2")
    push!(instrs, "store $typ %1, $typ* %ptr, align $alignment, !noalias !3")
    push!(instrs, "ret void")
    quote
        $(Expr(:meta, :inline))
        Base.llvmcall($(decl,join(instrs, "\n")), Cvoid, Tuple{Ptr{$T}, $T, $I}, ptr, v, i)
    end
end

# Fall back definitions
@inline vload(ptr::Ptr) = Base.unsafe_load(ptr)
# @inline vload(ptr::Ptr, i::Int) = vload(gep(ptr, i))
@inline vstore!(ptr::Ptr{T}, v::T) where {T} = Base.unsafe_store!(ptr, v)
@inline vload(::Type{T1}, ptr::Ptr{T2}) where {T1, T2} = vload(Base.unsafe_convert(Ptr{T1}, ptr))
@inline vstore!(ptr::Ptr{T1}, v::T2) where {T1,T2} = vstore!(ptr, convert(T1, v))
@inline vstore!(ptr::Ptr{T1}, v::T2) where {T1<:Integer,T2<:Integer} = vstore!(ptr, v % T1)
# @inline vstore!(ptr::Ptr{T1}, v::T2) where {T1,T2} = vstore!(ptr, convert(T1, v))

@inline tdot(a::Tuple{I1}, b::Tuple{I2}) where {I1,I2} = @inbounds a[1] * b[1]
@inline tdot(a::Tuple{I1,I3}, b::Tuple{I2,I4}) where {I1,I2,I3,I4} = @inbounds a[1]*b[1] + a[2]*b[2]
@inline tdot(a::Tuple{I1,I3,I5}, b::Tuple{I2,I4,I6}) where {I1,I2,I3,I4,I5,I6} = @inbounds a[1]*b[1] + a[2]*b[2] + a[3]*b[3]
# @inline tdot(a::NTuple{N,Int}, b::NTuple{N,Int}) where {N} = @inbounds a[1]*b[1] + tdot(Base.tail(a), Base.tail(b))

@inline tdot(a::Tuple{I}, b::Tuple{}) where {I} = @inbounds a[1]
@inline tdot(a::Tuple{}, b::Tuple{I}) where {I} = @inbounds b[1]
@inline tdot(a::Tuple{I1,Vararg}, b::Tuple{I2}) where {I1,I2} = @inbounds a[1]*b[1]
@inline tdot(a::Tuple{I1}, b::Tuple{I2,Vararg}) where {I1,I2} = @inbounds a[1]*b[1]
@inline tdot(a::Tuple{I1,Vararg}, b::Tuple{I2,Vararg}) where {I1,I2} = @inbounds a[1]*b[1] + tdot(Base.tail(a), Base.tail(b))

"""
A wrapper to the base pointer type, that supports pointer arithmetic.
Note that `VectorizationBase.vload` and `VectorizationBase.vstore!` are 0-indexed,
while `Base.unsafe_load` and `Base.unsafe_store!` are 1-indexed.
x = [1, 2, 3, 4, 5, 6, 7, 8];
ptrx = Pointer(x);
vload(ptrx)
# 1
vload(ptrx + 1)
# 2
ptrx[]
# 1
(ptrx+1)[]
# 2
ptrx[1]
# 1
ptrx[2]
# 2
"""
abstract type AbstractPointer{T} end

@generated function gepbyte(ptr::Ptr, i::I) where {T, I <: Integer}
    ptyp = JuliaPointerType
    ityp = llvmtype(I)
    instrs = String[]
    push!(instrs, "%ptr = inttoptr $ptyp %0 to i8*")
    push!(instrs, "%offsetptr = getelementptr inbounds i8, i8* %ptr, $ityp %1")
    push!(instrs, "%iptr = ptrtoint i8* %offsetptr to $ptyp")
    push!(instrs, "ret $ptyp %iptr")
    quote
        Base.llvmcall(
            $(join(instrs, "\n")),
            Ptr{$T}, Tuple{Ptr{$T}, $I},
            ptr, i
        )
    end
end
@generated function gep(ptr::Ptr{T}, i::I) where {T, I <: Integer}
    ptyp = JuliaPointerType
    typ = llvmtype(T)
    ityp = llvmtype(I)
    instrs = String[]
    push!(instrs, "%ptr = inttoptr $ptyp %0 to $typ*")
    push!(instrs, "%offsetptr = getelementptr inbounds $typ, $typ* %ptr, $ityp %1")
    push!(instrs, "%iptr = ptrtoint $typ* %offsetptr to $ptyp")
    push!(instrs, "ret $ptyp %iptr")
    quote
        Base.llvmcall(
            $(join(instrs, "\n")),
            Ptr{$T}, Tuple{Ptr{$T}, $I},
            ptr, i
        )
    end
end
@generated function gep(ptr::Ptr{T}, i::NTuple{W,Core.VecElement{I}}) where {W, T, I <: Integer}
    ptyp = JuliaPointerType
    typ = llvmtype(T)
    ityp = llvmtype(I)
    vityp = "<$W x $ityp>"
    vptyp = "<$W x $ptyp>"
    instrs = String[]
    push!(instrs, "%ptr = inttoptr $ptyp %0 to $typ*")
    push!(instrs, "%offsetptr = getelementptr inbounds $typ, $typ* %ptr, $vityp %1")
    push!(instrs, "%iptr = ptrtoint <$W x $typ*> %offsetptr to $vptyp")
    push!(instrs, "ret $vptyp %iptr")
    quote
        Base.llvmcall(
            $(join(instrs, "\n")),
            NTuple{$W,Core.VecElement{Ptr{$T}}}, Tuple{Ptr{$T}, NTuple{W,Core.VecElement{$I}}},
            ptr, i
        )
    end
end
@inline gep(ptr::Ptr, v::SVec) = gep(ptr, extract_data(v))
@inline gep(ptr::Ptr{Cvoid}, i::Integer) where {T} = ptr.ptr + i

struct Reference{T} <: AbstractPointer{T}
    ptr::Ptr{T}
    @inline Reference(ptr::Ptr{T}) where {T} = new{T}(ptr)
end
@inline Base.eltype(::AbstractPointer{T}) where {T} = T
# @inline gep(ptr::Pointer, i::Tuple{<:Integer}) = gep(ptr, first(i))
# @inline vstore!(ptr::AbstractPointer{T1}, v::T2, args...) where {T1,T2} = vstore!(ptr, convert(T1, v), args...)

@inline vload(ptr::AbstractPointer, i::Tuple) = vload(ptr.ptr, offset(ptr, staticm1(i)))
@inline vload(ptr::AbstractPointer, i::Tuple, u::Union{AbstractMask,Unsigned}) = vload(ptr.ptr, offset(ptr, staticm1(i)), u)
@inline vstore!(ptr::AbstractPointer, v, i::Tuple) = vstore!(ptr.ptr, v, offset(ptr, staticm1(i)))
@inline vstore!(ptr::AbstractPointer, v, i::Tuple, u::Union{AbstractMask,Unsigned}) = vstore!(ptr.ptr, v, offset(ptr, staticm1(i)), u)
@inline vnoaliasstore!(ptr::AbstractPointer, v, i::Tuple) = vnoaliasstore!(ptr.ptr, v, offset(ptr, staticm1(i)))
@inline vnoaliasstore!(ptr::AbstractPointer, v, i::Tuple, u::Union{AbstractMask,Unsigned}) = vnoaliasstore!(ptr.ptr, v, offset(ptr, staticm1(i)), u)
@inline vnoaliasstore!(args...) = vstore!(args...) # generic fallback

@inline vstore!(ptr::Ptr{T}, v::Number, i::Integer) where {T <: Number} = vstore!(ptr, convert(T, v), i)
@inline vstore!(ptr::Ptr{T}, v::Integer, i::Integer) where {T <: Integer} = vstore!(ptr, v % T, i)
# @inline vstore!(ptr::AbstractPointer{T}, v::T, i...) where {T<:Integer} = vstore!(ptr.ptr, v, offset(ptr, i...))
# @inline vstore!(ptr::AbstractPointer{T1}, v::T2, i...) where {T1<:Integer,T2<:Integer} = vstore!(ptr.ptr, v % T1, offset(ptr, i...))
# @inline vstore!(ptr::AbstractPointer{T1}, v::T2, i...) where {T1<:Number,T2<:Number} = vstore!(ptr.ptr, convert(T1,v), offset(ptr, i...))
# @inline vstore!(ptr::AbstractPointer{T}, v::T, i...) where {T<:Number} = vstore!(ptr.ptr, v, offset(ptr, i...))
# @inline vstore!(ptr::AbstractPointer{T}, v::T, i...) where {T<:Integer} = vstore!(ptr.ptr, v, offset(ptr, i...))
# @inline vstore!(ptr::AbstractPointer{T1}, v::T2, i...) where {T1<:Integer,T2<:Integer} = vstore!(ptr.ptr, v % T1, offset(ptr, i...))
# @inline vstore!(ptr::AbstractPointer{T1}, v::T2, i...) where {T1<:Number,T2<:Number} = vstore!(ptr.ptr, convert(T1,v), offset(ptr, i...))

abstract type AbstractStridedPointer{T} <: AbstractPointer{T} end
abstract type AbstractColumnMajorStridedPointer{T,N} <: AbstractStridedPointer{T} end
struct PackedStridedPointer{T,N} <: AbstractColumnMajorStridedPointer{T,N}
    ptr::Ptr{T}
    strides::NTuple{N,Int}
end
struct PackedStridedBitPointer{N} <: AbstractColumnMajorStridedPointer{Bool,N}
    ptr::Ptr{UInt}
    strides::NTuple{N,Int}
end

abstract type AbstractRowMajorStridedPointer{T,N} <: AbstractStridedPointer{T} end
struct RowMajorStridedPointer{T,N} <: AbstractRowMajorStridedPointer{T,N}
    ptr::Ptr{T}
    strides::NTuple{N,Int}
end
struct RowMajorStridedBitPointer{N} <: AbstractRowMajorStridedPointer{Bool,N}
    ptr::Ptr{UInt}
    strides::NTuple{N,Int}
end

abstract type AbstractSparseStridedPointer{T,N} <: AbstractStridedPointer{T} end
struct SparseStridedPointer{T,N} <: AbstractSparseStridedPointer{T,N}
    ptr::Ptr{T}
    strides::NTuple{N,Int}
end
struct SparseStridedBitPointer{N} <: AbstractSparseStridedPointer{Bool,N}
    ptr::Ptr{UInt}
    strides::NTuple{N,Int}
end

abstract type AbstractStaticStridedPointer{T,X} <: AbstractStridedPointer{T} end
struct StaticStridedPointer{T,X} <: AbstractStaticStridedPointer{T,X}
    ptr::Ptr{T}
end
struct StaticStridedBitPointer{X} <: AbstractStaticStridedPointer{Bool,X}
    ptr::Ptr{UInt}
end
const AbstractBitPointer = Union{PackedStridedBitPointer, RowMajorStridedBitPointer, SparseStridedBitPointer, StaticStridedBitPointer}

@inline offset(::AbstractColumnMajorStridedPointer, ::Tuple{}) = 0
@inline offset(::AbstractColumnMajorStridedPointer, i::Tuple{I}) where {I} = @inbounds i[1]
@inline offset(ptr::AbstractColumnMajorStridedPointer, i::Tuple{I,Vararg}) where {I} = @inbounds i[1] + tdot(Base.tail(i), ptr.strides)

@inline offset(ptr::AbstractColumnMajorStridedPointer{T,0}, i::Tuple{I,Vararg}) where {T,I} = @inbounds i[1]
@inline offset(ptr::AbstractColumnMajorStridedPointer{T,0}, i::Tuple{I}) where {T,I} = @inbounds i[1]
@inline offset(ptr::AbstractColumnMajorStridedPointer, i::Integer) = i
@inline offset(ptr::AbstractRowMajorStridedPointer, i::Integer) = i
# @inline offset(ptr::AbstractSparseStridedPointer, i::Integer) = i * @inbounds ptr.strides[1]
# @inline offset(ptr::AbstractStaticStridedPointer{<:Any,<:Tuple{1,Vararg}}, i::Integer) = i
# @inline offset(ptr::AbstractStaticStridedPointer{<:Any,<:Tuple{M,Vararg}}, i::Integer) where {M} = M*i
@inline gep(ptr::AbstractStridedPointer, i::Tuple) = gep(ptr.ptr, offset(ptr, i))
@inline gep(ptr::AbstractStridedPointer, i::Tuple{I}) where {I} = gep(ptr.ptr, first(offset(ptr, i)))
@inline gesp(ptr::AbstractStridedPointer, i) = similar(ptr, gep(ptr, i))

@inline Base.similar(p::PackedStridedPointer, ptr::Ptr) = PackedStridedPointer(ptr, p.strides)
@inline Base.similar(p::PackedStridedBitPointer, ptr::Ptr) = PackedStridedBitPointer(ptr, p.strides)
@inline Base.similar(p::RowMajorStridedPointer, ptr::Ptr) = RowMajorStridedPointer(ptr, p.strides)
@inline Base.similar(p::RowMajorStridedBitPointer, ptr::Ptr) = RowMajorStridedBitPointer(ptr, p.strides)
@inline Base.similar(p::SparseStridedPointer, ptr::Ptr) = SparseStridedPointer(ptr, p.strides)
@inline Base.similar(p::SparseStridedBitPointer, ptr::Ptr) = SparseStridedBitPointer(ptr, p.strides)
@inline Base.similar(::StaticStridedPointer{T,X}, ptr::Ptr) where {T,X} = StaticStridedPointer{T,X}(ptr)
@inline Base.similar(::StaticStridedBitPointer{X}, ptr::Ptr) where {X} = StaticStridedBitPointer{X}(ptr)

# @inline gesp(ptr::PackedStridedPointer, i) = PackedStridedPointer(gep(ptr, i), ptr.strides)
# @inline gesp(ptr::PackedStridedBitPointer, i) = PackedStridedBitPointer(gep(ptr, i), ptr.strides)
# @inline gesp(ptr::RowMajorStridedPointer, i) = RowMajorStridedPointer(gep(ptr, i), ptr.strides)
# @inline gesp(ptr::RowMajorStridedBitPointer, i) = RowMajorStridedBitPointer(gep(ptr, i), ptr.strides)
# @inline gesp(ptr::SparseStridedPointer, i) = SparseStridedPointer(gep(ptr, i), ptr.strides)
# @inline gesp(ptr::SparseStridedBitPointer, i) = SparseStridedBitPointer(gep(ptr, i), ptr.strides)
# @inline gesp(ptr::StaticStridedPointer{T,X}, i) where {T,X} = StaticStridedPointer{T,X}(gep(ptr, i))
# @inline gesp(ptr::StaticStridedBitPointer{X}, i) where {X} = StaticStridedBitPointer{X}(gep(ptr, i))

@inline LinearAlgebra.transpose(ptr::RowMajorStridedPointer) = PackedStridedPointer(ptr.ptr, ptr.strides)
@inline offset(ptr::AbstractRowMajorStridedPointer{T,0}, i::Tuple{I}) where {T,I} = @inbounds i[1]
@inline offset(ptr::AbstractRowMajorStridedPointer{T,1}, i::Tuple{I1,I2}) where {T,I1,I2} = @inbounds i[1]*ptr.strides[1] + i[2]
@inline offset(ptr::AbstractRowMajorStridedPointer{T,2}, i::Tuple{I1,I2,I3}) where {T,I1,I2,I3} = @inbounds i[1]*ptr.strides[2] + i[2]*ptr.strides[1] + i[3]
@inline offset(ptr::AbstractRowMajorStridedPointer{T}, i::Tuple) where {T,N} = (ri = reverse(i); @inbounds ri[1] + tdot(ptr.strides, Base.tail(ri)))
# @inline function offset(ptr::AbstractRowMajorStridedPointer{Cvoid,0}, i::Tuple{Int})
    # ptr.ptr + first(i)
# end
# @inline offset(ptr::AbstractRowMajorStridedPointer{T}, i::Tuple) where {T} = offset(PackedStridedPointer(ptr.ptr, reverse(ptr.strides)), i)

@inline offset(ptr::AbstractSparseStridedPointer{T}, i::Integer) where {T} = @inbounds ptr.strides[1]*i
@inline offset(ptr::AbstractSparseStridedPointer{T}, i::Tuple) where {T} = @inbounds tdot(i, ptr.strides)
# struct ZeroInitializedStaticStridedPointer{T,X} <: AbstractStaticStridedPointer{T,X}
    # ptr::Ptr{T}
# end

@generated function LinearAlgebra.transpose(ptr::StaticStridedPointer{T,X}) where {T,X}
    tup = Expr(:curly, :Tuple)
    N = length(X.parameters)
    for n ∈ N:-1:1
        push!(tup.args, (X.parameters[n])::Int)
    end
    Expr(:block, Expr(:meta, :inline), Expr(:call, Expr(:curly, :StaticStridedPointer, T, tup), Expr(:(.), :ptr, QuoteNode(:ptr))))
end

@inline offset(ptr::AbstractStaticStridedPointer{T,<:Tuple{1,Vararg}}, i::Integer) where {T} = i
@inline offset(ptr::AbstractStaticStridedPointer{T,<:Tuple{N,Vararg}}, i::Integer) where {N,T} = i*N
@inline offset(ptr::AbstractStaticStridedPointer{T,<:Tuple{1,Vararg}}, i::Tuple{I}) where {T,I<:Integer} = first(i)
@inline offset(ptr::AbstractStaticStridedPointer{T,<:Tuple{N,Vararg}}, i::Tuple{I}) where {N,T,I<:Integer} = first(i)*N
function indprod(X::Core.SimpleVector, i)
    Xᵢ = (X[i])::Int
    iᵢ = Expr(:ref, :i, i)
    Xᵢ == 1 ? iᵢ : Expr(:call, :*, Xᵢ, iᵢ)
end
@generated function offset(ptr::AbstractStaticStridedPointer{T,X}, i::I) where {T,X,I<:Tuple}
    N = length(I.parameters)
    Xv = X.parameters
    M = min(N, length(X.parameters))
    if M == 1
        ind = indprod(Xv, 1)
    else
        ind = Expr(:call, :+)
        for m ∈ 1:M
            push!(ind.args, indprod(Xv, m))
        end
    end
    Expr(
        :block,
        Expr(:meta,:inline),
        Expr(
            :macrocall,
            Symbol("@inbounds"),
            LineNumberNode(@__LINE__, Symbol(@__FILE__)), ind
            # Expr(:call, :offset, Expr(:(.), :ptr, QuoteNode(:ptr)), ind)
        )
    )
end

struct StaticStridedStruct{T,X,S} <: AbstractStaticStridedPointer{T,X}
    ptr::S
    offset::Int # keeps track of offset, incase of nested offset calls
end

const AbstractInitializedStridedPointer{T} = Union{
    PackedStridedPointer{T},
    RowMajorStridedPointer{T},
    SparseStridedPointer{T},
    StaticStridedPointer{T},
    StaticStridedStruct{T}
}
const AbstractInitializedPointer{T} = Union{
    # Pointer{T},
    PackedStridedPointer{T},
    RowMajorStridedPointer{T},
    SparseStridedPointer{T},
    StaticStridedPointer{T},
    StaticStridedStruct{T}
}

@inline Base.stride(ptr::AbstractColumnMajorStridedPointer, i) = isone(i) ? 1 : @inbounds ptr.strides[i-1]
@inline Base.stride(ptr::AbstractSparseStridedPointer, i) = @inbounds ptr.strides[i]
@generated function Base.stride(::AbstractStaticStridedPointer{T,X}, i) where {T,X}
    Expr(:block, Expr(:meta, :inline), Expr(:getindex, Expr(:tuple, X.parameters...), :i))
end
@inline LinearAlgebra.stride1(ptr::AbstractColumnMajorStridedPointer) = 1
@inline LinearAlgebra.stride1(ptr::AbstractSparseStridedPointer) = @inbounds first(ptr.strides)
@inline LinearAlgebra.stride1(::AbstractStaticStridedPointer{T,<:Tuple{X,Vararg}}) where {T,X} = X

@inline offset(ptr::AbstractPointer, i::CartesianIndex) = offset(ptr, i.I)

@inline Base.:+(ptr::AbstractPointer{T}, i) where {T} = similar(ptr, gep(ptr.ptr, i))
@inline Base.:+(i, ptr::AbstractPointer{T}) where {T} = similar(ptr, gep(ptr.ptr, i))
@inline Base.:-(ptr::AbstractPointer{T}, i) where {T} = similar(ptr, gep(ptr.ptr, - i))

@inline vload(ptr::AbstractPointer) = vload(ptr.ptr)
@inline Base.unsafe_load(ptr::AbstractPointer) = vload(ptr.ptr)
@inline Base.getindex(ptr::AbstractPointer) = vload(ptr.ptr)

@inline vload(ptr::AbstractPointer{T}, ::Tuple{}) where {T} = vload(ptr.ptr)
# @inline vload(ptr::AbstractPointer, i::Tuple) = vload(ptr.ptr, offset(ptr, i))
# @inline vload(ptr::AbstractPointer, i::Tuple, u::Unsigned) = vload(ptr.ptr, offset(ptr, i), u)
@inline Base.unsafe_load(ptr::AbstractPointer, i) = vload(ptr.ptr, offset(ptr, i - 1))
@inline Base.getindex(ptr::AbstractPointer, i) = vload(ptr.ptr, offset(ptr, i))

@inline vstore!(ptr::AbstractPointer{T}, v::T) where {T} = vstore!(ptr.ptr, v)
@inline Base.unsafe_store!(ptr::AbstractPointer{T}, v::T) where {T} = vstore!(ptr.ptr, v)
@inline Base.setindex!(ptr::AbstractPointer{T}, v::T) where {T} = vstore!(ptr.ptr, v)

# @inline vstore!(ptr::AbstractPointer{T}, v, i::Tuple) where {T} = vstore!(ptr.ptr, v, offset(ptr, i))
# @inline vstore!(ptr::AbstractPointer{T}, v, i::Tuple, u::Unsigned) where {T} = vstore!(ptr.ptr, v, offset(ptr, i), u)
@inline Base.unsafe_store!(ptr::AbstractPointer{T}, v::T, i) where {T} = vstore!(ptr.ptr, v, offset(ptr, i - 1))
@inline Base.setindex!(ptr::AbstractPointer{T}, v::T, i) where {T} = vstore!(ptr.ptr, v, offset(ptr, i))

# @inline Pointer(A) = Pointer(pointer(A))
@inline Base.pointer(ptr::AbstractPointer) = ptr.ptr
@inline Base.unsafe_convert(::Type{Ptr{T}}, ptr::AbstractPointer{T}) where {T} = ptr.ptr

@inline stridedpointer(x) = x#Pointer(x)
@inline stridedpointer(x::Ptr) = PackedStridedPointer(x, tuple())
@inline stridedpointer(x::Union{LowerTriangular,UpperTriangular}) = stridedpointer(parent(x))
# @inline stridedpointer(x::AbstractArray) = stridedpointer(parent(x))
@inline stridedpointer(A::AbstractArray) = PackedStridedPointer(pointer(A), Base.tail(strides(A)))
@inline tailstrides(A::AbstractArray) = Base.tail(strides(A))
@inline tailstrides(A::BitArray{1}) = tuple()
@inline tailstrides(A::BitArray{2}) = (size(A,1),)
@inline tailstrides(A::BitArray{3}) = (size(A,1),size(A,1)*size(A,2))
@generated function tailstrides(A::BitArray{N}) where {N}
    quote
        (Base.Cartesian.@ntuple $(N-1) s) = size(A)
        Base.Cartesian.@nexprs $(N-2) n -> s_{n+1} *= s_n
        (Base.Cartesian.@ntuple $(N-1) s)
    end
end
@inline stridedpointer(A::BitArray) = PackedStridedBitPointer(pointer(A.chunks), tailstrides(A))
@inline stridedpointer(A::AbstractArray{T,0}) where {T} = pointer(A)
@inline stridedpointer(A::SubArray{T,0,P,S}) where {T,P,S <: Tuple{Int,Vararg}} = pointer(A)
@inline stridedpointer(A::SubArray{T,N,P,S}) where {T,N,P,S <: Tuple{Int,Vararg}} = SparseStridedPointer(pointer(A), strides(A))
@inline stridedpointer(A::SubArray{T,N,P,S}) where {T,N,P,S} = PackedStridedPointer(pointer(A), Base.tail(strides(A)))
@inline stridedpointer(B::Union{Adjoint{T,A},Transpose{T,A}}) where {T,A <: AbstractVector{T}} = stridedpointer(parent(B))
@inline stridedpointer(A::PermutedDimsArray) = SparseStridedPointer(pointer(A), strides(A))

@inline function stridedpointer(B::Union{Adjoint{T,A},Transpose{T,A}}) where {T,N,A <: AbstractArray{T,N}}
    pB = parent(B)
    RowMajorStridedPointer(pointer(pB), Base.tail(strides(pB)))
end
@inline function stridedpointer(B::Union{Adjoint{Bool,A},Transpose{Bool,A}}) where {T,N,A <: BitArray{N}}
    pB = parent(B)
    RowMajorStridedBitPointer(pointer(pB.chunks), tailstrides(pB))
end
@inline function stridedpointer(C::Union{Adjoint{T,A},Transpose{T,A}}) where {T, P, B, A <: SubArray{T,2,P,Tuple{Int,Vararg},B}}
    pC = parent(C)
    SparseStridedPointer(pointer(pC), reverse(strides(pC)))
end

@inline stridedpointer(x::Number) = x
@inline stridedpointer(x::AbstractRange) = x
# @inline stridedpointer(ptr::Pointer) = PackedStridedPointer(pointer(ptr), tuple())
@inline stridedpointer(ptr::AbstractPointer) = ptr

@generated function noalias!(ptr::Ptr{T}) where {T}
    ptyp = JuliaPointerType
    typ = llvmtype(T)
    funcname = "noalias" * typ
    decls = "define noalias $typ* @$(funcname)($typ *%a) willreturn noinline { ret $typ* %a }"
    instrs = [
        "%ptr = inttoptr $ptyp %0 to $typ*",
        "%naptr = call $typ* @$(funcname)($typ* %ptr)",
        "%jptr = ptrtoint $typ* %naptr to $ptyp",
        "ret $ptyp %jptr"
    ]
    quote
        $(Expr(:meta,:inline))
        Base.llvmcall(
            $((decls, join(instrs, "\n"))),
            Ptr{$T}, Tuple{Ptr{$T}}, ptr
        )
    end
end
@inline noaliasstridedpointer(x) = stridedpointer(x)
@inline noaliasstridedpointer(x::AbstractRange) = stridedpointer(x)
@inline function noaliasstridedpointer(x::AbstractArray)
    sptr = stridedpointer(x)
    similar(sptr, noalias!(pointer(sptr)))
end
# @inline StaticStridedStruct{T,X}(s::S) where {T,X,S} = StaticStridedStruct{T,X,S}(s, 0)
# @inline StaticStridedStruct{T,X}(s::S, i::Int) where {T,X,S} = StaticStridedStruct{T,X,S}(s, i)
# @inline offset(ptr::StaticStridedStruct{T,X,S}, i::Integer) where {T,X,S} = StaticStridedStruct{T,X,S}(ptr.ptr, ptr.offset + i)
# # Trying to avoid generated functions
# @inline offset(ptr::StaticStridedStruct{T,X,S}, i::Tuple{<:Integer}) where {T,X,S} = @inbounds StaticStridedStruct{T,X,S}(ptr.ptr, ptr.offset + i[1])
# @inline function offset(ptr::StaticStridedStruct{T,Tuple{A},S}, i::Tuple{<:Integer,<:Integer}) where {T,A,S}
#     @inbounds StaticStridedStruct{T,Tuple{A},S}(ptr.ptr, ptr.offset + i[1] + A * i[2])
# end
# @inline function offset(ptr::StaticStridedStruct{T,Tuple{A,B},S}, i::Tuple{<:Integer,<:Integer,<:Integer}) where {T,A,B,S}
#     @inbounds StaticStridedStruct{T,Tuple{A,B},S}(ptr.ptr, ptr.offset + i[1] + A*i[2] + B*i[3])
# end
# @inline function offset(ptr::StaticStridedStruct{T,Tuple{A,B,C},S}, i::Tuple{<:Integer,<:Integer,<:Integer,<:Integer}) where {T,A,B,C,S}
#     @inbounds StaticStridedStruct{T,Tuple{A,B,C},S}(ptr.ptr, ptr.offset + i[1] + A*i[2] + B*i[3] + C*i[4] )
# end

# @generated tupletype_to_tuple(::Type{T}) where {T<:Tuple} = Expr(:block, Expr(:meta,:inline), Expr(:tuple, T.parameters...))
# @inline function offset(ptr::StaticStridedStruct{T,X,S}, i::NTuple{N}) where {T,X,S,N}
#     strides = tupletype_to_tuple(X)
#     StaticStridedStruct{T,X,S}(ptr.ptr, ptr.offset + first(i) + tdot(strides, Base.tail(i)))
# end

@inline vload(r::AbstractRange, i::Tuple{<:Integer}) = @inbounds r[i[1]]

@inline subsetview(ptr::PackedStridedPointer, ::Val{1}, i::Integer) = SparseStridedPointer(gep(ptr.ptr, i), ptr.strides)
@generated function subsetview(ptr::PackedStridedPointer{T, N}, ::Val{I}, i::Integer) where {I, T, N}
    I > N + 1 && return :ptr
    strides = Expr(:tuple, [Expr(:ref, :s, n) for n ∈ 1:N if n != I-1]...)
    offset = Expr(:call, :*, :i, Expr(:macrocall, Symbol("@inbounds"), LineNumberNode(@__LINE__, Symbol(@__FILE__)), Expr(:ref, :s, I - 1)))
    Expr(
        :block,
        Expr(:meta, :inline),
        Expr(:(=), :p, Expr(:(.), :ptr, QuoteNode(:ptr))),
        Expr(:(=), :s, Expr(:(.), :ptr, QuoteNode(:strides))),
        Expr(:(=), :strides, Expr(:macrocall, Symbol("@inbounds"), LineNumberNode(@__LINE__, Symbol(@__FILE__)), strides)),
        Expr(:(=), :gp, Expr(:call, :gep, :p, offset)),
        Expr(:call, :PackedStridedPointer, :gp, :strides)
    )
end

@generated function subsetview(ptr::RowMajorStridedPointer{T, N}, ::Val{I}, i::Integer) where {I, T, N}
    if N + 1 == I
        return Expr(
            :block, Expr(:meta, :inline),
            Expr(:(=), :p, Expr(:call, :gep, Expr(:(.), :ptr, QuoteNode(:ptr)), :i)),
            Expr(:call, :SparseStridedPointer, :p, Expr(:call, :reverse, Expr(:(.), :ptr, QuoteNode(:strides))))
        )
    elseif N + 1 < I
        return :ptr
    end
    strideind = N + 1 - I
    strides = Expr(:tuple, [Expr(:ref, :s, n) for n ∈ 1:N if n != strideind]...)
    offset = Expr(:call, :*, :i, Expr(:macrocall, Symbol("@inbounds"), LineNumberNode(@__LINE__, Symbol(@__FILE__)), Expr(:ref, :s, strideind)))
    Expr(
        :block,
        Expr(:meta, :inline),
        Expr(:(=), :p, Expr(:(.), :ptr, QuoteNode(:ptr))),
        Expr(:(=), :s, Expr(:(.), :ptr, QuoteNode(:strides))),
        Expr(:(=), :strides, Expr(:macrocall, Symbol("@inbounds"), LineNumberNode(@__LINE__, Symbol(@__FILE__)), strides)),
        Expr(:(=), :gp, Expr(:call, :gep, :p, offset)),
        Expr(:call, :RowMajorStridedPointer, :gp, :strides)
    )
end

@generated function subsetview(ptr::SparseStridedPointer{T, N}, ::Val{I}, i::Integer) where {I, T, N}
    I > N && return :ptr
    strides = Expr(:tuple, [Expr(:ref, :s, n) for n ∈ 1:N if n != I]...)
    offset = Expr(:call, :*, :i, Expr(:macrocall, Symbol("@inbounds"), LineNumberNode(@__LINE__, Symbol(@__FILE__)), Expr(:ref, :s, I)))
    Expr(
        :block,
        Expr(:meta, :inline),
        Expr(:(=), :p, Expr(:(.), :ptr, QuoteNode(:ptr))),
        Expr(:(=), :s, Expr(:(.), :ptr, QuoteNode(:strides))),
        Expr(:(=), :strides, Expr(:macrocall, Symbol("@inbounds"), LineNumberNode(@__LINE__, Symbol(@__FILE__)), strides)),
        Expr(:(=), :gp, Expr(:call, :gep, :p, offset)),
        Expr(:call, :SparseStridedPointer, :gp, :strides)
    )
end

@generated function subsetview(ptr::StaticStridedPointer{T, X}, ::Val{I}, i::Integer) where {I, T, X}
    I > length(X.parameters) && return :ptr
    Xa = Expr(:curly, :Tuple)
    Xparam = X.parameters
    for n ∈ 1:length(Xparam)
        n == I && continue
        push!(Xa.args, Xparam[n])
    end
    offset = Expr(:call, :*, :i, Xparam[I])
    Expr(
        :block,
        Expr(:meta, :inline),
        Expr(:(=), :p, Expr(:(.), :ptr, QuoteNode(:ptr))),
        Expr(:(=), :gp, Expr(:call, :gep, :p, offset)),
        Expr(:call, Expr(:curly, :StaticStridedPointer, T, Xa), :gp)
    )
end

@generated function subsetview(ptr::StaticStridedStruct{T, X}, ::Val{I}, i::Integer) where {I, T, X}
    I > length(X.parameters) && return :ptr
    Xa = Expr(:curly, :Tuple)
    Xparam = X.parameters
    for n ∈ 1:length(Xparam)
        n == I && continue
        push!(Xa.args, Xparam[n])
    end
    offset = Expr(:call, :*, :i, Xparam[I])
    Expr(
        :block,
        Expr(:meta, :inline),
        Expr(:(=), :p, Expr(:(.), :ptr, QuoteNode(:ptr))),
        Expr(:(=), :offset, Expr(:call, :+, Expr(:(.), :ptr, QuoteNote(:offset)), offset)),
        Expr(:call, Expr(:curly, :StaticStridedStruct, T, Xa), :p, :offset)
    )
end

@inline stridedpointer(A::AbstractArray, indices::Tuple) = stridedpointer(view(A, indices...))
@inline stridedpointer(A::AbstractArray, ::Type{Transpose}) = stridedpointer(transpose(A))
@inline stridedpointer(A::AbstractArray, ::Type{Adjoint}) = stridedpointer(adjoint(A))
@inline stridedpointer(A::AbstractArray, ::Nothing) = stridedpointer(A)

@inline filter_strides_by_dimequal1(sz::NTuple{N,Int}, st::NTuple{N,Int}) where {N} = @inbounds ntuple(n -> isone(sz[n]) ? 0 : st[n], Val{N}())

@inline stridedpointer_for_broadcast(A::AbstractRange) = A
@inline function stridedpointer_for_broadcast(A::AbstractArray{T,N}) where {T,N}
    PackedStridedPointer(pointer(A), filter_strides_by_dimequal1(Base.tail(size(A)), Base.tail(strides(A))))
end
@inline stridedpointer_for_broadcast(B::Union{Adjoint{T,A},Transpose{T,A}}) where {T,A <: AbstractVector{T}} = stridedpointer_for_broadcast(parent(B))
@inline stridedpointer_for_broadcast(A::SubArray{T,0,P,S}) where {T,P,S <: Tuple{Int,Vararg}} = pointer(A)
@inline function stridedpointer_for_broadcast(A::SubArray{T,N,P,S}) where {T,N,P,S <: Tuple{Int,Vararg}}
    SparseStridedPointer(pointer(A), filter_strides_by_dimequal1(size(A), strides(A)))
end
@inline function stridedpointer_for_broadcast(A::SubArray{T,N,P,S}) where {T,N,P,S}
    PackedStridedPointer(pointer(A), filter_strides_by_dimequal1(Base.tail(size(A)), Base.tail(strides(A))))
end
@inline stridedpointer_for_broadcast(A::BitArray) = stridedpointer(A)

struct MappedStridedPointer{F, T, P <: AbstractPointer{T}}
    f::F
    ptr::P
end
@inline vload(ptr::MappedStridedPointer) = ptr.f(vload(ptr.ptr))


text:
" Coloring scheme for the plot, in RGBA "
const colorMap = [
  sfColor(0,0,0,255),
  sfColor(66,46,15,255),
  sfColor(26,8,26,255),
  sfColor(10,0,46,255),
  sfColor(5,5,74,255),
  sfColor(0,8,99,255),
  sfColor(13,43,138,255),
  sfColor(56,125,209,255),
  sfColor(133,181,230,255),
  sfColor(209,235,247,255),
  sfColor(240,232,191,255),
  sfColor(247,201,94,255),
  sfColor(255,171,0,255),
  sfColor(204,128,0,255),
  sfColor(153,87,0,255),
  sfColor(105,51,3,255)
];

"""Linearly interpolate between two colors"""
function lerp(colorA, colorB, scale)
  r = UInt8(floor((1 - scale) * colorA.r + scale * colorB.r))
  g = UInt8(floor((1 - scale) * colorA.g + scale * colorB.g))
  b = UInt8(floor((1 - scale) * colorA.b + scale * colorB.b))
  a = UInt8(floor((1 - scale) * colorA.a + scale * colorB.a))

  return sfColor(r, g, b, a)
end

"""
getColor(val::Complex{Float64}, n::Integer)

Each color is colored exponentially according to its value after the simulation, and the number of steps needed for it to exceed that max value.

The color is chosen according to a logarithmic scale and linearly interpolated with the next closest color to provide a smooth transition between colors.
"""
function getColor(val, n, nmax)
  v = abs(val)
  if v < 2
    n_log = 0
  else
    n_log = log2(log2(v))
  end

  n = n - n_log

  n = (n / nmax) * (length(colorMap) - 1)

  colorA = colorMap[Int(floor(n)) + 1]
  colorB = colorMap[Int(floor(n + 1)) % length(colorMap) + 1]

  return lerp(colorA, colorB, n % 1)
end

max_stars_repo_name: benchislett/Skunkworks

max_issues_repo_name: benchislett/Skunkworks

max_forks_repo_name: benchislett/Skunkworks


text:
using Statistics
using StatsFuns: xlogx, xlogy

"Cache pairwise / marginal distribution for all variables in one dataset"
mutable struct DisCache
    pairwise::Array{Float64}
    marginal::Array{Float64}
end

@inline dimension(discache::DisCache) = size(discache.marginal)[1]
DisCache(num) = DisCache(Array{Float64}(undef, num, num, 4), Array{Float64}(undef, num, 2))

#####################
# Methods for pairwise and marginal distribution
#####################
@inline get_parameters(bm::AbstractMatrix{<:Bool}, α, w=nothing) = size(bm)[2], issomething(w) ? sum(w) : size(bm)[1], @. Float64(bm), @. Float64(!bm)

function cache_distributions(bm::AbstractMatrix{<:Bool}, w::Union{Nothing, AbstractVector{<:AbstractFloat}}=nothing; α, flag=(pairwise=true, marginal=true))
    # parameters
    D, N, (m, notm) = get_parameters(bm, α, w)
    dis_cache = DisCache(D)
    base = N + 4 * α
    w = isnothing(w) ? ones(Float64, N) : w

    # pairwise distribution
    if flag.pairwise
        dis_cache.pairwise[:,:,1] = (notm' * (notm .* w) .+ α) / base   # p00
        dis_cache.pairwise[:,:,2] = (notm' * (m .* w) .+ α) / base      # p01
        dis_cache.pairwise[:,:,3] = (m' * (notm .* w) .+ α) / base      # p10
        dis_cache.pairwise[:,:,4] = (m' * (m .* w) .+ α) / base         # p11
    end
    # marginal distribution

    if flag.marginal
        dis_cache.marginal[:, 1] = (sum(notm .* w, dims=1) .+ 2 * α) / base
        dis_cache.marginal[:, 2] = (sum(m .* w, dims=1).+ 2 * α) / base
    end
    dis_cache
end

#####################
# Mutual Information
#####################

function mutual_information(dis_cache::DisCache)
    D = dimension(dis_cache)
    p0 = @view dis_cache.marginal[:, 1]
    p1 = @view dis_cache.marginal[:, 2]
    pxpy = Array{Float64}(undef, D, D, 4)
    pxpy[:,:,1] = p0 * p0'
    pxpy[:,:,2] = p0 * p1'
    pxpy[:,:,3] = p1 * p0'
    pxpy[:,:,4] = p1 * p1'
    pxy_log_pxy = @. xlogx(dis_cache.pairwise)
    pxy_log_pxpy = @. xlogy(dis_cache.pairwise, pxpy)
    dropdims(sum(pxy_log_pxy - pxy_log_pxpy,dims=3), dims=3)
end

"Calculate mutual information of given bit matrix `bm`, example weights `w`, and smoothing pseudocount `α`"
function mutual_information(bm::AbstractMatrix{<:Bool}, w::Union{Nothing, AbstractVector{<:AbstractFloat}}=nothing; α)
    dis_cache = cache_distributions(bm, w; α=α)
    mi = mutual_information(dis_cache)
    return (dis_cache, mi)
end

function mutual_information(train_x::PlainXData, w::Union{Nothing, AbstractVector{<:AbstractFloat}}=nothing; α)
    mutual_information(feature_matrix(train_x), w; α=α)
end

"Calculate set mutual information"
function set_mutual_information(mi::Matrix, sets::Vector{Vector{Var}})::Matrix
    len = length(sets)
    if len == size(mi)[1]
        return mi
    end

    pmi = zeros(len, len)
    for i in 1 : len, j in i + 1 : len
        pmi[i, j] = pmi[j, i] = mean(mi[sets[i], sets[j]])
    end
    return pmi
end

#####################
# Entropy
#####################
function entropy(dis_cache::DisCache)
    D = dimension(dis_cache)
    px_log_px = @. xlogx(dis_cache.marginal)
    - dropdims(sum(px_log_px; dims=2); dims=2)
end

function entropy(bm::AbstractMatrix{<:Bool}, w::Union{Nothing, AbstractVector{<:AbstractFloat}}=nothing; α)
    dis_cache = cache_distributions(bm, w; α=α, flag=(pairwise=true, marginal=true))
    return (dis_cache, entropy(dis_cache))
end

function sum_entropy_given_x(bm::AbstractMatrix{<:Bool}, x::Var, w::Union{Nothing, AbstractVector{<:AbstractFloat}}=nothing; α)::Float64
    @assert x <= size(bm)[2]
    vars = [1 : x-1; x+1 : size(bm)[2]]
    indexes_left = bm[:,x].== 0
    indexes_right = bm[:,x] .== 1
    w1 = sum(Float64.(indexes_left))
    w2 = sum(Float64.(indexes_right))
    w1, w2 = w1 / (w1 + w2), w2 / (w1 + w2)
    subm_left = @view bm[indexes_left, vars]
    subm_right = @view bm[indexes_right, vars]

    w_left = issomething(w) ? (@view w[indexes_left]) : nothing
    w_right = issomething(w) ? (@view w[indexes_right]) : nothing

    w1 * sum(entropy(subm_left, w_left; α=α)[2]) + w2 * sum(entropy(subm_right, w_right; α=α)[2])
end

function conditional_entropy(dis_cache::DisCache)
    D = dimension(dis_cache)
    pxy_log_pxy = @. xlogx(dis_cache.pairwise)
    pxy = dis_cache.pairwise
    log_px = log.(dis_cache.marginal)

    pxy_log_px = Array{Float64}(undef, D, D, 4)
    pxy_log_px[:, :, 1] = view(pxy, :, :, 1) .* view(log_px, :, 1)
    pxy_log_px[:, :, 2] = view(pxy, :, :, 2) .* view(log_px, :, 1)
    pxy_log_px[:, :, 3] = view(pxy, :, :, 3) .* view(log_px, :, 2)
    pxy_log_px[:, :, 4] = view(pxy, :, :, 4) .* view(log_px, :, 2)
    h_y_given_x = - dropdims(sum(pxy_log_pxy - pxy_log_px; dims=3); dims=3)
end

function conditional_entropy(bm::AbstractMatrix{<:Bool}, w::Union{Nothing, AbstractVector{<:AbstractFloat}}=nothing; α)
    dis_cache = cache_distributions(bm, w; α=α)
    return (dis_cache, conditional_entropy(dis_cache))
end

function set_conditional_entropy(ce::Matrix{<:Float64})
    dropdims(sum(ce; dims=2); dims=2)
end


text:
# Functions for querying network properties.

######### Accessors: #########

"""
    species(network)

Given a [`ReactionSystem`](@ref), return a vector of species `Variable`s.

Notes:
- If `network.systems` is not empty, may allocate. Otherwise returns
  `network.states`.
"""
function species(network)
    isempty(network.systems) ? network.states : states(network)
end

"""
    params(network)

Given a [`ReactionSystem`](@ref), return a vector of parameter `Variable`s.

Notes:
- If `network.systems` is not empty, may allocate. Otherwise returns
  `network.ps`.
"""
function params(network)
    isempty(network.systems) ? network.ps : parameters(network)
end

"""
    reactions(network)

Given a [`ReactionSystem`](@ref), return a vector of all `Reactions` in the system.

Notes:
- If `network.systems` is not empty, may allocate. Otherwise returns
  `network.eqs`.
"""
function reactions(network)
    isempty(network.systems) ? network.eqs : equations(network)
end

"""
    speciesmap(network)

Given a [`ReactionSystem`](@ref), return a Dictionary mapping from species
`Variable`s to species indices. (Allocates)
"""
function speciesmap(network)
    Dict(S => i for (i,S) in enumerate(species(network)))
end

"""
    paramsmap(network)

Given a [`ReactionSystem`](@ref), return a Dictionary mapping from parameter
`Variable`s to parameter indices. (Allocates)
"""
function paramsmap(network)
    Dict(p => i for (i,p) in enumerate(params(network)))
end

"""
    numspecies(network)

Return the number of species within the given [`ReactionSystem`](@ref).
"""
function numspecies(network)
    ns = length(network.states)
    for sys in network.systems
        ns += numspecies(ns)
    end
    ns
end

"""
    numreactions(network)

Return the number of reactions within the given [`ReactionSystem`](@ref).
"""
function numreactions(network)
    nr = length(network.eqs)
    for sys in network.systems
        nr += numreactions(sys)
    end
    nr
end

"""
    numparams(network)

Return the number of parameters within the given [`ReactionSystem`](@ref).
"""
function numparams(network)
    np = length(network.ps)
    for sys in network.systems
        np += numparams(ns)
    end
    np
end

"""
    dependents(rx, network)

Given a [`Reaction`](@ref) and a [`ReactionSystem`](@ref), return a vector of
`ModelingToolkit.Operation`s corresponding to species the *reaction rate
law* depends on. E.g., for

`k*W, 2X + 3Y --> 5Z + W`

the returned vector would be `[W(t),X(t),Y(t)]`.

Notes:
- Allocates
- Does not check for dependents within any subsystems.
"""
function dependents(rx, network)
    if rx.rate isa Operation
        rvars = ModelingToolkit.get_variables(rx.rate, species(network))
        return union!(rvars, rx.substrates)
    end

    rx.substrates
end

"""
    dependents(rx, network)

See documentation for [`dependents`](@ref).
"""
function dependants(network, rxidx)
    dependents(network, rxidx)
end

######################## reaction network operators #######################

"""
    ==(rn1::Reaction, rn2::Reaction)

Tests whether two [`Reaction`](@ref)s are identical.

Notes:
- Ignores the order in which stoichiometry components are listed.
- *Does not* currently simplify rates, so a rate of `A^2+2*A+1` would be
    considered different than `(A+1)^2`.
"""
function (==)(rn1::Reaction, rn2::Reaction)
    isequal(rn1.rate, rn2.rate) || return false
    opit = sv -> (s.op for s in sv)
    issetequal(zip(opit(rn1.substrates),rn1.substoich), zip(opit(rn2.substrates),rn2.substoich)) || return false
    issetequal(zip(opit(rn1.products),rn1.prodstoich), zip(opit(rn2.products),rn2.prodstoich)) || return false
    issetequal(rn1.netstoich, rn2.netstoich)
end

"""
    ==(rn1::ReactionSystem, rn2::ReactionSystem)

Tests whether the underlying species `Variables`s, parameter `Variables`s
and reactions are the same in the two [`ReactionSystem`](@ref)s. Ignores
order network components were defined.

Notes:
- *Does not* currently simplify rates, so a rate of `A^2+2*A+1` would be
    considered different than `(A+1)^2`.
- Flattens subsystems, and hence may allocate, when checking equality.
"""
function (==)(rn1::ReactionSystem, rn2::ReactionSystem)
    issetequal(species(rn1), species(rn2)) || return false
    issetequal(params(rn1), params(rn2)) || return false
    isequal(rn1.iv, rn2.iv) || return false
    (numreactions(rn1) == numreactions(rn2)) || return false

    # the following fails for some reason, so need to use issubset
    #issetequal(equations(rn1), equations(rn2)) || return false
    (issubset(reactions(rn1),reactions(rn2)) && issubset(reactions(rn2),reactions(rn1))) || return false

    # BELOW SHOULD NOT BE NEEDED as species, params and equations flatten
    #issetequal(rn1.systems, rn2.systems) || return false
    # sys1 = rn1.systems; sys2 = rn2.systems
    # (issubset(sys1,sys2) && issubset(sys2,sys1)) || return false
    true
end

######################## functions to extend a network ####################

"""
    make_empty_network(; iv=Variable(:t))

Construct an empty [`ReactionSystem`](@ref). `iv` is the independent variable, usually time.
"""
function make_empty_network(; iv=Variable(:t))
    ReactionSystem(Reaction[], iv, Operation[], Operation[], Variable[], Equation[], gensym(:ReactionSystem), ReactionSystem[])
end

"""
    addspecies!(network::ReactionSystem, s::Variable; disablechecks=false)

Given a [`ReactionSystem`](@ref), add the species corresponding to the variable
`s` to the network (if it is not already defined). Returns the integer id of the
species within the system.

Notes:
- `disablechecks` will disable checking for whether the passed in variable is
  already defined, which is useful when adding many new variables to the system.
  *Do not disable checks* unless you are sure the passed in variable is a new
  variable, as this will potentially leave the system in an undefined state.
"""
function addspecies!(network::ReactionSystem, s::Variable; disablechecks=false)

    # we don't check subsystems since we will add it to the top-level system...
    curidx = disablechecks ? nothing : findfirst(S -> isequal(S, s), network.states)
    if curidx === nothing
        push!(network.states, s)
        return length(network.states)
    else
        return curidx
    end
end

"""
    addspecies!(network::ReactionSystem, s::Operation; disablechecks=false)

Given a [`ReactionSystem`](@ref), add the species corresponding to the
variable `s` to the network (if it is not already defined). Returns the
integer id of the species within the system.

- `disablechecks` will disable checking for whether the passed in variable is
  already defined, which is useful when adding many new variables to the system.
  *Do not disable checks* unless you are sure the passed in variable is a new
  variable, as this will potentially leave the system in an undefined state.
"""
function addspecies!(network::ReactionSystem, s::Operation; disablechecks=false)
    !(s.op isa Variable) && error("If the passed in species is an Operation, it must correspond to an underlying Variable.")
    addspecies!(network, convert(Variable,s); disablechecks=disablechecks)
end

"""
    addparam!(network::ReactionSystem, p::Variable; disablechecks=false)

Given a [`ReactionSystem`](@ref), add the parameter corresponding to the
variable `p` to the network (if it is not already defined). Returns the integer
id of the parameter within the system.

- `disablechecks` will disable checking for whether the passed in variable is
  already defined, which is useful when adding many new variables to the system.
  *Do not disable checks* unless you are sure the passed in variable is a new
  variable, as this will potentially leave the system in an undefined state.
"""
function addparam!(network::ReactionSystem, p::Variable; disablechecks=false)

    # we don't check subsystems since we will add it to the top-level system...
    curidx = disablechecks ? nothing : findfirst(S -> isequal(S, p), network.ps)
    if curidx === nothing
        push!(network.ps, p)
        return length(network.ps)
    else
        return curidx
    end
end

"""
    addparam!(network::ReactionSystem, p::Operation; disablechecks=false)

Given a [`ReactionSystem`](@ref), add the parameter corresponding to the
variable `p` to the network (if it is not already defined). Returns the
integer id of the parameter within the system.

- `disablechecks` will disable checking for whether the passed in variable is
  already defined, which is useful when adding many new variables to the system.
  *Do not disable checks* unless you are sure the passed in variable is a new
  variable, as this will potentially leave the system in an undefined state.
"""
function addparam!(network::ReactionSystem, p::Operation; disablechecks=false)
    !(p.op isa Variable) && error("If the passed in parameter is an Operation, it must correspond to an underlying Variable.")
    addparam!(network, convert(Variable,p); disablechecks=disablechecks)
end

"""
    addreaction!(network::ReactionSystem, rx::Reaction)

Add the passed in reaction to the [`ReactionSystem`](@ref). Returns the
integer id of `rx` in the list of `Reaction`s within `network`.

Notes:
- Any new species or parameters used in `rx` should be separately added to
    `network` using [`addspecies!`](@ref) and [`addparam!`](@ref).
"""
function addreaction!(network::ReactionSystem, rx::Reaction)
    push!(network.eqs, rx)
    length(network.eqs)
end

"""
    merge!(network1::ReactionSystem, network2::ReactionSystem)

Merge `network2` into `network1`.

Notes:
- Duplicate reactions between the two networks are not filtered out.
- [`Reaction`](@ref)s are not deepcopied to minimize allocations, so both networks will share underlying data arrays.
- Subsystems are not deepcopied between the two networks and will hence be shared.
- Returns `network1`.
"""
function merge!(network1::ReactionSystem, network2::ReactionSystem)
    isequal(network1.iv, network2.iv) || error("Reaction networks must have the same independent variable to be mergable.")
    specs = network1.states
    foreach(spec -> !(spec in specs) && push!(specs, spec), network2.states)
    ps = network1.ps
    foreach(p -> !(p in ps) && push!(ps, p), network2.ps)
    append!(network1.eqs, network2.eqs)
    append!(network1.systems, network2.systems)
    network1
end

"""
    merge(network1::ReactionSystem, network2::ReactionSystem)

Create a new network merging `network1` and `network2`.

Notes:
- Duplicate reactions between the two networks are not filtered out.
- [`Reaction`](@ref)s are not deepcopied to minimize allocations, so the new network will share underlying data arrays.
- Subsystems are not deepcopied between the two networks and will hence be shared.
- Returns the merged network.
"""
function merge(network1::ReactionSystem, network2::ReactionSystem)
    network = make_empty_network()
    merge!(network, network1)
    merge!(network, network2)
    network
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule librdkafka_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("librdkafka")
JLLWrappers.@generate_main_file("librdkafka", UUID("7943bfb0-7437-5acd-a008-22777931c7aa"))
end  # module librdkafka_jll


text:
module RandomLinearAlgebraSolvers

using LinearAlgebra, LLSModels, SparseArrays, Stopping
using StatsBase # for sample!

@inline krylov_dot(
  n::Integer,
  x::Vector{T},
  dx::Integer,
  y::Vector{T},
  dy::Integer,
) where {T <: BLAS.BlasReal} = BLAS.dot(n, x, dx, y, dy)
@inline krylov_dot(
  n::Integer,
  x::Vector{T},
  dx::Integer,
  y::Vector{T},
  dy::Integer,
) where {T <: BLAS.BlasComplex} = BLAS.dotc(n, x, dx, y, dy)
@inline krylov_dot(
  n::Integer,
  x::AbstractVector{T},
  dx::Integer,
  y::AbstractVector{T},
  dy::Integer,
) where {T <: Number} = dot(x, y)

macro kdot(n, x, y)
  return esc(:(krylov_dot($n, $x, 1, $y, 1)))
end

function get_matrix(lls::LLSModel)
  return sparse(lls.Arows, lls.Acols, lls.Avals, length(lls.b), lls.meta.nvar)
end

function get_vector(lls::LLSModel)
  return lls.b
end

function get_matrix(pb::LinearSystem)
  return pb.A
end

function get_vector(pb::LinearSystem)
  return pb.b
end

export RLAStopping

"""
    stp = RLAStopping(A, b::S; n_listofstates::Int = 0, kwargs...)

Creator of a `LAStopping` tailored for this package.
The problem, stp.pb, is a `Stopping.LinearSystem` if `A` is dense, and a `LLSModels.LLSModel` otherwise.
The state allocates space for the residual, `stp.current_state.res`, of length `|b|`.
This stopping uses the infinity norm of `stp.current_state.res` to declare optimality.
"""
function RLAStopping(A, b::S; n_listofstates::Int = 0, kwargs...) where {S}
  T = eltype(S)
  sparse = issparse(A)
  pb = sparse ? LLSModel(A, b) : LinearSystem(A, b)
  x0 = zeros(T, size(A, 2))
  state = GenericState(similar(x0), similar(b), res = similar(b))

  mcntrs = sparse ? init_max_counters_NLS() : init_max_counters_linear_operators()

  if n_listofstates > 0 && :list ∉ keys(kwargs)
    list = ListofStates(n_listofstates, Val{typeof(state)}())
    return LAStopping(
      pb,
      state,
      max_cntrs = mcntrs,
      list = list,
      optimality_check = (pb, state) -> state.res;
      kwargs...,
    )
  end

  return LAStopping(
    pb,
    state,
    max_cntrs = mcntrs,
    optimality_check = (pb, state) -> state.res;
    kwargs...,
  )
end

include("random_Kaczmarcz.jl")
include("random_coordinate_descent.jl")
include("random_projector.jl")
include("random_Newton.jl")
include("random_methods.jl")

end