text:
#
# This file is a part of MolecularGraph.jl
# Licensed under the MIT License http://opensource.org/licenses/MIT
#

export
    QueryFormula, QueryMol, querymol,
    tidyformula, findformula, findallformula,
    convertnotquery, convertnotquery!,
    inferatomaromaticity!, removehydrogens,
    recursiveatommatch,
    query_relationship, filter_queries

struct QueryFormula
    key::Symbol
    value::Any
end

Base.:(==)(a::QueryFormula, b::QueryFormula) = a.key == b.key && a.value == b.value

function Base.in(a::QueryFormula, b::Set{QueryFormula})
    # Necessary for set operations
    for i in b
        a == i && return true
    end
    return false
end

"""
    issubset(a::QueryFormula, b::QueryFormula) -> Bool

Check if fml1 contains fml2 (that is, all the query results of fml1 is included in the results of fml2)
"""
function Base.issubset(a::QueryFormula, b::QueryFormula)
    b == QueryFormula(:any, true) && return true
    if a.key === :and
        amap = Dict(i => v for (i, v) in enumerate(a.value))
        bset = b.key === :and ? b.value : [b]
        bmap = Dict(i => v for (i, v) in enumerate(bset))
        func = (x, y) -> issubset(amap[x], bmap[y])
        return maxcard(keys(amap), keys(bmap), func) == length(bmap)
    elseif b.key === :or
        bmap = Dict(i => v for (i, v) in enumerate(b.value))
        aset = a.key === :or ? a.value : [a]
        amap = Dict(i => v for (i, v) in enumerate(aset))
        func = (x, y) -> issubset(amap[x], bmap[y])
        return maxcard(keys(amap), keys(bmap), func) == length(amap)
    elseif a.key === :or && b.key === :not
        common = Set([c.key for c in a.value])
        length(common) == 1 || return false
        collect(common)[1] == b.value.key || return false
        return !(b.value.value in Set([c.value for c in a.value]))
    elseif a.key === :recursive
        if b.key === :recursive
            return hassubstructmatch(
                smartstomol(a.value), smartstomol(b.value),
                mandatory=Dict(1 => 1)
            )
        else
            amol = smartstomol(a.value)
            return issubset(nodeattr(amol, 1).query, b)
        end
    else
        return a == b
    end
end

struct QueryMol{A<:QueryAtom,B<:QueryBond} <: OrderedGraph
    neighbormap::Vector{Dict{Int,Int}}
    edges::Vector{Tuple{Int,Int}}
    nodeattrs::Vector{A}
    edgeattrs::Vector{B}
    cache::Dict{Symbol,Any}
    attributes::Dict{Symbol,Any}
    connectivity::Vector{Vector{Int}}
end

"""
    querymol() -> QueryMol

Generate empty `QueryMol`.
"""
querymol(
    ::Type{A}, ::Type{B}
) where {A<:QueryAtom,B<:QueryBond} = QueryMol{A,B}(
    [], [], [], [], Dict(), Dict(), [])

"""
    querymol(atoms::Vector{Atom}, bonds::Vector{Bond}) -> GraphMol

Generate `QueryMol` that has the given atom objects and edge objects.
"""
function querymol(edges, atoms::Vector{A}, bonds::Vector{B},
        connectivity::Vector{Vector{Int}}) where {A<:QueryAtom,B<:QueryBond}
    nbrmap = [Dict{Int,Int}() for i in 1:length(atoms)]
    edges = collect(edges)
    for (i, (u, v)) in enumerate(edges)
        nbrmap[u][i] = v
        nbrmap[v][i] = u
    end
    return QueryMol(
        nbrmap, edges, atoms, bonds,
        Dict{Symbol,Any}(), Dict{Symbol,Any}(), connectivity)
end

"""
    querymol(mol::SubgraphView{QueryMol}) -> QueryMol

Generate a new `QueryMol` from a substructure view.

Graph property caches and attributes are not inherited.
"""
function querymol(view::SubgraphView)
    newg = querymol(nodeattrtype(view), edgeattrtype(view))
    nkeys = sort(collect(nodeset(view)))
    ekeys = sort(collect(edgeset(view)))
    nmap = Dict{Int,Int}()
    for (i, n) in enumerate(nkeys)
        nmap[n] = i
        push!(newg.nodeattrs, nodeattr(view, n))
        push!(newg.neighbormap, Dict())
    end
    for (i, e) in enumerate(ekeys)
        (oldu, oldv) = getedge(view, e)
        u = nmap[oldu]
        v = nmap[oldv]
        push!(newg.edges, (u, v))
        push!(newg.edgeattrs, edgeattr(view, e))
        newg.neighbormap[u][i] = v
        newg.neighbormap[v][i] = u
    end
    return newg
end

"""
tidyformula(fml::QueryFormula) -> QueryFormula

Return tidy formulae.

- associative formulae will be juxtaposed
  (ex. :and => (A, :and => (B, C)) -> :and => (A, B, C))
- distributive formulae will be factored out
  (ex. :or => (:and => (A, B), :and => (A, C)) -> :and => (A, :or => (B, C)))
- Absorption
  (ex. :and => (A, :or => (A, B)) -> A
- `:any` absorbs everything
  (ex. :and => (:any => true, A) -> A, :or => (:any => true, A) -> :any => true
- `:not` would be inverted if possible
  (ex. :not => (:A => true) -> :A => false)
- `:not` has the highest precedence in SMARTS, but only in the case like [!C],
  De Morgan's law will be applied to remove `:and` under `:not`.
  (ex. :not => (:and => (:atomsymbol => :C, :isaromatic => false)
   -> :or => (:not => (:atomsymbol => :C), isaromatic => true)
"""
function tidyformula(fml::QueryFormula)
    # not
    if fml.key === :not
        child = fml.value
        if child.key === :and  # only the cases like [!C]
            return QueryFormula(:or, Set([
                tidyformula(QueryFormula(:not, c)) for c in child.value
            ]))
        elseif typeof(child.value) === Bool
            return QueryFormula(child.key, !child.value)
        end
    end
    fml.key in (:and, :or) || return fml
    childs = Set{QueryFormula}()
    # Association
    for child in fml.value
        cfml = tidyformula(child)
        if cfml.key === :any
            # Absorption
            (fml.key === :and) == cfml.value && continue
            return QueryFormula(:any, cfml.value)
        elseif cfml.key === fml.key
            union!(childs, cfml.value)
        else
            union!(childs, [cfml])
        end
    end
    length(childs) == 1 && return collect(childs)[1]
    # Distribution
    ckey = fml.key === :and ? :or : :and
    cc = collect(childs)
    bin = cc[1].key == ckey ? copy(cc[1].value) : Set([cc[1]])
    mono = Set([b.key for b in bin])
    for child in cc[2:end]
        @assert child.key !== fml.key  # already associated
        elems = child.key === ckey ? child.value : [child]
        intersect!(bin, elems)
        union!(mono, Set([e.key for e in elems]))
    end
    if isempty(bin)
        if fml.key === :and && length(mono) == 1
            # SMARTS primitives are disjoint, so the intersection should be an empty set.
            return QueryFormula(:any, false)
        else
            return QueryFormula(fml.key, childs)
        end
    end
    updated = Set{QueryFormula}()
    for child in childs
        cdiff = setdiff(child.key == ckey ? child.value : [child], bin)
        if isempty(cdiff)
            # Absorption
            return length(bin) == 1 ? collect(bin)[1] : QueryFormula(ckey, bin)
        elseif length(cdiff) == 1
            push!(updated, collect(cdiff)[1])
        else
            push!(updated, QueryFormula(ckey, cdiff))
        end
    end
    return tidyformula(QueryFormula(ckey, Set([bin..., QueryFormula(fml.key, updated)])))
end

"""
    findformula(q::QueryFormula, key::Symbol) -> Any

Find value of the given key.

If `deep=true` is set, it will also search under `:or` sets to match multiple formulae and
return the value after applying the `aggregate` function to them.
"""
function findformula(q::QueryFormula, key::Symbol; deep=false, aggregate=minimum)
    if q.key == key
        return aggregate([q.value])
    elseif q.key === :and
        for child in q.value
            if child.key == key
                return aggregate([child.value])
            elseif deep && child.key === :or
                gc = findformula(child, key, deep=deep, aggregate=aggregate)
                gc === nothing || return gc
            end
        end
    elseif deep && q.key === :or
        values = []
        for child in q.value
            child.key == key || return
            push!(values, child.value)
        end
        return aggregate(values)
    end
    return
end

"""
    findallformula(q::QueryFormula, key::Symbol) -> Any

Find all formula values that matched the given key regardless of logical operator roles.
"""
function findallformula(q::QueryFormula, key::Symbol)
    if q.key == key
        return [q.value]
    elseif q.key === :not
        return findallformula(q.value, key)
    elseif q.key in (:and, :or)
        values = []
        for child in q.value
            append!(values, findallformula(child, key))
        end
        return values
    end
    return []
end

NOTCONV = Dict(
    :isaromatic => Set([true, false]),
    :hydrogenconnected => Set([0, 1, 2, 3, 4]),
    :sssrcount => Set([0, 1, 2, 3]),
    :isringbond => Set([true, false]),
    :isaromaticbond => Set([true, false]),
    :stereo => Set([:up, :down, :unspecified])
)

"""
    convertnotquery(mol::QueryMol) -> QueryMol

Convert :not query to :or if possible.

- [!#1] would be *. Typical substructure mining tasks do not care stereochemistry,
so all hydrogen nodes can be removed.)
- following `tidyformula` is necessary.
"""
function convertnotquery(fml::QueryFormula)
    if fml.key === :not
        cfml = fml.value
        cfml == QueryFormula(:atomsymbol, :H) && return QueryFormula(:any, true)
        if haskey(NOTCONV, cfml.key)
            rest = setdiff(NOTCONV[cfml.key], [cfml.value])
            isempty(rest) && return QueryFormula(:any, false)
            length(rest) == 1 && return QueryFormula(cfml.key, collect(rest)[1])
            return QueryFormula(:or, Set([QueryFormula(cfml.key, v) for v in rest]))
        else
            return fml
        end
    end
    fml.key in (:and, :or) || return fml
    childs = Set{QueryFormula}(convertnotquery.(fml.value))
    return QueryFormula(fml.key, childs)
end

function convertnotquery!(qmol::QueryMol)
    nodes = []
    for n in 1:nodecount(qmol)
        na = nodeattr(qmol, n).query
        push!(nodes, SmartsAtom(tidyformula(convertnotquery(na))))
    end
    edges = []
    for e in 1:edgecount(qmol)
        ea = edgeattr(qmol, e).query
        push!(edges, SmartsBond(tidyformula(convertnotquery(ea))))
    end
    empty!(qmol.nodeattrs)
    append!(qmol.nodeattrs, nodes)
    empty!(qmol.edgeattrs)
    append!(qmol.edgeattrs, edges)
end

"""
    inferatomaromaticity!(qmol::QueryMol)

Infer aromaticity of atoms from connected bonds.

This directly add :isaromatic formulas to query atoms in the query molecule.
"""
function inferatomaromaticity!(qmol::QueryMol)
    for n in 1:nodecount(qmol)
        arom = false
        hcnt = 0
        heavyorder = 0
        for (inc, adj) in neighbors(qmol, n)
            iq = edgeattr(qmol, inc).query
            if findformula(iq, :isaromaticbond) === true
                arom = true
                break
            end
            aq = nodeattr(qmol, adj).query
            if findformula(aq, :atomsymbol) === :H
                hcnt += 1
            else
                heavyorder += something(findformula(iq, :bondorder, deep=true), 1)
            end
        end
        nq = nodeattr(qmol,n).query
        if arom
            newq = QueryFormula(:and, Set([nq, QueryFormula(:isaromatic, true)]))
            setnodeattr!(qmol, n, SmartsAtom(tidyformula(newq)))
            continue
        end
        hqcnt = findformula(nq, :hydrogenconnected, deep=true)
        if hqcnt !== nothing
            hcnt = hqcnt  # ignore hydrogen nodes if H query is already set
        end
        acceptable = Dict(:O => 0, :N => 1, :C => 1, :S => 0, :P => 1, :B => 1)
        if heavyorder + hcnt > get(acceptable, findformula(nq, :atomsymbol), 0)
            newq = QueryFormula(:and, Set([nq, QueryFormula(:isaromatic, false)]))
            setnodeattr!(qmol, n, SmartsAtom(tidyformula(newq)))
        end
    end
end

"""
    removehydrogens(mol::QueryMol) -> QueryMol

Return the molecular query with hydrogen nodes removed.

- Hydrogen nodes would be represented as `[*Hx]` in the connected heavy atoms.
(e.g. CN[H] -> C[N;H1,H2,H3,H4])
"""
function removehydrogens(mol::QueryMol)
    MAX_H_COUNT = 4
    hs = Set{Int}()
    hcntmap = Dict()
    for (i, atom) in enumerate(nodeattrs(mol))
        h = atom.query == QueryFormula(:atomsymbol, :H)
        h2 = atom.query == QueryFormula(
            :and, Set([QueryFormula(:atomsymbol, :H), QueryFormula(:isaromatic, false)]))
        (h || h2) || continue
        @assert degree(mol, i) == 1
        adj = iterate(adjacencies(mol, i))[1]
        if !haskey(hcntmap, adj)
            hcntmap[adj] = 0
        end
        hcntmap[adj] += 1
        push!(hs, i)
    end
    newmol = deepcopy(mol)
    for (i, hcnt) in hcntmap
        q = nodeattr(mol, i).query
        findformula(q, :hydrogenconnected, deep=true) === nothing || continue  # ignore if H is already set
        heavycnt = degree(mol, i) - hcnt
        maxh = MAX_H_COUNT - heavycnt
        if maxh <= hcnt
            hf = QueryFormula(:hydrogenconnected, hcnt)
        else
            hf = QueryFormula(
                :or, Set([QueryFormula(:hydrogenconnected, n) for n in hcnt:maxh]))
        end
        newq = QueryFormula(:and, Set([q, hf]))
        setnodeattr!(newmol, i, SmartsAtom(tidyformula(newq)))
    end
    ns = setdiff(nodeset(mol), hs)
    return querymol(nodesubgraph(newmol, ns))
end

function recursiveatommatch(qmol1::QueryMol, qmol2::QueryMol)
    return function (a1, a2)
        a1q = nodeattr(qmol1, a1).query
        a2q = nodeattr(qmol2, a2).query
        issubset(a1q, a2q) && return true
        recfmls= findformula(a2q, :recursive, deep=true, aggregate=collect)
        recfmls === nothing && return false
        for r in recfmls
            hassubstructmatch(
                qmol1, smartstomol(r), mandatory=Dict(a1 => 1)) && return true
        end
        return false
    end
end

const DEFAULT_QUERY_RELATIONS = let
    qrfile = joinpath(dirname(@__FILE__), "../../assets/const/default_query_relations.yaml")
    include_dependency(qrfile)
    qrfile
end

struct DictDiGraph <: OrderedDiGraph
    # TODO: should be moved to Graph module
    outneighbormap::Vector{Dict{Int,Int}}
    inneighbormap::Vector{Dict{Int,Int}}
    edges::Vector{Tuple{Int,Int}}
    nodeattrs::Vector{Dict}
    edgeattrs::Vector{Dict}
end

"""
    dictdigraph(view::DiSubgraphView{DictDiGraph}) -> DictDiGraph

Generate a new `DictDiGraph` from a substructure view.

Graph property caches and attributes are not inherited.
"""
function dictdigraph(view::DiSubgraphView{DictDiGraph})
    newg = DictDiGraph([], [], [], [], [])
    nkeys = sort(collect(nodeset(view)))
    ekeys = sort(collect(edgeset(view)))
    nmap = Dict{Int,Int}()
    for (i, n) in enumerate(nkeys)
        nmap[n] = i
        push!(newg.nodeattrs, nodeattr(view, n))
        push!(newg.outneighbormap, Dict())
        push!(newg.inneighbormap, Dict())
    end
    for (i, e) in enumerate(ekeys)
        (oldu, oldv) = getedge(view, e)
        u = nmap[oldu]
        v = nmap[oldv]
        push!(newg.edges, (u, v))
        push!(newg.edgeattrs, edgeattr(view, e))
        newg.outneighbormap[u][i] = v
        newg.inneighbormap[v][i] = u
    end
    return newg
end

"""
    query_relationship(;sourcefile=DEFAULT_QUERY_RELATIONS) -> DictDiGraph

Generate query relationship diagram.
"""
function query_relationship(;sourcefile=DEFAULT_QUERY_RELATIONS)
    graph = DictDiGraph([], [], [], [], [])
    keys = Dict()
    for (i, rcd) in enumerate(YAML.load(open(sourcefile)))
        addnode!(graph, rcd)
        keys[rcd["key"]] = i
    end
    for rcd in nodeattrs(graph)
        if haskey(rcd, "isa")
            for e in rcd["isa"]
                addedge!(graph, keys[rcd["key"]], keys[e], Dict("relation" => "isa"))
            end
        end
        if haskey(rcd, "has")
            for e in rcd["has"]
                addedge!(graph, keys[rcd["key"]], keys[e], Dict("relation" => "has"))
            end
        end
    end
    return graph
end

"""
    filter_queries(qr::DictDiGraph, mol::GraphMol) -> DictDiGraph

Filter query relationship diagram by the given molecule.
The filtered diagram represents query relationship that the molecule have.
"""
function filter_queries(qr::DictDiGraph, mol::GraphMol)
    matched = Set{Int}()
    qrc = deepcopy(qr)
    for n in reversetopologicalsort(qrc)
        rcd = nodeattr(qrc, n)
        if !issubset(successors(qrc, n), matched)
            continue
        end
        matches = collect(substructmatches(mol, smartstomol(rcd["query"])))
        if !isempty(matches)
            push!(matched, n)
            rcd["matched"] = Set([sort(collect(keys(m))) for m in matches])
        end
    end
    return dictdigraph(nodesubgraph(qrc, matched))
end


text:
module TestTikZWiringDiagrams
using Test

using Catlab.Theories, Catlab.WiringDiagrams
using Catlab: Graphics
using Catlab.Graphics: TikZ

to_tikz(f; kw...) = Graphics.to_tikz(f; outer_ports_layout=:fixed, kw...)

function stmts(expr::TikZ.Expression, type::Type)
  children = hasfield(typeof(expr), :stmts) ? expr.stmts : []
  vcat(expr isa type ? [expr] : [],
       mapreduce(expr -> stmts(expr, type), vcat, children; init=[]))
end
stmts(doc::TikZ.Document, type::Type) = stmts(doc.picture, type)

A, B = Ob(FreeSymmetricMonoidalCategory, :A, :B)
f, g = Hom(:f, A, B), Hom(:g, B, A)

pic = to_tikz(f)
@test length(stmts(pic, TikZ.Node)) == 2 # Includes outer box
@test length(stmts(pic, TikZ.Edge)) == 2

pic = to_tikz(id(A))
@test length(stmts(pic, TikZ.Node)) == 1
@test length(stmts(pic, TikZ.Edge)) == 1

pic = to_tikz(compose(f,g))
@test length(stmts(pic, TikZ.Node)) == 3
@test length(stmts(pic, TikZ.Edge)) == 3

pic = to_tikz(otimes(f,g))
@test length(stmts(pic, TikZ.Node)) == 3
@test length(stmts(pic, TikZ.Edge)) == 4

pic = to_tikz(braid(A,B))
@test length(stmts(pic, TikZ.Node)) == 1
@test length(stmts(pic, TikZ.Edge)) == 2

A, B = Ob(FreeBiproductCategory, :A, :B)

for expr in (mcopy(A), delete(A), mmerge(A), create(A))
  pic = to_tikz(expr)
  @test length(stmts(pic, TikZ.Node)) == 2
end

end

ext: jl


text:
using SigmaRidgeRegression
using Test
using Random
import StatsBase

ar1_design = BlockCovarianceDesign([
    AR1Design(;ρ=0.8),
    AR1Design(;ρ=0.5)])

id_design = BlockCovarianceDesign([
        IdentityCovarianceDesign(),
        IdentityCovarianceDesign()])

grp1 = GroupedFeatures([2000,4000])
grp2 = GroupedFeatures([800,500])

for grp in [grp1; grp2]
    for _design in [ar1_design; id_design]
        _design = set_groups(_design, grp)
        _n = 2000
        _γs = grp.ps ./ _n
        _λs = [2.0; 0.4]

        _αs = [1.0; 7.0]

        theory_risk = @show SigmaRidgeRegression.risk_formula(spectrum.(_design.blocks), _γs, _αs, _λs)

        _ridge_sim = GroupRidgeSimulationSettings(;
            groups=grp,
            ntrain=_n,
            ntest=20_000,
            Σ=_design,
            response_model=RandomLinearResponseModel(; αs=_αs, grp=grp),
        )
        Random.seed!(1)

        _sim_res = simulate(_ridge_sim)
        _X_train = _sim_res.X[_sim_res.resampling_idx[1][1],:]
        _Y_train = _sim_res.Y[_sim_res.resampling_idx[1][1]]

        _X_test = _sim_res.X[_sim_res.resampling_idx[1][2],:]
        _Y_test = _sim_res.Y[_sim_res.resampling_idx[1][2]]

        ridge_risk = mse_ridge(
                StatsBase.fit(
                    MultiGroupRidgeRegressor(; groups=grp, λs=_λs, center=false, scale=false),
                    _X_train,
                    _Y_train,
                    grp,
                ),
                _X_test,
                _Y_test,
            )
        @test ridge_risk ≈ theory_risk atol =1.0
    end
end


text:
using POMDPGym, POMDPs, Test

env = ZWrap(GymPOMDP(:CartPole), [1.0, 0.5])

s = rand(initialstate(env))
os = [rand(initialobs(env, s))]
i=1
while !isterminal(env, s) && i<100
	s, o, r = gen(env, s, rand([1,2]))
	push!(os, o)
	i += 1
end

@test all([o[1:2] == env.z for o in os])

env = ZWrap(PendulumMDP(), [1.5, 0.5])

s = rand(initialstate(env))
esses = [s]
i=1
while !isterminal(env, s) && i<100
	s, r = gen(env, s, rand(Float32))
	push!(esses, s)
	i += 1
end

@test all([s[1:2] == env.z for s in esses])


text:
import ...AMDGPU: libhsaruntime

const HSA_EXPORT_DECORATOR = identity
const HSA_CALL = nothing


text:
function [Grid] = CFPToGrid(Points,Faces,Cells,OrientFace)
Grid.Type='Hex';
NumNodes=size(Points,2);
Grid.Nodes(1:NumNodes)=Node([0 0 0],0);
NodeNumber=1;
for iN=1:NumNodes
  Grid.Nodes(NodeNumber)=Node([Points(iN).x;Points(iN).y;0],NodeNumber);
  NodeNumber=NodeNumber+1;
end
Grid.NumNodes=NumNodes;

NumEdges=size(Faces,2);
Grid.Edges(1:NumEdges)=Edge([1 2],Grid,0,0);

EdgeNumberI=0;
EdgeNumberB=0;
for iE=1:NumEdges
  if Faces(iE).Boundary==1
    EdgeNumberB=EdgeNumberB+1;
  else
    EdgeNumberI=EdgeNumberI+1;
  end
end
Grid.NumEdgesI=EdgeNumberI;
Grid.NumEdgesB=EdgeNumberB;

EdgeNumber=1;
EdgeNumberI=1;
EdgeNumberB=1;
EdgeOldToNew=zeros(1,NumEdges);
for iE=1:NumEdges
  if Faces(iE).Boundary==1
    Grid.Edges(EdgeNumber)=Edge(Faces(iE).P',Grid,EdgeNumber,...
      EdgeNumberB+Grid.NumEdgesI);
    EdgeOldToNew(iE)=EdgeNumberB+Grid.NumEdgesI;
    EdgeNumberB=EdgeNumberB+1;
  else
    Grid.Edges(EdgeNumber)=Edge(Faces(iE).P',Grid,EdgeNumber,...
      EdgeNumberI);
    EdgeOldToNew(iE)=EdgeNumberI;
    EdgeNumberI=EdgeNumberI+1;
  end
  EdgeNumber=EdgeNumber+1;
end
Grid.NumEdges=NumEdges;

NumFaces=size(Cells,2);
Faces1(1:NumFaces)=Face([0 0 0 0],Grid,0,'x');
FaceNumber=1;
Type='o';
for iF=1:NumFaces
  [Faces1(FaceNumber),Grid]=Face(Cells(iF).F,Grid,...
    FaceNumber,Type,OrientFace);
  FaceNumber=FaceNumber+1;
end
Grid.NumFaces=NumFaces;
Grid.Faces=Faces1;
end


text:
using Test
using AutoViz
using NGSIM
using NBInclude

@nbinclude(joinpath(dirname(pathof(NGSIM)), "..", "jnotebooks", "Demo.ipynb"))

# the following test requires downloading the dataset
# td1 = load_trajdata(1)
# scene = get!(Scene(500), td1, 1000)
# render(scene, ROADWAY_101)


text:
import Limbdark: Mn_series_coeff!

# Randomizer seed
Random.seed!(42)

#@testset "transit_poly" begin

# Now, integrate by hand:
function transit_poly_int(r,b,u,ns)
n = length(u)
if b < (1.0+r)
  s_1 = maximum([0.0,b-r]); s_2=minimum([1.0,b+r])
  ds=(s_2-s_1)/ns
  s = linearspace(s_1+.5*ds,s_2-.5*ds,ns)
  fobs = zero(r)
  for j=1:ns
    if s[j] < r-b
      dphi = 2pi
    else
      dphi = 2*acos((s[j]^2+b^2-r^2)/(2*b*s[j]))
    end
    z = sqrt(1-s[j]^2)
    imu = 1.0
    for i=1:n
      imu -= u[i]*(1-z)^i
    end
    fobs += s[j]*ds*dphi*imu
  end
  norm = 1.
  for i=1:n
    norm -= 2*u[i]/(2+3*i+i^2)
  end
  fobs /= pi*norm
  fobs = 1-fobs
  return fobs
else
  return 1.0
end
end

r0=[0.01,0.01,0.01,0.01,0.01,0.1,0.1,0.1,0.1,0.1,1.0,1.0,1.0,10.,10.,10.,100.,100.,100.]
b0=[0.,0.01,0.99,1.0,1.01, 0.,0.1,0.9,1.0,1.1, 0.,1.0,2.0, 9.,10.,11., 99.,100.,101.]

n = 2+ceil(Int64,rand()*20); u = rand(n); u *= rand()/sum(u)

ns = 5000
trans = transit_init(r0[1],b0[1],u,false)
for i=1:length(r0)
  trans.r = r0[i]; trans.b = b0[i]
  flux = transit_poly!(trans)
  @time flux = transit_poly!(trans)
  f_num = transit_poly_int(r0[i],b0[i],u,ns)
  @time f_num = transit_poly_int(r0[i],b0[i],u,ns)
  println("r: ",r0[i]," b: ",b0[i]," f_an: ",flux," f_num: ",f_num," diff: ",flux-f_num)
  @test isapprox(f_num,flux,atol=1e-5)
end

# Now, compute light curves and plot:

npts = 1000
iu = 30
u_n = ones(iu) / iu
b0 = zeros(npts)
for i=1:npts
    b0[i] = 3.0 * (i / float(npts) - 0.5)
end

lc_ana = zeros(npts); lc_num=zeros(npts)
rtest = 0.1
trans = transit_init(rtest,b0[1],u_n,false)

for i=1:length(b0)
#  if mod(i,10) == 0
#    println("i: ",i," b: ",b)
#  end
#  lc_ana[i] = transit_poly(r,abs(b),u_n)
  trans.b = abs(b0[i])
  lc_ana[i] = transit_poly!(trans)
  lc_num[i] = transit_poly_int(rtest,abs(b0[i]),u_n,ns)
end

using PyPlot

plot(b0,lc_ana, linewidth=2, label="Analytic")
plot(b0,lc_num, linewidth=1, label="Numeric")

println("Maximum difference of lightcurve for N=",iu,": ",maximum(abs,lc_ana-lc_num))
@test isapprox(lc_ana,lc_num,atol=1e-5)

#end


text:
using AVLTrees
using Test

@testset "AVLTrees.jl" begin
    include("node.jl")
    include("tree.jl")
    include("set.jl")
end


text:
using TreeKnit
using TreeKnit.SplitGraph
using Test
using TreeTools

println("##### splitgraph #####")

include("$(dirname(pathof(TreeKnit)))/..//test/splitgraph/basic/test.jl")
include("$(dirname(pathof(TreeKnit)))/..//test/splitgraph/3solutions/test.jl")
include("$(dirname(pathof(TreeKnit)))/..//test/splitgraph/w_resolution/test.jl")

max_stars_repo_name: PierreBarrat/RecombTools

max_issues_repo_name: PierreBarrat/TreeKnit

max_forks_repo_name: PierreBarrat/TreeKnit


text:
using GFlops
using Test

function my_axpy!(a, x, y)
    @inbounds @simd for i in eachindex(x)
        y[i] += a*x[i]
    end
end

function my_prod(m, v)
    res = similar(v, size(m, 1))

    for i in 1:size(m, 1)
        acc = zero(eltype(m))
        for j in 1:length(v)
            acc += m[i,j]*v[j]
        end
        res[i] = acc
    end
    res
end

import BenchmarkTools: @benchmark
struct FakeResults
    times
    allocs
    memory
end
macro benchmark(e)
    quote
        FakeResults(#= times  =# [2.0, 3.0],
                    #= allocs =# 1,
                    #= memory =# 1042)
    end
end

@testset "GFlops" begin
    @testset "Counter" begin
        @testset "display empty" begin
            let
                cnt = GFlops.Counter()
                str = string(cnt)
                @test str == "Flop Counter: 0 flop"
            end
        end

        @testset "display non-empty" begin
            let
                cnt = GFlops.Counter()
                cnt.add32 = 1
                str = string(cnt)
                @test  occursin("Float32", str)
                @test !occursin("Float64", str)
                @test  occursin("add",     str)
                @test !occursin("mul",     str)
                @test  occursin(" 1 ",     str)
            end
        end
    end

    @testset "@count_ops" begin
        @testset "mul+add 64" begin
            let N = 100
                a = 2.5
                x = rand(N)
                y = similar(x)

                cnt = @count_ops my_axpy!(a, x, y)
                @test cnt.add64 == N
                @test cnt.mul64 == N
                @test GFlops.flop(cnt) == 2*N

                m = rand(N, N)
                v = rand(N)

                cnt = @count_ops(my_prod(m, v))
                @test cnt.add64 == N*N
                @test cnt.mul64 == N*N
                @test GFlops.flop(cnt) == 2*N*N
            end
        end

        @testset "mul+add 32" begin
            let N = 100
                a = 2.5f0
                x = rand(Float32, N)
                y = similar(x)

                cnt = @count_ops my_axpy!(a, x, y)
                @test cnt.add32 == N
                @test cnt.mul32 == N
                @test GFlops.flop(cnt) == 2*N

                m = rand(Float32, N, N)
                v = rand(Float32, N)

                cnt = @count_ops(my_prod(m, v))
                @test cnt.add32 == N*N
                @test cnt.mul32 == N*N
                @test GFlops.flop(cnt) == 2*N*N
            end
        end

        @testset "neg" begin
            let cnt = @count_ops -(4.2)
                @test cnt.neg64 == 1
                @test GFlops.flop(cnt) == 1
            end

            let cnt = @count_ops -(4.2f0)
                @test cnt.neg32 == 1
                @test GFlops.flop(cnt) == 1
            end
        end

        @testset "abs" begin
            let cnt = @count_ops abs(-4.2)
                @test cnt.abs64 == 1
                @test GFlops.flop(cnt) == 1
            end

            let cnt = @count_ops abs(-4.2f0)
                @test cnt.abs32 == 1
                @test GFlops.flop(cnt) == 1
            end
        end

        @testset "sqrt" begin
            let cnt = @count_ops sqrt(4.2)
                @test cnt.sqrt64 == 1
                @test GFlops.flop(cnt) == 1
            end

            let cnt = @count_ops sqrt(4.2f0)
                @test cnt.sqrt32 == 1
                @test GFlops.flop(cnt) == 1
            end
        end

        @testset "rem" begin
            let cnt = @count_ops rem(12.0, 5.0)
                @test cnt.rem64 == 1
                @test GFlops.flop(cnt) == 1
            end

            let cnt = @count_ops rem(12.0f0, 5.0f0)
                @test cnt.rem32 == 1
                @test GFlops.flop(cnt) == 1
            end
        end

        @testset "fma" begin
            let cnt = @count_ops fma(1.0, 2.0, 3.0)
                @test cnt.fma64 == 1
                @test GFlops.flop(cnt) == 2
            end

            let cnt = @count_ops fma(1.0f0, 2.0f0, 3.0f0)
                @test cnt.fma32 == 1
                @test GFlops.flop(cnt) == 2
            end
        end

        @testset "muladd" begin
            let cnt = @count_ops muladd(1.0, 2.0, 3.0)
                @test cnt.muladd64 == 1
                @test GFlops.flop(cnt) == 2
            end

            let cnt = @count_ops muladd(1.0f0, 2.0f0, 3.0f0)
                @test cnt.muladd32 == 1
                @test GFlops.flop(cnt) == 2
            end
        end

        @testset "interpolated arguments" begin
            let N = 100

                T = Float64
                cnt = @count_ops my_axpy!(pi, $(rand(T, N)), $(rand(T, N)))
                @test cnt.add64 == N
                @test cnt.mul64 == N
                @test GFlops.flop(cnt) == 2*N

                T = Float32
                cnt = @count_ops my_axpy!(pi, $(rand(T, N)), $(rand(T, N)))
                @test cnt.add32 == N
                @test cnt.mul32 == N
                @test GFlops.flop(cnt) == 2*N
            end
        end

        @testset "broadcast" begin
            let N = 100

                x = 42.0
                cnt1 = @count_ops sin(x)
                cnt2 = @count_ops sin.($(fill(x, N)))
                @test GFlops.flop(cnt1) != 0
                @test cnt2 == N*cnt1

                x = 42.0f0
                cnt1 = @count_ops sin(x)
                cnt2 = @count_ops sin.($(fill(x, N)))
                @test GFlops.flop(cnt1) != 0
                @test cnt2 == N*cnt1
            end
        end
    end

    @testset "@gflops" begin
        let
            N = 100
            a = 2.5
            x = rand(N)
            y = Vector{Float64}(undef, N)

            @test @gflops(my_axpy!(a, x, y))          == N
            @test @gflops(my_axpy!(π, $(rand(N)), y)) == N
        end

        let
            N = 100
            m = rand(N, N)
            v = rand(N)
            @test @gflops(my_prod(m, v)) == N*N
        end
    end
end


text:
using Test
using ParallelStencil
import ParallelStencil: @reset_parallel_stencil, @is_initialized, SUPPORTED_PACKAGES, PKG_CUDA, PKG_THREADS
import ParallelStencil: @require
using ParallelStencil.FiniteDifferences2D
@static if (PKG_CUDA in SUPPORTED_PACKAGES) import CUDA end

@static for package in SUPPORTED_PACKAGES  eval(:(
    @testset "$(basename(@__FILE__)) (package: $(nameof($package)))" begin
        @init_parallel_stencil($package, Float64, 2)
        @require @is_initialized()
        nx, ny = 7, 5
        A      =  @rand(nx,   ny  );
        Ax     =  @rand(nx+1, ny  );
        Ay     =  @rand(nx,   ny+1);
        Axy    =  @rand(nx+1, ny+1);
        Axx    =  @rand(nx+2, ny  );
        Ayy    =  @rand(nx,   ny+2);
        Axyy   =  @rand(nx+1, ny+2);
        Axxy   =  @rand(nx+2, ny+1);
        Axxyy  =  @rand(nx+2, ny+2);
        R      = @zeros(nx,   ny  );
        Rxxyy  = @zeros(nx+2, ny+2);
        @testset "1. compute macros" begin
            @testset "differences" begin
                @require @is_initialized()
                @parallel  d_xa!(R, Ax)    = (@all(R) = @d_xa(Ax); return)
                @parallel  d_ya!(R, Ay)    = (@all(R) = @d_ya(Ay); return)
                @parallel  d_xi!(R, Axyy)  = (@all(R) = @d_xi(Axyy); return)
                @parallel  d_yi!(R, Axxy)  = (@all(R) = @d_yi(Axxy); return)
                @parallel d2_xi!(R, Axxyy) = (@all(R) = @d2_xi(Axxyy); return)
                @parallel d2_yi!(R, Axxyy) = (@all(R) = @d2_yi(Axxyy); return)
                R.=0; @parallel  d_xa!(R, Ax);     @test all(R .== Ax[2:end,    :].-Ax[1:end-1,      :])
                R.=0; @parallel  d_ya!(R, Ay);     @test all(R .== Ay[    :,2:end].-Ay[      :,1:end-1])
                R.=0; @parallel  d_xi!(R, Axyy);   @test all(R .== Axyy[2:end  ,2:end-1].-Axyy[1:end-1,2:end-1])
                R.=0; @parallel  d_yi!(R, Axxy);   @test all(R .== Axxy[2:end-1,2:end  ].-Axxy[2:end-1,1:end-1])
                R.=0; @parallel d2_xi!(R, Axxyy);  @test all(R .== (Axxyy[3:end,2:end-1].-Axxyy[2:end-1,2:end-1]).-(Axxyy[2:end-1,2:end-1].-Axxyy[1:end-2,2:end-1]))
                R.=0; @parallel d2_yi!(R, Axxyy);  @test all(R .== (Axxyy[2:end-1,3:end].-Axxyy[2:end-1,2:end-1]).-(Axxyy[2:end-1,2:end-1].-Axxyy[2:end-1,1:end-2]))
            end;
            @testset "selection" begin
                @require @is_initialized()
                @parallel all!(R, A)     = (@all(R) = @all(A); return)
                @parallel inn!(R, Axxyy) = (@all(R) = @inn(Axxyy); return)
                @parallel inn_x!(R, Axx) = (@all(R) = @inn_x(Axx); return)
                @parallel inn_y!(R, Ayy) = (@all(R) = @inn_y(Ayy); return)
                R.=0; @parallel all!(R, A);      @test all(R .== A)
                R.=0; @parallel inn!(R, Axxyy);  @test all(R .== Axxyy[2:end-1,2:end-1])
                R.=0; @parallel inn_x!(R, Axx);  @test all(R .== Axx[2:end-1,      :])
                R.=0; @parallel inn_y!(R, Ayy);  @test all(R .== Ayy[      :,2:end-1])
            end;
            @testset "averages" begin
                @require @is_initialized()
                @parallel av!(R, Axy)     = (@all(R) = @av(Axy); return)
                @parallel av_xa!(R, Ax)   = (@all(R) = @av_xa(Ax); return)
                @parallel av_ya!(R, Ay)   = (@all(R) = @av_ya(Ay); return)
                @parallel av_xi!(R, Axyy) = (@all(R) = @av_xi(Axyy); return)
                @parallel av_yi!(R, Axxy) = (@all(R) = @av_yi(Axxy); return)
                R.=0; @parallel av!(R, Axy);      @test all(R .== (Axy[1:end-1,1:end-1].+Axy[2:end,1:end-1].+Axy[1:end-1,2:end].+Axy[2:end,2:end])*0.25)
                R.=0; @parallel av_xa!(R, Ax);    @test all(R .== (Ax[2:end,    :].+Ax[1:end-1,      :]).*0.5)
                R.=0; @parallel av_ya!(R, Ay);    @test all(R .== (Ay[    :,2:end].+Ay[      :,1:end-1]).*0.5)
                R.=0; @parallel av_xi!(R, Axyy);  @test all(R .== (Axyy[2:end  ,2:end-1].+Axyy[1:end-1,2:end-1]).*0.5)
                R.=0; @parallel av_yi!(R, Axxy);  @test all(R .== (Axxy[2:end-1,2:end  ].+Axxy[2:end-1,1:end-1]).*0.5)
            end;
            @testset "others" begin
                @require @is_initialized()
                @parallel maxloc!(R, Axxyy) = (@all(R) = @maxloc(Axxyy); return)
                R.=0; @parallel maxloc!(R, Axxyy); @test all(R .== max.(max.(max.(max.(Axxyy[1:end-2,2:end-1],Axxyy[3:end,2:end-1]),Axxyy[2:end-1,2:end-1]),Axxyy[2:end-1,1:end-2]),Axxyy[2:end-1,3:end]))
            end;
        end;
        @testset "2. apply masks" begin
            @testset "selection" begin
                @require @is_initialized()
                @parallel inn_all!(Rxxyy, A)     = (@inn(Rxxyy) = @all(A); return)
                @parallel inn_inn!(Rxxyy, Axxyy) = (@inn(Rxxyy) = @inn(Axxyy); return)
                Rxxyy.=0; @parallel inn_all!(Rxxyy, A);      @test all(Rxxyy[2:end-1,2:end-1] .== A)
                Rxxyy[2:end-1,2:end-1].=0;                   @test all(Rxxyy .== 0)  # Test that boundary values remained zero.
                Rxxyy.=0; @parallel inn_inn!(Rxxyy, Axxyy);  @test all(Rxxyy[2:end-1,2:end-1] .== Axxyy[2:end-1,2:end-1])
                Rxxyy[2:end-1,2:end-1].=0;                   @test all(Rxxyy .== 0)  # Test that boundary values remained zero.
            end;
            @testset "differences" begin
                @require @is_initialized()
                @parallel  inn_d_xa!(Rxxyy, Ax)     = (@inn(Rxxyy) = @d_xa(Ax); return)
                @parallel  inn_d_yi!(Rxxyy, Axxy)   = (@inn(Rxxyy) = @d_yi(Axxy); return)
                @parallel  inn_d2_yi!(Rxxyy, Axxyy) = (@inn(Rxxyy) = @d2_yi(Axxyy); return)
                Rxxyy.=0; @parallel  inn_d_xa!(Rxxyy, Ax);     @test all(Rxxyy[2:end-1,2:end-1] .== Ax[2:end,    :].-Ax[1:end-1,      :])
                Rxxyy[2:end-1,2:end-1].=0;                     @test all(Rxxyy .== 0)  # Test that boundary values remained zero.
                Rxxyy.=0; @parallel  inn_d_yi!(Rxxyy, Axxy);   @test all(Rxxyy[2:end-1,2:end-1] .== Axxy[2:end-1,2:end  ].-Axxy[2:end-1,1:end-1])
                Rxxyy[2:end-1,2:end-1].=0;                     @test all(Rxxyy .== 0)  # Test that boundary values remained zero.
                Rxxyy.=0; @parallel inn_d2_yi!(Rxxyy, Axxyy);  @test all(Rxxyy[2:end-1,2:end-1] .== (Axxyy[2:end-1,3:end].-Axxyy[2:end-1,2:end-1]).-(Axxyy[2:end-1,2:end-1].-Axxyy[2:end-1,1:end-2]))
                Rxxyy[2:end-1,2:end-1].=0;                     @test all(Rxxyy .== 0)  # Test that boundary values remained zero.
            end;
        end;
        @reset_parallel_stencil()
    end;
)) end == nothing || true;


text:
"""
---
Quadrature Amplitude Modulation (QAM) hard decoding function
Apply symbol hard demapping to a input symbol sequence (of size 1xN) with constellation size M.
Output is a binary (1xL) with N = L / log2(M)
Conventional gray demapping is used.
Input constellation is Array{Complex{Float64}}
Output bitsream is Array{Int8}
# --- Syntax
	  bitDemappingQAM!(hardBits,M,qamVect)
# --- Input parameters
- hardBits	: Vector of bits to populate [Array{UInt8}, length(qamVect)/log2(M)]
- M			: Constellation size (i.e from 4 to 256)
- qamVect	: Complex observation vector to decode.
# --- Output parameters
- []
"""
function bitDemappingQAM!(hardBits,M, qamVect)
	# ----------------------------------------------------
	# --- Overall parameters
	# ----------------------------------------------------
	# --- Bit sequence and initialisation
	nbBits	= Int( length(qamVect) * log2(M));
	# ----------------------------------------------------
	# --- Switch on modulation order
	# ----------------------------------------------------
	# Here all Voronoi region are defined with normalized constellation
	# Symbol demapper are recall from bitMappingQAM (without rescaling factor)
	if M == 2
		## BPSK demodulator
		hardBits .= Int.( real(qamVect) .> 0);
	elseif M == 4
		## --- QPSK demodulator
		# --- Gray encoding
		#  -1  1
		#   0  1
		# --- Voronoi region based on sign
		# < 0 --> -1
		# > 0 --> 1
		# --- Iterative hard demapping
		for n = 1 : 1 : length(qamVect)
			hardBits[(n-1)*2+1] = Int8(real(qamVect[n]) > 0);
			hardBits[(n-1)*2+2] = Int8(imag(qamVect[n]) > 0);
		end
	elseif M == 16
		## 16-QAM hard demapper
		# --- Gray encoding
		#  -3     -1     1     3
		#  10     11    01    00
		# --- Voronoi region definition
		bounds	= [-2 0 2]/(sqrt(2/3*15));
		#bounds	= [-sqrt(2)/2 0 sqrt(2)/2];
		# --- Iterative hard demapping
		for n = 1 : 1 : length(qamVect)
			# --- Decision for real part
			# Getting only real part
			realI	= real(qamVect[n]);
			# Getting area of interest: from area 1 to 4
			decI	= 1+ (realI > bounds[1]) + (realI>bounds[2]) + (realI>bounds[3]);
			# Convert into binary sequence based on Gray encoding scheme.
			if decI == 1
				hardBits[4(n-1)+1] = 1;
				hardBits[4(n-1)+3] = 0;
			elseif decI == 2
				hardBits[4(n-1)+1] = 1;
				hardBits[4(n-1)+3] = 1;
			elseif decI == 3
				hardBits[4(n-1)+1] = 0;
				hardBits[4(n-1)+3] = 1;
			elseif decI == 4
				hardBits[4(n-1)+1] = 0;
				hardBits[4(n-1)+3] = 0;
			end
			# --- Decision for imag part
			# Getting only imag part
			imagP	= imag(qamVect[n]);
			# Getting area of interest: from area 1 to 4
			decQ	= 1+ (imagP > bounds[1]) + (imagP>bounds[2]) + (imagP>bounds[3]);
			# Convert into binary sequence based on Gray encoding scheme.
			if decQ == 1
				hardBits[4(n-1)+2] = 1;
				hardBits[4(n-1)+4] = 0;
			elseif decQ == 2
				hardBits[4(n-1)+2] = 1;
				hardBits[4(n-1)+4] = 1;
			elseif decQ == 3
				hardBits[4(n-1)+2] = 0;
				hardBits[4(n-1)+4] = 1;
			elseif decQ == 4
				hardBits[4(n-1)+2] = 0;
				hardBits[4(n-1)+4] = 0;
			end
		end
	elseif M == 64
		## 64-QAM hard demapper
		# --- Gray encoding
		#   -7     -5    -3    -1    +1    +3     +5     +7
		#  100    101   111   110   010   011    001    000
		# --- Voronoi region definition
		bounds = [-6 -4 -2 0 2 4 6]/(sqrt(2/3*63));
		# --- Iterative hard demapping
		for n = 1 : 1 : length(qamVect)
			# --- Decision for real part
			# Getting only real part
			realI	= real(qamVect[n]);
			# Getting area of interest: from area 1 to 4
			decI	= 1+ (realI > bounds[1]) + (realI>bounds[2]) + (realI>bounds[3])+ (realI>bounds[4])+ (realI>bounds[5])+ (realI>bounds[6])+ (realI>bounds[7]);
			# Convert into binary sequence based on Gray encoding scheme.
			if decI == 1
				hardBits[6(n-1)+1] = 1;
				hardBits[6(n-1)+3] = 0;
				hardBits[6(n-1)+5] = 0;
			elseif decI == 2
				hardBits[6(n-1)+1] = 1;
				hardBits[6(n-1)+3] = 0;
				hardBits[6(n-1)+5] = 1;
			elseif decI == 3
				hardBits[6(n-1)+1] = 1;
				hardBits[6(n-1)+3] = 1;
				hardBits[6(n-1)+5] = 1;
			elseif decI == 4
				hardBits[6(n-1)+1] = 1;
				hardBits[6(n-1)+3] = 1;
				hardBits[6(n-1)+5] = 0;
			elseif decI == 5
				hardBits[6(n-1)+1] = 0;
				hardBits[6(n-1)+3] = 1;
				hardBits[6(n-1)+5] = 0;
			elseif decI == 6
				hardBits[6(n-1)+1] = 0;
				hardBits[6(n-1)+3] = 1;
				hardBits[6(n-1)+5] = 1;
			elseif decI == 7
				hardBits[6(n-1)+1] = 0;
				hardBits[6(n-1)+3] = 0;
				hardBits[6(n-1)+5] = 1;
			elseif decI == 8
				hardBits[6(n-1)+1] = 0;
				hardBits[6(n-1)+3] = 0;
				hardBits[6(n-1)+5] = 0;
			end
			# --- Decision for imag part
			# Getting only imag part
			realQ	= imag(qamVect[n]);
			# Getting area of interest: from area 1 to 8
			decQ	= 1+ (realQ > bounds[1]) + (realQ>bounds[2]) + (realQ>bounds[3])+ (realQ>bounds[4])+ (realQ>bounds[5])+ (realQ>bounds[6])+(realQ>bounds[7]);
			# Convert into binary sequence based on Gray encoding scheme.
			if decQ == 1
				hardBits[6(n-1)+2] = 1;
				hardBits[6(n-1)+4] = 0;
				hardBits[6(n-1)+6] = 0;
			elseif decQ == 2
				hardBits[6(n-1)+2] = 1;
				hardBits[6(n-1)+4] = 0;
				hardBits[6(n-1)+6] = 1;
			elseif decQ == 3
				hardBits[6(n-1)+2] = 1;
				hardBits[6(n-1)+4] = 1;
				hardBits[6(n-1)+6] = 1;
			elseif decQ == 4
				hardBits[6(n-1)+2] = 1;
				hardBits[6(n-1)+4] = 1;
				hardBits[6(n-1)+6] = 0;
			elseif decQ == 5
				hardBits[6(n-1)+2] = 0;
				hardBits[6(n-1)+4] = 1;
				hardBits[6(n-1)+6] = 0;
			elseif decQ == 6
				hardBits[6(n-1)+2] = 0;
				hardBits[6(n-1)+4] = 1;
				hardBits[6(n-1)+6] = 1;
			elseif decQ == 7
				hardBits[6(n-1)+2] = 0;
				hardBits[6(n-1)+4] = 0;
				hardBits[6(n-1)+6] = 1;
			elseif decQ == 8
				hardBits[6(n-1)+2] = 0;
				hardBits[6(n-1)+4] = 0;
				hardBits[6(n-1)+6] = 0;
			end
		end
	elseif M == 256
		## 64-QAM hard demapper
		# --- Gray encoding
		# ------------------------------------------------#
		#	   |       |         |        |        |        |        |        |
		#	 -15     -13 .     -11 .     -9 .     -7 .     -5 .     -3 .     -1
		#	0100    0101 .    0111 .   0110 .   0010 .   0011 .   0001 .   0000
		#	   |       |         |        |        |        |        |        |
		#     +1 .     +3 .     +5 .     +7 .     +9 .    +11 .    +13      +15.
		#   1000 .   1001 .   1011 .   1010 .   1110 .   1111 .   1101 .   1100
		# --- Voronoi region definition
		bounds = [-14 -12 -10 -8 -6 -4 -2 0 2 4 6 8 10 12 14]/(sqrt(2/3*255));
		# --- Iterative hard demapping
		for n = 1 : 1 : length(qamVect)
			# --- Decision for real part
			# Getting only real part
			realI	= real(qamVect[n]);
			# Getting area of interest: from area 1 to 15
			decI	= 1+ (realI > bounds[1]) + (realI>bounds[2]) + (realI>bounds[3])+ (realI>bounds[4])+ (realI>bounds[5])+ (realI>bounds[6])+ (realI>bounds[7])+ (realI>bounds[8])+(realI>bounds[9])+(realI>bounds[10])+(realI>bounds[11])+(realI>bounds[12])+(realI>bounds[13])+ (realI>bounds[14])+ (realI>bounds[15]) ;
			# Convert into binary sequence based on Gray encoding scheme.
			if decI == 1
				hardBits[8(n-1)+1] = 0;
				hardBits[8(n-1)+3] = 1;
				hardBits[8(n-1)+5] = 0;
				hardBits[8(n-1)+7] = 0;
			elseif decI == 2
				hardBits[8(n-1)+1] = 0;
				hardBits[8(n-1)+3] = 1;
				hardBits[8(n-1)+5] = 0;
				hardBits[8(n-1)+7] = 1;
			elseif decI == 3
				hardBits[8(n-1)+1] = 0;
				hardBits[8(n-1)+3] = 1;
				hardBits[8(n-1)+5] = 1;
				hardBits[8(n-1)+7] = 1;
			elseif decI == 4
				hardBits[8(n-1)+1] = 0;
				hardBits[8(n-1)+3] = 1;
				hardBits[8(n-1)+5] = 1;
				hardBits[8(n-1)+7] = 0;
			elseif decI == 5
				hardBits[8(n-1)+1] = 0;
				hardBits[8(n-1)+3] = 0;
				hardBits[8(n-1)+5] = 1;
				hardBits[8(n-1)+7] = 0;
			elseif decI == 6
				hardBits[8(n-1)+1] = 0;
				hardBits[8(n-1)+3] = 0;
				hardBits[8(n-1)+5] = 1;
				hardBits[8(n-1)+7] = 1;
			elseif decI == 7
				hardBits[8(n-1)+1] = 0;
				hardBits[8(n-1)+3] = 0;
				hardBits[8(n-1)+5] = 0;
				hardBits[8(n-1)+7] = 1;
			elseif decI == 8
				hardBits[8(n-1)+1] = 0;
				hardBits[8(n-1)+3] = 0;
				hardBits[8(n-1)+5] = 0;
				hardBits[8(n-1)+7] = 0;
			elseif decI == 9
				hardBits[8(n-1)+1] = 1;
				hardBits[8(n-1)+3] = 0;
				hardBits[8(n-1)+5] = 0;
				hardBits[8(n-1)+7] = 0;
			elseif decI == 10
				hardBits[8(n-1)+1] = 1;
				hardBits[8(n-1)+3] = 0;
				hardBits[8(n-1)+5] = 0;
				hardBits[8(n-1)+7] = 1;
			elseif decI == 11
				hardBits[8(n-1)+1] = 1;
				hardBits[8(n-1)+3] = 0;
				hardBits[8(n-1)+5] = 1;
				hardBits[8(n-1)+7] = 1;
			elseif decI == 12
				hardBits[8(n-1)+1] = 1;
				hardBits[8(n-1)+3] = 0;
				hardBits[8(n-1)+5] = 1;
				hardBits[8(n-1)+7] = 0;
			elseif decI == 13
				hardBits[8(n-1)+1] = 1;
				hardBits[8(n-1)+3] = 1;
				hardBits[8(n-1)+5] = 1;
				hardBits[8(n-1)+7] = 0;
			elseif decI == 14
				hardBits[8(n-1)+1] = 1;
				hardBits[8(n-1)+3] = 1;
				hardBits[8(n-1)+5] = 1;
				hardBits[8(n-1)+7] = 1;
			elseif decI == 15
				hardBits[8(n-1)+1] = 1;
				hardBits[8(n-1)+3] = 1;
				hardBits[8(n-1)+5] = 0;
				hardBits[8(n-1)+7] = 1;
			elseif decI == 16
				hardBits[8(n-1)+1] = 1;
				hardBits[8(n-1)+3] = 1;
				hardBits[8(n-1)+5] = 0;
				hardBits[8(n-1)+7] = 0;
			end
			# --- Decision for imag part
			# Getting only imag part
			realQ	= imag(qamVect[n]);
			# Getting area of interest: from area 1 to 15
			decQ	= 1+ (realQ > bounds[1]) + (realQ>bounds[2]) + (realQ>bounds[3])+ (realQ>bounds[4])+ (realQ>bounds[5])+ (realQ>bounds[6])+ (realQ>bounds[7])+ (realQ>bounds[8])+(realQ>bounds[9])+(realQ>bounds[10])+(realQ>bounds[11])+(realQ>bounds[12])+(realQ>bounds[13])+ (realQ>bounds[14])+ (realQ>bounds[15]) ;
			# Convert into binary sequence based on Gray encoding scheme.
			if decQ == 1
				hardBits[8(n-1)+2] = 0;
				hardBits[8(n-1)+4] = 1;
				hardBits[8(n-1)+6] = 0;
				hardBits[8(n-1)+8] = 0;
			elseif decQ == 2
				hardBits[8(n-1)+2] = 0;
				hardBits[8(n-1)+4] = 1;
				hardBits[8(n-1)+6] = 0;
				hardBits[8(n-1)+8] = 1;
			elseif decQ == 3
				hardBits[8(n-1)+2] = 0;
				hardBits[8(n-1)+4] = 1;
				hardBits[8(n-1)+6] = 1;
				hardBits[8(n-1)+8] = 1;
			elseif decQ == 4
				hardBits[8(n-1)+2] = 0;
				hardBits[8(n-1)+4] = 1;
				hardBits[8(n-1)+6] = 1;
				hardBits[8(n-1)+8] = 0;
			elseif decQ == 5
				hardBits[8(n-1)+2] = 0;
				hardBits[8(n-1)+4] = 0;
				hardBits[8(n-1)+6] = 1;
				hardBits[8(n-1)+8] = 0;
			elseif decQ == 6
				hardBits[8(n-1)+2] = 0;
				hardBits[8(n-1)+4] = 0;
				hardBits[8(n-1)+6] = 1;
				hardBits[8(n-1)+8] = 1;
			elseif decQ == 7
				hardBits[8(n-1)+2] = 0;
				hardBits[8(n-1)+4] = 0;
				hardBits[8(n-1)+6] = 0;
				hardBits[8(n-1)+8] = 1;
			elseif decQ == 8
				hardBits[8(n-1)+2] = 0;
				hardBits[8(n-1)+4] = 0;
				hardBits[8(n-1)+6] = 0;
				hardBits[8(n-1)+8] = 0;
			elseif decQ == 9
				hardBits[8(n-1)+2] = 1;
				hardBits[8(n-1)+4] = 0;
				hardBits[8(n-1)+6] = 0;
				hardBits[8(n-1)+8] = 0;
			elseif decQ == 10
				hardBits[8(n-1)+2] = 1;
				hardBits[8(n-1)+4] = 0;
				hardBits[8(n-1)+6] = 0;
				hardBits[8(n-1)+8] = 1;
			elseif decQ == 11
				hardBits[8(n-1)+2] = 1;
				hardBits[8(n-1)+4] = 0;
				hardBits[8(n-1)+6] = 1;
				hardBits[8(n-1)+8] = 1;
			elseif decQ == 12
				hardBits[8(n-1)+2] = 1;
				hardBits[8(n-1)+4] = 0;
				hardBits[8(n-1)+6] = 1;
				hardBits[8(n-1)+8] = 0;
			elseif decQ == 13
				hardBits[8(n-1)+2] = 1;
				hardBits[8(n-1)+4] = 1;
				hardBits[8(n-1)+6] = 1;
				hardBits[8(n-1)+8] = 0;
			elseif decQ == 14
				hardBits[8(n-1)+2] = 1;
				hardBits[8(n-1)+4] = 1;
				hardBits[8(n-1)+6] = 1;
				hardBits[8(n-1)+8] = 1;
			elseif decQ == 15
				hardBits[8(n-1)+2] = 1;
				hardBits[8(n-1)+4] = 1;
				hardBits[8(n-1)+6] = 0;
				hardBits[8(n-1)+8] = 1;
			elseif decQ == 16
				hardBits[8(n-1)+2] = 1;
				hardBits[8(n-1)+4] = 1;
				hardBits[8(n-1)+6] = 0;
				hardBits[8(n-1)+8] = 0;
			end
		end
	end
	# ---  Output on decoded sequence
	return hardBits;
end

"""
---
Quadrature Amplitude Modulation (QAM) hard decoding function
Apply symbol hard demapping to a input symbol sequence (of size 1xN) with constellation size M.
Output is a binary (1xL) with N = L / log2(M)
Conventional gray demapping is used.
Input constellation is Array{Complex{Float64}}
Output bitsream is Array{Int8}
# --- Syntax
hardBits = 	 bitDemappingQAM!(hardBits,M,qamVect)
# --- Input parameters
- M			: Constellation size (i.e from 4 to 256)
- qamVect	: Complex observation vector to decode.
# --- Output parameters
- hardBits	: Vector of bits to populate [Array{UInt8}, length(qamVect)/log2(M)]
"""
function bitDemappingQAM(M,qamVect);
	# --- Create receive bit vector
	hardBits = zeros(UInt8, Int(length(qamVect)*log2(M)));
	# --- Call bang method
	bitDemappingQAM!(hardBits,M,qamVect);
	return hardBits;
end


text:
module Decoder

using Base.Iterators: take

using ...Client

include("field.jl")
include("process.jl")
include("ticktype.jl")

# Make a shortcut
const pop = popfirst!

function decode(msg, w, ver)

  it = Iterators.Stateful(Field.(msg))

  # The first field is the message ID
  id::Int = pop(it)

  # The second field (version) is ignored for id < 75 and != 3, 5, 11, 17
  if id  < 75 && id ∉ [3, 5, 11, 17]             ||
     id ==  3 && ver < Client.MARKET_CAP_PRICE   ||
     id ==  5 && ver < Client.ORDER_CONTAINER    ||
     id == 11 && ver < Client.LAST_LIQUIDITY     ||
     id == 17 && ver < Client.SYNT_REALTIME_BARS

    pop(it)
  end

  f = get(process, id, nothing)

  if isnothing(f)
    @warn "Decoder: unknown message" ID=id

  else
    try
      f(it, w, ver)
    catch e
      if e isa EOFError
        @warn "Decoder: reached end of message" ID=id
      else
        rethrow()
      end
    end

    isempty(it) || @warn "Decoder: messsage not fully parsed" ID=id
  end

end

end


text:
using mna
using Test

@testset "SYMBOLIC: Butterworth Low Pass Filter (cuttoff = 200MHz)" begin
    Z₀ = 50;
    Vₛ = 1;

    testbench = Dict([ "Vs" => (:V, 1, 0, Vₛ)
                     , "Rs" => (:R, 1, 2, Z₀)
                     , "Rl" => (:R, 3, 0, Z₀)
                     ]);

    lowPass = Dict([ "C1" => (:C, 2, 0, 15e-12)
                   , "C2" => (:C, 3, 0, 15e-12)
                   , "L1" => (:L, 2, 3, 82e-9)
                   ]);

    netlist = merge(testbench, lowPass);

    x = mna.analyze(netlist);
    H = mna.transferFunction(x);

    samples = 201;
    f = range(1, stop=2e9, length=samples);
    ω = 2π .* f;
    tf = hcat(H.(ω)...);

    S₁₁ = ((tf[2,:] ./ -tf[4,:]) .- Z₀) ./ ((tf[2,:] ./ -tf[4,:]) .+ Z₀);
    S₂₁ = (tf[3,:] ./ tf[2,:]) .* (1 .+ S₁₁);

    S₁₁dB = 20 .* log10.(abs.(S₁₁));
    S₂₁dB = 20 .* log10.(abs.(S₂₁));

    fc = 200e6;     # Target: 200MHz
    ftol = 1;       # Tolerance: can be off by 1Hz (depends on sample rate)
    fc_ = first(f[S₂₁dB .< -3]);
    @test (fc_ - 1) <= fc <= (fc_ + 1)
end

@testset "NUMERIC: Butterworth Low Pass Filter (cuttoff = 200MHz)" begin
    Z₀ = 50;
    Vₛ = 1;

    testbench = Dict([ "Vs" => (:V, 1, 0, Vₛ)
                     , "Rs" => (:R, 1, 2, Z₀)
                     , "Rl" => (:R, 3, 0, Z₀)
                     ]);

    lowPass = Dict([ "C1" => (:C, 2, 0, 15e-12)
                   , "C2" => (:C, 3, 0, 15e-12)
                   , "L1" => (:L, 2, 3, 82e-9)
                   ]);

    netlist = merge(testbench, lowPass);

    samples = 201;
    f = range(1, stop=2e9, length=samples);
    tf = hcat([mna.analyze(netlist, ω) for ω in  (2π .* f)]...);

    S₁₁ = ((tf[2,:] ./ -tf[4,:]) .- Z₀) ./ ((tf[2,:] ./ -tf[4,:]) .+ Z₀);
    S₂₁ = (tf[3,:] ./ tf[2,:]) .* (1 .+ S₁₁);

    S₁₁dB = 20 .* log10.(abs.(S₁₁));
    S₂₁dB = 20 .* log10.(abs.(S₂₁));

    fc = 200e6;     # Target: 200MHz
    ftol = 1;       # Tolerance: can be off by 1Hz (depends on sample rate)
    fc_ = first(f[S₂₁dB .< -3]);
    @test (fc_ - 1) <= fc <= (fc_ + 1)
end


text:
# ALL TESTS THAT CONCERN MODFICIATION

# include
using LatticePhysics
using Base.Test

# begin the Modification testblock
path_testset = @testset "Interaction Matrices" begin

################################################################################
#
#   METHODS FOR CONSTRUCTION INTERACTION MATRICES FOR LATTICES
#
################################################################################

################################################################################
#
#   1) INTERACTION MATRICES IN REAL SPACE
#
################################################################################

# begin the testset
@testset "Interaction matrices in real space" begin

    # begin the testset
    @testset "2D (honeycomb)" begin

        # get a unitcell
        unitcell = getUnitcellHoneycomb()
        # get a lattice
        lattice = getLattice(unitcell, [-4,-2])

        # test the function for unitcell
        @test typeof(getInteractionMatrixRealSpace(unitcell)) == Array{Complex, 2}
        @test ishermitian(getInteractionMatrixRealSpace(unitcell, enforce_hermitian=true))
        @test size(getInteractionMatrixRealSpace(unitcell)) == (2,2)

        # test the function for lattice
        @test typeof(getInteractionMatrixRealSpace(lattice)) == Array{Complex, 2}
        @test ishermitian(getInteractionMatrixRealSpace(lattice, enforce_hermitian=true))
        @test size(getInteractionMatrixRealSpace(lattice)) == (16,16)

    # end the testset here
    end;

    # begin the testset
    @testset "2D (square)" begin

        # get a unitcell
        unitcell = getUnitcellSquare()
        # get a lattice
        lattice = getLattice(unitcell, [-4,-2])

        # test the function for unitcell
        @test typeof(getInteractionMatrixRealSpace(unitcell)) == Array{Complex, 2}
        @test ishermitian(getInteractionMatrixRealSpace(unitcell, enforce_hermitian=true))
        @test size(getInteractionMatrixRealSpace(unitcell)) == (1,1)

        # test the function for lattice
        @test typeof(getInteractionMatrixRealSpace(lattice)) == Array{Complex, 2}
        @test ishermitian(getInteractionMatrixRealSpace(lattice, enforce_hermitian=true))
        @test size(getInteractionMatrixRealSpace(lattice)) == (8,8)

    # end the testset here
    end;

    # begin the testset
    @testset "3D (diamond)" begin

        # get a unitcell
        unitcell = getUnitcellDiamond()
        # get a lattice
        lattice = getLattice(unitcell, [-4,-2,-2])

        # test the function for unitcell
        @test typeof(getInteractionMatrixRealSpace(unitcell)) == Array{Complex, 2}
        @test ishermitian(getInteractionMatrixRealSpace(unitcell, enforce_hermitian=true))
        @test size(getInteractionMatrixRealSpace(unitcell)) == (2,2)

        # test the function for lattice
        @test typeof(getInteractionMatrixRealSpace(lattice)) == Array{Complex, 2}
        @test ishermitian(getInteractionMatrixRealSpace(lattice, enforce_hermitian=true))
        @test size(getInteractionMatrixRealSpace(lattice)) == (32,32)

    # end the testset here
    end;

# end the testset here
end;

################################################################################
#
#   2) INTERACTION MATRICES IN MOMENTUM SPACE
#
################################################################################

# begin the testset
@testset "Interaction matrices in momentum space" begin

    # begin the testset
    @testset "2D (honeycomb)" begin

        # get a unitcell
        unitcell = getUnitcellHoneycomb()
        # get a lattice
        lattice = getLattice(unitcell, [-4,-2])

        # k vector (random)
        k = [0.1, 0.151235]

        # test the function for unitcell
        @test typeof(getInteractionMatrixKSpace(unitcell, k)) == Array{Complex, 2}
        @test ishermitian(getInteractionMatrixKSpace(unitcell, k, enforce_hermitian=true))
        @test size(getInteractionMatrixKSpace(unitcell, k)) == (2,2)

        # test the function for lattice
        @test typeof(getInteractionMatrixKSpace(lattice, k)) == Array{Complex, 2}
        @test ishermitian(getInteractionMatrixKSpace(lattice, k, enforce_hermitian=true))
        @test size(getInteractionMatrixKSpace(lattice, k)) == (16,16)

    # end the testset here
    end;

    # begin the testset
    @testset "2D (square)" begin

        # get a unitcell
        unitcell = getUnitcellSquare()
        # get a lattice
        lattice = getLattice(unitcell, [-4,-2])

        # k vector (random)
        k = [0.1, 0.151235]

        # test the function for unitcell
        @test typeof(getInteractionMatrixKSpace(unitcell, k)) == Array{Complex, 2}
        @test ishermitian(getInteractionMatrixKSpace(unitcell, k, enforce_hermitian=true))
        @test size(getInteractionMatrixKSpace(unitcell, k)) == (1,1)

        # test the function for lattice
        @test typeof(getInteractionMatrixKSpace(lattice, k)) == Array{Complex, 2}
        @test ishermitian(getInteractionMatrixKSpace(lattice, k, enforce_hermitian=true))
        @test size(getInteractionMatrixKSpace(lattice, k)) == (8,8)

    # end the testset here
    end;

    # begin the testset
    @testset "3D (diamond)" begin

        # get a unitcell
        unitcell = getUnitcellDiamond()
        # get a lattice
        lattice = getLattice(unitcell, [-4,-2,-2])

        # k vector (random)
        k = [0.1, 0.151235, -0.21]

        # test the function for unitcell
        @test typeof(getInteractionMatrixKSpace(unitcell, k)) == Array{Complex, 2}
        @test ishermitian(getInteractionMatrixKSpace(unitcell, k, enforce_hermitian=true))
        @test size(getInteractionMatrixKSpace(unitcell, k)) == (2,2)

        # test the function for lattice
        @test typeof(getInteractionMatrixKSpace(lattice, k)) == Array{Complex, 2}
        @test ishermitian(getInteractionMatrixKSpace(lattice, k, enforce_hermitian=true))
        @test size(getInteractionMatrixKSpace(lattice, k)) == (32,32)

    # end the testset here
    end;

# end the testset here
end;

# end the path test block here
end;


text:
julia> args = [ "first", (1,2,17), "last" ]
3-element Array{Any,1}:
 "first"
 (1,2,17)
 "last

julia> print_each(args...)
first
(1,2,17)
last

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
abstract CoordinateSystem

type OneD <: CoordinateSystem
	L
end

type Cartesian <: CoordinateSystem
	dimensions
	limits
end

type TwoDPolar <: CoordinateSystem
	a
end

function λ(c::OneD, j::Integer)
	((j * π / (2c.L)) ^ 2)::Float64
end

function φ(c::OneD, x::Array{Float64, 1}, j::Integer)
	x = x[1]
	(sin(π * j * (x .+ c.L) / (2c.L)) / sqrt(c.L))::Float64
end

function λ(c::Cartesian, js::Array{Int64, 1})
	ans = 0
	for d = 1:c.dimensions
		d
		L = c.limits[d]
		ans += ((js[d] * π / (2L)) ^ 2)::Float64
	end
	ans
end

function φ(c::Cartesian, x::Array{Float64, 1}, js::Array{Int64, 1})
	ans = 1
	for d = 1:c.dimensions
		L = c.limits[d]
		ans *= (sin(π * js[d] * (x[d] .+ L) / (2L)) / sqrt(L))::Float64
	end
	ans
end


text:
using BinDeps

@BinDeps.setup

@static if is_apple()
  ostype = "apple"
  elseif is_windows()
    ostype = "windows"
  elseif is_linux()
    ostype = "linux"
  else
    error("No wkhtmltox library available for this OS")
  end

ostype *= Int==Int64 ? "64" : "32"

ver = "2.1.1"

archivemap = Dict(
 "apple32"   => "-macosx.zip",
 "windows32" => "-windows.zip",
 "linux32"   => "-linux-i686.tar.bz2",
 "apple64"   => "-macosx.zip",
 "windows64" => "-windows.zip",
 "linux64"   => "-linux-x86_64.tar.bz2",
 )

url = "https://cnpmjs.org/downloads/phantomjs-$ver" * archivemap[ostype]
downloadname = basename(url)

exemap = Dict(
 "apple32"   => "phantomjs-$ver-macosx/bin/phantomjs",
 "windows32" => "phantomjs-$ver-windows/bin/phantomjs.exe",
 "linux32"   => "phantomjs-$ver-linux-i686/bin/phantomjs",
 "apple64"   => "phantomjs-$ver-macosx/bin/phantomjs",
 "windows64" => "phantomjs-$ver-windows/bin/phantomjs.exe",
 "linux64"   => "phantomjs-$ver-linux-x86_64/bin/phantomjs",
 )

exepath = exemap[ostype]
exefile = basename(exepath)

destdir      = joinpath(dirname(@__FILE__), "usr/bin")
unzipdir     = joinpath(dirname(@__FILE__), "src")
downloadsdir = joinpath(dirname(@__FILE__), "downloads")

type FileCopyRule <: BinDeps.BuildStep
    src::AbstractString
    dest::AbstractString
end
Base.run(fc::FileCopyRule) = isfile(fc.dest) || cp(fc.src, fc.dest)

type MakeExeRule <: BinDeps.BuildStep
    path::AbstractString
end
Base.run(rule::MakeExeRule) = chmod(rule.path, 0o755)

run(@build_steps begin
      CreateDirectory(downloadsdir, true)
      CreateDirectory(unzipdir, true)
      CreateDirectory(destdir, true)
      @build_steps FileRule( [joinpath(destdir, exefile)],
          @build_steps begin
            FileDownloader(url, joinpath(downloadsdir, downloadname))
            FileUnpacker(joinpath(downloadsdir, downloadname), unzipdir, exepath)
            FileCopyRule(joinpath(unzipdir,exepath), joinpath(destdir,exefile))
            MakeExeRule(joinpath(destdir,exefile))
          end )
      end)


text:
using FiniteFloats
using Test

@test isfinite(Finite64(Inf))
@test isfinite(Finite32(-Inf32))
@test isfinite(Finite16(Inf16))


text:
using BinaryBuilder

# Collection of sources required to build Nettle
name = "Nettle"
version = v"3.4.1"
sources = [
    ArchiveSource("https://ftp.gnu.org/gnu/nettle/nettle-$(version).tar.gz",
                  "f941cf1535cd5d1819be5ccae5babef01f6db611f9b5a777bae9c7604b8a92ad"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/nettle-*/
atomic_patch -p1 ../patches/alloca.patch

# Force c99 mode
export CFLAGS="${CFLAGS} -std=c99"

update_configure_scripts
./configure --prefix=$prefix --host=$target --with-include-path=$prefix/include
make -j${nproc}
make install
install_license COPYING*
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    LibraryProduct("libnettle", :libnettle),
    LibraryProduct("libhogweed", :libhogweed),
    ExecutableProduct("nettle-hash", :nettle_hash)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("GMP_jll", v"6.1.2"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: blenessy/Yggdrasil

max_issues_repo_name: blenessy/Yggdrasil

max_forks_repo_name: blenessy/Yggdrasil


text:
using Glob, CSV, DataFrames, Dates

# Change this to the directory which stores EMI and Cliflo data CSVs
# This does not have to be stored in the same location as this script
dir = "C:\\Users\\calro\\Documents\\Academic\\JuDGE project\\LDC_New"
# dir = "C:\\Users\\calro\\Documents\\Academic\\JuDGE project\\LDC_Jamie\\Solar\\Raw"

validation = true   # Run the validation script
only2019data = true     # Multiple years of data (Rep days) or just 2019 (LDCs)

function windfarm_dict()
    # Dictionary of wind farm gen_codes => actual names
    # and the existing capacity at each wind farm
    # Note that the MW amounts are EXISTING capacity, not investments.

    # 1. Names
    # keys: Gen_Code, vals: Location

    # 2. Capcity
    # keys: Location, vals: capacity

    Gen_Code = ["te_apiti"
    , "te_rere_hau"
    , "te_uku"
    , "twf_12"
    , "twf_3"
    , "west_wind"
    , "white_hill"]

    Location = ["TeApiti"
    , "TeRereHau"
    , "TeUku"
    , "Tararua_Embed"
    , "Tararua"
    , "Westwind"
    , "Whitehill"]

    Capacity = [90.75
    , 48.5
    , 0.0
    , 0.0
    , 93.0
    , 143.0
    , 0.0]

    # Gen_Code = ["te_apiti"
    # , "te_rere_hau"
    # , "te_uku"
    # , "twf_12"
    # , "twf_3"
    # , "west_wind"
    # , "white_hill"]
    #
    # Location = ["TeApiti"
    # , "TeRereHau"
    # , "TeUku"
    # , "Tararua"
    # , "Tararua"
    # , "Westwind"
    # , "Whitehill"]
    #
    # Capacity = [90.75
    # , 48.5
    # , 64.4
    # , 161.0
    # , 143
    # , 58]

    # return Dict(zip(Gen_Code, Location)), Dict(zip(unique(Location), Capacity))
    return Dict(zip(Gen_Code, Location)), Dict(zip(Location, Capacity))
end

function solar_dict()
    # Dictionary of solar cells that we have collected data for.
    # Note that the MW amounts are EXISTING capacity, not investments.

    Station = ["Pukekohe Ews"
    , "Lauder Ews"
    , "Blenheim Research Ews"
    , "Christchurch Aero"]

    Region = ["Auckland"
    , "Otago"
    , "Blenheim"
    , "Christchurch"]

    Island = ["North"
    , "South"
    , "South"
    , "South"]

    # North = 77.72 MW
    # South = 36.68 MW
    # MW = [77.72
    # , 36.68]

    MW = [0.0
    , 0.0]

    # First returns a dictionary of station name and island
    # Second returns the MW production in each island
    return Dict(zip(Station, Island)), Dict(zip(unique(Island), MW))
end

function investments()
    # Input an array of all individual investments, this script will then
    # iterate over the power set of all investments.

    # This is a tad tricky. I would recommend initially running all functions up
    # to and including "Merge()" to see what the column names are for each gen
    # type.
    # CURRENTLY:
        # The WIND sites can be referred to by their LOCATION in
        # windfarm_dict(),
        # SOLAR sites are to be referred to as: "Solar_NI" and "Solar_SI" for
        # the amount in either of these islands. This will probably change

    # Probably worth adding some assert statements that validate that the user
    # input is correct, but for now, you're on your own.

    Investment_Type = ["Wind"
                        , "Wind"
                        , "Wind"
                        , "Wind"
                        , "Wind"
                        ]

    Investment_Location = ["TeApiti"
                            , "Westwind"
                            , "TeUku"
                            , "Tararua"
                            , "Whitehill"
                            ]

    Investment_Amount = [600.0
                        , 700.0
                        , 800.0
                        , 900.0
                        , 1000.0
                        ]

    return Dict(zip(Investment_Location, Investment_Type)),
         Dict(zip(Investment_Location, Investment_Amount))
end

function validate_capacity()
    # This function prints out the current existing capacity of windsites,
    # solarsites, and any proposed combination of investments.
    WindLocs, Windcap = windfarm_dict()
    SolarLocs, Solarcap = solar_dict()
    InvLocs, Invs = investments()

    cnt = length(InvLocs)

    println("\n\nCurrent windfarm capacity:")
    println(stack(DataFrame(Windcap)))
    println("\n\nCurrent solarfarm capacity:")
    println(stack(DataFrame(Solarcap)))

    # Convert investment df into a more readable form
    Invs = stack(DataFrame(Invs))
    rename!(Invs, Dict(:variable => "Investment_Location", :value => "Investment_Amount"))
    Invs[!, :Generation_Type] .= ""
    InvLocs = DataFrame(InvLocs)
    for i in names(InvLocs)
        Invs[Invs.Investment_Location .== i, :Generation_Type] .= InvLocs[1, i]
    end

    println("\n\nAll investments (individual):")
    println(Invs)

    println("\n\nIndividual investments: ", cnt, ".")
    println("Total portfolio size: ", 2^cnt, " investment options.")
end

function powerset()

    println("\n\nGenerating portfolio of investments ...")
    InvLocs, _ = investments()
    length(InvLocs)
    bin_rep = []

    #Power set of all investments
    for i in 0:2^length(InvLocs) - 1
     # Calculate the binary representation for the 'i'th combination of
     # investments. Pad the MSBs with zeros for binary representations that
     # are shorter than the number of investments
     push!(bin_rep, string(i, base = 2, pad = length(InvLocs)))
    end
    return bin_rep
end

function Format_TP_Demand(dir)
    println("\n\nFormatting demand data ...")
    # Read in raw demand data
    demand = reduce(vcat, CSV.read.(glob("*_import.csv", dir)))
    # Reshape the demand df such that each trading period becomes a row
    select!(demand, Not([:POC, :NWK_Code, :GENERATION_TYPE, :TRADER, :UNIT_MEASURE, :FLOW_DIRECTION, :STATUS]))
    # This is the new version of "melt" : Produces the same output as: demand = melt(demand, [:TRADING_DATE])
    demand = stack(demand, names(demand)[2:end])
    # Sum all TP on the same day, convert load from kWh to MW
    demand = combine(groupby(demand, [:TRADING_DATE, :variable]), :value => sum)
    demand[!, :value_sum] = demand[!, :value_sum] * 2 / 1000
    # Create a new column and input the numerical value of which trading period it is
    demand[!, :TP] .= parse.(Int64, reduce(hcat, split.(String.(demand[!, :variable]), "P"))[2, :])
    # Drop unnecessary columns and reorder the remaining columns
    # select!(demand, Not(:variable))
    select!(demand, [:TRADING_DATE, :TP, :value_sum])
    rename!(demand,  Dict(:TRADING_DATE => "Trading_date", :TP => "Trading_period", :value_sum => "Load"))
    # Drop TP's > 48 (Daylight savings)
    delete!(demand, demand[!, :Trading_period] .> 48)
    # Parse the trading date and add seasons
    demand[!, :Trading_date] = Date.(demand[!, :Trading_date],"d/m/y")
    seasons = ["Summer", "Autumn", "Winter", "Spring", "Summer"]
    insertcols!(demand, 4, :Season => seasons[convert.(Int64,floor.(Dates.month.(demand[!, :Trading_date]) / 3) .+ 1)])
    # Sort chronologically and return
    sort!(demand, [:Trading_date, :Trading_period], rev = false)
    demand[!, :Load_raw] = demand[!, :Load]                         # Store the gross load
    # CSV.write(joinpath(@__DIR__, "demanddata.csv"),  DataFrame(demand), header=true)
    return demand
end

function Format_TP_Price(dir)
    println("Formatting price data ...")
    # Read in price data
    price = reduce(vcat, CSV.read.(glob("*_prices.csv", dir)))
    # Prices at PEN0221
    price = price[price[!, :Node] .== "PEN0221",:]
    # Convert from KWh to MW - Not sure about this
    # price[!, :Price] = price[!, :Price] * 2 / 1000
    # Sort chronologically and return
    select!(price, Not(:Node))
    sort!(price, [:Trading_date, :Trading_period], rev = false)
    # CSV.write(joinpath(@__DIR__, "pricedata.csv"),  DataFrame(price), header=true)
    return price
end

function Format_TP_Wind(dir)
    println("Formatting wind data ...")
    # Extract wind wind data
    wind = reduce(vcat, CSV.read.(glob("*Generation_MD.csv", dir)))
    wind = wind[wind[!, :Tech_Code] .== "Wind" ,:]
    select!(wind, Not([:Site_Code, :Nwk_Code, :Fuel_Code, :Tech_Code]))
    # Melt data
    wind = stack(wind, names(wind)[4:end])
    rename!(wind, :value => :Load)
    # Extract the number from TP and assign to Trading_period and drop TP's > 48
    wind[!, :Trading_period] .= parse.(Int64, reduce(hcat, split.(String.(wind[!, :variable]), "P"))[2, :])
    delete!(wind, wind[!, :Trading_period] .> 48)
    # Parse the trading date and add seasons
    # Only run the commented line below if Trading_date is not parsed already
    # wind[!, :Trading_date] = Date.(wind[!, :Trading_date],"d/m/y")
    seasons = ["Summer", "Autumn", "Winter", "Spring", "Summer"]
    insertcols!(wind, 7, :Season => seasons[convert.(Int64,floor.(Dates.month.(wind[!, :Trading_date]) / 3) .+ 1)])
    # Change names of locations
    locationdict, _ = windfarm_dict()
    wind[:, :Location] .= ""
    [wind[i, :Location] = locationdict[wind[i, :Gen_Code]] for i in 1:size(wind)[1]]
    # Sum over all TP on the same day, convert load from kWh to MW
    select!(wind, [:Trading_period, :Trading_date, :Season, :Load, :Location])
    dropmissing!(wind)
    wind = combine(groupby(wind, [:Trading_date, :Trading_period, :Season, :Location]), :Load => sum)
    wind[!, :Load_sum] = wind[!, :Load_sum] * 2 / 1000
    # Transform rows back to wide format and remove the two empty daylight
    # savings rows in September. If more than two TPs disappear, check what's up
    wind = unstack(wind, [:Trading_period, :Trading_date, :Season], :Location, :Load_sum)
    dropmissing!(wind)
    # Divide each element by the max of each column to get the wind scale factors
    [wind[!, i+3] = wind[!, i+3] / maximum(wind[!, i+3]) for i in 1:7]
    sort!(wind, [:Trading_date, :Trading_period], rev = false)
    # CSV.write(joinpath(@__DIR__, "windfactordata.csv"),  DataFrame(wind), header=true)
    return wind
end

function Format_TP_Solar(dir)
    println("Formatting solar data ...")
    # Read in solar data
    solar = reduce(vcat, CSV.read.(glob("NIWA_Solar_*.csv", dir)))
    rename!(solar, ["aa", "Station", "date", "Irradiance", "bb", "cc", "dd"])
    # Extract date and hour
    solar[!, :Trading_date] = Date.(solar[!, :date],"yyyymmdd:HM")
    transform!(solar, :date => x -> chop.(x, head=9, tail=0))
    select!(solar, [:Station, :Trading_date, :date_function, :Irradiance])
    rename!(solar, :date_function => :Hour)
    # Format the trading_periods (hourly data)
    solar[!, :Hour] .= parse.(Int64, string.(solar[!, :Hour]))
    solar[!, :Hour] .= convert.(Int64, floor.(solar[!, :Hour] / 100))
    # Format island column by calling the dictionary
    locationdict, _ = solar_dict()
    solar[:, :Island] .= ""
    [solar[i, :Island] = locationdict[solar[i, :Station]] for i in 1:size(solar)[1]]
    # Calculate the W produced from the MJ/(hr m^2)
    # # # Commenting this out for the moment seeings as we are working with
    # "solar-factors" anyway?
    # solar[!, :Irradiance] .= solar[!, :Irradiance].* 1000 .*area.*efficiency./3600
    # Halve it for the half hourly data
    solar[!, :Irradiance] .= solar[!, :Irradiance] ./ 2
    solar2 = DataFrame(Station = String[], Trading_date = Date[], Hour = Float64[], Irradiance = Float64[], Island = String[])
    # Duplicate each row such that
    for i in 1:size(solar)[1]
        push!(solar2, solar[i, :])  # Push the original, hourly data
        push!(solar2, solar[i, :])   # Add an additional row to account for the half hour TP
        solar2[end, :Hour] = solar2[end, :Hour] + 0.5
    end
    # Assign the trading_period labels
    solar2[!, :Hour] .= convert.(Int64, 2 .* (solar2[!, :Hour] .+ 0.5))
    rename!(solar2, :Hour => :Trading_period)
    # Melt the data back from long to wide format
    unique!(solar2)
    solar = unstack(solar2, [:Trading_date, :Trading_period], :Island, :Irradiance)
    solar[!, :Solar] = solar[!, :North] + solar[!, :South]
    select!(solar, [:Trading_date, :Trading_period, :Solar])
    dropmissing!(solar)
    # 2019 data
    delete!(solar, Dates.year.(solar[!, :Trading_date]) .<2019)
    delete!(solar,  ((Dates.year.(solar[!, :Trading_date]) .<=2019) .* (Dates.month.(solar[!, :Trading_date]) .<=2)))
    # Divide each element by the max of each column to get the solar scale factors
    # [solar[!, i+2] = round.(solar[!, i+2] / maximum(solar[!, i+2]); digits = 6) for i in 1:2]
    solar[!, :Solar] = solar[!, :Solar] ./ maximum(solar[!, :Solar])
    # CSV.write(joinpath(@__DIR__, "CSV_files", "solardata.csv"),  DataFrame(solar), header=true)
    return solar
end

function Merge(demand, price, wind, solar, file)
    println("Merging raw data ...")
    # Merge price, demand, wind, and solar data
    merged = innerjoin(demand, price, on = [:Trading_date, :Trading_period])
    merged = innerjoin(merged, wind, on = [:Trading_date, :Trading_period, :Season])
    merged = innerjoin(merged, solar, on = [:Trading_date, :Trading_period])
    # rename!(merged, [:North => "Solar_NI", :South => "Solar_SI"])
    # Return the df chronologically sorted
    # select!(merged, [:Trading_period, :Load, :Price, :Trading_date, :Season])
    sort!(merged, [:Trading_date, :Trading_period], rev = false)
    CSV.write(joinpath(@__DIR__, file),  DataFrame(merged), header=true)
    println("Raw data exported to CSV: ", file)
    return merged
end

function Capacity_Processing(filename, portfolio)
    println("\n\nSubtracting preexisting generation ...")
    _, Windcap = windfarm_dict()
    # _, Solarcap = solar_dict()
    # existingcap = hcat(DataFrame(Windcap), DataFrame(Solarcap))
    # Call a function to subtract existing generation
    # df = Capacity_Subtraction(filename[1], existingcap)
    df = Capacity_Subtraction(filename[1], DataFrame(Windcap))
    # Store this as a CSV too
    CSV.write(joinpath(@__DIR__, filename[2]),  DataFrame(df), header=true)
    # Prepare investment
    _, Invs = investments()
    # Iterate over portfolio
    println("Subtracting investment portfolio ...")
    for set in portfolio
        invscap = deepcopy(DataFrame(Invs))
        for (j, inv) in enumerate(set)
            invscap[1, j] = invscap[1, j] * parse(Int, inv, base = 10)
        end
        # Call a function to subtract what is required
        netdf = Capacity_Subtraction(filename[2], invscap)
        # Save the CSV based on the binary representation in the name
        # Note that this saves into the Portfolio_Files dir
        CSV.write(joinpath(dirname(@__DIR__), "Portfolio_Files", filename[3]*set*".csv"),  DataFrame(netdf), header=true)
    end
    println("Files are stored in: ", joinpath(dirname(@__DIR__), "Portfolio_Files"))
end

function Capacity_Subtraction(file, capacity)
    # Read in the dataframe from a csv
    df = CSV.read(joinpath(@__DIR__, file), DataFrame)
    df[!, :Load] = collect(skipmissing(df[!, :Load]))
    # Iterate over all locs in capacity and subtract each's generation from Load
    for loc in names(capacity)
        df[!, :Load] = df[!, :Load] .- (df[!, loc] .* capacity[1, loc])
    end

    # Any negative values are just set to 0
    df[df[!, :Load] .< 0, :Load] .= 0

    return df
end

function Portfolio_CSV(filename)
    println("Dividing each portfolio into seasons ...")

    dir = joinpath(dirname(@__DIR__), "Portfolio_Files")

    # Extract a list of all files in the Portfolio database
    files = []
    for f in split.(glob(filename[3]*"*", dir), "\\")
        push!(files, split(f[end], ".")[1])
    end

    for file in files
        df = CSV.read(joinpath(dir, file*".csv"))
        select!(df, [:Trading_period, :Load, :Price, :Trading_date, :Season, :Solar, :Load_raw])

        # Create a column YearSeasonStarts, such that we can correctly assign
        # the summer season which crosses years
        df[!, :YearSeasonStarts] .= Dates.year.(df[!, :Trading_date])
        df[Int64.(Dates.month.(df[!, :Trading_date])) .<= 2, :YearSeasonStarts] .-= 1
        df[!, :YearSeasonStarts] .= string.(df[!, :YearSeasonStarts])

        for seas in unique(df[!, :Season])

            output = deepcopy(df[df[!, :Season] .== seas, :])
            output = output[output[!, :YearSeasonStarts] .== "2019", :]

            sort!(output, [:Load], rev = true)
            select!(output, [:Trading_period, :Load, :Price, :Trading_date, :Season, :Solar, :Load_raw])
            if size(output)[1] > 0
                CSV.write(joinpath(dirname(@__DIR__), "LDC_CSV_Files"
                        , file * "_" * seas * ".csv")
                        , DataFrame(output), header=true)
            end
        end
    end
    println("CSVs stored in: ", joinpath(dirname(@__DIR__), "LDC_CSV_Files"))
end

function solar_dict_Adapted_For_Jamie()
    # Dictionary of solar cells that we have collected data for.
    # Note that the MW amounts are EXISTING capacity, not investments.

    Station = ["Pukekohe Ews"
    , "Baring Head"
    , "Lauder Ews"
    , "Blenheim Research Ews"
    , "Lincoln, Broadfield Ews"]

    Region = ["Auckland"
    , "Wellington"
    , "Queenstown"
    , "Blenheim"
    , "Christchurch"]

    # First returns a dictionary of station name and region
    return Dict(zip(Station, Region))
end

function Portfolio_CSV_Demand_Adapted_For_Jamie(filename)
    println("Dividing each portfolio into seasons ...")

    # dir = joinpath(dirname(@__DIR__), "Portfolio_Files")
    dir = joinpath(dirname(@__DIR__), "temp")

    file = "demanddata"
    df = CSV.read(joinpath(dir, file*".csv"))
    # select!(df, [:Trading_period, :Load, :Price, :Trading_date, :Season])
    select!(df, [:Trading_period, :Load, :Trading_date, :Season])

    # Create a column YearSeasonStarts, such that we can correctly assign
    # the summer season which crosses years
    df[!, :YearSeasonStarts] .= Dates.year.(df[!, :Trading_date])
    df[Int64.(Dates.month.(df[!, :Trading_date])) .<= 2, :YearSeasonStarts] .-= 1
    df[!, :YearSeasonStarts] .= string.(df[!, :YearSeasonStarts])

    for yrs in unique(df[!, :YearSeasonStarts])
        for seas in unique(df[!, :Season])

            output = deepcopy(df[df[!, :YearSeasonStarts] .== yrs, :])
            output = deepcopy(output[output[!, :Season] .== seas, :])

            sort!(demand, [:Trading_date, :Trading_period], rev = false)
            select!(output, [:Trading_period, :Load, :Trading_date, :Season])
            if size(output)[1] > 0
                CSV.write(joinpath(dirname(@__DIR__), "temp"
                        , "Demand_" * seas * "_" * yrs * ".csv")
                        , DataFrame(output), header=true)
            end
        end
    end
    println("CSVs stored in: ", joinpath(dirname(@__DIR__), "LDC_CSV_Files"))
end

function Portfolio_CSV_Solar_Adapted_For_Jamie(filename)
    println("Dividing each portfolio into seasons ...")

    # dir = joinpath(dirname(@__DIR__), "Portfolio_Files")
    dir = joinpath(dirname(@__DIR__), "temp")

    file = "solardata"
    df = CSV.read(joinpath(dir, file*".csv"))
    # select!(df, [:Trading_period, :Load, :Price, :Trading_date, :Season])
    select!(df, [:Trading_period, :Irradiance, :Trading_date, :Region])

    seasons = ["Summer", "Autumn", "Winter", "Spring", "Summer"]
    insertcols!(df, 4, :Season => seasons[convert.(Int64,floor.(Dates.month.(df[!, :Trading_date]) / 3) .+ 1)])

    # Create a column YearSeasonStarts, such that we can correctly assign
    # the summer season which crosses years
    df[!, :YearSeasonStarts] .= Dates.year.(df[!, :Trading_date])
    df[Int64.(Dates.month.(df[!, :Trading_date])) .<= 2, :YearSeasonStarts] .-= 1
    df[!, :YearSeasonStarts] .= string.(df[!, :YearSeasonStarts])

    for site in unique(df[!, :Region])
        for seas in unique(df[!, :Season])

            output = deepcopy(df[df[!, :Season] .== seas, :])
            output = deepcopy(output[output[!, :Region] .== site, :])

            sort!(output, [:Trading_date, :Trading_period], rev = false)
            select!(output, [:Trading_period, :Irradiance, :Trading_date, :Season, :Region])
            if size(output)[1] > 0
                CSV.write(joinpath(dirname(@__DIR__), "temp"
                        , "Solar_" * seas * "_" * site * ".csv")
                        , DataFrame(output), header=true)
            end
        end
    end
    println("CSVs stored in: ", joinpath(dirname(@__DIR__), "temp"))
end

function Format_TP_Solar_Adapted_For_Jamie(dir)
    println("Formatting solar data ...")
    # Read in solar data
    solar = reduce(vcat, CSV.read.(glob("NIWA_Solar_*.csv", dir)))
    rename!(solar, ["Station", "date", "Irradiance", "bb", "cc", "dd"])
    # Extract date and hour
    solar[!, :Trading_date] = Date.(solar[!, :date],"yyyymmdd:HM")
    transform!(solar, :date => x -> chop.(x, head=9, tail=0))
    select!(solar, [:Station, :Trading_date, :date_function, :Irradiance])
    rename!(solar, :date_function => :Hour)
    # Format the trading_periods (hourly data)
    solar[!, :Hour] .= parse.(Int64, string.(solar[!, :Hour]))
    solar[!, :Hour] .= convert.(Int64, floor.(solar[!, :Hour] / 100))
    # Format island column by calling the dictionary
    locationdict = solar_dict_Adapted_For_Jamie()
    solar[:, :Region] .= ""
    [solar[i, :Region] = locationdict[solar[i, :Station]] for i in 1:size(solar)[1]]
    # Calculate the W produced from the MJ/(hr m^2)
    # # # Commenting this out for the moment seeings as we are working with
    # "solar-factors" anyway?
    # solar[!, :Irradiance] .= solar[!, :Irradiance].* 1000 .*area.*efficiency./3600
    # Halve it for the half hourly data
    solar[!, :Irradiance] .= solar[!, :Irradiance] ./ 2
    solar2 = DataFrame(Station = String[], Trading_date = Date[], Hour = Float64[], Irradiance = Float64[], Region = String[])
    # Duplicate each row such that
    for i in 1:size(solar)[1]
        push!(solar2, solar[i, :])  # Push the original, hourly data
        push!(solar2, solar[i, :])   # Add an additional row to account for the half hour TP
        solar2[end, :Hour] = solar2[end, :Hour] + 0.5
    end
    # Assign the trading_period labels
    solar2[!, :Hour] .= convert.(Int64, 2 .* (solar2[!, :Hour] .+ 0.5))
    rename!(solar2, :Hour => :Trading_period)
    # Melt the data back from long to wide format
    unique!(solar2)

    temp = unstack(solar2, [:Station, :Trading_date, :Region], :Trading_period, :Irradiance)
    dropmissing!(temp)
    solar3 = stack(temp, names(temp)[4:end])
    rename!(solar3, Dict(:variable => "Trading_period", :value => "Irradiance"))

    # solar = unstack(solar2, [:Trading_date, :Trading_period], :Region, :Irradiance)
    # dropmissing!(solar)
    # # Divide each element by the max of each column to get the wind scale factors
    # [solar[!, i+2] = round.(solar[!, i+2] / maximum(solar[!, i+2]); digits = 6) for i in 1:2]
    CSV.write(joinpath(@__DIR__, "solardata.csv"),  DataFrame(solar3), header=true)
    return solar2
end

# Raw data processing
@time begin
println("\n\nProcessing ...") # Doesn't always print this line. Def a Julia bug.

filename = ["df_Raw.csv", "df_Gen_Sub.csv", "df_net_"]
demand   = Format_TP_Demand(dir)
price    = Format_TP_Price(dir)
wind     = Format_TP_Wind(dir)
solar    = Format_TP_Solar(dir)
merged   = Merge(demand, price, wind, solar, filename[1])

# Validate the set of investments and existing capacities
if validation
    validate_capacity()
end

# Power set of all investments
portfolio = powerset()
Capacity_Processing(filename, portfolio)
Portfolio_CSV(filename)
println("\nFinished ...\n")
end

max_stars_repo_name: calroughan/EPSILON

max_issues_repo_name: calroughan/EPSILON

max_forks_repo_name: calroughan/EPSILON


text:
#
# map-tests.jl -
#
# Tests for basic mappings.
#
module TestingLazyAlgebraMappings

using Test
using LazyAlgebra
using LazyAlgebra.Foundations
using LinearAlgebra: ⋅, UniformScaling

include("common.jl")
Scaled = LazyAlgebra.Scaled
Sum = LazyAlgebra.Sum
Composition = LazyAlgebra.Composition
Endomorphism = LazyAlgebra.Endomorphism
MorphismType = LazyAlgebra.MorphismType
identical = LazyAlgebra.identical
ALPHAS = (0, 1, -1,  2.71, π)
BETAS = (0, 1, -1, -1.33, Base.MathConstants.φ)

@testset "Rules for uniform scaling             " begin
    # Check + operator.
    @test Id + UniformScaling(1) === 2Id
    @test Id + UniformScaling(2) === 3Id
    @test UniformScaling(1) + Id === 2Id
    @test UniformScaling(2) + Id === 3Id
    # Check - operator.
    @test Id - UniformScaling(1) === 0Id
    @test Id - UniformScaling(2) === -Id
    @test UniformScaling(1) - Id === 0Id
    @test UniformScaling(2) - Id === Id
    # Check * operator.
    @test Id*UniformScaling(1) === Id
    @test Id*UniformScaling(2) === 2Id
    @test UniformScaling(1)*Id === Id
    @test UniformScaling(2)*Id === 2Id
    # Check \circ operator.
    @test Id∘UniformScaling(1) === Id
    @test Id∘UniformScaling(2) === 2Id
    @test UniformScaling(1)∘Id === Id
    @test UniformScaling(2)∘Id === 2Id
    # \cdot is specific.
    @test Id⋅UniformScaling(1) === Id
    @test Id⋅UniformScaling(2) === 2Id
    @test UniformScaling(1)⋅Id === Id
    @test UniformScaling(2)⋅Id === 2Id
    # Check / operator.
    @test Id/UniformScaling(1) === Id
    @test Id/UniformScaling(2) === (1/2)*Id
    @test UniformScaling(1)/Id === Id
    @test UniformScaling(2)/Id === 2Id
    # Check \ operator.
    @test Id\UniformScaling(1) === Id
    @test Id\UniformScaling(2) === 2Id
    @test UniformScaling(1)\Id === Id
    @test UniformScaling(2)\Id === (1/2)*Id
end # testset

@testset "Uniform scaling              ($T)" for T in (Float32, Float64)
    dims = (3,4,5)
    x = randn(T, dims)
    y = randn(T, dims)
    λ = sqrt(2)
    U = λ*Id
    atol, rtol = zero(T), sqrt(eps(T))
    @test U*x  ≈ λ*x       atol=atol rtol=rtol
    @test U'*x ≈ λ*x       atol=atol rtol=rtol
    @test U\x  ≈ (1/λ)*x   atol=atol rtol=rtol
    @test U'\x ≈ (1/λ)*x   atol=atol rtol=rtol
    for α in ALPHAS,
        β in BETAS
        for P in (Direct, Adjoint)
            @test apply!(α, P, U, x, β, vcopy(y)) ≈
                T(α*λ)*x + T(β)*y atol=atol rtol=rtol norm=vnorm2
        end
        for P in (Inverse, InverseAdjoint)
            @test apply!(α, P, U, x, β, vcopy(y)) ≈
                T(α/λ)*x + T(β)*y atol=atol rtol=rtol norm=vnorm2
        end
    end
end # testset

@testset "Rank 1 operators             ($T)" for T in (Float32, Float64)
    dims = (3,4,5)
    n = prod(dims)
    w = rand(T, dims)
    x = rand(T, dims)
    y = rand(T, dims)
    A = RankOneOperator(w, w)
    B = RankOneOperator(w, y)
    C = SymmetricRankOneOperator(w)
    atol, rtol = zero(T), sqrt(eps(T))
    @test LinearType(A) === Linear()
    @test LinearType(C) === Linear()
    @test MorphismType(C) === Endomorphism()
    @test A*Id === A
    @test Id*A === A
    @test A*x  ≈ sum(w.*x)*w atol=atol rtol=rtol norm=vnorm2
    @test A'*x ≈ sum(w.*x)*w atol=atol rtol=rtol norm=vnorm2
    @test B*x  ≈ sum(y.*x)*w atol=atol rtol=rtol norm=vnorm2
    @test B'*x ≈ sum(w.*x)*y atol=atol rtol=rtol norm=vnorm2
    @test C*x  ≈ sum(w.*x)*w atol=atol rtol=rtol norm=vnorm2
    @test C'*x ≈ sum(w.*x)*w atol=atol rtol=rtol norm=vnorm2
    z = vmul!(vcreate(y), B, x)
    @test z == vmul(B, x)
    @test z == B*x
    for α in ALPHAS,
        β in BETAS
        for P in (Direct, Adjoint)
            @test apply!(α, P, C, x, β, vcopy(y)) ≈
                T(α*vdot(w,x))*w + T(β)*y atol=atol rtol=rtol norm=vnorm2
        end
    end
end # testset

@testset "Non-uniform scaling          ($T)" for T in (Float32, Float64)
    dims = (3,4,5)
    n = prod(dims)
    w = randn(T, dims)
    for i in eachindex(w)
        while w[i] == 0
            w[i] = randn(T)
        end
    end
    x = randn(T, dims)
    y = randn(T, dims)
    z = vcreate(y)
    S = NonuniformScaling(w)
    @test diag(S) === w
    @test Diag(w) === S
    atol, rtol = zero(T), sqrt(eps(T))
    @test S*x  ≈ w.*x atol=atol rtol=rtol norm=vnorm2
    @test S'*x ≈ w.*x atol=atol rtol=rtol norm=vnorm2
    @test S\x ≈ x./w atol=atol rtol=rtol norm=vnorm2
    @test S'\x ≈ x./w atol=atol rtol=rtol norm=vnorm2
    for α in ALPHAS,
        β in BETAS
        for P in (Direct, Adjoint)
            @test apply!(α, P, S, x, β, vcopy(y)) ≈
                T(α)*w.*x + T(β)*y atol=atol rtol=rtol norm=vnorm2
        end
        for P in (Inverse, InverseAdjoint)
            @test apply!(α, P, S, x, β, vcopy(y)) ≈
                T(α)*x./w + T(β)*y atol=atol rtol=rtol norm=vnorm2
        end
    end
end # testset

@testset "Non-uniform scaling (Complex{$T})" for T in (Float32, Float64)
    dims = (3,4,5)
    n = prod(dims)
    w = complex.(randn(T, dims), randn(T, dims))
    for i in eachindex(w)
        while w[i] == 0
            w[i] = complex(randn(T), randn(T))
        end
    end
    x = complex.(randn(T, dims), randn(T, dims))
    y = complex.(randn(T, dims), randn(T, dims))
    wx = w.*x
    qx = x./w
    z = vcreate(y)
    S = NonuniformScaling(w)
    atol, rtol = zero(T), sqrt(eps(T))
    @test S*x  ≈ w.*x atol=atol rtol=rtol norm=vnorm2
    @test S'*x ≈ conj.(w).*x atol=atol rtol=rtol norm=vnorm2
    @test S\x ≈ x./w atol=atol rtol=rtol norm=vnorm2
    @test S'\x ≈ x./conj.(w) atol=atol rtol=rtol norm=vnorm2
    for α in ALPHAS,
        β in BETAS
        @test apply!(α, Direct, S, x, β, vcopy(y)) ≈
            T(α)*w.*x + T(β)*y atol=atol rtol=rtol norm=vnorm2
        @test apply!(α, Adjoint, S, x, β, vcopy(y)) ≈
            T(α)*conj.(w).*x + T(β)*y atol=atol rtol=rtol norm=vnorm2
        @test apply!(α, Inverse, S, x, β, vcopy(y)) ≈
            T(α)*x./w + T(β)*y atol=atol rtol=rtol norm=vnorm2
        @test apply!(α, InverseAdjoint, S, x, β, vcopy(y)) ≈
            T(α)*x./conj.(w) + T(β)*y atol=atol rtol=rtol norm=vnorm2
    end
end # testset

@testset "Generalized matrices         ($T)" for T in (Float32, Float64)
    rows, cols = (2,3,4), (5,6)
    nrows, ncols = prod(rows), prod(cols)
    A = randn(T, rows..., cols...)
    x = randn(T, cols)
    y = randn(T, rows)
    G = GeneralMatrix(A)
    atol, rtol = zero(T), sqrt(eps(T))
    mA = reshape(A, nrows, ncols)
    vx = reshape(x, ncols)
    vy = reshape(y, nrows)
    Gx = G*x
    Gty = G'*y
    @test Gx  ≈ reshape(mA*vx,  rows) atol=atol rtol=rtol norm=vnorm2
    @test Gty ≈ reshape(mA'*vy, cols) atol=atol rtol=rtol norm=vnorm2
    test_api(Direct, G, x, y)
    test_api(Adjoint, G, x, y)
end

nothing # testset

end # mdoule


text:
using PyCall
pyimport_conda("wntr", "wntr")


text:
"""
    RSOCtoSOCBridge{T, F, G}

The `RotatedSecondOrderCone` is `SecondOrderCone` representable; see [BN01, p. 104].
Indeed, we have ``2tu = (t/√2 + u/√2)^2 - (t/√2 - u/√2)^2`` hence
```math
2tu \\ge \\lVert x \\rVert_2^2
```
is equivalent to
```math
(t/√2 + u/√2)^2 \\ge \\lVert x \\rVert_2^2 + (t/√2 - u/√2)^2.
```
We can therefore use the transformation ``(t, u, x) \\mapsto (t/√2+u/√2, t/√2-u/√2, x)``.
Note that the linear transformation is a symmetric involution (i.e. it is its own transpose and its own inverse).
That means in particular that the norm of constraint primal and dual values are preserved by the tranformation.

[BN01] Ben-Tal, Aharon, and Nemirovski, Arkadi. *Lectures on modern convex optimization: analysis, algorithms, and engineering applications*. Society for Industrial and Applied Mathematics, 2001.
"""
struct RSOCtoSOCBridge{T,F,G} <:
       SetMapBridge{T,MOI.SecondOrderCone,MOI.RotatedSecondOrderCone,F,G}
    constraint::CI{F,MOI.SecondOrderCone}
end

function rotate_function_type(G::Type{<:MOI.AbstractVectorFunction}, T::Type)
    S = MOIU.promote_operation(/, T, MOIU.scalar_type(G), T)
    Y = MOIU.promote_operation(-, T, S, S)
    Z = MOIU.promote_operation(+, T, S, S)
    return MOIU.promote_operation(vcat, T, Z, Y, G)
end

function concrete_bridge_type(
    ::Type{<:RSOCtoSOCBridge{T}},
    G::Type{<:MOI.AbstractVectorFunction},
    ::Type{MOI.RotatedSecondOrderCone},
) where {T}
    return RSOCtoSOCBridge{T,rotate_function_type(G, T),G}
end

"""
    SOCtoRSOCBridge{T, F, G}

We simply do the inverse transformation of [`RSOCtoSOCBridge`](@ref). In fact, as the
transformation is an involution, we do the same transformation.
"""
struct SOCtoRSOCBridge{T,F,G} <:
       SetMapBridge{T,MOI.RotatedSecondOrderCone,MOI.SecondOrderCone,F,G}
    constraint::CI{F,MOI.RotatedSecondOrderCone}
end

function concrete_bridge_type(
    ::Type{<:SOCtoRSOCBridge{T}},
    G::Type{<:MOI.AbstractVectorFunction},
    ::Type{MOI.SecondOrderCone},
) where {T}
    return SOCtoRSOCBridge{T,rotate_function_type(G, T),G}
end


text:
import Pkg; Pkg.activate("../../")

import CytofRepFAM.Model.ImportanceSampling: exp_f_rfam, min_pairwise_dist_ge_u, rep_fn
using PyPlot
PyPlot.matplotlib.use("Agg")

function get_exp_f(phis, us; N, J, K)
  [[let
      println("phi: $phi, u: $u")
      exp_f_rfam(N=N, J=J, K=K,
                 repulsive_fn=Z->rep_fn(Z, phi),
                 f=Z->min_pairwise_dist_ge_u(Z, u))
    end for u in us] for phi in phis]
end

function plot_results(phis, us, exp_f, savepath)
  for z in zip(phis, exp_f)
    phi, f = z
    plt.plot(us, f, label="ϕ = $phi", marker="o")
  end
  plt.xticks(fontsize=15)
  plt.yticks(fontsize=15)
  plt.xlabel("threshold " * L"d̲", fontsize=18)
  plt.ylabel(L"\Pr\left(\min_{1 ≤ k₁ < k₂ ≤ K}~d(\mathbf{z}_{k_1}, \mathbf{z}_{k_2}) ≥ d̲ \mid \phi \right)",
             fontsize=18)
  plt.axhline(0.95, ls=":", color="grey")
  plt.legend(loc="lower left")
  plt.savefig(savepath, bbox_inches="tight")
  plt.close()
end

# phis = 0:5:50
phis = 0:5:30
us = 0:5

exp_f_data = get_exp_f(phis, us, N=10000, J=15, K=25)
exp_f_sim = get_exp_f(phis, us, N=10000, J=21, K=15)

plt.figure(figsize=(6, 6))
plot_results(phis, us, exp_f_data, "img/data_prior_analysis.pdf")
plt.figure(figsize=(6, 6))
plot_results(phis, us, exp_f_sim, "img/sim_prior_analysis.pdf")