text:
"""
    function growGeneratorProductions(gpExperimentInfo::GEInfo, rng::Random.AbstractRNG,
                    initInd::Array{UInt8}=Array{UInt8}(undef, 0))::GEGenotype

Generates a Grammatical Evolution tree ensuring that the individual is valid having
productions in GEInfo as limit.

# Self-provided Arguments
- `gpExperimentInfo::GEInfo`: information about grammatical evolution.
- `rng::Random.AbstractRNG`: random number generator for random number consistency along an experiment.

!!! note
    `Selfprovided Arguments` are already provided by the library, so only `User Arguments` must be provided.
"""
function growGeneratorProductions(gpExperimentInfo::GEInfo, rng::Random.AbstractRNG,
                initInd::Array{UInt8}=Array{UInt8}(undef, 0))::GEGenotype

    productions::UInt16 = gpExperimentInfo._maxProductions
    gramm = gpExperimentInfo._grammar

    ind = initInd
    if isempty(ind)
        ind = rand(rng, UInt8, productions)
    end

    Root = GERep(getInitialSymbol(gramm))
    stack = [Root]

    # A production is always made due to initial symbol
    actualProductions =  productions - 0x2

    actualDepth::UInt8 = 0x0
    symbolsPerDepth = zeros(UInt8, 2)
    symbolsPerDepth[1] = 0x1

    selections = UInt8[]
    i = 1
    while !isempty(stack)
        if symbolsPerDepth[actualDepth + 0x1] == 0x0
            actualDepth += 0x1
        end

        symbolsPerDepth[actualDepth + 0x1] -= 0x1

        node = popfirst!(stack)
        symbol = getSymbol(node)
        if !isTerminalSymbol(gramm, symbol)

            actualProductions = actualProductions + getDepth(gramm, symbol) - 0x1

            indexes, nIndexes = getPosibleConsequentProductions(gramm, symbol, actualProductions)

            selected = ind[i] % nIndexes + 0x1
            selected = indexes[selected]

            nChoices = getNConsequent(gramm, symbol)
            ind[i] = ind[i] % (div(256,nChoices)) * nChoices + selected - 0x1

            push!(selections, selected)

            consequent = getConsequent(gramm, symbol, selected)
            if length(symbolsPerDepth) < actualDepth + 0x2
                push!(symbolsPerDepth, getNSymbols(consequent))
            else
                symbolsPerDepth[actualDepth + 0x2] += getNSymbols(consequent)
            end
            actualProductions = actualProductions - getMinProductions(consequent)
            symbols = consequent._symbols
            for j=length(symbols):-1:1
                child = GERep(symbols[j], actualDepth + 0x1)
                pushChild!(node, child)
                push!(stack, child)
            end
            i = i + 1
        end
    end

    return GEGenotype(Root, IntegerGenotype{UInt8}(ind), selections, i-1, true)
end # function

"""
    function growGeneratorDepth(gpExperimentInfo::GEInfo, rng::Random.AbstractRNG,
                    initInd::Array{UInt8}=Array{UInt8}(undef, 0))::GEGenotype

Generates a Grammatical Evolution tree ensuring that the individual is valid having
productions and depth in GEInfo as limit.

# Self-provided Arguments
- `gpExperimentInfo::GEInfo`: information about grammatical evolution.
- `rng::Random.AbstractRNG`: random number generator for random number consistency along an experiment.

!!! note
    `Selfprovided Arguments` are already provided by the library, so only `User Arguments` must be provided.
"""
function growGeneratorDepth(gpExperimentInfo::GEInfo, rng::Random.AbstractRNG,
                        initInd::Array{UInt8}=Array{UInt8}(undef, 0))::GEGenotype

    maxDepth::UInt8 = gpExperimentInfo._maxDepth
    maxProductions::UInt16 = gpExperimentInfo._maxProductions
    gramm = gpExperimentInfo._grammar

    ind = initInd

    if isempty(ind)
        ind = rand(rng, UInt8, maxProductions)
    end

    Root = GERep(getInitialSymbol(gramm))
    stack = [Root]

    actualProductions =  maxProductions - 0x2

    actualDepth::UInt8 = 0x0
    symbolsPerDepth = zeros(UInt8, maxDepth+1)
    symbolsPerDepth[1] = 0x1
    selections = UInt8[]
    i = 1
    while !isempty(stack)
        if symbolsPerDepth[actualDepth + 0x1] == 0x0
            actualDepth += 0x1
        end

        symbolsPerDepth[actualDepth + 0x1] -= 0x1
        node = popfirst!(stack)
        symbol = getSymbol(node)

        if !isTerminalSymbol(gramm, symbol)
            actualProductions = actualProductions + getDepth(gramm, symbol) - 0x1
            indexes, nIndexes = getPosibleConsequentDepth(gramm, symbol, maxDepth-actualDepth, actualProductions)

            selected = ind[i] % nIndexes + 0x1
            selected = indexes[selected]

            nChoices = getNConsequent(gramm, symbol)
            ind[i] = ind[i] % (div(256,nChoices)) * nChoices + selected - 0x1

            push!(selections, selected)

            consequent = getConsequent(gramm, symbol, selected)
            symbolsPerDepth[actualDepth + 0x2] += getNSymbols(consequent)
            actualProductions = actualProductions - getMinProductions(consequent)
            symbols = consequent._symbols
            for j=length(symbols):-1:1
                child = GERep(symbols[j], actualDepth+0x1)
                pushChild!(node, child)
                push!(stack, child)
            end
            i = i + 1
        end
    end

    return GEGenotype(Root, IntegerGenotype(ind), selections, i-1, true)
end # function

"""
    function growGenerator(gpExperimentInfo::GEInfo, rng::Random.AbstractRNG,
                    initInd::Array{UInt8}=Array{UInt8}(undef, 0))::GEGenotype

Generates a Grammatical Evolution tree ensuring that the individual is valid having
productions or productions and depth in GEInfo as limit.

# Self-provided Arguments
- `gpExperimentInfo::GEInfo`: information about grammatical evolution.
- `rng::Random.AbstractRNG`: random number generator for random number consistency along an experiment.

!!! note
    `Selfprovided Arguments` are already provided by the library, so only `User Arguments` must be provided.
"""
function growGenerator(gpExperimentInfo::GEInfo, rng::Random.AbstractRNG,
                        initInd::Array{UInt8}=Array{UInt8}(undef, 0))::GEGenotype

    if gpExperimentInfo._maxDepth < 0
        return growGeneratorProductions(gpExperimentInfo, rng, initInd)
    else
        return growGeneratorDepth(gpExperimentInfo, rng, initInd)
    end
end

"""
    function fullGeneratorProductions(gpExperimentInfo::GEInfo, rng::Random.AbstractRNG,
                    initInd::Array{UInt8}=Array{UInt8}(undef, 0))::GEGenotype

Generates a Grammatical Evolution tree ensuring that the individual is valid having
productions in GEInfo as limit and trying to reach that limit.

# Self-provided Arguments
- `gpExperimentInfo::GEInfo`: information about grammatical evolution.
- `rng::Random.AbstractRNG`: random number generator for random number consistency along an experiment.

!!! note
    `Selfprovided Arguments` are already provided by the library, so only `User Arguments` must be provided.
"""
function fullGeneratorProductions(gpExperimentInfo::GEInfo, rng::Random.AbstractRNG,
                                initInd::Array{UInt8}=Array{UInt8}(undef, 0))::GEGenotype

    productions::UInt16 = gpExperimentInfo._maxProductions
    gramm = gpExperimentInfo._grammar

    ind = initInd
    if isempty(ind)
        ind = rand(rng, UInt8, productions)
    end

    Root = GERep(getInitialSymbol(gramm))
    stack = [Root]

    # A production is always made due to initial symbol
    actualProductions =  productions - 0x2

    actualDepth::UInt8 = 0x0
    symbolsPerDepth = zeros(UInt8, 2)
    symbolsPerDepth[1] = 0x1

    selections = UInt8[]
    i = 1
    while !isempty(stack)
        if symbolsPerDepth[actualDepth + 0x1] == 0x0
            actualDepth += 0x1
        end

        symbolsPerDepth[actualDepth + 0x1] -= 0x1

        node = popfirst!(stack)
        symbol = getSymbol(node)
        if !isTerminalSymbol(gramm, symbol)
            actualProductions = actualProductions + getDepth(gramm, symbol) - 0x1
            indexes, nIndexes = getPosibleRecursiveProductions(gramm, symbol, actualProductions)
            """
            if actualProductions <= 0x1
                actualProductions = actualProductions + getDepth(gramm, symbol) - 0x1
                indexes, nIndexes = getPosibleConsequentProductions(gramm, symbol, actualProductions)
            else
                actualProductions = actualProductions + getDepth(gramm, symbol) - 0x1
                indexes, nIndexes = getPosibleRecursiveProductions(gramm, symbol, actualProductions)
                if nIndexes == 0x0
                    indexes, nIndexes = getPosibleConsequentProductions(gramm, symbol, actualProductions)

                end
            end
            """

            selected = ind[i] % nIndexes + 0x1
            selected = indexes[selected]

            nChoices = getNConsequent(gramm, symbol)
            ind[i] = ind[i] % (div(256,nChoices)) * nChoices + selected - 0x1

            push!(selections, selected)

            consequent = getConsequent(gramm, symbol, selected)
            if length(symbolsPerDepth) < actualDepth + 0x2
                push!(symbolsPerDepth, getNSymbols(consequent))
            else
                symbolsPerDepth[actualDepth + 0x2] += getNSymbols(consequent)
            end
            actualProductions = actualProductions - getMinProductions(consequent)
            symbols = consequent._symbols
            for j=length(symbols):-1:1
                child = GERep(symbols[j], actualDepth + 0x1)
                pushChild!(node, child)
                push!(stack, child)
            end
            i = i + 1
        end
    end

    return GEGenotype(Root, IntegerGenotype(ind), selections, i-1, true)
end # function

"""
    function fullGeneratorDepth(gpExperimentInfo::GEInfo, rng::Random.AbstractRNG,
                    initInd::Array{UInt8}=Array{UInt8}(undef, 0))::GEGenotype

Generates a Grammatical Evolution tree ensuring that the individual is valid having
productions and depth in GEInfo as limit and trying to reach that limit.

# Self-provided Arguments
- `gpExperimentInfo::GEInfo`: information about grammatical evolution.
- `rng::Random.AbstractRNG`: random number generator for random number consistency along an experiment.

!!! note
    `Selfprovided Arguments` are already provided by the library, so only `User Arguments` must be provided.
"""
function fullGeneratorDepth(gpExperimentInfo::GEInfo, rng::Random.AbstractRNG,
                            initInd::Array{UInt8}=Array{UInt8}(undef, 0))::GEGenotype

    maxDepth::UInt8 = gpExperimentInfo._maxDepth
    maxProductions::UInt16 = gpExperimentInfo._maxProductions
    gramm = gpExperimentInfo._grammar

    ind = initInd
    if isempty(ind)
        ind = rand(rng, UInt8, maxProductions)
    end

    Root = GERep(getInitialSymbol(gramm))
    stack = [Root]

    actualProductions =  maxProductions - 0x2
    # A production is always made due to initial symbol
    actualDepth::UInt8 = 0x0
    symbolsPerDepth = zeros(UInt8, maxDepth+1)
    symbolsPerDepth[1] = 0x1
    selections = UInt8[]
    i = 1
    while !isempty(stack)
        if symbolsPerDepth[actualDepth+0x1] == 0x0
            actualDepth += 0x1
        end

        symbolsPerDepth[actualDepth+0x1] -= 0x1
        node = popfirst!(stack)
        symbol = getSymbol(node)

        if !isTerminalSymbol(gramm, symbol)
            leftDepth = maxDepth-actualDepth
            actualProductions = actualProductions + getDepth(gramm, symbol) - 0x1

            indexes, nIndexes = getPosibleRecursiveDepth(gramm, symbol, maxDepth-actualDepth, actualProductions)
            """
            if leftDepth == 0x1
                indexes, nIndexes = getPosibleConsequentDepth(gramm, symbol, leftDepth)
            else
                indexes, nIndexes = getPosibleRecursiveDepth(gramm, symbol, leftDepth)
                if nIndexes == 0x0
                    indexes, nIndexes = getPosibleConsequentDepth(gramm, symbol, leftDepth)
                end
            end
            """

            selected = ind[i] % nIndexes + 0x1
            selected = indexes[selected]

            nChoices = getNConsequent(gramm, symbol)
            ind[i] = ind[i] % (div(256,nChoices)) * nChoices + selected - 0x1

            push!(selections, selected)

            consequent = getConsequent(gramm, symbol, selected)
            symbolsPerDepth[actualDepth + 0x2] += getNSymbols(consequent)
            actualProductions = actualProductions - getMinProductions(consequent)
            symbols = consequent._symbols
            for j=length(symbols):-1:1
                child = GERep(symbols[j], actualDepth+0x1)
                pushChild!(node, child)
                push!(stack, child)
            end
            i = i + 1
        end
    end

    return GEGenotype(Root, IntegerGenotype(ind), selections, i-1, true)
end # function

"""
    function fullGenerator(gpExperimentInfo::GEInfo, rng::Random.AbstractRNG,
                    initInd::Array{UInt8}=Array{UInt8}(undef, 0))::GEGenotype

Generates a Grammatical Evolution tree ensuring that the individual is valid having
productions or productions and depth in GEInfo as limit and trying to reach that limit.

# Self-provided Arguments
- `gpExperimentInfo::GEInfo`: information about grammatical evolution.
- `rng::Random.AbstractRNG`: random number generator for random number consistency along an experiment.

!!! note
    `Selfprovided Arguments` are already provided by the library, so only `User Arguments` must be provided.
"""
function fullGenerator(gpExperimentInfo::GEInfo, rng::Random.AbstractRNG,
                        initInd::Array{UInt8}=Array{UInt8}(undef, 0))::GEGenotype

    if gpExperimentInfo._maxDepth < 0
        return fullGeneratorProductions(gpExperimentInfo, rng, initInd)
    else
        return fullGeneratorDepth(gpExperimentInfo, rng, initInd)
    end
end

"""
    function generateTreeProductions(gpExperimentInfo::GEInfo, rng::Random.AbstractRNG,
                    codons::Array{UInt8})::GEGenotype

Generates a Grammatical Evolution tree from codons having productions in GEInfo as limit, if
the limit is overcome, then the individual is not valid.

# Self-provided Arguments
- `gpExperimentInfo::GEInfo`: information about grammatical evolution.
- `rng::Random.AbstractRNG`: random number generator for random number consistency along an experiment.

# User Arguments
- `codons::Array{UInt8}`: array with the proteins that represents a number related with the selections to be done.

!!! note
    `Selfprovided Arguments` are already provided by the library, so only `User Arguments` must be provided.
"""
function generateTreeProductions(gpExperimentInfo::GEInfo, rng::Random.AbstractRNG,
                                codons::Array{UInt8})::GEGenotype

    gramm = gpExperimentInfo._grammar
    maxProductions = gpExperimentInfo._maxProductions

    Root = GERep(getInitialSymbol(gramm))
    stack = [Root]
    selections = UInt8[]

    actualDepth::UInt8 = 0x0
    symbolsPerDepth = zeros(UInt8, 2)
    symbolsPerDepth[1] = 0x1

    i = 1
    while !isempty(stack) && i<=maxProductions
        if symbolsPerDepth[actualDepth + 0x1] == 0
            actualDepth += 0x1
        end

        symbolsPerDepth[actualDepth + 0x1] -= 0x1

        node = popfirst!(stack)
        symbol = getSymbol(node)
        if !isTerminalSymbol(gramm, symbol)

            selected = codons[i] % getNConsequent(gramm, symbol) + 0x1

            push!(selections, selected)
            consequent = getConsequent(gramm, symbol, selected)
            if length(symbolsPerDepth) < actualDepth + 0x2
                push!(symbolsPerDepth, getNSymbols(consequent))
            else
                symbolsPerDepth[actualDepth + 0x2] += getNSymbols(consequent)
            end
            symbols = consequent._symbols
            for j=length(symbols):-1:1
                child = GERep(symbols[j], actualDepth + 0x1)
                pushChild!(node, child)
                push!(stack, child)
            end
            i = i + 1
        end
    end
    if i > maxProductions
        return GEGenotype(Root, IntegerGenotype{UInt8}(codons), selections, i-1, false)
    else
        return GEGenotype(Root, IntegerGenotype{UInt8}(codons), selections, i-1, true)
    end
end

"""
    function generateTreeDepth(gpExperimentInfo::GEInfo, rng::Random.AbstractRNG,
                    codons::Array{UInt8})::GEGenotype

Generates a Grammatical Evolution tree from codons having productions and depth
in GEInfo as limit, if the limit is overcome, then the individual is not valid.

# Self-provided Arguments
- `gpExperimentInfo::GEInfo`: information about grammatical evolution.
- `rng::Random.AbstractRNG`: random number generator for random number consistency along an experiment.

# User Arguments
- `codons::Array{UInt8}`: array with the proteins that represents a number related with the selections to be done.

!!! note
    `Selfprovided Arguments` are already provided by the library, so only `User Arguments` must be provided.
"""
function generateTreeDepth(gpExperimentInfo::GEInfo, rng::Random.AbstractRNG,
                                codons::Array{UInt8})::GEGenotype

    gramm = gpExperimentInfo._grammar
    maxProductions = gpExperimentInfo._maxProductions
    maxDepth = gpExperimentInfo._maxDepth

    Root = GERep(getInitialSymbol(gramm))
    stack = [Root]
    selections = UInt8[]

    actualDepth::UInt8 = 0x0
    symbolsPerDepth = zeros(UInt8, 2)
    symbolsPerDepth[1] = 0x1

    i = 1
    while !isempty(stack) && i<=maxProductions && actualDepth <= maxDepth
        if symbolsPerDepth[actualDepth + 0x1] == 0
            actualDepth += 0x1
        end

        symbolsPerDepth[actualDepth + 0x1] -= 0x1

        node = popfirst!(stack)
        symbol = getSymbol(node)
        if !isTerminalSymbol(gramm, symbol)

            selected = codons[i] % getNConsequent(gramm, symbol) + 0x1

            push!(selections, selected)
            consequent = getConsequent(gramm, symbol, selected)
            if length(symbolsPerDepth) < actualDepth + 0x2
                push!(symbolsPerDepth, getNSymbols(consequent))
            else
                symbolsPerDepth[actualDepth + 0x2] += getNSymbols(consequent)
            end
            symbols = consequent._symbols
            for j=length(symbols):-1:1
                child = GERep(symbols[j], actualDepth + 0x1)
                pushChild!(node, child)
                push!(stack, child)
            end
            i = i + 1
        end
    end
    if i > maxProductions || actualDepth > maxDepth
        return GEGenotype(Root, IntegerGenotype{UInt8}(codons), selections, i-1, false)
    else
        return GEGenotype(Root, IntegerGenotype{UInt8}(codons), selections, i-1, true)
    end
end

"""
    function generateTree(gpExperimentInfo::GEInfo, rng::Random.AbstractRNG,
                                    codons::Array{UInt8})::GEGenotype

Generates a Grammatical Evolution tree ensuring that the individual is valid having
productions or productions and depth in GEInfo as limit and trying to reach that limit.

# Self-provided Arguments
- `gpExperimentInfo::GEInfo`: information about grammatical evolution.
- `rng::Random.AbstractRNG`: random number generator for random number consistency along an experiment.

# User Arguments
- `codons::Array{UInt8}`: array with the proteins that represents a number related with the selections to be done.

!!! note
    `Selfprovided Arguments` are already provided by the library, so only `User Arguments` must be provided.
"""
function generateTree(gpExperimentInfo::GEInfo, rng::Random.AbstractRNG,
                                codons::Array{UInt8})::GEGenotype

    if gpExperimentInfo._maxDepth < 0
        return generateTreeProductions(gpExperimentInfo, rng, codons)
    else
        return generateTreeDepth(gpExperimentInfo, rng, codons)
    end
end

# TODO Cruce basado en probabilidad.
"""

function geGenerator(gpExperimentInfo::GEInfo, rng::Random.AbstractRNG, initInd::Array{Integer}=Array{Integer}(undef, 0))

    productions = gpExperimentInfo._maxProductions
    gramm = gpExperimentInfo._grammar
    growRatio = gpExperimentInfo._growRatio

    ind = initInd
    if length(ind) == 0
        ind = [rand(rng, UInt16) for x=1:productions]
    end

    Root = GERep(gramm._S)
    stack = [Root]
    auxStack = [getSymbol(Root)]
    # A production is always made due to initial symbol
    actualProductions = productions - 2

    actualDepth = 0
    symbolsPerDepth = zeros(UInt8, 2)
    symbolsPerDepth[1] = 1

    selections = []
    i = 1
    while !isempty(stack)
        if symbolsPerDepth[actualDepth + 1] == 0
            actualDepth += 1
        end

        symbolsPerDepth[actualDepth + 1] -= 1

        node = popfirst!(stack)
        symbol = getSymbol(node)
        if !isTerminalSymbol(gramm, symbol)
            if actualProductions <= 1
                actualProductions = actualProductions + getDepth(gramm, symbol) - 1
                indexes = getPosibleConsequentProductions(gramm, symbol, actualProductions)
            else
                actualProductions = actualProductions + getDepth(gramm, symbol) - 1

                if rand(rng) < ratio

                end
            end

            indexes = getPosibleConsequentProductions(gramm, symbol, actualProductions)

            nIndexes = length(indexes)
            selected = ind[i] % nIndexes + 1
            selected = indexes[selected]

            nChoices = getNConsequent(gramm, symbol)
            ind[i] = ind[i] % (256/nChoices) * nChoices + selected

            push!(selections, selected)
            consequent = getConsequent(gramm, symbol, selected)
            if length(symbolsPerDepth) < actualDepth + 2
                push!(symbolsPerDepth, getNSymbols(consequent))
            else
                symbolsPerDepth[actualDepth + 2] += getNSymbols(consequent)
            end
            actualProductions = actualProductions - getMinProductions(consequent)
            symbols = consequent._symbols
            for j=1:length(symbols)
                child = GERep(symbols[j], actualDepth + 1)
                pushChild!(node, child)
                push!(stack, child)
                push!(auxStack, symbols[j])
            end
            i = i + 1
        end
    end

    return GEGenotype(Root, IntegerGenotype(ind), selections, i-1)
end # function
"""

max_stars_repo_name: SergioGmezM/GeneticJulia

max_issues_repo_name: SergioGmezM/GeneticJulia

max_forks_repo_name: SergioGmezM/GeneticJulia


text:
#MSCOCO 2014
using JSON
#file = open("ids_coco_bm4");
file = open("candidate_ids.txt");
candidate_ids = readlines(file);
candidate_ids = map(p->parse(Int,p), candidate_ids)

file = open("../data/MsCoCo/captions_val2014.json");
captions = JSON.parse(file);
annotations = captions["annotations"];

caps = Dict();
for item in annotations
         arr = get!(caps, item["image_id"], [])
         if length(arr) == 5
           continue;
         end
         cap = strip(item["caption"]);
         cap = strip(cap,['.']);
         cap = cap * " ."
         push!(arr, lowercase(cap))
end

refs = Any[]
for i=0:4
  push!(refs, open("./coco_refs/ref$(i)","w"))
end

for id in candidate_ids
           arr = get(caps,id,nothing)
           for i=1:length(refs)
             println(refs[i],strip(arr[i]))
           end
end

map(p->close(p),refs);
println("MSCOCO Scores")
run(pipeline(`perl multi-bleu.perl ./coco_refs/ref`, stdin="candidates.txt"))
#Flickr30k
file = open("candidate_ids_flickr");
candidate_ids = readlines(file);
candidate_ids = map(p->parse(Int,p), candidate_ids)

file = open("../data/Flickr30k/results_20130124.token");
captions = readlines(file);

caps = Dict();
for line in captions
          info = split(line,"#")
          id = info[1]
          id = split(id,".")[1]
          id = parse(Int,id)
          cap = info[2]
          cap = split(cap,"\t")[2]
          cap = lowercase(strip(cap))
          arr = get!(caps,id, [])
          push!(arr,cap)
end

refs = Any[]
for i=0:4
  push!(refs, open("./flickr_refs/f_ref$(i)","w"))
end

for id in candidate_ids
           arr = get(caps,id,nothing)
           if arr==nothing
             println("id is missing in reference: ", id)
             break;
           end
           for i=1:length(refs)
             println(refs[i],strip(arr[i]))
           end
end
map(p->close(p),refs);

println("Flickr30k Scores")
run(pipeline(`perl multi-bleu.perl ./flickr_refs/f_ref`, stdin="candidates_flickr"))

#run(`perl multi-bleu.perl f_ref < candidates_f.txt`)
#run(`perl multi-bleu.perl ref < candidates.txt`)

ext: jl

max_issues_repo_name: ekinakyurek/LRCN

max_forks_repo_name: ekinakyurek/LRCN


text:
include(joinpath(dirname(@__DIR__), "src", "BestApproximation.jl")); using .BestApproximation
using Test

@testset "BestApproximation.jl" begin
    @test best_approx(factorial(big(26)), 1:4) == (2, 88)
    @test best_approx(factorial(big(26)), 4) == (2, 88)
    @test best_approx(factorial(big(26)), 10) == (5, 38)
    @test best_approx(factorial(big(26)), 100) == (46, 16)
    @test best_approx(factorial(big(26)), 1000) == (904, 9)
    @test best_approx(factorial(big(26)), 10000) == (904, 9)
    @test best_approx(factorial(big(26)), 100000) == (27181, 6)
    @test best_approx(factorial(big(26)), 1000000) == (209470, 5)
    @test best_approx(1234567890, 94) == (5, 13)
    @test best_approx(12345678901234567890, 94) == (3, 40)
    @test best_approx(123456789012345678901234567890, 94878) == (70564, 6)
    @test best_approx(factorial(big(26)), 1:4, 1:100) == (2, 88)
    @test best_approx(factorial(big(26)), 4, 100) == (2, 88)
    @test best_approx(factorial(big(26)), 10, 100) == (5, 38)
end

max_stars_repo_name: jakewilliami/BestApproximation

max_issues_repo_name: jakewilliami/BestApproximation

max_forks_repo_name: jakewilliami/BestApproximation


text:
struct DistributionFreeTest{E<:SKCE,B,V} <: CalibrationTest
    """Calibration estimator."""
    estimator::E
    """Upper bound of estimator."""
    bound::B
    """Number of data points."""
    n::Int
    """SKCE estimate."""
    estimate::V
end

function DistributionFreeTest(estimator::SKCE,
                              data::Tuple{<:AbstractMatrix{<:Real},<:AbstractVector{<:Integer}};
                              bound = maximum(estimator))
    # check if the number of predictions and labels is equal
    predictions, labels = CalibrationErrors.get_predictions_labels(data)

    DistributionFreeTest(estimator, bound, size(predictions, 2),
                         calibrationerror(estimator, data))
end

HypothesisTests.testname(::DistributionFreeTest) = "Distribution-free calibration test"

function HypothesisTests.pvalue(test::DistributionFreeTest{<:BiasedSKCE})
    @unpack bound, n, estimate = test

    s = sqrt(n * estimate / bound) - 1

    # if estimate is below threshold B₂
    s < zero(s) && (s = zero(s))

    exp(- s^2 / 2)
end

function HypothesisTests.pvalue(test::DistributionFreeTest{<:Union{QuadraticUnbiasedSKCE,
                                                                   LinearUnbiasedSKCE}})
    @unpack bound, n, estimate = test

    exp(- div(n, 2) * estimate^2 / (2 * bound ^ 2))
end


text:
function NDMFile()
    return NDMFile(())
end


text:
using YaoLang
using YaoLang.Compiler
using YaoLang.Compiler: device_m, enable_timings, timings

# enable_timings()

ex = :(function qft(n::Int)
    1 => H
    for k in 2:n
        @ctrl k 1 => shift(2π / 2^k)
    end

    if n > 1
        2:n => qft(n - 1)
    end
end)

@time device_m(@__MODULE__, ex);

# timings()


text:
module PlaneProblems

using ..Input: CellData
using ..TypeSynonyms
using ..Symbolic: NormedPolynomial, add!, z_conj_diff, re_conjugate!

using ..General

abstract type PlaneProblem <: ProblemType end

struct PlaneLayer
    E       ::Float
    ν       ::Float
    r_inner ::Float
    r_outer ::Float
    ϕ       ::NormedPolynomial{Complex{Float}, Float}
    z_bar_Φ ::NormedPolynomial{Complex{Float}, Float}
    bar_ψ   ::NormedPolynomial{Complex{Float}, Float}
end

function couple!(dest::PlaneLayer, source::PlaneLayer)
    if dest.r_inner != source.r_outer
        error("Layers don't touch")
    end

    E1 = source.E
    ν1 = source.ν
    G1 = E1/(2*(1+ν1))
    κ1 = 3 - 4*ν1

    E2 = dest.E
    ν2 = dest.ν
    κ2 = 3 - 4*ν2
    G2 = E2/(2*(1+ν2))

    empty!(dest.ϕ)
    empty!(dest.z_bar_Φ)
    empty!(dest.bar_ψ)

    add!(dest.ϕ, source.ϕ,       (1+κ1*G2/G1)/(1+κ2))
    add!(dest.ϕ, source.z_bar_Φ, (1-G2/G1)   /(1+κ2))
    add!(dest.ϕ, source.bar_ψ,   (1-G2/G1)   /(1+κ2))

    dest.z_bar_Φ = z_conj_diff(dest.ϕ, dest.r_inner)

    add!(dest.bar_ψ, source.ϕ)
    add!(dest.bar_ψ, source.z_bar_Φ)
    add!(dest.bar_ψ, source.bar_ψ)
    add!(dest.bar_ψ, dest.ϕ, -1)
    add!(dest.bar_ψ, dest.z_bar_Φ, -1)

    re_conjugate!(dest.bar_ψ,   dest.r_inner, dest.r_outer)
    re_conjugate!(dest.z_bar_Φ, dest.r_inner, dest.r_outer)
end

function displacements(layer::PlaneLayer)
    E = layer.E
    ν = layer.ν
    G = E/(2*(1+ν))
    κ = 3 - 4*ν

    displ = NormedPolynomial{Complex{Float}, Float}()
    add!(displ, layer.ϕ,        κ/(2G))
    add!(displ, layer.z_bar_Φ, -1/(2G))
    add!(displ, layer.bar_ψ,   -1/(2G))

    return displ
end

function forces(layer::PlaneLayer)
    f = NormedPolynomial{Complex{Float}, Float}()

    add!(f, layer.ϕ)
    add!(f, layer.z_bar_Φ, -1)
    add!(f, layer.bar_ψ,   -1)

    return f
end

#######################################

struct PlaneFiber
    layers ::Vector{PlaneLayer}
end

struct PlaneFiber_Displacements <: Fiber{PlaneProblem}
    fiber ::PlaneFiber
end

struct PlaneFiber_Displacements <: Fiber{PlaneProblem}
    fiber ::PlaneFiber
end

#=
using ..Symbolic: Polynomial, max_abs_index, add!

struct Layer
    E::Float64
    ν::Float64

    ϕ         ::Polynomial
    z_bar_Phi ::Polynomial
    bar_ψ     ::Polynomial
end

function fill_slae!(layer::Layer, slae)
    ϕ = layer.ϕ
    ψ = layer.ψ

    E = layer.E
    ν = layer.ν
    G = E / (2.0*(1.0+ν))
    κ = 3.0-4.0*ν

    max_index = maximum(max_abs_index(ϕ), max_abs_index(ψ))
    displacements = Polynomial(-max_index : max_index)
    stresses      = Polynomial(-max_index : max_index)

    add_plain!(displacements, ϕ, κ/G)
    add_z_conj_diff!(displacements, ϕ, -1.0/G)
    add_conj!(displacements, ψ, -1.0/G)

    add_plain!(stresses, ϕ)
    add_z_conj_diff!(stresses, ϕ)
    add_conj!(stresses, ψ)

    # TODO fill slae
end

struct Fiber
    layers ::Vector{Layer}
end

function fill_slae!(fiber::Fiber, slae)
    fill_slae!(fiber.layers[end], slae)
end

struct Cohesive

end

struct PlaneProblem
    fibers    ::Dict{Pole, Fiber}
    cohesive  ::Cohesive

    slae_rows ::Dict{Pole, UnitRange{Int64}}
    slae_left ::Matrix{Float64}
    slae_right::Vector{Float64}
end

function PlaneProblem(cell::CellData)
    # TODO
end

function fill_slae!(cohesive::Cohesive, pole::Pole, slae)
    # TODO
end

function fill_slae!(plane_problem::PlaneProblem)
    slae_left = plane_problem.slae_left
    rows = plane_problem.slae_rows

    for pole in keys(plane_problem.fibers)
        slae_view = view(slae_left[rows[pole]])
        fill_slae!(plane_problem.fibers[pole], slae_view)
        fill_slae!(plane_problem.cohesive, pole, slae_view)
    end

    # TODO averaging
end

=#

end # module PlaneProblems


text:
function JSif(ex::Expr)
    # The if statement is outside of a macro, so it gets turned into a mux to verilog
    # The wires are not turned into regs
    jexpr = Expr(:Block)

    # Convert Juliog conditions to JIR if necessary
    ex = JSifconds(ex)

    # Recursively solve for nested if assignments
    jargs = JSifblock(ex)

    l = length(jargs)
    for i = 1:l
        jargs[i] = JSblockhelper(jarg[i])
    end

    jexpr.args = jargs
    return jexpr
end

function JSifconds(ex::Expr)
    if !isa(ex.args[1], Expr)
        error("Condition for if statement is not an expr:\n$(ex)")
    end

    # Convert top level condition
    ex.args[1] = JSifcond(ex.args[1])

    # Plunge into if block, convert any nested conditions
    for i = 1:length(ex.args[2].args)
        arg = ex.args[2].args[i]
        if isa(arg, Expr) && arg.head == :if
            ex.args[2].args[i] = JSifconds(arg)
        end
    end

    # Plunge into else block, if it exists, convert any nested conditions
    if length(ex.args) == 3
        for i = 1:length(ex.args[3].args)
            arg = ex.args[3].args[i]
            if isa(arg, Expr) && arg.head == :if
                ex.args[3].args[i] = JSifconds(arg)
            end
        end
    end

    return ex
end

function JSifcond(ex::Expr)
    global arrowDict

    # Check for warnings in condition syntax
    JSrhscall(ex)

    if     ex.args[3] == 1 || ex.args[3] == 0
        if     isa(ex.args[2], Symbol)
            bc = arrowDict[ex.args[2]].BitCount
        elseif isa(ex.args[2], Expr)
            # MOON Right here solve matrix notation and swap ex.args[2] to it
            bc,endi = JSreadrhs(ex.args[2])
        end

        if bc == 1
            if     ex.args[3] == 1
                return ex.args[2]
            elseif ex.args[3] == 0
                return [:~ ; ex.args[2]]
            end
        else
            error("Condition of if statement does not equal 1 bit:\n$(ex)")
        end
    elseif ex.args[2] == 1 || ex.args[2] == 0
        if     isa(ex.args[3], Symbol)
            bc = arrowDict[ex.args[3]].BitCount
        elseif isa(ex.args[3], Expr)
            # MOON Right here solve matrix notation and swap ex.args[2] to it
            bc,endi = JSreadrhs(ex.args[3])
        end

        if bc == 1
            if     ex.args[2] == 1
                return ex.args[3]
            elseif ex.args[2] == 0
                return [:~ ; ex.args[3]]
            end
        else
            error("Condition of if statement does not equal 1 bit:\n$(ex)")
        end
    end

    # Can't be reduced to a single wire, so just return as is
    return ex
end

function JSifblock(ex::Expr)
    # Need to recuresively search for nested if statements
    if length(ex.args) == 3
        # Has both an if statement and an else statement
        # We may need to merge if assignments and else statements
        return JSifblockhelper(ex)
    else
        # No else statement so just raise the if block
        elseargs = Array{Any}(0)
        l = length(ex.args[2].args)
        for i = 1:l
            arg = ex.args[2].args[i]
            h = arg.head
            if     h == :(=)
                JSifequals(arg)
                push!(elseargs, arg)
            elseif h == :if
                append!(elseargs, JSifblock(arg))
            else
                error("Found incorrect expr head $(h) in if expr:\n$(ex)")
            end
        end
        l = length(elseargs)
        for i = 1:l
            arg = elseargs[i]
            elseargs[i].args[2] = Expr(:if, ex.args[1], arg.args[2], arg.args[1])
            warn("Latch of expr $(elseargs[i].args[2]) created in if expr:\n$(ex)")
        end
        return elseargs
    end
end

function JSifblockhelper(ex::Expr)
    ifargs = Array{Any}(0)
    elseargs = Array{Any}(0)

    # Iterate through if statements and create if-assignments
    l = length(ex.args[2].args)
    for i = 1:l
        arg = ex.args[2].args[i]
        h = arg.head
        if     h == :(=)
            JSifblockequals(arg)
            push!(ifargs, arg)
        elseif h == :if
            args = JSifblock(arg)
            append!(ifargs, args)
        else
            error("Found incorrect expr head $(h) in if expr:\n$(ex)")
        end
    end

    # Iterate through else statements and create if-assignments
    l = length(ex.args[3].args)
    for i = 1:l
        arg = ex.args[3].args[i]
        h = arg.head
        if     h == :(=)
            JSifblockequals(arg)
            push!(elseargs, arg)
        elseif h == :if
            args = JSifblock(arg)
            append!(elseargs, args)
        else
            error("Found incorrect expr head $(h) in if expr:\n$(ex)")
        end
    end

    # loop through to find redundancy in elseargs
    l = length(elseargs)
    elseargslhs = Array{Any}(l)
    for i = 1:l
        lhs = elseargs[i].args[1]
        if find(elseargslhs .== lhs) != []
            error("Attempting to write to lhs of $(elseargs[i]) in:\n$(ex)")
        end
        elseargslhs[i] = lhs
    end
    # loop through to find redundancy in ifargs
    l = length(ifargs)
    ifargslhs = Array{Any}(l)
    for i = 1:l
        lhs = ifargs[i].args[1]
        if find(ifargslhs .== lhs) != []
            error("Attempting to write to lhs of $(ifargs[i]) in:\n$(ex)")
        end
        ifargslhs[i] = lhs
    end

    # TODO compare ifargs against elseargs
    jargs = Array{Any}(0)
    for i = 1:l

        jarg = jargs[j]
        jargs[j].args[2] = Expr(:if, condition, jarg.args[2], jarg.args[1])
    end

    # finish up with elseargs
    l = length(elseargs)
    for i = 1:l

    end

    append!(jargs, ifargs)
    ex.args[2].args[i] = jexpr

    return jargs
end

function JSifblockequals(ex::Expr)
    # Check rhs correctness
    bcrhs, endirhs = JSrhs(ex.args[2])

    # check lhs correctness
    if isa(ex.args[1], Expr)
        # lhs is a reference
        JSifblockequalsref(ex, bcrhs, endirhs)
    else
        # lhs is just a symbol for a wire
        JSifblockequalssymbol(ex, bcrhs, endirhs)
    end
end

function JSifblockequalsref(ex::Expr, bcrhs::Int, endirhs::Int)
    # LHS is a reference
    #   Maybe a vcat or tuple, can't do anything about that yet
    name = ex.args[1].args[1]

    # Check to see if the rhs and lhs bitcounts agree
    #  Throw a warning if that is not the case
    arg = ex.args[1].args[2]
    if isa(arg, Expr)
        # Multiple bit or single bit
        bl = arg.args[1]
        l = length(arg.args)
        if l == 3
            br = arg.args[3]
        else
            br = arg.args[2]
        end

        # Determine the endianness and bitcount of the lhs
        if bl == br
            bclhs = 1
            endilhs = 0
        elseif bl > br
            bclhs = bl - br + 1
            endilhs = 1
        else
            bclhs = br - bl + 1
            endilhs = 2
        end

        if     bclhs > bcrhs
            warn("Zero Bit Padding necessitated: $(bclhs) bits on lhs, and $(bcrhs) bits on rhs:\n$(ex)")
        elseif bclhs < bcrhs
            warn("Truncation of Bits necessitated: $(bclhs) bits on lhs, and $(bcrhs) bits on rhs:\n$(ex)")
        end

        if     endilhs == 1 && endirhs == 2
            error("Endian Mismatch: LHS is Big Endian, RHS is Little Endian:\n$(ex)")
        elseif endilhs == 2 && endilhs == 1
            error("Endian Mismatch: LHS is Little Endian, RHS is Big Endian:\n$(ex)")
        end
    else
        # Single bit
        bclhs = 1
        bl = arg
        br = arg
        endilhs = 0
        if     1 > bcrhs
            warn("Zero Bit Padding necessitated: $(1) bit on lhs, and $(bcrhs) bits on rhs:\n$(ex)")
        elseif 1 < bcrhs
            warn("Truncation of Bits necessitated: $(1) bit on lh, and $(bcrhs) bits on rhs:\n$(ex)")
        end
    end
end

function JSifblockequalssymbol(ex::Expr, bcrhs::Int, endirhs::Int)
    # Must be a wire name
    name = ex.args[1]

    # If the lhs and rhs bits don't agree, throw out warnings
    if haskey(ifAD, name)
        bclhs = ifAD[name].BitCount
        if     bclhs > bcrhs
            warn("Zero Bit Padding necessitated: $(bclhs) bits on lhs, and $(bcrhs) bits on rhs:\n$(ex)")
        elseif bclhs < bcrhs
            warn("Truncation of Bits necessitated: $(bclhs) bit on lh, and $(bcrhs) bits on rhs:\n$(ex)")
        end
    end
end

max_stars_repo_name: trevortheblack/Juliog

max_issues_repo_name: trevortheblack/Juliog

max_forks_repo_name: trevortheblack/Juliog


text:
"""
    caxs = centeraxes(axs)

Return a set of axes centered on zero. Specifically, if `axs[i]` is a range,
then `caxs[i]` is a UnitRange of the same length that is approximately symmetric
around 0.
"""
centeraxes(axs) = map(centeraxis, axs)

function centeraxis(ax)
    f, l = first(ax), last(ax)
    n = l - f + 1
    nhalf = (n+1) ÷ 2
    return 1-nhalf:n-nhalf
end

function extract1(u::AbstractArray{V}, N, ssz) where V<:SVector
    if ndims(u) == N+1
        ϕ = GridDeformation(reshape(u, size(u)[1:end-1]), ssz)
    else
        ϕ = GridDeformation(u, ssz)
    end
    ϕ
end
extract1(ϕs::Vector{D}, N, ssz) where {D<:AbstractDeformation} = ϕs[1]

function extracti(u::AbstractArray{V}, i, ssz) where V<:SVector
    colons = [Colon() for d = 1:ndims(u)-1]
    GridDeformation(u[colons..., i], ssz)
end
extracti(ϕs::Vector{D}, i, _) where {D<:AbstractDeformation} = ϕs[i]

function checkϕdims(u::AbstractArray{V}, N, n) where V<:SVector
    ndims(u) == N+1 || error("u's dimensionality $(ndims(u)) is inconsistent with the number of spatial dimensions $N of the image")
    if size(u)[end] != n
        error("Must have one `u` slice per image")
    end
    nothing
end
checkϕdims(ϕs::Vector{D}, N, n) where {D<:AbstractDeformation} = length(ϕs) == n || error("Must have one `ϕ` per image")

# TODO?: do we need to return real values beyond-the-edge for a SubArray?

floattype(::Type{T}) where T<:AbstractFloat = T
floattype(::Type{<:Integer}) = Float64
floattype(::Type{SVector{N,T}}) where {N,T} = floattype(T)

nanvalue(::Type{T}) where T<:Real = convert(promote_type(T, Float32), NaN)
nanvalue(::Type{C}) where C<:AbstractGray = Gray(nanvalue(eltype(C)))
nanvalue(::Type{C}) where C<:AbstractRGB  = (x = nanvalue(eltype(C)); RGB(x, x, x))

to_etp(img) = extrapolate(interpolate(img, BSpline(Linear())), nanvalue(eltype(img)))

to_etp(itp::AbstractInterpolation) = extrapolate(itp, convert(promote_type(eltype(itp), Float32), NaN))

to_etp(etp::AbstractExtrapolation) = etp

to_etp(img, A::AffineMap) = TransformedArray(to_etp(img), A)

getcoefs(ϕ::GridDeformation) = getcoefs(ϕ.u)
getcoefs(itp::AbstractInterpolation) = getcoefs(itp.itp)
getcoefs(itp::Interpolations.BSplineInterpolation) = itp.coefs
getcoefs(u::AbstractArray{<:SVector}) = u

# Extensions to Interpolations and StaticArrays
@generated function vecindex(A::AbstractArray, x::SVector{N}) where N
    args = [:(x[$d]) for d = 1:N]
    meta = Expr(:meta, :inline)
    quote
        $meta
        getindex(A, $(args...))
    end
end

@generated function vecindex(A::AbstractInterpolation, x::SVector{N}) where N
    args = [:(x[$d]) for d = 1:N]
    meta = Expr(:meta, :inline)
    quote
        $meta
        A($(args...))
    end
end

@generated function vecgradient!(g, itp::AbstractArray, x::SVector{N}) where N
    args = [:(x[$d]) for d = 1:N]
    meta = Expr(:meta, :inline)
    quote
        $meta
        Interpolations.gradient!(g, itp, $(args...))
    end
end

function convert_to_fixed(u::Array{T}, sz=size(u)) where T
    N = sz[1]
    convert_to_fixed(SVector{N, T}, u, tail(sz))
end

# Unlike the one above, this is type-stable
function convert_to_fixed(::Type{SVector{N,T}}, u::AbstractArray{T}, sz=tail(size(u))) where {T,N}
    reshape(reinterpret(SVector{N,T}, vec(u)), sz)
end

@generated function copy_ctf!(dest::Array{SVector{N,T}}, src::Array) where {N,T}
    exvec = [:(src[offset+$d]) for d=1:N]
    quote
        for i = 1:length(dest)
            offset = (i-1)*N
            dest[i] = SVector($(exvec...))
        end
        dest
    end
end

function convert_from_fixed(uf::AbstractArray{SVector{N,T}}, sz=size(uf)) where {N,T}
    if isbitstype(T) && isa(uf, Array)
        u = reshape(reinterpret(T, vec(uf)), (N, sz...))
    else
        u = Array{T}(undef, N, sz...)
        for i = 1:length(uf)
            for d = 1:N
                u[d,i] = uf[i][d]
            end
        end
    end
    u
end

# # Note this is a bit unsafe as it requires the user to specify C correctly
# @generated function Base.convert{R,C,T}(::Type{SMatrix{Tuple{R,C},T}}, v::Vector{SVector{R,T}})
#     args = [:(Tuple(v[$d])) for d = 1:C]
#     :(SMatrix{Tuple{R,C},T}(($(args...),)))
# end

# Wrapping functions to interface with CoordinateTransfromations instead of AffineTransfroms module
tformeye(m::Int) = AffineMap(Matrix{Float64}(I,m,m), zeros(m))
tformtranslate(trans::Vector) = (m = length(trans); AffineMap(Matrix{Float64}(I,m,m), trans))

rotation2(angle) = RotMatrix(angle)
function tformrotate(angle)
    A = RotMatrix(angle)
    AffineMap(A, zeros(eltype(A),2))
end

function rotationparameters(R::Matrix)
    size(R, 1) == size(R, 2) || error("Matrix must be square")
    if size(R, 1) == 2
        return [atan(-R[1,2],R[1,1])]
    elseif size(R, 1) == 3
        aa = AngleAxis(R)
        return rotation_angle(aa)*rotation_axis(aa)
    else
        error("Rotations in $(size(R, 1)) dimensions not supported")
    end
end

function rotation3(axis::Vector{T}, angle) where T
    n = norm(axis)
    axisn = n>0 ? axis/n : (tmp = zeros(T,length(axis)); tmp[1] = 1; tmp)
    AngleAxis(angle, axisn...)
end

function rotation3(axis::Vector{T}) where T
    n = norm(axis)
    axisn = n>0 ? axis/n : (tmp = zeros(typeof(one(T)/1),length(axis)); tmp[1] = 1; tmp)
    AngleAxis(n, axisn...)
end

function tformrotate(axis::Vector, angle)
    if length(axis) == 3
        return AffineMap(rotation3(axis, angle), zeros(eltype(axis),3))
    else
        error("Dimensionality ", length(axis), " not supported")
    end
end

function tformrotate(x::Vector)
    if length(x) == 3
        return AffineMap(rotation3(x), zeros(eltype(x),3))
    else
        error("Dimensionality ", length(x), " not supported")
    end
end

#=
# The following assumes uaxis is normalized
function _rotation3(uaxis::Vector, angle)
    if length(uaxis) != 3
        error("3d rotations only")
    end
    ux, uy, uz = uaxis[1], uaxis[2], uaxis[3]
    c = cos(angle)
    s = sin(angle)
    cm = one(typeof(c)) - c
    R = [c+ux*ux*cm       ux*uy*cm-uz*s      ux*uz*cm+uy*s;
         uy*ux*cm+uz*s      c+uy*uy*cm       uy*uz*cm-ux*s;
         uz*ux*cm-uy*s    uz*uy*cm+ux*s      c+uz*uz*cm]
end

function rotation3{T}(axis::Vector{T}, angle)
    n = norm(axis)
    axisn = n>0 ? axis/n : (tmp = zeros(T,length(axis)); tmp[1] = 1)
    _rotation3(axisn, angle)
end

# Angle/axis representation where the angle is the norm of the vector (so axis is not normalized)
function rotation3{T}(axis::Vector{T})
    n = norm(axis)
    axisn = n>0 ? axis/n : (tmp = zeros(typeof(one(T)/1),length(axis)); tmp[1] = 1; tmp)
    _rotation3(axisn, n)
end
=#


text:
"""
Constructs a Differentiable Optimizer model from a MOI Optimizer.
Supports `forward` and `backward` methods for solving and differentiating the model respectectively.

## Note
Currently supports differentiating linear and quadratic programs only.
"""

"""
    diff_optimizer(optimizer_constructor)::Optimizer

Creates a `DiffOpt.Optimizer`, which is an MOI layer with an internal optimizer
and other utility methods. Results (primal, dual and slack values) are obtained
by querying the internal optimizer instantiated using the
`optimizer_constructor`. These values are required for find jacobians with respect to problem data.

One define a differentiable model by using any solver of choice. Example:

```julia
julia> import DiffOpt, GLPK

julia> model = DiffOpt.diff_optimizer(GLPK.Optimizer)
julia> model.add_variable(x)
julia> model.add_constraint(...)

julia> _backward_quad(model)  # for convex quadratic models

julia> _backward_quad(model)  # for convex conic models
```
"""
function diff_optimizer(optimizer_constructor)::Optimizer
    return Optimizer(MOI.instantiate(optimizer_constructor, with_bridge_type=Float64))
end

Base.@kwdef struct QPCache
    problem_data::Tuple{
        SparseArrays.SparseMatrixCSC{Float64,Int64}, Vector{Float64}, # Q, q
        SparseArrays.SparseMatrixCSC{Float64,Int64}, Vector{Float64}, # G, h
        SparseArrays.SparseMatrixCSC{Float64,Int64}, Vector{Float64}, # A, b
        Int, Vector{VI}, # nz, var_list
        Int, Vector{MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}}}, # nineq_le, le_con_idx
        Int, Vector{MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}}}, # nineq_ge, ge_con_idx
        Int, Vector{MOI.ConstraintIndex{MOI.VariableIndex, MOI.LessThan{Float64}}}, # nineq_sv_le, le_con_sv_idx
        Int, Vector{MOI.ConstraintIndex{MOI.VariableIndex, MOI.GreaterThan{Float64}}}, # nineq_sv_ge, ge_con_sv_idx
        Int, Vector{MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}}}, # neq, eq_con_idx
        Int, Vector{MOI.ConstraintIndex{MOI.VariableIndex, MOI.EqualTo{Float64}}}, # neq_sv, eq_con_sv_idx
    }
    inequality_duals::Vector{Float64}
    equality_duals::Vector{Float64}
    var_primals::Vector{Float64}
    lhs::SparseMatrixCSC{Float64, Int}
    index_map::MOIU.IndexMap
end

Base.@kwdef struct ConicCache
    M::SparseMatrixCSC{Float64, Int}
    vp::Vector
    Dπv::BlockDiagonals.BlockDiagonal{Float64, Matrix{Float64}}
    xys::NTuple{3, Vector{Float64}}
    A::SparseMatrixCSC{Float64, Int}
    b::Vector{Float64}
    c::Vector{Float64}
    index_map::MOIU.IndexMap
    cones::ProductOfSets{Float64}
end

const CACHE_TYPE = Union{
    Nothing,
    QPCache,
    ConicCache,
}

Base.@kwdef struct QPForwBackCache
    dz::Vector{Float64}
    dλ::Vector{Float64}
    dν::Vector{Float64}
end
Base.@kwdef struct ConicForwCache
    du::Vector{Float64}
    dv::Vector{Float64}
    dw::Vector{Float64}
end
Base.@kwdef struct ConicBackCache
    g::Vector{Float64}
    πz::Vector{Float64}
end

const CACHE_FORW_TYPE = Union{
    Nothing,
    QPForwBackCache,
    ConicForwCache,
}
const CACHE_BACK_TYPE = Union{
    Nothing,
    QPForwBackCache,
    ConicBackCache,
}

const MOIDD = MOI.Utilities.DoubleDicts

Base.@kwdef mutable struct DiffInputCache
    dx::Dict{VI, Float64} = Dict{VI, Float64}()# dz for QP
    # ds
    # dy #= [d\lambda, d\nu] for QP
    # FIXME Would it be possible to have a DoubleDict where the value depends
    #       on the function type ? Here, we need to create two dicts to have
    #       concrete value types.
    # `scalar_constraints` and `vector_constraints` includes `A` and `b` for CPs
    # or `G` and `h` for QPs
    scalar_constraints::MOIDD.DoubleDict{MOI.ScalarAffineFunction{Float64}} = MOIDD.DoubleDict{MOI.ScalarAffineFunction{Float64}}() # also includes G for QPs
    vector_constraints::MOIDD.DoubleDict{MOI.VectorAffineFunction{Float64}} = MOIDD.DoubleDict{MOI.VectorAffineFunction{Float64}}() # also includes G for QPs
    objective::Union{Nothing,MOI.AbstractScalarFunction} = nothing
end

"""
    ForwardInObjective <: MOI.AbstractModelAttribute

A `MOI.AbstractModelAttribute` to set input data to forward differentiation, that
is, problem input data.
The possible values are any `MOI.AbstractScalarFunction`.
A `MOI.ScalarQuadraticFunction` can only be used in linearly constrained
quadratic models.

For instance, if the objective contains `θ * (x + 2y)`, for the purpose of
computing the derivative with respect to `θ`, the following should be set:
```julia
MOI.set(model, DiffOpt.ForwardInObjective(), 1.0 * x + 2.0 * y)
```
where `x` and `y` are the relevant `MOI.VariableIndex`.
"""
struct ForwardInObjective <: MOI.AbstractModelAttribute end

"""
    ForwardInConstraint <: MOI.AbstractConstraintAttribute

A `MOI.AbstractConstraintAttribute` to set input data to forward differentiation, that
is, problem input data.

For instance, if the scalar constraint of index `ci` contains `θ * (x + 2y) <= 5θ`,
for the purpose of computing the derivative with respect to `θ`, the following
should be set:
```julia
MOI.set(model, DiffOpt.ForwardInConstraint(), ci, 1.0 * x + 2.0 * y - 5.0)
```
Note that we use `-5` as the `ForwardInConstraint` sets the tangent of the
ConstraintFunction so we consider the expression `θ * (x + 2y - 5)`.
"""
struct ForwardInConstraint <: MOI.AbstractConstraintAttribute end

"""
    ForwardOutVariablePrimal <: MOI.AbstractVariableAttribute

A `MOI.AbstractVariableAttribute` to get output data from forward
differentiation, that is, problem solution.

For instance, to get the tangent of the variable of index `vi` corresponding to
the tangents given to `ForwardInObjective` and `ForwardInConstraint`, do the
following:
```julia
MOI.get(model, DiffOpt.ForwardOutVariablePrimal(), vi)
```
"""
struct ForwardOutVariablePrimal <: MOI.AbstractVariableAttribute end
MOI.is_set_by_optimize(::ForwardOutVariablePrimal) = true

"""
    BackwardInVariablePrimal <: MOI.AbstractVariableAttribute

A `MOI.AbstractVariableAttribute` to set input data to backward
differentiation, that is, problem solution.

For instance, to set the tangent of the variable of index `vi`, do the
following:
```julia
MOI.set(model, DiffOpt.BackwardInVariablePrimal(), x)
```
"""
struct BackwardInVariablePrimal <: MOI.AbstractVariableAttribute end

"""
    BackwardOutObjective <: MOI.AbstractModelAttribute

A `MOI.AbstractModelAttribute` to get output data to backward differentiation,
that is, problem input data.

For instance, to get the tangent of the objective function corresponding to
the tangent given to `BackwardInVariablePrimal`, do the
following:
```julia
func = MOI.get(model, DiffOpt.BackwardOutObjective)
```
Then, to get the sensitivity of the linear term with variable `x`, do
```julia
JuMP.coefficient(func, x)
```
To get the sensitivity with respect to the quadratic term with variables `x`
and `y`, do either
```julia
JuMP.coefficient(func, x, y)
```
or
```julia
DiffOpt.quad_sym_half(func, x, y)
```

!!! warning
    These two lines are **not** equivalent in case `x == y`, see
    [`quad_sym_half`](@ref) for the details on the difference between these two
    functions.
"""
struct BackwardOutObjective <: MOI.AbstractModelAttribute end
MOI.is_set_by_optimize(::BackwardOutObjective) = true

"""
    BackwardOutConstraint

An `MOI.AbstractConstraintAttribute` to get output data to backward differentiation, that
is, problem input data.

For instance, if the following returns `x + 2y + 5`, it means that the tangent
has coordinate `1` for the coefficient of `x`, coordinate `2` for the
coefficient of `y` and `5` for the function constant.
If the constraint is of the form `func == constant` or `func <= constant`,
the tangent for the constant on the right-hand side is `-5`.
```julia
MOI.get(model, DiffOpt.BackwardOutConstraint(), ci)
```
"""
struct BackwardOutConstraint <: MOI.AbstractConstraintAttribute end
MOI.is_set_by_optimize(::BackwardOutConstraint) = true

"""
    @enum ProgramClassCode QUADRATIC CONIC AUTOMATIC

Program class used by DiffOpt. DiffOpt implements differentiation of two
different program class:
1) Quadratic Program (QP): quadratic objective and linear constraints and
2) Conic Program (CP): linear objective and conic constraints.

`AUTOMATIC` which means that the class will be automatically selected given the
problem data: if any constraint is conic, CP is used and QP is used otherwise.
See [`ProgramClass`](@ref).
"""
@enum ProgramClassCode QUADRATIC CONIC AUTOMATIC

mutable struct Optimizer{OT <: MOI.ModelLike} <: MOI.AbstractOptimizer
    optimizer::OT

    program_class::ProgramClassCode

    # storage for problem data in matrix form
    # includes maps from matrix indices to problem data held in `optimizer`
    # also includes KKT matrices
    # also includes the solution
    gradient_cache::CACHE_TYPE

    # caches for sensitivity output
    # result from solving KKT/residualmap linear systems
    # this allows keeping the same `gradient_cache`
    # if only sensitivy input changes
    forw_grad_cache::CACHE_FORW_TYPE
    back_grad_cache::CACHE_BACK_TYPE

    # sensitivity input cache using MOI like sparse format
    input_cache::DiffInputCache

    function Optimizer(optimizer::OT) where {OT <: MOI.ModelLike}
        new{OT}(
            optimizer,
            AUTOMATIC,
            nothing,
            nothing,
            nothing,
            DiffInputCache(),
        )
    end
end

"""
    ProgramClass <: MOI.AbstractOptimizerAttribute

Determines which program class to used from [`ProgramClassCode`](@ref). The
default is `AUTOMATIC`.

One important advantage of setting the class explicitly is that it will allow
necessary bridges to be used. If the class is `AUTOMATIC` then
`DiffOpt.Optimizer` will report that it supports both objective and constraints
of the QP and CP classes. For instance, it will reports that is supports both
quadratic objective and conic constraints. However, at the differentiation
stage, we won't be able to differentiate since QP does not support conic
constraints and CP does not support quadratic objective. On the other hand, if
the `ProgramClass` is set to `CONIC` then `DiffOpt.Optimizer` will report that
it does not support quadratic objective hence it will be bridged to second-order
cone constraints and we will be able to use CP to differentiate.
"""
struct ProgramClass <: MOI.AbstractOptimizerAttribute end

MOI.supports(::Optimizer, ::ProgramClass) = true
MOI.get(model::Optimizer, ::ProgramClass) = model.program_class
function MOI.set(model::Optimizer, ::ProgramClass, class::ProgramClassCode)
    model.program_class = class
end

"""
    ProgramClassUsed <: MOI.AbstractOptimizerAttribute

Program class actually used, same as [`ProgramClass`](@ref) except that it does
not return `AUTOMATIC` but the class automatically chosen instead. This attribute
is read-only, it cannot be set, set [`ProgramClass`](@ref) instead.
"""
struct ProgramClassUsed <: MOI.AbstractOptimizerAttribute end

function MOI.get(model::Optimizer, ::ProgramClassUsed)
    if model.program_class == AUTOMATIC
        if _qp_supported(model.optimizer)
            return QUADRATIC
        else
            return CONIC
        end
    else
        return model.program_class
    end
end

function MOI.get(model::Optimizer, ::ForwardOutVariablePrimal, vi::VI)
    return _get_dx(model, vi)
end
_get_dx(model::Optimizer, vi) = _get_dx(model.forw_grad_cache, model.gradient_cache, vi)
function _get_dx(cache::QPForwBackCache, g_cache::QPCache, vi)
    i = g_cache.index_map[vi].value
    return cache.dz[i]
end
function _get_dx(f_cache::ConicForwCache, g_cache::ConicCache, vi)
    i = g_cache.index_map[vi].value
    du = f_cache.du
    dw = f_cache.dw
    x = g_cache.xys[1]
    return - (du[i] - x[i] * dw[])
end

function MOI.get(model::Optimizer, ::BackwardInVariablePrimal, vi::VI)
    return get(model.input_cache.dx, vi, 0.0)
end
function MOI.set(model::Optimizer, ::BackwardInVariablePrimal, vi::VI, val)
    model.input_cache.dx[vi] = val
    return
end

function lazy_combination(op::F, α, a, β, b) where {F<:Function}
    return LazyArrays.ApplyArray(
        op,
        LazyArrays.@~(α .* a),
        LazyArrays.@~(β .* b),
    )
end
# Workaround for Julia v1.0
@static if VERSION < v"1.6"
    _view(x, I) = x[I]
else
    _view(x, I) = view(x, I)
end
function lazy_combination(op::F, α, a, β, b, I::UnitRange) where {F<:Function}
    return lazy_combination(op, α, _view(a, I), β, _view(b, I))
end
function lazy_combination(op::F, a, b, i::Integer, args::Vararg{Any,N}) where {F<:Function,N}
    return lazy_combination(op, a[i], b, b[i], a, args...)
end
function lazy_combination(op::F, a, b, i::UnitRange, I::UnitRange) where {F<:Function}
    return lazy_combination(op, _view(a, i), b, _view(b, i), a, I)
end

function MOI.get(model::Optimizer, ::BackwardOutObjective)
    return IndexMappedFunction(
        _back_obj(model.back_grad_cache, model.gradient_cache),
        model.gradient_cache.index_map,
    )
end
function _back_obj(b_cache::ConicBackCache, g_cache::ConicCache)
    g = b_cache.g
    πz = b_cache.πz
    dc = lazy_combination(-, πz, g, length(g))
    return VectorScalarAffineFunction(dc, 0.0)
end
function _back_obj(b_cache::QPForwBackCache, g_cache::QPCache)
    ∇z = b_cache.dz
    z = g_cache.var_primals
    # `∇z * z' + z * ∇z'` doesn't work, see
    # https://github.com/JuliaArrays/LazyArrays.jl/issues/178
    dQ = LazyArrays.@~ (∇z .* z' + z .* ∇z') / 2.0
    return MatrixScalarQuadraticFunction(
        VectorScalarAffineFunction(b_cache.dz, 0.0),
        dQ,
    )
end

function MOI.get(model::Optimizer, ::ForwardInObjective)
    return model.input_cache.objective
end
function MOI.set(model::Optimizer, ::ForwardInObjective, objective)
    model.input_cache.objective = objective
    return
end

_lazy_affine(vector, constant::Number) = VectorScalarAffineFunction(vector, constant)
_lazy_affine(matrix, vector) = MatrixVectorAffineFunction(matrix, vector)
function MOI.get(model::Optimizer, ::BackwardOutConstraint, ci::CI)
    return IndexMappedFunction(
        _lazy_affine(_get_dA(model, ci), _get_db(model, ci)),
        model.gradient_cache.index_map,
    )
end
_get_db(model::Optimizer, ci) = _get_db(model.back_grad_cache, model.gradient_cache, ci)
function _get_db(b_cache::ConicBackCache, g_cache::ConicCache, ci::CI{F,S}
) where {F<:MOI.AbstractVectorFunction,S}
    cf = g_cache.cones
    _ci = g_cache.index_map[ci]
    i = MOI.Utilities.rows(cf, _ci) # vector
    # i = g_cache.index_map[ci].value
    (x, _, _) = g_cache.xys
    n = length(x) # columns in A
    # db = - dQ[n+1:n+m, end] + dQ[end, n+1:n+m]'
    g = b_cache.g
    πz = b_cache.πz
    return lazy_combination(-, πz, g, length(g), n .+ i)
end
function _get_db(b_cache::ConicBackCache, g_cache::ConicCache, ci::CI{F,S}
) where {F<:MOI.AbstractScalarFunction,S}
    i = g_cache.index_map[ci].value
    (x, _, _) = g_cache.xys
    n = length(x) # columns in A
    # db = - dQ[n+1:n+m, end] + dQ[end, n+1:n+m]'
    g = b_cache.g
    πz = b_cache.πz
    dQ_ni_end = - g[n+i] * πz[end]
    dQ_end_ni = - g[end] * πz[n+i]
    return - dQ_ni_end + dQ_end_ni
end
_neg_if_gt(x, ::Type{<:MOI.LessThan}) = x
_neg_if_gt(x, ::Type{<:MOI.GreaterThan}) = -x
function _get_db(b_cache::QPForwBackCache, g_cache::QPCache, ci::CI{F,S}
) where {F,S}
    i = g_cache.index_map[ci].value
    # dh = -Diagonal(λ) * dλ
    dλ = b_cache.dλ
    λ = g_cache.inequality_duals
    return _neg_if_gt(λ[i] * dλ[i], S)
end
function _get_db(b_cache::QPForwBackCache, g_cache::QPCache, ci::CI{F,S}
) where {F,S<:MOI.EqualTo}
    i = g_cache.index_map[ci].value
    dν = b_cache.dν
    return dν[i]
end

function MOI.get(model::Optimizer,
    ::ForwardInConstraint, ci::CI{MOI.ScalarAffineFunction{T},S}
) where {T,S}
    return get(model.input_cache.scalar_constraints, ci, zero(MOI.ScalarAffineFunction{T}))
end
function MOI.get(model::Optimizer,
    ::ForwardInConstraint, ci::CI{MOI.VectorAffineFunction{T},S}
) where {T,S}
    func = get(model.input_cache.vector_constraints, ci, nothing)
    if func === nothing
        set = MOI.get(model, MOI.ConstraintSet(), ci)
        dim = MOI.dimension(set)
        return MOI.Utilities.zero_with_output_dimension(MOI.VectorAffineFunction{T}, dim)
    else
        return func
    end
end
function MOI.set(model::Optimizer,
    ::ForwardInConstraint,
    ci::CI{MOI.ScalarAffineFunction{T},S},
    func::MOI.ScalarAffineFunction{T},
) where {T,S}
    model.input_cache.scalar_constraints[ci] = func
    return
end
function MOI.set(model::Optimizer,
    ::ForwardInConstraint,
    ci::CI{MOI.VectorAffineFunction{T},S},
    func::MOI.VectorAffineFunction{T},
) where {T,S}
    model.input_cache.vector_constraints[ci] = func
    return
end

_get_dA(model::Optimizer, ci) = _get_dA(model.back_grad_cache, model.gradient_cache, ci)
function _get_dA(b_cache::ConicBackCache, g_cache::ConicCache, ci::CI{F,S}
) where {F<:MOI.AbstractScalarFunction,S}
    j = g_cache.index_map[vi].value
    i = g_cache.index_map[ci].value
    (x, y, _) = g_cache.xys
    n = length(x) # columns in A
    m = length(y) # lines in A
    # dA = - dQ[1:n, n+1:n+m]' + dQ[n+1:n+m, 1:n]
    g = b_cache.g
    πz = b_cache.πz
    return lazy_combination(-, g, πz, i, n .+ (1:n))
end
function _get_dA(b_cache::ConicBackCache, g_cache::ConicCache, ci::CI{F,S}
) where {F<:MOI.AbstractVectorFunction,S}
    cf = g_cache.cones
    _ci = g_cache.index_map[ci]
    i = MOI.Utilities.rows(cf, _ci) # vector
    # i = g_cache.index_map[ci].value
    (x, y, _) = g_cache.xys
    n = length(x) # columns in A
    m = length(y) # lines in A
    # dA = - dQ[1:n, n+1:n+m]' + dQ[n+1:n+m, 1:n]
    g = b_cache.g
    πz = b_cache.πz
    return lazy_combination(-, g, πz, i, n .+ (1:n))
end
# quadratic matrix indexes are split by type either == or (<=/>=)
function _get_dA(b_cache::QPForwBackCache, g_cache::QPCache, ci::CI{F,S}
) where {F, S<:MOI.EqualTo}
    i = g_cache.index_map[ci].value
    z = g_cache.var_primals
    dz = b_cache.dz
    ν = g_cache.equality_duals
    dν = b_cache.dν
    return lazy_combination(+, dν[i], z, ν[i], dz)
end
function _get_dA(b_cache::QPForwBackCache, g_cache::QPCache, ci::CI{F,S}
) where {F,S}
    i = g_cache.index_map[ci].value
    z = g_cache.var_primals
    dz = b_cache.dz
    λ = g_cache.inequality_duals
    dλ = b_cache.dλ
    l = _neg_if_gt(λ[i], S)
    return lazy_combination(+, l * dλ[i], z, l * λ[i], dz)
end

function MOI.optimize!(model::Optimizer)
    model.gradient_cache = nothing
    MOI.optimize!(model.optimizer)

    # do not fail. interferes with MOI.Tests.linear12test
    if !in(MOI.get(model.optimizer, MOI.TerminationStatus()),  (MOI.LOCALLY_SOLVED, MOI.OPTIMAL))
        @warn "problem status: $(MOI.get(model.optimizer, MOI.TerminationStatus()))"
        return
    end

    return
end

"""
    backward(model::Optimizer)

Wrapper method for the backward pass.
This method will consider as input a currently solved problem and differentials
with respect to the solution set with the [`BackwardInVariablePrimal`](@ref) attribute.
The output problem data differentials can be queried with the
attributes [`BackwardOutObjective`](@ref) and [`BackwardOutConstraint`](@ref).
"""
function backward(model::Optimizer)
    if MOI.get(model, ProgramClassUsed()) == QUADRATIC
        return _backward_quad(model)
    else
        return _backward_conic(model)
    end
end

"""
    forward(model::Optimizer)

Wrapper method for the forward pass.
This method will consider as input a currently solved problem and
differentials with respect to problem data set with
the [`ForwardInObjective`](@ref) and  [`ForwardInConstraint`](@ref) attributes.
The output solution differentials can be queried with the attribute
[`ForwardOutVariablePrimal`](@ref).
"""
function forward(model::Optimizer)
    if MOI.get(model, ProgramClassUsed()) == QUADRATIC
        return _forward_quad(model)
    else
        return _forward_conic(model)
    end
end

"""
    _backward_quad(model::Optimizer)

Method to differentiate optimal solution `z` and return
product of jacobian matrices (`dz / dQ`, `dz / dq`, etc) with
the backward pass vector `dl / dz`

The method computes the product of
1. jacobian of problem solution `z*` with respect to
    problem parameters set with the [`BackwardInVariablePrimal`](@ref)
2. a backward pass vector `dl / dz`, where `l` can be a loss function

Note that this method *does not returns* the actual jacobians.

For more info refer eqn(7) and eqn(8) of https://arxiv.org/pdf/1703.00443.pdf
"""
function _backward_quad(model::Optimizer)

    if model.gradient_cache === nothing
        build_quad_diff_cache!(model)
    end
    (
        Q, q, G, h, A, b, nz, var_list,
        nineq_le, le_con_idx,
        nineq_ge, ge_con_idx,
        nineq_sv_le, le_con_sv_idx,
        nineq_sv_ge, ge_con_sv_idx,
        neq, eq_con_idx,
        neq_sv, eq_con_sv_idx,
    ) = model.gradient_cache.problem_data
    z = model.gradient_cache.var_primals
    λ = model.gradient_cache.inequality_duals
    ν = model.gradient_cache.equality_duals
    LHS = model.gradient_cache.lhs

    index_map = model.gradient_cache.index_map
    dl_dz = zeros(length(z))
    for (vi, val) in model.input_cache.dx
        inner_index = index_map[vi].value
        dl_dz[inner_index] = val
    end

    nineq_total = nineq_le + nineq_ge + nineq_sv_le + nineq_sv_ge
    RHS = [dl_dz; zeros(neq + neq_sv + nineq_total)]

    partial_grads = if norm(Q) ≈ 0
        -lsqr(LHS, RHS)
    else
        -LHS \ RHS
    end

    dz = partial_grads[1:nz]
    dλ = partial_grads[nz+1:nz+nineq_total]
    dν = partial_grads[nz+nineq_total+1:nz+nineq_total+neq+neq_sv]

    model.back_grad_cache = QPForwBackCache(dz, dλ, dν)
    return nothing
    # dQ = 0.5 * (dz * z' + z * dz')
    # dq = dz
    # dG = Diagonal(λ) * (dλ * z' + λ * dz') # was: Diagonal(λ) * dλ * z' - λ * dz')
    # dh = -Diagonal(λ) * dλ
    # dA = dν * z'+ ν * dz' # was: dν * z' - ν * dz'
    # db = -dν
    # todo, check MOI signs for dA and dG
end

"""
    _forward_quad(model::Optimizer)
"""
function _forward_quad(model::Optimizer)
    if model.gradient_cache === nothing
        build_quad_diff_cache!(model)
    end
    (
        Q, q, G, h, A, b, nz, var_list,
        nineq_le, le_con_idx,
        nineq_ge, ge_con_idx,
        nineq_sv_le, le_con_sv_idx,
        nineq_sv_ge, ge_con_sv_idx,
        neq, eq_con_idx,
        neq_sv, eq_con_sv_idx,
    ) = model.gradient_cache.problem_data
    z = model.gradient_cache.var_primals
    λ = model.gradient_cache.inequality_duals
    ν = model.gradient_cache.equality_duals
    LHS = model.gradient_cache.lhs
    index_map = model.gradient_cache.index_map

    objective_function = _convert(MOI.ScalarQuadraticFunction{Float64}, model.input_cache.objective)
    sparse_array_obj = sparse_array_representation(objective_function, LinearAlgebra.checksquare(Q), index_map)
    dQ = sparse_array_obj.quadratic_terms
    dq = sparse_array_obj.affine_terms

    # The user sets the constraint function in the sense `func`-in-`set` while
    # `db` and `dh` corresponds to the tangents of the set constants. Therefore,
    # we should multiply the constant by `-1`. For `GreaterThan`, we needed to
    # multiply by `-1` to transform it to `LessThan` so it cancels out.
    db = zeros(length(b))
    _fill(isequal(MOI.EqualTo{Float64}), (::Type{MOI.EqualTo{Float64}}) -> true, model, _QPSets(), db)
    dh = zeros(length(h))
    _fill(!isequal(MOI.EqualTo{Float64}), !isequal(MOI.GreaterThan{Float64}), model, _QPSets(), dh)

    nz = nnz(A)
    (lines, cols) = size(A)
    dAi = zeros(Int, 0)
    dAj = zeros(Int, 0)
    dAv = zeros(Float64, 0)
    sizehint!(dAi, nz)
    sizehint!(dAj, nz)
    sizehint!(dAv, nz)
    _fill(isequal(MOI.EqualTo{Float64}), isequal(MOI.GreaterThan{Float64}), model, _QPSets(), dAi, dAj, dAv)
    dA = sparse(dAi, dAj, dAv, lines, cols)

    nz = nnz(G)
    (lines, cols) = size(G)
    dGi = zeros(Int, 0)
    dGj = zeros(Int, 0)
    dGv = zeros(Float64, 0)
    sizehint!(dGi, nz)
    sizehint!(dGj, nz)
    sizehint!(dGv, nz)
    _fill(!isequal(MOI.EqualTo{Float64}), isequal(MOI.GreaterThan{Float64}), model, _QPSets(), dGi, dGj, dGv)
    dG = sparse(dGi, dGj, dGv, lines, cols)

    RHS = [
        dQ * z + dq + dG' * λ + dA' * ν
        λ .* (dG * z) - λ .* dh
        dA * z - db
    ]

    partial_grads = if norm(Q) ≈ 0
        -lsqr(LHS', RHS)
    else
        -_linsolve(LHS', RHS)
    end

    nv = length(z)
    nineq_total = nineq_le + nineq_ge + nineq_sv_le + nineq_sv_ge
    dz = partial_grads[1:nv]
    dλ = partial_grads[nv+1:nv+nineq_total]
    dν = partial_grads[nv+nineq_total+1:nv+nineq_total+neq+neq_sv]

    model.forw_grad_cache = QPForwBackCache(dz, dλ, dν)
    return nothing
end

_linsolve(A, b) = A \ b
# See https://github.com/JuliaLang/julia/issues/32668
_linsolve(A, b::SparseVector) = A \ Vector(b)

"""
    π(v::Vector{Float64}, model::MOI.ModelLike, cones::ProductOfSets, index_map::MOIU.IndexMap)

Given a `model`, its `cones` and the `index_map` from the indices of
`model` to the indices of `cones`, find the projection of the vectors `v`
of length equal to the number of rows in the conic form onto the cartesian
product of the cones corresponding to these rows.
For more info, refer to https://github.com/matbesancon/MathOptSetDistances.jl
"""
function π(v::Vector{T}, model::MOI.ModelLike, cones::ProductOfSets, index_map::MOIU.IndexMap) where T
    return map_rows(model, cones, index_map, Flattened{T}()) do ci, r
        MOSD.projection_on_set(
            MOSD.DefaultDistance(),
            v[r],
            MOI.dual_set(MOI.get(model, MOI.ConstraintSet(), ci))
        )
    end
end

"""
    Dπ(v::Vector{Float64}, model, cones::ProductOfSets, index_map::MOIU.IndexMap)

Given a `model`, its `cones` and the `index_map` from the indices of
`model` to the indices of `cones`, find the gradient of the projection of
the vectors `v` of length equal to the number of rows in the conic form onto the
cartesian product of the cones corresponding to these rows.
For more info, refer to https://github.com/matbesancon/MathOptSetDistances.jl
"""
function Dπ(v::Vector{T}, model::MOI.ModelLike, cones::ProductOfSets, index_map::MOIU.IndexMap) where T
    return BlockDiagonals.BlockDiagonal(
        map_rows(model, cones, index_map, Nested{Matrix{T}}()) do ci, r
            MOSD.projection_gradient_on_set(
                MOSD.DefaultDistance(),
                v[r],
                MOI.dual_set(MOI.get(model, MOI.ConstraintSet(), ci)),
            )
        end
    )
end

# See the docstring of `map_rows`.
struct Nested{T} end
struct Flattened{T} end

# Store in `x` the values `y` corresponding to the rows `r` and the `k`th
# constraint.
function _assign_mapped!(x, y, r, k, ::Nested)
    x[k] = y
end
function _assign_mapped!(x, y, r, k, ::Flattened)
    x[r] = y
end

# Map the rows corresponding to `F`-in-`S` constraints and store it in `x`.
function _map_rows!(f::Function, x::Vector, model, cones::ProductOfSets, index_map::MOIU.DoubleDicts.IndexDoubleDictInner{F, S}, map_mode, k) where {F, S}
    for ci in MOI.get(model, MOI.ListOfConstraintIndices{F, S}())
        r = MOI.Utilities.rows(cones, index_map[ci])
        k += 1
        _assign_mapped!(x, f(ci, r), r, k, map_mode)
    end
    return k
end

# Allocate a vector for storing the output of `map_rows`.
_allocate_rows(cones, ::Nested{T}) where {T} = Vector{T}(undef, length(cones.dimension))
_allocate_rows(cones, ::Flattened{T}) where {T} = Vector{T}(undef, MOI.dimension(cones))

"""
    map_rows(f::Function, model, cones::ProductOfSets, index_map::MOIU.IndexMap, map_mode::Union{Nested{T}, Flattened{T}})

Given a `model`, its `cones`, the `index_map` from the indices of `model`
to the indices of `cones` and `map_mode` of type `Nested` (resp.
`Flattened`), return a `Vector{T}` of length equal to the number of cones (resp.
rows) in the conic form where the value for the index (resp. rows) corresponding
to each cone is equal to `f(ci, r)` where `ci` is the corresponding constraint
index in `model` and `r` is a `UnitRange` of the corresponding rows in the conic
form.
"""
function map_rows(f::Function, model, cones::ProductOfSets, index_map::MOIU.IndexMap, map_mode::Union{Nested, Flattened})
    x = _allocate_rows(cones, map_mode)
    k = 0
    for (F, S) in MOI.get(model, MOI.ListOfConstraintTypesPresent())
        # Function barrier for type unstability of `F` and `S`
        # `con_map` is a `MOIU.DoubleDicts.MainIndexDoubleDict`, we index it at `F, S`
        # which returns a `MOIU.DoubleDicts.IndexWithType{F, S}` which is type stable.
        # If we have a small number of different constraint types and many
        # constraint of each type, this mostly removes type unstabilities
        # as most the time is in `_map_rows!` which is type stable.
        k = _map_rows!(f, x, model, cones, index_map.con_map[F, S], map_mode, k)
    end
    return x
end

function _check_termination_status(model::Optimizer)
    if !in(
        MOI.get(model, MOI.TerminationStatus()), (MOI.LOCALLY_SOLVED, MOI.OPTIMAL)
        )
        error("problem status: ", MOI.get(model.optimizer, MOI.TerminationStatus()))
    end
end

"""
    _forward_conic(model::Optimizer)

Method to compute the product of the derivative (Jacobian) at the
conic program parameters `A`, `b`, `c`  to the perturbations `dA`, `db`, `dc`.
This is similar to [`forward`](@ref).

For theoretical background, refer Section 3 of Differentiating Through a Cone Program, https://arxiv.org/abs/1904.09043
"""
function _forward_conic(model::Optimizer)
    _check_termination_status(model)

    if model.gradient_cache === nothing
        build_conic_diff_cache!(model)
    end

    M = model.gradient_cache.M
    vp = model.gradient_cache.vp
    Dπv = model.gradient_cache.Dπv
    (x, y, s) = model.gradient_cache.xys
    A = model.gradient_cache.A
    b = model.gradient_cache.b
    c = model.gradient_cache.c
    index_map = model.gradient_cache.index_map

    objective_function = _convert(MOI.ScalarAffineFunction{Float64}, model.input_cache.objective)
    sparse_array_obj = sparse_array_representation(objective_function, length(c), index_map)
    dc = sparse_array_obj.terms

    db = zeros(length(b))
    _fill(S -> false, model, model.gradient_cache.cones, db)
    (lines, cols) = size(A)
    nz = nnz(A)
    dAi = zeros(Int, 0)
    dAj = zeros(Int, 0)
    dAv = zeros(Float64, 0)
    sizehint!(dAi, nz)
    sizehint!(dAj, nz)
    sizehint!(dAv, nz)
    _fill(S -> false, model, model.gradient_cache.cones, dAi, dAj, dAv)
    dA = sparse(dAi, dAj, dAv, lines, cols)

    m = size(A, 1)
    n = size(A, 2)
    N = m + n + 1
    # NOTE: w = 1 systematically since we asserted the primal-dual pair is optimal
    (u, v, w) = (x, y - s, 1.0)

    # g = dQ * Π(z/|w|) = dQ * [u, vp, 1.0]
    RHS = [dA' * vp + dc; -dA * u + db; -dc ⋅ u - db ⋅ vp]

    dz = if norm(RHS) <= 1e-400 # TODO: parametrize or remove
        RHS .= 0 # because M is square
    else
        lsqr(M, RHS)
    end

    du, dv, dw = dz[1:n], dz[n+1:n+m], dz[n+m+1]
    model.forw_grad_cache = ConicForwCache(du, dv, [dw])
    return nothing
    # dx = du - x * dw
    # dy = Dπv * dv - y * dw
    # ds = Dπv * dv - dv - s * dw
    # return -dx, -dy, -ds
end

# Just a hack that will be removed once we use `MOIU.MatrixOfConstraints`
struct _QPSets end
MOI.Utilities.rows(::_QPSets, ci::MOI.ConstraintIndex) = ci.value

function _fill(neg::Function, model::Optimizer, cones, args...)
    _fill(S -> true, neg, model, cones, args...)
end
function _fill(filter::Function, neg::Function, model::Optimizer, cones, args...)
    con_map = model.gradient_cache.index_map.con_map
    var_map = model.gradient_cache.index_map.var_map
    for (F, S) in MOI.get(model, MOI.ListOfConstraintTypesPresent())
        filter(S) || continue
        if F == MOI.ScalarAffineFunction{Float64}
            _fill(args..., neg(S), cones, con_map[F,S], var_map, model.input_cache.scalar_constraints[F,S])
        elseif F == MOI.VectorAffineFunction{Float64}
            _fill(args..., neg(S), cones, con_map[F,S], var_map, model.input_cache.vector_constraints[F,S])
        end
    end
    return
end

function _fill(vector::Vector, neg::Bool, cones, constraint_map, variable_map, dict)
    for (ci, func) in dict
        r = MOI.Utilities.rows(cones, constraint_map[ci])
        vector[r] = neg ? -MOI.constant(func) : MOI.constant(func)
    end
end
function _fill(I::Vector, J::Vector, V::Vector, neg::Bool, cones, constraint_map, variable_map, dict)
    for (ci, func) in dict
        r = MOI.Utilities.rows(cones, constraint_map[ci])
        for term in func.terms
            _push_term(I, J, V, neg, r, term, variable_map)
        end
    end
end
function _push_term(I::Vector, J::Vector, V::Vector, neg::Bool, r::Integer, term::MOI.ScalarAffineTerm, variable_map)
    push!(I, r)
    push!(J, variable_map[term.variable].value)
    push!(V, neg ? -term.coefficient : term.coefficient)
end
function _push_term(I::Vector, J::Vector, V::Vector, neg::Bool, r::UnitRange, term::MOI.VectorAffineTerm, variable_map)
    _push_term(I, J, V, neg, r[term.output_index], term.scalar_term, variable_map)
end

"""
    _backward_conic(model::Optimizer, dx::Vector{Float64}, dy::Vector{Float64}, ds::Vector{Float64})

Method to compute the product of the transpose of the derivative (Jacobian) at the
conic program parameters `A`, `b`, `c`  to the perturbations `dx`, `dy`, `ds`.
This is similar to [`backward`](@ref).

For theoretical background, refer Section 3 of Differentiating Through a Cone Program, https://arxiv.org/abs/1904.09043
"""
function _backward_conic(model::Optimizer)
    _check_termination_status(model)

    if model.gradient_cache === nothing
        build_conic_diff_cache!(model)
    end

    M = model.gradient_cache.M
    vp = model.gradient_cache.vp
    Dπv = model.gradient_cache.Dπv
    (x, y, s) = model.gradient_cache.xys
    A = model.gradient_cache.A
    b = model.gradient_cache.b
    c = model.gradient_cache.c

    index_map = model.gradient_cache.index_map
    dx = zeros(length(c))
    for (vi, val) in model.input_cache.dx
        inner_index = index_map[vi].value
        dx[inner_index] = val
    end
    dy = zeros(length(b))
    ds = zeros(length(b))

    m = size(A, 1)
    n = size(A, 2)
    N = m + n + 1
    # NOTE: w = 1 systematically since we asserted the primal-dual pair is optimal
    (u, v, w) = (x, y - s, 1.0)

    # dz = D \phi (z)^T (dx,dy,dz)
    dz = [
        dx
        Dπv' * (dy + ds) - ds
        - x' * dx - y' * dy - s' * ds
    ]

    g = if norm(dz) <= 1e-4 # TODO: parametrize or remove
        dz .= 0 # because M is square
    else
        lsqr(M, dz)
    end

    πz = [
        u
        vp
        1.0
    ]

    # TODO: very important
    # contrast with:
    # http://reports-archive.adm.cs.cmu.edu/anon/2019/CMU-CS-19-109.pdf
    # pg 97, cap 7.4.2

    model.back_grad_cache = ConicBackCache(g, πz)
    return nothing
    # dQ = - g * πz'
    # dA = - dQ[1:n, n+1:n+m]' + dQ[n+1:n+m, 1:n]
    # db = - dQ[n+1:n+m, end] + dQ[end, n+1:n+m]'
    # dc = - dQ[1:n, end] + dQ[end, 1:n]'
    # return dA, db, dc
end


text:
@testset "Elastic net ($T)" for T in [Float32, Float64, ComplexF32, ComplexF64]

    using ProximalOperators
    using ProximalAlgorithms
    using LinearAlgebra

    A = T[  1.0  -2.0   3.0  -4.0  5.0;
            2.0  -1.0   0.0  -1.0  3.0;
           -1.0   0.0   4.0  -3.0  2.0;
           -1.0  -1.0  -1.0   1.0  3.0]
    b = T[1.0, 2.0, 3.0, 4.0]

    m, n = size(A)

    R = real(T)

    reg = ElasticNet(R(1), R(1))
    reg1 = NormL1(R(1))
    reg2 = SqrNormL2(R(1))
    loss = Translate(SqrNormL2(R(1)), -b)

    L = opnorm(A)^2

    x_star = T[-0.6004983388704322, 0.0, 0.0, 0.195182724252491, 0.764119601328903]

    @testset "DavisYin" begin

        # with known initial iterate

        x0 = zeros(T, n)
        solver = ProximalAlgorithms.DavisYin{R}(tol=R(1e-6))
        xf_dys, xg_dys, it_dys = solver(
            x0, f=reg1, g=reg2, h=loss, A=A, L=L
        )
        @test eltype(xf_dys) == T
        @test eltype(xg_dys) == T
        @test norm(xf_dys - x_star, Inf) <= 1e-3
        @test norm(xg_dys - x_star, Inf) <= 1e-3
        @test it_dys <= 140

        # with random initial iterate

        x0 = randn(T, n)
        solver = ProximalAlgorithms.DavisYin{R}(tol=R(1e-6))
        xf_dys, xg_dys, it_dys = solver(
            x0, f=reg1, g=reg2, h=loss, A=A, L=L
        )
        @test eltype(xf_dys) == T
        @test eltype(xg_dys) == T
        @test norm(xf_dys - x_star, Inf) <= 1e-3
        @test norm(xg_dys - x_star, Inf) <= 1e-3

    end

    afba_test_params = [
        (R(2), R(0), 130),
        (R(1), R(1), 1890),
        (R(0), R(1), 320),
        (R(0), R(0), 194),
        (R(1), R(0), 130)
    ]

    @testset "AFBA" for (theta, mu, maxit) in afba_test_params

        # with known initial iterates

        x0 = zeros(T, n)
        y0 = zeros(T, m)

        solver = ProximalAlgorithms.AFBA{R}(theta=theta, mu=mu, tol=R(1e-6))
        x_afba, y_afba, it_afba = solver(
            x0, y0, f=reg2, g=reg1, h=loss, L=A, beta_f=R(1),
        )
        @test eltype(x_afba) == T
        @test eltype(y_afba) == T
        @test norm(x_afba - x_star, Inf) <= 1e-4
        @test it_afba <= maxit

        # with random initial iterates

        x0 = randn(T, n)
        y0 = randn(T, m)

        solver = ProximalAlgorithms.AFBA{R}(theta=theta, mu=mu, tol=R(1e-6))
        x_afba, y_afba, it_afba = solver(
            x0, y0, f=reg2, g=reg1, h=loss, L=A, beta_f=R(1),
        )
        @test eltype(x_afba) == T
        @test eltype(y_afba) == T
        @test norm(x_afba - x_star, Inf) <= 1e-4

    end
end


text:
using ACE1pack, Test

@info("Test constructing ace basis and parameters")
ace_basis = basis_params(type="ace", species = :Si, N = 3, maxdeg = 10)
ACE1pack.generate_basis(ace_basis)

@info("Test constructing degree and degree params")
degree = degree_params()
D = ACE1pack.generate_degree(degree)
degreeM = degree_params(type = "sparseM", Dd = Dict("default" => 10, 1 => 7, (2, "H") => 3))
D = ACE1pack.generate_degree(degreeM)

@info("Test constructing tranforms and transform parameters")

# polytransform
polytransform_params = transform_params(r0 = 1.1)
transform_poly = ACE1pack.generate_transform(polytransform_params)

#multitransform
cutoffs = Dict(
    (:C, :C) => (1.5, 2.3),
    (:H, :H) => (0.4, 2.1),
    (:C, :H) => (1.2, 3.4))

transforms = Dict(
    (:C, :C) => transform_params(r0 = 2),
    (:H, :H) => transform_params(r0 = 1.1),
    (:C, :H) => transform_params(r0 = 2))

multitransform_params = transform_params(
    type = "multitransform",
    transforms = transforms,
    cutoffs = cutoffs)

transform_mult = ACE1pack.generate_transform(multitransform_params)

@info("Test constructing pair basis and parameters")
# with polytransform transform
pair_basis = basis_params(type="pair", species = [:Ti, :Al], maxdeg = 6, r0 = 1.2)
ACE1pack.generate_basis(pair_basis)

# with multitransform
pair_basis = basis_params(type="pair", species = [:Ti, :Al], maxdeg = 6, r0 = 1.2, transform = multitransform_params)
ACE1pack.generate_basis(pair_basis)

@info("Test constructing radial basis parameters and ace radial basis")
rad_basis = basis_params(type="radial", rin = 0.0, pin = 0)
ACE1pack.generate_radial_basis(rad_basis, D, 6, :Si, transform_poly)
ACE1pack.generate_radial_basis(rad_basis, D, 6, :Si, transform_mult)


text:
test_gage(GeometricCenter,
    [1.111, 0, 0],
    [1.111, 0, 0],
    "test 1")


text:
using Logging, Random
using JuMP, Gurobi
using DecisionProgramming

rng = MersenneTwister(111)

C, D, V = random_diagram(rng, 5, 3, 2, 3, 3)
S = States(rng, [2, 3], length(C) + length(D))
X = [Probabilities(rng, c, S; n_inactive=0) for c in C]
Y = [Consequences(rng, v, S, low=-1.0, high=1.5) for v in V]

validate_influence_diagram(S, C, D, V)
sort!.((C, D, V, X, Y), by = x -> x.j)

P = DefaultPathProbability(C, X)
U = DefaultPathUtility(V, Y)

U⁺ = PositivePathUtility(S, U)
model = Model()
z = DecisionVariables(model, S, D)
x_s = PathCompatibilityVariables(model, z, S, P)

α = 0.1
w = 0.5
EV = expected_value(model, x_s, U⁺, P)
CVaR = conditional_value_at_risk(model, x_s, U⁺, P, α)
@objective(model, Max, w * EV + (1 - w) * CVaR)

optimizer = optimizer_with_attributes(
    () -> Gurobi.Optimizer(Gurobi.Env()),
    "IntFeasTol"      => 1e-9,
    "LazyConstraints" => 1,
)
set_optimizer(model, optimizer)
optimize!(model)

Z = DecisionStrategy(z)

@info("Printing decision strategy:")
print_decision_strategy(S, Z)

@info("Printing state probabilities:")
sprobs = StateProbabilities(S, P, Z)
print_state_probabilities(sprobs, [c.j for c in C])
print_state_probabilities(sprobs, [d.j for d in D])

@info("Computing utility distribution.")
udist = UtilityDistribution(S, P, U, Z)

@info("Printing utility distribution.")
print_utility_distribution(udist)

@info("Printing statistics")
print_statistics(udist)

@info("Printing risk measures")
print_risk_measures(udist, [α])


text:
using Smile
using Base.Test

function test_vec_approx_equal{T<:Real, S<:Real}(a::Vector{T}, b::Vector{S})
	@test length(a) == length(b)
	for i = 1 : length(a)
		@test isapprox(a[i], b[i])
	end
	true
end
function test_vec_approx_equal{T<:Real, S<:Real}(a::Vector{T}, b::Vector{S}, tol::Float64)
	@test length(a) == length(b)
	for i = 1 : length(a)
		@test abs(a[i] - b[i]) < tol
	end
	true
end

# write your own tests here
include("test_intarray.jl")
include("test_pattern.jl")
include("test_dataset.jl")
include("test_learning.jl")
include("test_sampling.jl")
include("tut1_create_a_bayesian_network.jl")
include("high_level_learning.jl")


text:
include("./src/plotter.jl")
include("./src/plot_hyperplane_dist.jl")
include("src/contact_structure.jl")

plotly()
plot_along_param_curve(t -> (5t-2.5,5t-2.5,0), standard_structure, 0.125;  midpoint_marker_color=:red, ε=0.7)
plot!(xlims=(-3, 3), ylims=(-3,3), zlims=(-3,3))
savefig("standard_contact_str_r3_1.html")

plot_along_param_curve(t -> (5t-2.5,-2.5,0), standard_structure, 0.125;  midpoint_marker_color=:red, ε=0.6)
plot_along_param_curve!(t -> (5t-2.5,-1.25,0), standard_structure, 0.125;  midpoint_marker_color=:red, ε=0.6)
plot_along_param_curve!(t -> (5t-2.5,0,0), standard_structure, 0.125;  midpoint_marker_color=:red, ε=0.6)
plot_along_param_curve!(t -> (5t-2.5,1.25,0), standard_structure, 0.125;  midpoint_marker_color=:red, ε=0.6)
plot_along_param_curve!(t -> (5t-2.5,2.5,0), standard_structure, 0.125;  midpoint_marker_color=:red, ε=0.6)
plot!(xlims=(-3, 3), ylims=(-3,3), zlims=(-3,3))
savefig("standard_contact_str_r3_2.html")


text:
hello()   = println("Hello, world!")
quine()   = println(src)
bottles() = for i = 99:-1:1 print("\n$i bottles of beer on the wall\n$i bottles of beer\nTake one down, pass it around\n$(i-1) bottles of beer on the wall\n") end
acc = 0
incr()    = global acc += 1

const dispatch = Dict(
'h' => hello,
'q' => quine,
'9' => bottles,
'+' => incr)

if length(ARGS) < 1
    println("Usage: julia ./HQ9+.jl file.hq9")
    exit(1)
else
    file = ARGS[1]
end

try
    open(file) do s
        global src = readstring(s)
    end
catch
    warning("can't open $file")
    exit(1)
end

for i in lowercase(src)
    if haskey(dispatch, i) dispatch[i]() end
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
    """
    `RacahAlgebra.integralRulesForOneYlm(rex::RacahAlgebra.RacahExpression)`
        ... attempts to find a simplification of the given Racah expression by using sum rules for one Wigner 3j symbol.
            Once a simplification is found, no attempt is made to find another simplifcation for this set of rules.
            A (istrue, rex)::Tuple{Bool, RacahExpression} is returned but where rex has no meaning for !istrue.
    """
    function integralRulesForOneYlm(rex::RacahAlgebra.RacahExpression)

        error("Not yet implemented.")
        # Loop through all Wigner 3j symbols
        for  (iaW3j, aW3j) in enumerate(rex.w3js)
            aRexList = RacahAlgebra.symmetricForms(aW3j)
            for  xaRex in aRexList
                ww = xaRex.w3js[1]
                #
                #   Rule:                  -m   ( j  j J )         -j
                #   ----        Sum(m) (-1)     (        )  =  (-1)    (2j + 1)  delta(J,0)  delta(M,0)
                #                               ( m -m M )
                #
                specialW3j = W3j(ww.ja, ww.ja, ww.jc, ww.ma, -ww.ma, ww.mc)
                if  ww == specialW3j
                    newPhase  = rex.phase  + xaRex.phase;     testPhase = newPhase   + ww.ma
                    newWeight = rex.weight * xaRex.weight;    testWeight = newWeight
                    newW3js   = W3j[];       for (ibW3j, bW3j) in enumerate(rex.w3js)  if  iaW3j != ibW3j   push!(newW3js, bW3j)   end   end
                    #
                    if  RacahAlgebra.hasIndexOld(ww.ma, rex.summations)          &&  RacahAlgebra.hasNoVars([ww.ma], testPhase)      &&
                        RacahAlgebra.hasNoVars([ww.ma], testWeight)           &&  RacahAlgebra.hasNoVars([ww.ma], rex.deltas)     &&
                        RacahAlgebra.hasNoVars([ww.ma], rex.triangles)        &&  RacahAlgebra.hasNoVars([ww.ma], newW3js)        &&
                        RacahAlgebra.hasNoVars([ww.ma], rex.w6js)             &&  RacahAlgebra.hasNoVars([ww.ma], rex.w9js)
                        newSummations = RacahAlgebra.removeIndex(ww.ma, rex.summations)
                        newDeltas = rex.deltas;      push!( newDeltas, Kronecker(ww.jc, 0));     push!( newDeltas, Kronecker(ww.mc, 0))
                        wa = RacahExpression( newSummations, rex.integrals, RacahAlgebra.purifyPhase(newPhase + ww.ma - ww.ja), newWeight * (2*ww.ja+1),
                                              newDeltas, rex.triangles, newW3js, rex.w6js, rex.w9js, rex.ylms, rex.djpqs )
                        println(">> Apply sum rule for one W3j -- Sum(m) (-1)^m ...")
                        return( (true, wa) )
                    end
                end
            end
        end

        return( (false, RacahExpression()) )
    end