text:
METADATA_MAP = Dict{String, String}("confidence" => "value.Conficence AS Confidence",
                 "type" => "value.Type AS Type",
                 "sequence" => "value.Sequence AS Sequence",
                 "comments" => "value.Comments AS Comments",
                 "previousIds" => "value.PreviousIds AS PreviousIds",
                 "url" => "value.URL AS Url",
                 "chromosome" => "value.Chromosome AS Chromosome",
                 "regionStart" => "value.RegionStart AS RegionStart",
                 "regionEnd" => "value.RegionEnd AS RegionEnd",
                 "strand" => "value.Strand AS Strand",
                 "communityAnnotation" => "value.CommunityAnnotation AS CommunityAnnotation",
                 "nExperiments" => "value.nExperiments AS nExperiments",
                 "reads" => "value.Reads AS Reads",
                 "evidenceType" => "value.EvidenceType AS EvidenceType");

"""
    convertToCurrentMirbaseVersion(queryId::String; species::String, metadata::Union{String, Vector{String}})
    convertToCurrentMirbaseVersion(queryId::Vector{String}; species::String, metadata::Union{String, Vector{String}})

Translate input IDs to the name in the current miRBase release version.

As input all kinds of miRNA identifiers are accepted, such as mature miRNA names, MI-accession, MIMAT-accession
and precursor-miRNA name. Moreover, metadata from the latest miRBase release version
can be returned with the name, such as the sequence, genomic location and annotations (
see `getValidMirnaMetadataValues` for further information on metadata values).

### Arguments
- `queryId` : either of type `String` or `Vector{String}` containing any miRNA identifier
- `species::String` : Three letter species value, such as `"hsa"` or `"mmu"`
- `metadata` : `String` or `Vector{String}` with metadata keywords (see `getValidMirnaMetadataValues()`)

**See also**
`getValidMirnaMetadataValues()`

### Examples
```julia-repl
> convertToCurrentMirbaseVersion("hsa-mir-29a")
1×4 DataFrames.DataFrame
│ Row │ InputId     │ Accession │ CurrentMirna │ CurrentVersion │
├─────┼─────────────┼───────────┼──────────────┼────────────────┤
│ 1   │ hsa-mir-29a │ MI0000087 │ hsa-mir-29a  │ 22.0           │

> convertToCurrentMirbaseVersion(["hsa-mir-29a", "hsa-mir-134a", "hsa-miR-29a"], metadata = ["sequence", "chromosome"])
10-Apr 16:23:41:WARNING:root:Some input identifiers were not found in the database.
Union{Missings.Missing, String}["hsa-mir-134a"]2×6 DataFrames.DataFrame
│ Row │ InputId     │ Accession    │ CurrentMirna   │ CurrentVersion │ Sequence                                                         │ Chromosome │
├─────┼─────────────┼──────────────┼────────────────┼────────────────┼──────────────────────────────────────────────────────────────────┼────────────┤
│ 1   │ hsa-miR-29a │ MIMAT0000086 │ hsa-miR-29a-3p │ 22.0           │ UAGCACCAUCUGAAAUCGGUUA                                           │ missing    │
│ 2   │ hsa-mir-29a │ MI0000087    │ hsa-mir-29a    │ 22.0           │ AUGACUGAUUUCUUUUGGUGUUCAGAGUCAAUAUAAUUUUCUAGCACCAUCUGAAAUCGGUUAU │ 7          │
```
"""
function convertToCurrentMirbaseVersion(
   queryId::S;
   species::S = "",
   metadata::U = [""]) where {S<:String, U<:Union{String, Vector{String}}}

   return convertToCurrentMirbaseVersion(
         [queryId],
         species = species,
         metadata = isa(metadata, String) ? [metadata] : metadata);

end

function convertToCurrentMirbaseVersion(
   queryId::T;
   species::S = "",
   metadata::U = [""]) where {S<:String, T<:Vector{String}, U<:Union{String, Vector{String}}}

   metadata = isa(metadata, String) ? String[metadata] : metadata;

   q = "CALL rcsi.mirna.convertToCurrentMirbaseVersion(" *
      "{queryId}, {species}" *
      (metadata[1] == "" ? "" : (", {metadata}")) *
      ") " *
      "YIELD value " *
      "RETURN DISTINCT value.InputId AS InputId, " *
      "value.Accession AS Accession, " *
      "value.CurrentMirna AS CurrentMirna, " *
      "value.CurrentVersion AS CurrentVersion";

   if (all(metadata .!= ""))
      q *= ", " *
         join(
            map(i->haskey(METADATA_MAP, i) ?
               METADATA_MAP[i] :
               error("provided metadata value '$i' not supported"), metadata), ", ");
   end

   x = Neo4j.cypherQuery(
         NEO4J_CONNECTION,
         q,
         "queryId" => queryId,
         "species" => species,
         "metadata" => metadata);

   # Check for multi matches and not found input Ids
   multiMatches = @> begin
      x[[:InputId, :Accession]]
      unique()
      @by(:InputId, nRows = length(:Accession))
      @where(:nRows .> 1)
   end

   if (DataFrames.nrow(multiMatches) != 0)
      warn("Some input identifiers match to more than one MIMAT accession!");
      show(multiMatches);
   end
   postCheckMirnas(x, :Accession, queryId);

   return x;

end

"""
    convertMatureMirnasToVersions(queryId::String; targetVersion::Union{Union{Float64, Int64}, Union{Vector{Float64}, Vector{Int64}}}, species::String, sequence::Boolean)
    convertMatureMirnasToVersions(queryId::Vector{String}; targetVersion::Union{Union{Float64, Int64}, Union{Vector{Float64}, Vector{Int64}}}, species::String, sequence::Boolean)

Translate mature miRNA names to different versions

Taking single or multiple `queryId`s return the name(s) for a single or
#' multiple `targetVersion`(s). Optionally the `species` can be provided and if
#' `sequence` == TRUE the sequence for each version respecively returned.

### Arguments
- `queryId` : either of type `String` or `Vector{String}` containing any miRNA identifier
- `targetVersion` : The target miRBase version(s), such as `[20, 21]` (default: highest supported version).
- `species::String` : Three letter species value, such as `"hsa"` or `"mmu"` (default: "")
- `sequence::Bool` : A boolean depicting if sequence for respective miRBase versions is returned or not (default: `false`).

**See also**
`convertToCurrentMirbaseVersion()` if only current version is needed.

### Examples
```julia-repl
> convertMatureMirnasToVersions("hsa-mir-29a")
1×4 DataFrames.DataFrame
│ Row │ InputId     │ Accession │ CurrentMirna │ CurrentVersion │
├─────┼─────────────┼───────────┼──────────────┼────────────────┤
│ 1   │ hsa-mir-29a │ MI0000087 │ hsa-mir-29a  │ 22.0           │

> convertMatureMirnasToVersions(["hsa-mir-29a", "hsa-mir-134a", "hsa-miR-29a"], metadata = ["sequence", "chromosome"])
10-Apr 16:23:41:WARNING:root:Some input identifiers were not found in the database.
Union{Missings.Missing, String}["hsa-mir-134a"]2×6 DataFrames.DataFrame
│ Row │ InputId     │ Accession    │ CurrentMirna   │ CurrentVersion │ Sequence                                                         │ Chromosome │
├─────┼─────────────┼──────────────┼────────────────┼────────────────┼──────────────────────────────────────────────────────────────────┼────────────┤
│ 1   │ hsa-miR-29a │ MIMAT0000086 │ hsa-miR-29a-3p │ 22.0           │ UAGCACCAUCUGAAAUCGGUUA                                           │ missing    │
│ 2   │ hsa-mir-29a │ MI0000087    │ hsa-mir-29a    │ 22.0           │ AUGACUGAUUUCUUUUGGUGUUCAGAGUCAAUAUAAUUUUCUAGCACCAUCUGAAAUCGGUUAU │ 7          │
```
"""
function convertMatureMirnasToVersions(
   queryId::S;
   targetVersion::F = [0.0],
   species::S = "",
   sequence::B = false) where {S<:String, F<:Union{Union{Float64, Int64}, Union{Vector{Float64}, Vector{Int64}}}, B<:Bool}

   # targetVersionF = normaliseTargetVersionParameter(targetVersion);

   return convertMatureMirnasToVersions(
         [queryId],
         targetVersion = targetVersion,
         species = species,
         sequence = sequence);

end

function convertMatureMirnasToVersions(
   queryId::T;
   targetVersion::F = [0.0],
   species::S = "",
   sequence::B = false) where {S<:String, T<:Vector{String}, F<:Union{Union{Float64, Int64}, Union{Vector{Float64}, Vector{Int64}}}, B<:Bool}

   targetVersionF = normaliseTargetVersionParameter(targetVersion);

   q = "CALL rcsi.mirna.convertMatureMirnas(" *
         "{queryId}, {targetVersion}, {species}, {sequence}" *
         ") " *
         "YIELD value " *
         "RETURN DISTINCT value.InputId AS InputId, " *
            "value.MatureAccession AS MatureAccession, " *
            "value.miRBaseVersion AS miRBaseVersion, " *
            "value.TargetMirna AS TargetMirna";

   if (sequence)
      q *= ", value.TargetSequence AS TargetSequence";
   end

   x = Neo4j.cypherQuery(
         NEO4J_CONNECTION,
         q,
         "queryId" => queryId,
         "targetVersion" => targetVersionF,
         "species" => species,
         "sequence" => sequence);

   # Check for multi matches and not found input Ids
   if DataFrames.nrow(x) > 0
      multiMatches = @> begin
         x[[:InputId, :MatureAccession]]
         unique()
         @by(:InputId, nRows = length(:MatureAccession))
         @where(:nRows .> 1)
      end
   else
      return x;
   end

   if (DataFrames.nrow(multiMatches) != 0)
      warn("Some input identifiers match to more than one MIMAT accession!");
      show(multiMatches);
   end
   postCheckMirnas(x, :MatureAccession, queryId);

   return x;

end

function normaliseTargetVersionParameter(targetVersion::Union{Union{Float64, Int64}, Union{Vector{Float64}, Vector{Int64}}})
   targetVersionF = [0.0];
   if isa(targetVersion, Int64) || isa(targetVersion, Float64)
      targetVersionF = Float64[targetVersion];
   elseif isa(targetVersion, Vector{Int64})
      targetVersionF = collect(Float64,targetVersion);
   elseif isa(targetVersion, Vector{Float64})
      targetVersionF = targetVersion;
   end
   return targetVersionF;
end

max_stars_repo_name: StefanHaunsberger/GeneNameGenieJ

max_issues_repo_name: StefanHaunsberger/GeneNameGenieJ

max_forks_repo_name: StefanHaunsberger/GeneNameGenieJ


text:
using TimeWarp
 #using TimeWarp.WarpPlots
using TimeWarp.Datasets
using PyPlot # use pyplot
plt = PyPlot

# UCI data repository must be downloaded first, run TimeWarp.Datasets.download_ucr()
data, labels = ucr_traindata("MedicalImages");

# c = [ cl==1 ? :blue : :red for cl in class ]'
# plot(y, linecolor=c, legend=false)
nclust = 3
init_centers = TimeWarp.dbaclust_initial_centers(data, nclust, FastDTW(10))
centers, clustids, result = dbaclust(data, nclust, FastDTW(10); init_centers=deepcopy(init_centers), iterations=1)

figure()
plot(data, color="0.75")
plot(centers, color="red")
plot(init_centers, color="green")

if is_linux()
  plt.show()
end

ext: jl


text:
# editing

function edit(file::String, line::Integer)
    if OS_NAME == :Windows || OS_NAME == :Darwin
        default_editor = "open"
    elseif isreadable("/etc/alternatives/editor")
        default_editor = "/etc/alternatives/editor"
    else
        default_editor = "emacs"
    end
    editor = get(ENV,"JULIA_EDITOR", get(ENV,"VISUAL", get(ENV,"EDITOR", default_editor)))
    if ispath(editor)
        if isreadable(editor)
            edpath = realpath(editor)
            edname = basename(edpath)
        else
            error("can't find \"$editor\"")
        end
    else
        edpath = edname = editor
    end
    issrc = length(file)>2 && file[end-2:end] == ".jl"
    if issrc
        file = find_source_file(file)
    end
    if beginswith(edname, "emacs")
        jmode = joinpath(JULIA_HOME, "..", "..", "contrib", "julia-mode.el")
        if issrc && isreadable(jmode)
            run(`$edpath $file --eval "(progn
                                     (require 'julia-mode \"$jmode\")
                                     (julia-mode)
                                     (goto-line $line))"`)
        else
            run(`$edpath $file --eval "(goto-line $line)"`)
        end
    elseif edname == "vim"
        run(`$edpath $file +$line`)
    elseif edname == "textmate" || edname == "mate"
        spawn(`$edpath $file -l $line`)
    elseif beginswith(edname, "subl")
        spawn(`$edpath $file:$line`)
    elseif OS_NAME == :Windows && (edname == "start" || edname == "open")
        spawn(`start /b $file`)
    elseif OS_NAME == :Darwin && (edname == "start" || edname == "open")
        spawn(`open -t $file`)
    elseif edname == "kate"
        spawn(`$edpath $file -l $line`)
    elseif edname == "nano"
        run(`$edpath +$line $file`)
    else
        run(`$(shell_split(edpath)) $file`)
    end
    nothing
end
edit(file::String) = edit(file, 1)

function less(file::String, line::Integer)
    pager = get(ENV, "PAGER", "less")
    run(`$pager +$(line)g $file`)
end
less(file::String) = less(file, 1)

edit(f::Union(Function,DataType))    = edit(functionloc(f)...)
edit(f::Union(Function,DataType), t) = edit(functionloc(f,t)...)
less(f::Union(Function,DataType))    = less(functionloc(f)...)
less(f::Union(Function,DataType), t) = less(functionloc(f,t)...)

function edit( m::Method )
    tv, decls, file, line = arg_decl_parts(m)
    edit( string(file), line )
end

# clipboard copy and paste

@osx_only begin
    function clipboard(x)
        w,p = writesto(`pbcopy`)
        print(w,x)
        close(w)
        wait(p)
    end
    clipboard() = readall(`pbpaste`)
end

@linux_only begin
    _clipboardcmd = nothing
    function clipboardcmd()
        global _clipboardcmd
        _clipboardcmd !== nothing && return _clipboardcmd
        for cmd in (:xclip, :xsel)
            success(`which $cmd` |> DevNull) && return _clipboardcmd = cmd
        end
        error("no clipboard command found, please install xsel or xclip")
    end
    function clipboard(x)
        c = clipboardcmd()
        cmd = c == :xsel  ? `xsel --nodetach --input --clipboard` :
              c == :xclip ? `xclip -quiet -in -selection clipboard` :
            error("unexpected clipboard command: $c")
        w,p = writesto(cmd)
        print(w,x)
        close(w)
        wait(p)
    end
    function clipboard()
        c = clipboardcmd()
        cmd = c == :xsel  ? `xsel --nodetach --output --clipboard` :
              c == :xclip ? `xclip -quiet -out -selection clipboard` :
            error("unexpected clipboard command: $c")
        readall(cmd)
    end
end

@windows_only begin
    function clipboard(x::String)
        ccall((:OpenClipboard, "user32"), stdcall, Bool, (Ptr{Void},), C_NULL)
        ccall((:EmptyClipboard, "user32"), stdcall, Bool, ())
        p = ccall((:GlobalAlloc, "kernel32"), stdcall, Ptr{Void}, (Uint16,Int32), 2, length(x)+1)
        p = ccall((:GlobalLock, "kernel32"), stdcall, Ptr{Void}, (Ptr{Void},), p)
        # write data to locked, allocated space
        ccall(:memcpy, Ptr{Void}, (Ptr{Void},Ptr{Uint8},Int32), p, x, length(x)+1)
        ccall((:GlobalUnlock, "kernel32"), stdcall, Void, (Ptr{Void},), p)
        # set clipboard data type to 13 for Unicode text/string
        p = ccall((:SetClipboardData, "user32"), stdcall, Ptr{Void}, (Uint32, Ptr{Void}), 1, p)
        ccall((:CloseClipboard, "user32"), stdcall, Void, ())
    end
    clipboard(x) = clipboard(sprint(io->print(io,x))::ByteString)

    function clipboard()
        ccall((:OpenClipboard, "user32"), stdcall, Bool, (Ptr{Void},), C_NULL)
        s = bytestring(ccall((:GetClipboardData, "user32"), stdcall, Ptr{Uint8}, (Uint32,), 1))
        ccall((:CloseClipboard, "user32"), stdcall, Void, ())
        return s
    end
end

if !isdefined(:clipboard)
    clipboard(x="") = error("clipboard functionality not implemented for $OS_NAME")
end

# system information

function versioninfo(io::IO=STDOUT, verbose::Bool=false)
    println(io,             "Julia Version $VERSION")
    if !isempty(GIT_VERSION_INFO.commit_short)
      println(io,             "Commit $(GIT_VERSION_INFO.commit_short) ($(GIT_VERSION_INFO.date_string))")
    end
    if ccall(:jl_is_debugbuild, Cint, ())!=0
        println(io, "DEBUG build")
    end
    println(io,             "Platform Info:")
    println(io,             "  System: ", Sys.OS_NAME, " (", Sys.MACHINE, ")")

    cpu = Sys.cpu_info()
    println(io,         "  CPU: ", cpu[1].model)
    println(io,             "  WORD_SIZE: ", Sys.WORD_SIZE)
    if verbose
        lsb = ""
        @linux_only try lsb = readchomp(`lsb_release -ds` .> DevNull) end
        @windows_only try lsb = strip(readall(`$(ENV["COMSPEC"]) /c ver`)) end
        if lsb != ""
            println(io,     "           ", lsb)
        end
        println(io,         "  uname: ",readchomp(`uname -mprsv`))
        println(io,         "Memory: $(Sys.total_memory()/2^30) GB ($(Sys.free_memory()/2^20) MB free)")
        try println(io,     "Uptime: $(Sys.uptime()) sec") catch end
        print(io,           "Load Avg: ")
        print_matrix(io,    Sys.loadavg()')
        println(io          )
        Sys.cpu_summary(io)
        println(io          )
    end
    if Base.libblas_name == "libopenblas" || blas_vendor() == :openblas
        openblas_config = openblas_get_config()
        println(io,         "  BLAS: libopenblas (", openblas_config, ")")
    else
        println(io,         "  BLAS: ",libblas_name)
    end
    println(io,             "  LAPACK: ",liblapack_name)
    println(io,             "  LIBM: ",libm_name)
    if verbose
        println(io,         "Environment:")
        for (k,v) in ENV
            if !is(match(r"JULIA|PATH|FLAG|^TERM$|HOME",bytestring(k)), nothing)
                println(io, "  $(k) = $(v)")
            end
        end
        println(io          )
        println(io,         "Package Directory: ", Pkg.dir())
        Pkg.status(io)
    end
end
versioninfo(verbose::Bool) = versioninfo(STDOUT,verbose)

# searching definitions

function which(f::Callable, args...)
    if !isgeneric(f)
        throw(ErrorException("not a generic function, no methods available"))
    end
    ms = methods(f, map(a->(isa(a,Type) ? Type{a} : typeof(a)), args))
    isempty(ms) && throw(MethodError(f, args))
    ms[1]
end

macro which(ex0)
    if isa(ex0,Expr) &&
        any(a->(Meta.isexpr(a,:kw) || Meta.isexpr(a,:parameters)), ex0.args)
        # keyword args not used in dispatch, so just remove them
        args = filter(a->!(Meta.isexpr(a,:kw) || Meta.isexpr(a,:parameters)), ex0.args)
        return Expr(:call, :which, map(esc, args)...)
    end
    if isa(ex0, Expr) && ex0.head == :call
        return Expr(:call, :which, map(esc, ex0.args)...)
    end
    ex = expand(ex0)
    exret = Expr(:call, :error, "expression is not a function call")
    if !isa(ex, Expr)
        # do nothing -> error
    elseif ex.head == :call
        if any(e->(isa(e,Expr) && e.head==:(...)), ex0.args) &&
            isa(ex.args[1],TopNode) && ex.args[1].name == :apply
            exret = Expr(:call, ex.args[1], :which,
                         Expr(:tuple, esc(ex.args[2])),
                         map(esc, ex.args[3:end])...)
        else
            exret = Expr(:call, :which, map(esc, ex.args)...)
        end
    elseif ex.head == :body
        a1 = ex.args[1]
        if isa(a1, Expr) && a1.head == :call
            a11 = a1.args[1]
            if a11 == :setindex!
                exret = Expr(:call, :which, a11, map(esc, a1.args[2:end])...)
            end
        end
    elseif ex.head == :thunk
        exret = Expr(:call, :error, "expression is not a function call, or is too complex for @which to analyze; "
                                  * "break it down to simpler parts if possible")
    end
    exret
end

# `methodswith` -- shows a list of methods using the type given

function methodswith(t::Type, m::Module, showparents::Bool=false)
    meths = Method[]
    for nm in names(m)
        try
           mt = eval(m, nm)
           d = mt.env.defs
           while !is(d,())
               if any(map(x -> x == t || (showparents && t <: x && x != Any && x != ANY && !isa(x, TypeVar)), d.sig))
                   push!(meths, d)
               end
               d = d.next
           end
        end
    end
    return meths
end

function methodswith(t::Type, showparents::Bool=false)
    meths = Method[]
    mainmod = current_module()
    # find modules in Main
    for nm in names(mainmod)
        if isdefined(mainmod,nm)
            mod = eval(mainmod, nm)
            if isa(mod, Module)
                meths = [meths, methodswith(t, mod, showparents)]
            end
        end
    end
    return meths
end

## file downloading ##

downloadcmd = nothing
@unix_only function download(url::String, filename::String)
    global downloadcmd
    if downloadcmd === nothing
        for checkcmd in (:curl, :wget, :fetch)
            if success(`which $checkcmd` |> DevNull)
                downloadcmd = checkcmd
                break
            end
        end
    end
    if downloadcmd == :wget
        run(`wget -O $filename $url`)
    elseif downloadcmd == :curl
        run(`curl -o $filename -L $url`)
    elseif downloadcmd == :fetch
        run(`fetch -f $filename $url`)
    else
        error("no download agent available; install curl, wget, or fetch")
    end
    filename
end

@windows_only function download(url::String, filename::String)
    b_dest = bytestring(filename)
    b_url = bytestring(url)
    res = ccall((:URLDownloadToFileA,:urlmon),stdcall,Cuint,
                (Ptr{Void},Ptr{Uint8},Ptr{Uint8},Cint,Ptr{Void}),0,b_url,b_dest,0,0)
    if res != 0
        error("automatic download failed (error: $res): $url")
    end
end

function download(url::String)
    filename = tempname()
    download(url, filename)
end

max_stars_repo_name: alyst/julia

max_issues_repo_name: alyst/julia

max_forks_repo_name: alyst/julia


text:
using Gurobi
using CSV
using DataFrames
using JuMP
using ProgressMeter
using ElectronDisplay

include("model/src/Sesam.jl")
using .Sesam

# Power plants
pp_df = CSV.read("data/test/power_plants.csv")
avail_con_df = CSV.read("data/test/avail_con.csv")
prices_df = CSV.read("data/test/prices.csv")

avail_pv = CSV.read("data/test/avail_pv.csv")
avail_windon = CSV.read("data/test/avail_windon.csv")
avail_windoff = CSV.read("data/test/avail_windoff.csv")

avail = Dict(:PV => avail_pv,
	:WindOnshore => avail_windon,
	:WindOffshore => avail_windoff,
	:global => avail_con_df)

pp = PowerPlants(pp_df, avail = avail_con_df, prices = prices_df)

res_df = CSV.read("data/test/res.csv")
res = RenewableEnergySource(res_df, avail)

# Load
load_df = CSV.read("data/test/load.csv") |>
    df -> sum.(eachrow(df))
load = load_df

T = collect(1:(24*7*10))
P = pp.id
R = res.id

# JuMP model
economic_dispatch = JuMP.Model(with_optimizer(Gurobi.Optimizer))

@variable(economic_dispatch, G[T, P] >= 0)
@variable(economic_dispatch, G_R[T, R] >= 0)

@constraint(economic_dispatch, MarketClearing[t = T],
    sum(G[t, p] for p in P) + sum(G_R[t, r] for r in R) == load[t]);

@constraint(economic_dispatch, MaxGenerationDisp[p = P, t = T],
    G[t, p] <= pp.capacity[p][t]);

@constraint(economic_dispatch, MaxGenerationDispR[r = R, t = T],
    G_R[t, r] <= res.infeed[r][t]);

@objective(economic_dispatch, Min,
    sum(pp.mc[p][t] * G[t, p] for p in P, t in T));

###
JuMP.optimize!(economic_dispatch)
JuMP.termination_status(economic_dispatch)
JuMP.primal_status(economic_dispatch)
JuMP.dual_status(economic_dispatch)

JuMP.objective_value(economic_dispatch)
JuMP.value.(G)
JuMP.value.(G_R)
value(G[1,"BNA0005"])

G[1,"LN"]

df = DataFrame(T=[],ID=[],P=[])
for t in T
	for p in P
		val = JuMP.value(G[t, p])
		push!(df, [t p val])
	end
end

CSV.write("output/dispatch.csv",df)

electrondisplay(df)

# Dual value of market clearing
has_duals(economic_dispatch)
JuMP.dual.(MarketClearing)
JuMP.shadow_price.(MarketClearing)


text:
# Loader of svmlight / liblinear format files

# Copyright (c) 2015 Ishita Takeshi
# License is MIT

module SVMLightLoader

export
    SVMLightFile,
    load_svmlight_file,
    InvalidFormatError,
    NoDataException

include("exception.jl")
include("loader.jl")
include("iterable.jl")

end # module


text:
# compute running crc for microchip RTDM
# pass crc = 0xffff to start
const polynomial = 0xa001
function rtdm_crc(crc::UInt16, data::UInt8)
  for i in 1:8
    if ((crc $ data)&0x0001) == 0x0001
      crc = (crc>>>1) $ polynomial
    else
      crc = crc>>>1
    end
    data >>>=1
  end
  return crc
end

pipecrc = PipeBuffer()
const buffercrc = Array(UInt8,1)
function rtdm_crc(crc::UInt16, data)
  write(pipecrc,data)
  while ~eof(pipecrc)
    buffercrc[1] = read(pipecrc, UInt8)
    crc = rtdm_crc(crc, buffercrc[1])
  end
  return crc
end

function rtdm_crc(crc::UInt16, data::Array{Any})
  for element in data
    crc = rtdm_crc(crc, element)
  end
  return crc
end

rtdm_crc(data) = rtdm_crc(0xffff,data)  # or call with one parameter to start


text:
using Base.Test

using OPCSPs

positions = Vector{Float64}[[0,0], [1,0], [1,1], [0,1], [0.5,0.5]]
r = [0.0, 10.0, 100.0, 1.0, 20.0]

op = SimpleOP(r, positions, 3.0, 1, 4)
path = solve_op(op)
@test path==[1,5,3,4]

op = SimpleOP(r, positions, 2.7, 2, 1)
path = solve_op(op)
@test path==[2,3,5,1]


text:
module DiffOperatorsTests

using Test
using InplaceArrays.TensorValues
using InplaceArrays.Fields
using InplaceArrays.Fields: MockField
using LinearAlgebra
using FillArrays

np = 4
p = Point(1,2)
x = fill(p,np)

v = 3.0
d = 2
_f = MockField{d}(v)

l = 10
_af = Fill(_f,l)

for f in (_f,_af)

  @test ∇(f) == gradient(f)

  @test divergence(f) == tr(gradient(f))

  @test curl(f) == grad2curl(gradient(f))

  @test ∇*f == divergence(f)

  @test cross(∇,f) == curl(f)

  @test outer(∇,f) == ∇(f)

  @test outer(f,∇) == transpose(∇(f))

  @test ε(f) == symmetic_part(gradient(f))

  @test Δ(f) == ∇*∇(f)

end

end # module


text:
export getAnalyticalConstGrad2D,getAnalyticalConstGrad3D,getAnalyticalConstGradInv2D,getAnalyticalConstGradInv3D,getSmoothGaussianMedium,getSmoothFactoredModel,getSmoothFactoredModel3D

function getAnalyticalConstGrad2D(n::Array{Int64,1},h::Array{Float64,1})
src = [1,div(n[2],2)];

source1 = (src[1]-1)*h[1];
source2 = (src[2]-1)*h[2];

(X1,X2) = ndgrid((0:(n[1]-1))*h[1],(0:(n[2]-1))*h[2]);

a  = -0.4;
s0 = 2.0;

kappaSquared = s0^2 .+ 2.0*a*(X1.-source1);

# matshow(kappaSquared); colorbar();
SBarSquared = s0^2 .+ a*(X1.-source1);
rSquared = ((X1.-source1).^2 .+ (X2.-source2).^2);
sigmaSquarred = 2*rSquared./(SBarSquared + sqrt.(SBarSquared.^2 .- (a^2).*rSquared));
sigma = sqrt.(sigmaSquarred);
T_exact =  SBarSquared.*sigma .- (a^2).*(sigma.^3)./6;

return kappaSquared,src,T_exact;
end

function getAnalyticalConstGrad3D(n::Array{Int64,1},h::Array{Float64,1})
src = [div(n[1],2),div(n[2],2),1];

source1 = (src[1]-1)*h[1];
source2 = (src[2]-1)*h[2];
source3 = (src[3]-1)*h[3];

(X1,X2,X3) = ndgrid((0:(n[1]-1))*h[1],(0:(n[2]-1))*h[2],(0:(n[3]-1))*h[3]);

# a  = -1.75;
a  = -1.65;
s0 = 2.0;

kappaSquared = s0^2 .+ 2*a*(X3.-source3);
SBarSquared = s0^2 .+ a*(X3.-source3);
rSquared = (X1.-source1).^2 .+ (X2.-source2).^2 .+ (X3.-source3).^2;
sigmaSquarred = 2*rSquared./(SBarSquared .+ sqrt.(SBarSquared.^2 .- (a^2).*rSquared));
sigma = sqrt.(sigmaSquarred);
T_exact =  SBarSquared.*sigma .- (a^2).*(sigma.^3)./6;

return kappaSquared,src,T_exact;
end

function getAnalyticalConstGradInv2D(n::Array{Int64,1},h::Array{Float64,1})

src = [1,div(n[2],2)];

source1 = (src[1]-1)*h[1];
source2 = (src[2]-1)*h[2];

(X1,X2) = ndgrid((0:(n[1]-1))*h[1],(0:(n[2]-1))*h[2]);

a  = 1.0;
ainv = (1.0./a);
s0 = 2.0;

Kappa = 1.0./(1.0./s0 .+ a*(X1.-source1));
T_exact = ainv*acosh.(1.0.+(0.5*s0*a*a).*Kappa.*((X1.-source1).^2 .+ (X2.-source2).^2));
G2_exact = ainv*(1.0./sqrt.((1.0.+(0.5*s0*a*a).*Kappa.*((X1.-source1).^2 .+ (X2 .- source2).^2)).^2 .- 1)).*(0.5*s0*a*a).*Kappa.*2.0.*(X2.-source2);
G1_exact = ainv*(1.0./sqrt.((1.0.+(0.5*s0*a*a).*Kappa.*((X1.-source1).^2 .+ (X2 .- source2).^2)).^2 .- 1)).*((0.5*s0*a*a).*Kappa.*2.0.*(X1.-source1)
			+ (0.5*s0*a*a).*((X1.-source1).^2 + (X2.-source2).^2).*(-(Kappa.^2)).*a);
kappaSquared = Kappa.^2;

return kappaSquared,src,T_exact;

end

function getAnalyticalConstGradInv3D(n::Array{Int64,1},h::Array{Float64,1})

src = [div(n[1],2),div(n[1],2),1];

source1 = (src[1]-1)*h[1];
source2 = (src[2]-1)*h[2];
source3 = (src[3]-1)*h[3];

(X1,X2,X3) = ndgrid((0:(n[1]-1))*h[1],(0:(n[2]-1))*h[2],(0:(n[3]-1))*h[3]);

a  = 1.0;
ainv = (1.0./a);
s0 = 2.0;
Kappa = 1.0./(1.0./s0 .+ a*(X3.-source3));
Rsquared = (X1.-source1).^2 .+ (X2.-source2).^2 .+ (X3.-source3).^2;
T_exact = ainv*acosh.(1.0.+(0.5*s0*a*a).*Kappa.*(Rsquared));
G1_exact = ainv*(1.0./sqrt.((1.0.+(0.5*s0*a*a).*Kappa.*Rsquared).^2 .- 1.0)).*(0.5*s0*a*a).*Kappa.*2.0.*(X1.-source1);
G2_exact = ainv*(1.0./sqrt.((1.0.+(0.5*s0*a*a).*Kappa.*Rsquared).^2 .- 1.0)).*(0.5*s0*a*a).*Kappa.*2.0.*(X2.-source2);
G3_exact = ainv*(1.0./sqrt.((1.0.+(0.5*s0*a*a).*Kappa.*Rsquared).^2 .- 1.0)).*((0.5*s0*a*a).*Kappa.*2.0.*(X3.-source3)
			.+ (0.5*s0*a*a).*(Rsquared).*(-(Kappa.^2)).*a);
kappaSquared = Kappa.^2;

return kappaSquared,src,T_exact;

end

function getSmoothGaussianMedium(n::Array{Int64,1},h::Array{Float64,1})
src = div.(n,4);
kappaSquared = [];
T_exact = [];
if length(n)==2
	(T1_exact,G11_exact,G12_exact) = getSmoothFactoredModel(n,h);
	(T0,G01,G02,L0) = getAnalytic2DeikonalSolutionAll(n,h,src);
	G1_exact = T0.*G11_exact + G01.*T1_exact;
	G2_exact = T0.*G12_exact + G02.*T1_exact;
	kappaSquared = G1_exact.*G1_exact + G2_exact.*G2_exact;
	# matshow(sqrt(kappaSquared)); colorbar();
	# xlabel("y");
	# ylabel("x");
	# xticks(0:div(n[2],8):n[2],0:1:8);
	# yticks(0:div(n[1],4):n[1],0:1:4);

	#figure()
	#CS = contour(kappaSquared,50);colorbar();title("kappaSquared");
	#clabel(CS, inline=1, fontsize=10)
	#figure()
	#contour(G1_exact);colorbar();title("G1_exact");
	T_exact = T0.*T1_exact;
	#figure()
	#contour(T_exact,100);colorbar();title("T exact");

else

	(T1_exact,G11_exact,G12_exact,G13_exact) = getSmoothFactoredModel3D(n,h);

	(T0,G01,G02,G03) = getAnalytic3DeikonalSolutionAll(n,h,src);

	G1_exact = T0.*G11_exact + G01.*T1_exact;
	G2_exact = T0.*G12_exact + G02.*T1_exact;
	G3_exact = T0.*G13_exact + G03.*T1_exact;
	kappaSquared = G1_exact.*G1_exact + G2_exact.*G2_exact + G3_exact.*G3_exact;

	# Er = kappaSquared;
	# for k=1:n[3]
		# matshow(reshape(Er[:,:,k],n[1],n[2]));colorbar();title("kappaSQ");
	# end

	T_exact = T0.*T1_exact;

end

return kappaSquared,src,T_exact;
end

function getSmoothFactoredModel(n::Array{Int64,1},h::Array{Float64,1})
xsrc = zeros(2);
src_kappa = zeros(Int64,2);
src_kappa[1] = div(n[1],3);
src_kappa[2] = div(n[2],4);

xsrc[1] = (src_kappa[1]-1)*h[1];
xsrc[2] = (src_kappa[2]-1)*h[2];

X1,X2 = ndgrid((0:(n[1]-1))*h[1],(0:(n[2]-1))*h[2]);

sigma = 0.1;
T1_exact = (exp.( - (sigma*((X1 .- xsrc[1]).^2) + 4*sigma*((X2.-xsrc[2]).^2))) .+ 1)/2;
G11_exact = -2*sigma*(X1 .- xsrc[1]).*exp.( - (sigma*((X1 .- xsrc[1]).^2) + 4*sigma*((X2.-xsrc[2]).^2)))/2;
G12_exact = -8*sigma*(X2 .- xsrc[2]).*exp.( - (sigma*((X1 .- xsrc[1]).^2) + 4*sigma*((X2.-xsrc[2]).^2)))/2;

return T1_exact,G11_exact,G12_exact;
end

function getSmoothFactoredModel3D(n::Array{Int64,1},h::Array{Float64,1})
xsrc = zeros(3);
src_kappa = zeros(Int64,3);
src_kappa[1] = div(n[1],3);
src_kappa[2] = div(n[2],4);
src_kappa[3] = div(n[3],2);

xsrc[1] = (src_kappa[1]-1)*h[1];
xsrc[2] = (src_kappa[2]-1)*h[2];
xsrc[3] = (src_kappa[3]-1)*h[3];

X1,X2,X3 = ndgrid((0:(n[1]-1))*h[1],(0:(n[2]-1))*h[2],(0:(n[3]-1))*h[3]);

Sigma = [0.1,0.4,0.2];
EXPRSkewedSquared = exp.(-Sigma[1]*((X1 .- xsrc[1]).^2) .- Sigma[2]*((X2.-xsrc[2]).^2) .- Sigma[3]*((X3.-xsrc[3]).^2));

T1_exact = EXPRSkewedSquared/2.0 .+ 0.5;
G11_exact = -Sigma[1]*(X1 .- xsrc[1]).*EXPRSkewedSquared;
G12_exact = -Sigma[2]*(X2 .- xsrc[2]).*EXPRSkewedSquared;
G13_exact = -Sigma[3]*(X3 .- xsrc[3]).*EXPRSkewedSquared;

return T1_exact,G11_exact,G12_exact,G13_exact;
end


text:
import Base: eachindex, rand
import Statistics: mean, cov
import Distributions: logpdf, ContinuousMultivariateDistribution
export AbstractGP, kernel, rand, logpdf, elbo, marginal_cov, marginal_std, marginals,
    mean_vec, xcov, GPC

# A collection of GPs (GPC == "GP Collection"). Used to keep track of internals.
mutable struct GPC
    n::Int
    GPC() = new(0)
end

"""
    AbstractGaussianProcess <: ContinuousMultivariateDistribution

Any subtype of `AbstractGaussianProcess` must define methods of `mean` and `kernel`
"""
abstract type AbstractGaussianProcess <: ContinuousMultivariateDistribution end
const AbstractGP = AbstractGaussianProcess

"""
    ==(f::AbstractGP, g::AbstractGP)

Two `AbstractGP`s said to be equal if their marginal distributions are equal. That is, if
the mean function and kernels of both are equal.
"""
==(f::AbstractGP, g::AbstractGP) = (mean(f) == mean(g)) && (kernel(f) == kernel(g))

"""
    length(f::AbstractGP)

The dimensionality of the space over which `f` distributes. May be finite (multivariate
Normal) or infinite (usually, distributes over a function).
"""
length(f::AbstractGP) = length(mean(f))

"""
    eachindex(f::AbstractGP)

Enumerate the dimensions of the space over which `f` distributes if such an enumeration
exists and is finite.
"""
eachindex(f::AbstractGP) = eachindex(mean(f))

"""
    mean_vec(f::AbstractGP)

The `Vector` representation of the mean function of `f`, if `isfinite(length(f))`.
"""
mean_vec(f::AbstractGP) = AbstractVector(mean(f))

"""
    cov(f::AbstractGP)

The covariance matrix of `f`, if `isfinite(length(f))`.
"""
cov(f::AbstractGP) = AbstractMatrix(kernel(f))

"""
    xcov(f::AbstractGP, g::AbstractGP)

The cross-covariance matrix between `f` and `g`, if
`isfinite(length(f)) && isfinite(length(g))`.
"""
xcov(f::AbstractGP, g::AbstractGP) = AbstractMatrix(kernel(f, g))

"""
    marginal_cov(f::AbstractGP)

Efficiently compute of `diag(cov(f))`.
"""
marginal_cov(f::AbstractGP) = map(kernel(f), eachindex(f))

"""
    marginal_std(f::AbstractGP)

Efficiently compute `sqrt.(diag(cov(f)))`.
"""
marginal_std(f::AbstractGP) = sqrt.(marginal_cov(f))

"""
    marginals(f::AbstractGP)

Sugar, equivalent to `(mean(f), marginal_std(f))`.
"""
marginals(f::AbstractGP) = (mean_vec(f), marginal_std(f))

"""
    rand(rng::AbstractRNG, f::AbstractGP, N::Int=1)

Obtain `N` independent samples from the GP `f` using `rng`, if `isfinite(length(f))`.
"""
function rand(rng::AbstractRNG, f::AbstractGP, N::Int)
    return mean_vec(f) .+ chol(cov(f))' * randn(rng, length(f), N)
end
rand(rng::AbstractRNG, f::AbstractGP) = vec(rand(rng, f, 1))
rand(f::AbstractGP, N::Int) = rand(Random.GLOBAL_RNG, f, N)
rand(f::AbstractGP) = vec(rand(f, 1))

"""
    logpdf(f::AbstractGP, y::AbstractVector{<:Real})

The log probability density of `y` under `f`.
"""
function logpdf(f::AbstractGP, y::AbstractVector{<:Real})
    μ, Σ = mean_vec(f), cov(f)
    return -0.5 * (length(y) * log(2π) + logdet(Σ) + Xt_invA_X(Σ, y - μ))
end

"""
    elbo(f::AbstractGP, y::AbstractVector{<:Real}, u::AbstractGP, σ::Real)

The saturated Titsias-ELBO. Requires a reasonable degree of care.
"""
function elbo(f::AbstractGP, y::AV{<:Real}, u::AbstractGP, σ::Real)
    Γ = (chol(cov(u))' \ xcov(u, f)) ./ σ
    Ω, δ = LazyPDMat(Γ * Γ' + I, 0), y - mean_vec(f)
    return -0.5 * (length(y) * log(2π * σ^2) + logdet(Ω) - sum(abs2, Γ) +
        (sum(abs2, δ) - sum(abs2, chol(Ω)' \ (Γ * δ)) + sum(marginal_cov(f))) / σ^2)
end


text:
g!(du::DNSField{T},u,p,t) where T = du .= p.F
g!(du::GSSField{T},u,p,t) where T = du.x[1] .= p.F.x[1]

function f!(du::DNSField{T},u::DNSField{T},p::NLParams{T},t) where T<:AbstractFloat
    nx,ny,A,B,Cp,Cm = p.nx,p.ny,p.A,p.B,p.C⁺,p.C⁻
    du .= zero(Complex{T})
    # zonal constant and linear terms
    @inbounds for n1=1:ny-1
        m1 = 0
        du[n1+ny,m1+1] += A[n1+ny]
        du[n1+ny,m1+1] += B[n1+ny,m1+1]*u[n1+ny,m1+1]
    end
    @inbounds for m1=1:nx-1
        @inbounds for n1=-ny+1:ny-1
            # non-zonal linear terms
            du[n1+ny,m1+1] += B[n1+ny,m1+1]*u[n1+ny,m1+1]
            # ++ interactions
            @inbounds for m2=0:min(m1,nx-1-m1)
                n2min = m2 == 0 ? 1 : -ny+1
                @inbounds for n2=max(n2min,-ny+1-n1):min(ny-1,ny-1-n1)
                    m = m1 + m2
                    n = n1 + n2
                    du[n+ny,m+1] += Cp[n2+ny,m2+1,n1+ny,m1+1]*u[n1+ny,m1+1]*u[n2+ny,m2+1]
                end
            end
            # +- interactions
            @inbounds for m2=0:m1
                n2min = m2 == 0 ? 1 : -ny+1
                n2max = m2 == m1 ? n1 - 1 : ny-1
                @inbounds for n2=max(n2min,n1-ny+1):min(n2max,n1+ny-1)
                    m = m1 - m2
                    n = n1 - n2
                    du[n+ny,m+1] += Cm[n2+ny,m2+1,n1+ny,m1+1]*u[n1+ny,m1+1]*conj(u[n2+ny,m2+1])
                end
            end
        end
    end
    nothing
end

function f!(du::DNSField{T},u::DNSField{T},p::GQLParams{T},t) where T<:AbstractFloat
    nx,ny,Λ,A,B,Cp,Cm = p.nx,p.ny,p.Λ,p.A,p.B,p.C⁺,p.C⁻
    du .= zero(Complex{T})
    # zonal constant and linear terms
    @inbounds for n1=1:ny-1
        m1 = 0
        du[n1+ny,m1+1] += A[n1+ny]
        du[n1+ny,m1+1] += B[n1+ny,m1+1]*u[n1+ny,m1+1]
    end
    @inbounds for m1=1:Λ
        @inbounds for n1=-ny+1:ny-1
            # low mode linears
            du[n1+ny,m1+1] += B[n1+ny,m1+1]*u[n1+ny,m1+1]
            # L + L = L
            @inbounds for m2=0:min(m1,Λ-m1)
                n2min = m2 == 0 ? 1 : -ny+1
                @inbounds for n2=max(n2min,-ny+1-n1):min(ny-1,ny-1-n1)
                    m = m1 + m2
                    n = n1 + n2
                    du[n+ny,m+1] += Cp[n2+ny,m2+1,n1+ny,m1+1]*u[n1+ny,m1+1]*u[n2+ny,m2+1]
                end
            end
            # L - L = L
            for m2=0:m1
                n2min = m2 == 0 ? 1 : -ny+1
                n2max = m2 == m1 ? n1 - 1 : ny-1
                for n2=max(n2min,n1-ny+1):min(n2max,n1+ny-1)
                    m = m1 - m2
                    n = n1 - n2
                    du[n+ny,m+1] += Cm[n2+ny,m2+1,n1+ny,m1+1]*u[n1+ny,m1+1]*conj(u[n2+ny,m2+1])
                end
            end
        end
    end
    @inbounds for m1=Λ+1:nx-1
        @inbounds for n1=-ny+1:ny-1
            # high mode linears
            du[n1+ny,m1+1] += B[n1+ny,m1+1]*u[n1+ny,m1+1]
            # H - H = L
            @inbounds for m2=max(Λ+1,m1-Λ):m1
                n2max = m2 == m1 ? n1 - 1 : ny-1
                @inbounds for n2=max(-ny+1,n1-ny+1):min(n2max,n1+ny-1)
                    m = m1 - m2
                    n = n1 - n2
                    du[n+ny,m+1] += Cm[n2+ny,m2+1,n1+ny,m1+1]*u[n1+ny,m1+1]*conj(u[n2+ny,m2+1])
                end
            end
            # H + L = H
            @inbounds for m2=0:min(nx-1-m1,Λ)
                n2min = m2 == 0 ? 1 : -ny+1
                @inbounds for n2=max(n2min,-ny+1-n1):min(ny-1,ny-1-n1)
                    m = m1 + m2
                    n = n1 + n2
                    du[n+ny,m+1] += Cp[n2+ny,m2+1,n1+ny,m1+1]*u[n1+ny,m1+1]*u[n2+ny,m2+1]
                end
            end
            # H - L = H
            @inbounds for m2=0:1:min(Λ,m1 - Λ - 1)
                n2min = m2 == 0 ? 1 : -ny+1
                @inbounds for n2=max(n2min,n1-ny+1):1:min(ny-1,n1+ny-1)
                    m = m1 - m2
                    n = n1 - n2
                    du[n+ny,m+1] += Cm[n2+ny,m2+1,n1+ny,m1+1]*u[n1+ny,m1+1]*conj(u[n2+ny,m2+1])
                end
            end
        end
    end
    nothing
end

function f!(du::GSSField{T},u::GSSField{T},p::GCE2Params,t) where T<:AbstractFloat
    nx,ny,Λ,A,B,Cp,Cm,F,dx,dy,temp = p.nx,p.ny,p.Λ,p.A,p.B,p.C⁺,p.C⁻,p.F,p.dx,p.dy,p.temp
    du .= zero(Complex{T})
    # zonal terms
    dx .= zero(Complex{T})
    @inbounds for n1=1:ny-1
        m1 = 0
        dx[n1+ny,1] += A[n1+ny]
        dx[n1+ny,m1+1] += B[n1+ny,m1+1]*u.x[1][n1+ny,m1+1]
    end
    # low modes
    @inbounds for m1=1:Λ
        @inbounds for n1=-ny+1:ny-1
            dx[n1+ny,m1+1] += B[n1+ny,m1+1]*u.x[1][n1+ny,m1+1]
            # L + L = L
            @inbounds for m2=0:min(m1,Λ-m1)
                n2min = m2 == 0 ? 1 : -ny+1
                @inbounds for n2=max(n2min,-ny+1-n1):min(ny-1,ny-1-n1)
                    m = m1 + m2
                    n = n1 + n2
                    dx[n+ny,m+1] += Cp[n2+ny,m2+1,n1+ny,m1+1]*u.x[1][n1+ny,m1+1]*u.x[1][n2+ny,m2+1]
                end
            end
            # L - L = L
            @inbounds for m2=0:m1
                n2min = m2 == 0 ? 1 : -ny+1
                n2max = m2 == m1 ? n1 - 1 : ny-1
                @inbounds for n2=max(n2min,n1-ny+1):min(n2max,n1+ny-1)
                    m = m1 - m2
                    n = n1 - n2
                    dx[n+ny,m+1] += Cm[n2+ny,m2+1,n1+ny,m1+1]*u.x[1][n1+ny,m1+1]*conj(u.x[1][n2+ny,m2+1])
                end
            end
        end
    end
    # high modes
    temp .= zero(Complex{T})
    @inbounds for m1=Λ+1:nx-1
        @inbounds for n1=-ny+1:ny-1
            # H - H = L
            @inbounds for m2=max(Λ+1,m1-Λ):m1
                n2max = m2 == m1 ? n1 - 1 : ny-1
                @inbounds for n2=max(-ny+1,n1-ny+1):min(n2max,n1+ny-1)
                    m = m1 - m2
                    n = n1 - n2
                    # note: u.x[2] contains H2*conj(H1) so H-H is conj(H2)*H1
                    dx[n+ny,m+1] += Cm[n2+ny,m2+1,n1+ny,m1+1]*conj(u.x[2][n2+ny,m2-Λ,n1+ny,m1-Λ])
                end
            end
            # H + L = H
            @inbounds for m2=0:min(nx-1-m1,Λ)
                n2min = m2 == 0 ? 1 : -ny+1
                @inbounds for n2=max(n2min,-ny+1-n1):min(ny-1,ny-1-n1)
                    m = m1 + m2
                    n = n1 + n2
                    temp[n1+ny,m1-Λ,n+ny,m-Λ] += Cp[n2+ny,m2+1,n1+ny,m1+1]*u.x[1][n2+ny,m2+1]
                end
            end
            # H - L = H
            @inbounds for m2=0:min(Λ,m1 - Λ - 1)
                n2min = m2 == 0 ? 1 : -ny+1
                @inbounds for n2=max(n2min,n1-ny+1):min(ny-1,n1+ny-1)
                    m = m1 - m2
                    n = n1 - n2
                    temp[n1+ny,m1-Λ,n+ny,m-Λ] += Cm[n2+ny,m2+1,n1+ny,m1+1]*conj(u.x[1][n2+ny,m2+1])
                end
            end
        end
    end
    copyto!(du.x[1],dx)
    # du.x[1] .= dx
    # H'*H
    @inbounds for m3=Λ+1:nx-1
        @inbounds for n3=-ny+1:ny-1
            @inbounds for m=Λ+1:nx-1
                @inbounds for n=-ny+1:ny-1
                    # accumulator = zero(Complex{T})
                    accumulator::Complex{T} = B[n+ny,m+1]*u.x[2][n+ny,m-Λ,n3+ny,m3-Λ]
                    accumulator += F.x[2][n+ny,m-Λ,n3+ny,m3-Λ]
                    # from H+L
                    @inbounds for m1=max(Λ+1,m-Λ):min(nx-1,m)
                        n2min = m1 == m ? 1 : -ny+1
                        @inbounds for n1=max(-ny+1,n-ny+1):min(n-n2min,ny-1)
                            accumulator += temp[n1+ny,m1-Λ,n+ny,m-Λ]*u.x[2][n1+ny,m1-Λ,n3+ny,m3-Λ]
                        end
                    end
                    # from H-L
                    @inbounds for m1=max(Λ+1,m):min(nx-1,m+Λ)
                        n2max = m1 == m ? -1 : ny-1
                        @inbounds for n1=max(-ny+1,n-n2max):min(n+ny-1,ny-1)
                            accumulator += temp[n1+ny,m1-Λ,n+ny,m-Λ]*u.x[2][n1+ny,m1-Λ,n3+ny,m3-Λ]
                        end
                    end
                    dy[n+ny,m-Λ,n3+ny,m3-Λ] = accumulator
                end
            end
        end
    end
    # permutedims!(temp,du.x[2],[3,4,1,2])
    @inbounds for m3=Λ+1:nx-1
        @inbounds for n3=-ny+1:ny-1
            @inbounds for m=Λ+1:nx-1
                @inbounds for n=-ny+1:ny-1
                    du.x[2][n+ny,m-Λ,n3+ny,m3-Λ] = dy[n+ny,m-Λ,n3+ny,m3-Λ] + conj(dy[n3+ny,m3-Λ,n+ny,m-Λ])
                end
            end
        end
    end
    nothing
end

function f!(du::DSSField{T},u::DSSField{T},p::CE2Params,t) where T<:AbstractFloat
    nx,ny,A,B,Cp,Cm,F,dy,temp = p.nx,p.ny,p.A,p.B,p.C⁺,p.C⁻,p.F,p.dy,p.temp
    du .= zero(Complex{T})
    @inbounds for n1=1:ny-1
        du.x[1][n1+ny] += A[n1+ny]
        du.x[1][n1+ny] += B[n1+ny,1]*u.x[1][n1+ny]
        # M + M = M
        @inbounds for n2=max(1,-ny+1-n1):min(ny-1,ny-1-n1)
            du.x[1][n1+n2+ny] += Cp[n2+ny,1,n1+ny,1]*u.x[1][n1+ny]*u.x[1][n2+ny]
        end
        # M - M = M
        @inbounds for n2=max(1,n1-ny+1):min(n1-1,n1+ny-1)
            du.x[1][n1-n2+ny] += Cm[n2+ny,1,n1+ny,1]*u.x[1][n1+ny]*conj(u.x[1][n2+ny])
        end
    end
    temp .= zero(Complex{T})
    @inbounds for m=1:nx-1
        @inbounds for n1=-ny+1:ny-1
            # E - E = M
            @inbounds for n2=max(-ny+1,n1-ny+1):min(n1-1,n1+ny-1)
                # note: u.x[2] contains H2*conj(H1) so H-H is conj(H2)*H1
                du.x[1][n1-n2+ny] += Cm[n2+ny,m+1,n1+ny,m+1]*conj(u.x[2][n2+ny,n1+ny,m])
            end
            # E + M = E
            @inbounds for n2=max(1,-ny+1-n1):min(ny-1,ny-1-n1)
                temp[n1+ny,n1+n2+ny,m] += Cp[n2+ny,1,n1+ny,m+1]*u.x[1][n2+ny]
            end
            # E - M = E
            @inbounds for n2=max(1,n1-ny+1):min(ny-1,n1+ny-1)
                temp[n1+ny,n1-n2+ny,m] += Cm[n2+ny,1,n1+ny,m+1]*conj(u.x[1][n2+ny])
            end
        end
    end
    # H'*H
    @inbounds for m3=1:nx-1
        @inbounds for n3=-ny+1:ny-1
            @inbounds for n=-ny+1:ny-1
                du.x[2][n+ny,n3+ny,m3] += B[n+ny,m3+1]*u.x[2][n+ny,n3+ny,m3]
                du.x[2][n+ny,n3+ny,m3] += F[n+ny,n3+ny,m3]

                accumulator = zero(Complex{T})
                @inbounds for n1=max(-ny+1,n-(ny-1)):min(n-1,ny-1)
                    accumulator += temp[n1+ny,n+ny,m3]*u.x[2][n1+ny,n3+ny,m3]
                end
                @inbounds for n1=max(-ny+1,n+1):min(n+ny-1,ny-1)
                    accumulator += temp[n1+ny,n+ny,m3]*u.x[2][n1+ny,n3+ny,m3]
                end
                du.x[2][n+ny,n3+ny,m3] += accumulator
                dy[n3+ny,n+ny,m3] = conj(du.x[2][n+ny,n3+ny,m3]) # build conjugate transpose
            end
        end
    end
    @inbounds for m3=1:nx-1
        @inbounds for n3=-ny+1:ny-1
            @inbounds for n=-ny+1:ny-1
                du.x[2][n+ny,n3+ny,m3] += dy[n+ny,n3+ny,m3]
            end
        end
    end
    nothing
end


text:
"""
`small(x)` converts `x` to `Int`-level precision if there is
no loss. Applies to integer, rational, Gaussian integer,
Gaussian rational, and `SimplePolynomial`. It is a way of
undoing `big` (if possible).
"""
function small(x::T) where {T<:Integer}
    try
        x = Int(x)
    catch
    end
    return x
end

function small(x::Complex{T}) where {T<:Integer}
    try
        x = Complex{Int}(x)
    catch
    end
    return x
end

function small(x::Rational{T}) where {T<:Integer}
    try
        x = Rational{Int}(x)
    catch
    end
    return x
end

function small(x::Complex{Rational{T}}) where {T<:Integer}
    try
        x = Complex{Rational{Int}}(x)
    catch
    end
    return x
end

small(x::AbstractMod) = x


text:
# Create the sky and the entire Obs_dataset
include("mapmaking.jl")

using FFTW
using ConjGrad
using LinearAlgebra, SparseArrays
using Test
using Plots

# GCP!!!!!!!!!!!!!!!!!!!!!!!, No λ functions !
# Costruisco FTZ
# Arrivo qui che ho il dataset delle baselines

output = Array{Float64, 2}(undef, 1024, 1024)


text:
using MomentTensors, Test

@testset "Conversion" begin
    @test m0(mw(1e18)) ≈ 1e18
    @test mw(m0(5)) ≈ 5
    @test Array(MT{Float16}(0, 45, 90, 1)) isa Array{Float16,2}
end


text:
#
# Block
#

struct Block{N}
  data::NTuple{N,UInt}
  hash::UInt
  function Block{N}(data::NTuple{N,UInt}) where {N}
    h = _hash(data)
    return new{N}(data, h)
  end
  function Block{0}(::Tuple{})
    h = _hash(())
    return new{0}((), h)
  end
end

#
# Constructors
#

Block{N}(t::Tuple{Vararg{<:Any,N}}) where {N} = Block{N}(UInt.(t))

Block{N}(I::CartesianIndex{N}) where {N} = Block{N}(I.I)

Block{N}(v::MVector{N}) where {N} = Block{N}(Tuple(v))

Block{N}(v::SVector{N}) where {N} = Block{N}(Tuple(v))

Block(b::Block) = b

Block(I::CartesianIndex{N}) where {N} = Block{N}(I)

Block(v::MVector{N}) where {N} = Block{N}(v)

Block(v::SVector{N}) where {N} = Block{N}(v)

Block(t::NTuple{N,UInt}) where {N} = Block{N}(t)

Block(t::Tuple{Vararg{<:Any,N}}) where {N} = Block{N}(t)

Block(::Tuple{}) where {N} = Block{0}(())

Block(I::Union{Integer,Block{1}}...) = Block(I)

#
# Conversions
#

CartesianIndex(b::Block) = CartesianIndex(Tuple(b))

Tuple(b::Block{N}) where {N} = NTuple{N,UInt}(b.data)

convert(::Type{Block}, I::CartesianIndex{N}) where {N} = Block{N}(I.I)

convert(::Type{Block{N}}, I::CartesianIndex{N}) where {N} = Block{N}(I.I)

convert(::Type{Block}, t::Tuple) where {N} = Block(t)

convert(::Type{Block{N}}, t::Tuple) where {N} = Block{N}(t)

(::Type{IntT})(b::Block{1}) where {IntT<:Integer} = IntT(only(b))

#
# Getting and setting fields
#

gethash(b::Block) = b.hash[]

sethash!(b::Block, h::UInt) = (b.hash[] = h; return b)

#
# Basic functions
#

length(::Block{N}) where {N} = N

isless(b1::Block, b2::Block) = isless(Tuple(b1), Tuple(b2))

iterate(b::Block, args...) = iterate(b.data, args...)

@propagate_inbounds function getindex(b::Block, i::Integer)
  return b.data[i]
end

@propagate_inbounds function setindex(b::Block{N}, val, i::Integer) where {N}
  return Block{N}(setindex(b.data, UInt(val), i))
end

ValLength(::Type{<:Block{N}}) where {N} = Val{N}

deleteat(b::Block, pos) = Block(deleteat(Tuple(b), pos))

insertafter(b::Block, val, pos) = Block(insertafter(Tuple(b), UInt.(val), pos))

getindices(b::Block, I) = getindices(Tuple(b), I)

#
# checkbounds
#

# XXX: define this properly
checkbounds(::Tensor, ::Block) = nothing

#
# Hashing
#

# Borrowed from:
# https://github.com/JuliaLang/julia/issues/37073
# This is the same as Julia's Base tuple hash, but is
# a bit faster.
_hash(t::Tuple) = _hash(t, zero(UInt))
_hash(::Tuple{}, h::UInt) = h + Base.tuplehash_seed
@generated function _hash(b::NTuple{N}, h::UInt) where {N}
  quote
    out = h + Base.tuplehash_seed
    @nexprs $N i -> out = hash(b[$N - i + 1], out)
  end
end

if VERSION < v"1.7.0-DEV.933"
  # Stop inlining after some number of arguments to avoid code blowup
  function _hash(t::Base.Any16, h::UInt)
    out = h + Base.tuplehash_seed
    for i in length(t):-1:1
      out = hash(t[i], out)
    end
    return out
  end
else
  # Stop inlining after some number of arguments to avoid code blowup
  function _hash(t::Base.Any32, h::UInt)
    out = h + Base.tuplehash_seed
    for i in length(t):-1:1
      out = hash(t[i], out)
    end
    return out
  end
end

hash(b::Block) = UInt(b.hash)
hash(b::Block, h::UInt) = h + hash(b)

#
# Custom NTuple{N, Int} hashes
# These are faster, but have a lot of collisions
#

# Borrowed from:
# https://stackoverflow.com/questions/20511347/a-good-hash-function-for-a-vector
# This seems to have a lot of clashes
#function hash(b::Block, seed::UInt)
#  h = UInt(0x9e3779b9)
#  for n in b
#    seed ⊻= n + h + (seed << 6) + (seed >> 2)
#  end
#  return seed
#end

# Borrowed from:
# http://www.docjar.com/html/api/java/util/Arrays.java.html
# Could also consider uring the CPython tuple hash:
# https://github.com/python/cpython/blob/0430dfac629b4eb0e899a09b899a494aa92145f6/Objects/tupleobject.c#L406
#function hash(b::Block, h::UInt)
#  h += Base.tuplehash_seed
#  for n in b
#    h = 31 * h + n ⊻ (n >> 32)
#  end
#  return h
#end

#
# Printing for Block type
#

show(io::IO, mime::MIME"text/plain", b::Block) = print(io, "Block$(Int.(Tuple(b)))")

show(io::IO, b::Block) = show(io, MIME("text/plain"), b)


text:
using Hankel
using Base.Test

using PyCall
pygui(:qt)
using PyPlot

function hankel_test(f1, f2, fun_latex,
                     p, N, R, Va, Na=1000,
                     plot_error=false,
                     fignum=1)
    r = linspace(0,R,Na)
    ν_a = linspace(0,Va,Na)
    figure(fignum,figsize=[12, 17])
    clf()
    for i = 1:length(p)
        subplot(length(p)+1,1,1)
        plot(r, f1(r,p[i]), label = latexstring("\$f_{1,$(p[i])}\$"))
        margins(0.05, 0.1)
        title(latexstring("\$f_{1,p}=$(fun_latex)\$"))
        legend()

        subplot(length(p)+1,1,i+1)
        ax1 = gca()
        ax2 = ax1[:twinx]()
        lines = []
        append!(lines, ax1[:plot](ν_a, abs(f2(ν_a,p[i])),
                                  "k", label="Exact", linewidth = 2.0))
        function test_discrete(trf, name, args...)
            ν_d,f2_d = trf(r -> f1(r,p[i]), p[i], args...)
            append!(lines, ax1[:plot](ν_d, abs(f2_d), "-", label=name))
            if plot_error
                append!(lines, ax2[:semilogy](ν_d, abs(abs(f2(ν_d,p[i]))-abs(f2_d))./abs(maximum(f2_d)), ":",
                                              label="$name error"))
            end
        end

        test_discrete(guizar, "Guizar", R, N)

        ax1[:margins](0.05, 0.1)
        ax1[:set_ylabel]("Function value")
        ax1[:set_title](latexstring("\$f_{2,$(p[i])}\$"))
        ax2[:margins](0.05, 0.1)
        ax2[:set_ylabel]("Dynamic error")
        legend(lines, [l[:get_label]() for l in lines])
    end
    tight_layout()
end

γ = 5
hankel_test((r,p) -> sinc(2π*γ*r),
            function(ν,p)
            function f(ν,p)
            s2 = γ^2 - ν.^2
            if 0<=ν && ν<=γ
            ν^p*cos(p*0.5π)/(2π*γ*√(s2)*(γ*√(s2))^p)
            elseif ν==γ
            Inf
            else
            sin(p*asin(γ./ν))/(2π*γ*√(-s2))
            end
            end
            map(x -> f(x,p), ν)
            end,
            "\\operatorname{sinc}(2\\pi\\gamma r),\\quad \\gamma=$γ",
                0:2, 256, 3.0, 20.0, 1000, true, 1)

a = 0.5
hankel_test((r,p) -> r.^p.*exp(-a*r.^2),
            (ν,p) -> (π/a).^(p+1) * ν.^p .* exp(-π^2*ν.^2/a),
            "r^p\\exp(-ar^2),\\quad a=$a",
            0:4, 256, 20, 7, 1000, true, 2)

hankel_test((r,p) -> ((0 .<= r) .== (r .<= 1))*1.0,
            (ν,p) -> besselj(1,2π*ν)./ν,
            "\\theta(r)",
            0, 128, 20, 15, 1000, true, 3)

# write your own tests here
@test 1 == 1


text:
using DataDeps: DataDep, register, unpack, @datadep_str

"""
   GraphDataset

Abstract graph dataset type.

## mandatory methods to implement
* [`loadgraph`](@ref) and/or [`loadgraphs`](@ref)
* [`dataset_root_directory_name`](@ref)
* [`dataset_url`](@ref)
* [`dataset_hash`](@ref)
* [`dataset_message`](@ref)

## optional methods to implement
* [`loadreadme`](@ref)
* [`dataset_path`](@ref)
"""
abstract type GraphDataset end

## ----------------------------------------
##    methods to implement
## ----------------------------------------

function dataset_root_directory_name end

function dataset_url end

dataset_hash(ds::GraphDataset) = nothing

function dataset_message(ds::GraphDataset)

    @warn "No message specified for dataset $ds"
    return ""
end

"""
    loadgraphs(ds::GraphDataset; resolve_categories::Bool=false)

Loads multiple graphs from a dataset `ds`.

Either this method or `loadgraph` must be implement for new subytpes of `GraphDataset`.

# Keywords
- `resolve_categories`: Some graph metadata might be of categorical form (i.e strings).
    If this argument is true, try to resolve that category instead of keeping a numerical
    value.

### See also
[`loadgraph`](@Ref), `GraphDataset`](@Ref)
"""
loadgraphs(ds::GraphDataset; resolve_categories::Bool)

## ----------------------------------------
##    optional methods to implement
## ----------------------------------------

function loadreadme end

dataset_path(ds::GraphDataset) = @datadep_str(dataset_root_directory_name(ds))

## ----------------------------------------
##    various methods
## ----------------------------------------

_registered_datasets = IdDict{Type{<:GraphDataset}, String}()

"""
    list_datasets()

List available graph datasets
"""
function list_datasets()

    foreach(println, sort(collect(values(_registered_datasets))))
    return nothing
end

function register_dataset(ds::GraphDataset, prefix_modules = String[])

    register(DataDep(
        dataset_root_directory_name(ds),
        dataset_message(ds::GraphDataset),
        dataset_url(ds),
        dataset_hash(ds),
        post_fetch_method = unpack
    ))

    name_in_list = join(map(s -> "$s.", prefix_modules)) * string(nameof(typeof(ds)))
    _registered_datasets[typeof(ds)] = name_in_list
end


text:
# Code used to use phi, a 1-parameter function induced by f and s
# Need to pass s as an explicit parameter
function alphatry{T}(alpha::T,
                     df,
                     x::Array,
                     s::Array,
                     xtmp::Array,
                     gtmp::Array,
                     lsr::LineSearchResults,
                     psi1::Real = convert(T,0.2),
                     psi2::Real = convert(T,2),
                     psi3::Real = convert(T,0.1),
                     iterfinitemax::Integer = ceil(Integer, -log2(eps(T))),
                     alphamax::Real = convert(T, Inf),
                     verbose::Bool = false)
    f_calls = 0
    g_calls = 0

    phi0 = lsr.value[1]
    dphi0 = lsr.slope[1]

    alphatest = psi1 * alpha
    alphatest = min(alphatest, alphamax)

    # Use xtmp here
    phitest = df.f(x + alphatest * s)
    f_calls += 1

    iterfinite = 1
    while !isfinite(phitest)
        alphatest = psi3 * alphatest
        # Use xtmp here
        phitest = df.f(x + alphatest * s)
        f_calls += 1
        lsr.nfailures += 1
        iterfinite += 1
        if iterfinite >= iterfinitemax
            return zero(T), true, f_calls, g_calls
#             error("Failed to achieve finite test value; alphatest = ", alphatest)
        end
    end
    a = ((phitest-phi0)/alphatest - dphi0)/alphatest  # quadratic fit
    if verbose == true
        println("quadfit: alphatest = ", alphatest,
                ", phi0 = ", phi0,
                ", phitest = ", phitest,
                ", quadcoef = ", a)
    end
    mayterminate = false
    if isfinite(a) && a > 0 && phitest <= phi0
        alpha = -dphi0 / 2 / a # if convex, choose minimum of quadratic
        if alpha == 0
            error("alpha is zero. dphi0 = ", dphi0, ", phi0 = ", phi0, ", phitest = ", phitest, ", alphatest = ", alphatest, ", a = ", a)
        end
        if alpha <= alphamax
            mayterminate = true
        else
            alpha = alphamax
            mayterminate = false
        end
        if verbose == true
            println("alpha guess (quadratic): ", alpha,
                    ",(mayterminate = ", mayterminate, ")")
        end
    else
        if phitest > phi0
            alpha = alphatest
        else
            alpha *= psi2 # if not convex, expand the interval
        end
    end
    alpha = min(alphamax, alpha)
    if verbose == true
        println("alpha guess (expand): ", alpha)
    end
    return alpha, mayterminate, f_calls, g_calls
end

# Generate initial guess for step size (HZ, stage I0)
function alphainit{T}(alpha::Real,
                      x::Array{T},
                      gr::Array,
                      f_x::Real,
                      psi0::T = convert(T,0.01))
    if isnan(alpha)
        alpha = one(T)
        gr_max = maximum(abs(gr))
        if gr_max != 0.0
            x_max = maximum(abs(x))
            if x_max != 0.0
                alpha = psi0 * x_max / gr_max
            elseif f_x != 0.0
                alpha = psi0 * abs(f_x) / vecnorm(gr)
            end
        end
    end
    return alpha
end


text:
##############################################################################
##
## Old one
##
##
##############################################################################

function oldparse_fixedeffect(df::AbstractDataFrame, feformula::FormulaTerm)
    fe = FixedEffect[]
    id = Symbol[]
    for term in eachterm(feformula.rhs)
        result = oldparse_fixedeffect(df, term, feformula)
        if result != nothing
            push!(fe, result[1])
            push!(id, result[2])
        end
    end
    return fe, id
end

# Constructors from dataframe + Term
function oldparse_fixedeffect(df::AbstractDataFrame, a::Term, feformula::FormulaTerm)
    v = df[!, Symbol(a)]
    if isa(v, CategoricalVector)
        return FixedEffect(v), Symbol(a)
    else
        # x from x*id -> x + id + x&id
        if !any(isa(term, InteractionTerm) & (a ∈ terms(term)) for term in eachterm(feformula.rhs))
               error("The term $(a) in fe= is a continuous variable. Convert it to a categorical variable using 'categorical'.")
        end
    end
end

# Constructors from dataframe + InteractionTerm
function oldparse_fixedeffect(df::AbstractDataFrame, a::InteractionTerm, feformula::FormulaTerm)
    factorvars, interactionvars = _split(df, a)
    if !isempty(factorvars)
        # x1&x2 from (x1&x2)*id
        fe = FixedEffect((df[!, v] for v in factorvars)...; interaction = old_multiply(df, interactionvars))
        id = old_name(Symbol.(terms(a)))
        return fe, id
    end
end

function _split(df::AbstractDataFrame, a::InteractionTerm)
    factorvars, interactionvars = Symbol[], Symbol[]
    for s in terms(a)
        s = Symbol(s)
        isa(df[!, s], CategoricalVector) ? push!(factorvars, s) : push!(interactionvars, s)
    end
    return factorvars, interactionvars
end

function old_multiply(df, ss::Vector{Symbol})
    if isempty(ss)
        out = Ones(size(df, 1))
    else
        out = ones(size(df, 1))
        for j in eachindex(ss)
            old_multiply!(out, df[!, ss[j]])
        end
    end
    return out
end

function old_multiply!(out, v)
    for i in eachindex(out)
        if v[i] === missing
            # may be missing when I remove singletons
            out[i] = 0.0
        else
            out[i] = out[i] * v[i]
        end
    end
end

function old_name(s::Vector{Symbol})
    if isempty(s)
        out = nothing
    else
        out = Symbol(reduce((x1, x2) -> string(x1)*"x"*string(x2), s))
    end
    return out
end

struct ModelTerm
    f::FormulaTerm
    dict::Dict{Symbol, Any}
end

ModelTerm(f::FormulaTerm; kwargs...) = ModelTerm(f, Dict(pairs(kwargs)...))
function Base.show(io::IO, m::ModelTerm)
    println(io, m.f)
    for (k, v) in m.dict
        println(io, k, ": ", v)
    end
end

import StatsModels: capture_call
macro model(ex, kws...)
    f = StatsModels.terms!(StatsModels.sort_terms!(StatsModels.parse!(ex)))
    d = Dict{Symbol, Any}()
    for kw in kws
       isa(kw, Expr) &&  kw.head== :(=) || throw("All arguments of @model, except the first one, should be keyboard arguments")
       if kw.args[1] == :fe
            @warn "The keyword argument fe is deprecated. Instead of @model(y ~ x, fe = state + year),  write @formula(y ~ x + fe(state) + fe(year))"
            d[:feformula] = kw.args[2]
        elseif kw.args[1] == :ife
                 @warn "The keyword argument ife is deprecated. Instead of @model(y ~ x, ife = (state + year, 2)),  write @formula(y ~ x + ife(state, year, 2))"
            d[:ifeformula] = kw.args[2]
        elseif kw.args[1] == :vcov
            d[:vcovformula] = kw.args[2]
            @warn "The keyword argument vcov is deprecated. Instead of reg(df, @model(y ~ x, vcov = cluster(State))),  write reg(df, @formula(y ~ x), Vcov.cluster(:State))"
        elseif kw.args[1] == :subset
            d[:subsetformula] = kw.args[2]
            @warn "The keyword argument subset is deprecated. Instead of reg(df, @model(y ~ x, subset = State .>= 30),  write reg(df, @formula(y ~ x), subset = df.State .>= 30))"
        elseif kw.args[1] == :weight
                 d[:weight] = kw.args[2]
                 @warn "The keyword argument weight is deprecated. Instead of reg(df, @model(y ~ x, weight = Pop),  write reg(df, @formula(y ~ x), weight = :Pop)"
        else
           d[kw.args[1]] = kw.args[2]
       end
    end
    :(ModelTerm($f, $d))
end

function evaluate_subset(df, ex::Expr)
    if ex.head == :call
        return Expr(ex.head, ex.args[1], (evaluate_subset(df, ex.args[i]) for i in 2:length(ex.args))...)
    else
        return Expr(ex.head, (evaluate_subset(df, ex.args[i]) for i in 1:length(ex.args))...)
    end
end
evaluate_subset(df, ex::Symbol) = df[!, ex]
evaluate_subset(df, ex)  = ex

function reg(df, m::ModelTerm;kwargs...)
    reg(df, m.f; m.dict..., kwargs...)
end

function partial_out(df, m::ModelTerm; kwargs...)
    partial_out(DataFrame(df), m.f; m.dict..., kwargs...)
end

function fes(args...)
    @warn "fes() is deprecated. Use fe()"
    fe(args...)
end


text:
@testset "integrate_gauss_legendre" begin
    # TODO: Test on randomly generated polynomials
    # TODO: Test with alternating endpoints
    # TODO: Test with wide endpoints
    # TODO: Test with non-finite endpoints

    @test Arblib.overlaps(
        ArbExtras.integrate_gauss_legendre(sin, Arb(0), Arb(1)),
        1 - cos(Arb(1)),
    )

    @test Arblib.overlaps(
        ArbExtras.integrate_gauss_legendre(x -> x^3, Arb(0), Arb(1)),
        Arb(1 // 4),
    )

    @test Arblib.overlaps(
        ArbExtras.integrate_gauss_legendre(x -> 1 + x + x^2 + x^3, Arb(1), Arb(3)),
        Arb(104 // 3),
    )

    @test Arblib.overlaps(
        ArbExtras.integrate_gauss_legendre(sin ∘ exp, Arb(0), Arb(4)),
        real(Arblib.integrate(sin ∘ exp, 0, 4)),
    )

    @test isnan(ArbExtras.integrate_gauss_legendre(inv, Arb(-1), Arb(1)))
end

@testset "integrate" begin
    # TODO: Add more tests

    @test Arblib.overlaps(ArbExtras.integrate(sin, Arb(0), Arb(1)), 1 - cos(Arb(1)))

    @test Arblib.overlaps(ArbExtras.integrate(x -> x^3, Arb(0), Arb(1)), Arb(1 // 4))

    @test Arblib.overlaps(
        ArbExtras.integrate(sin ∘ exp, Arb(0), Arb(4)),
        real(Arblib.integrate(sin ∘ exp, 0, 4)),
    )

    @test isnan(ArbExtras.integrate(inv, Arb(-1), Arb(1)))

    # depth_start > 0
    @test Arblib.overlaps(
        ArbExtras.integrate(sin, Arb(0), Arb(1), depth_start = 4),
        1 - cos(Arb(1)),
    )

    @test Arblib.overlaps(
        ArbExtras.integrate(x -> x^3, Arb(0), Arb(1), depth_start = 5),
        Arb(1 // 4),
    )
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder
using Pkg

julia_version = v"1.5.3"

name = "FastJet_Julia_Wrapper"
version = v"0.8.5"

# Collection of sources required to build FastJet_Julia_Wrapper
sources = [
	GitSource("https://github.com/jstrube/FastJet_Julia_Wrapper.git", "3fab5f91e88f2b348ec700c27733e0fcc661bd13"),
]

# Bash recipe for building across all platforms
script = raw"""
cd ${WORKSPACE}/srcdir/FastJet_Julia_Wrapper
mkdir build && cd build
cmake -DJulia_PREFIX=${prefix} -DCMAKE_INSTALL_PREFIX=${prefix} -DCMAKE_FIND_ROOT_PATH=${prefix} -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} -DCMAKE_BUILD_TYPE=Release ..
make -j${nproc}
make install
install_license $WORKSPACE/srcdir/FastJet_Julia_Wrapper/LICENSE.md
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = Platform[
    Platform("x86_64", "linux"; libc="glibc"),
    Platform("x86_64", "macos")
]
platforms = expand_cxxstring_abis(platforms)
# The products that we will ensure are always built
products = [
    LibraryProduct("libfastjetwrap", :libfastjetwrap)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="libcxxwrap_julia_jll",version=v"0.8.5")),
    Dependency("FastJet_jll"),
    BuildDependency(PackageSpec(name="libjulia_jll", version=julia_version))
]

build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies;
    preferred_gcc_version=v"7", julia_compat = "^$(julia_version.major).$(julia_version.minor)")

max_stars_repo_name: Djoop/Yggdrasil

max_issues_repo_name: Djoop/Yggdrasil

max_forks_repo_name: Djoop/Yggdrasil


text:
# Functions to manage pedigree models, primarily utilise JWAS package - Utilities for pedigree_run.jl
function get_namesA(markerInfo)
     return markerInfo[:,3] .* "_" .* markerInfo[:,1] .* "_" .* markerInfo[:,2] .* "_A"
end

function make_genotype_table(genoprefix)
    geno = read_binary_genotypes(genoprefix)
    markerInfo = geno.markerInfo
    namesA = get_namesA(markerInfo)
    g_tbl_names = vcat(["anml_key"], namesA)
    g_tbl_names = map(x -> Symbol(x), g_tbl_names)
    return table(map(x -> parse(Int64,x), geno.IID), bedvals[geno.binarygenotypes][:]; names = g_tbl_names), markerInfo
end
#
function join_tables(y_tbl, g_tbl)
    return join(y_tbl, g_tbl; how = :inner, lkey = :anml_key, rkey = :anml_key)
end

function center_phenotype(tbl)
    meanY = mean(JuliaDB.select(tbl, :yd_trait))
    tbl = IndexedTables.transform(tbl, :yd_trait => :yd_trait => x -> x - meanY)
    return tbl
end
#
function generate_genotypeclass_columns(tbl, markerInfo)
    m = size(markerInfo, 1)
    namesA = get_namesA(markerInfo)
    namesG11 = markerInfo[:,3] .* "_" .* markerInfo[:,1] .* "_" .* markerInfo[:,2] .* "_G11"
    namesG12 = markerInfo[:,3] .* "_" .* markerInfo[:,1] .* "_" .* markerInfo[:,2] .* "_G12"
    namesG22 = markerInfo[:,3] .* "_" .* markerInfo[:,1] .* "_" .* markerInfo[:,2] .* "_G22"
    tbl = IndexedTables.transform(tbl, Symbol(namesG11[1]) => Symbol(namesA[1]) => x -> x == 0 ? 1 : 0)
    tbl = IndexedTables.transform(tbl, Symbol(namesG12[1]) => Symbol(namesA[1]) => x -> x == 1 ? 1 : 0)
    tbl = IndexedTables.transform(tbl, Symbol(namesG22[1]) => Symbol(namesA[1]) => x -> x == 2 ? 1 : 0)
    tbl = JuliaDB.select(tbl, (:anml_key, Symbol(namesG11[1]), Symbol(namesG12[1]), Symbol(namesG22[1])))
    return tbl
end

function calc_pvalue_norm(t)
    dist = Normal()
    return 2*ccdf.(dist, abs.(t))
end

function summarise(pheno_header, basedir, markerlist; burnin=200)
    m = size(markerlist, 1)
    for pheno in pheno_header
        dir = basedir .* "/" .* pheno
        summary = [repeat([pheno],m) markerlist zeros(Float64, m, 4)]
        for i in 1:m
            file = dir .* "/" .* markerlist[i] .* "_1:" .* markerlist[i] .* ".txt"
            input = readdlm(file, '\t'; header = true)
            chain = input[1][201:end]# remove burnin
            summary[i,3] = mean(chain)
            summary[i,4] = std(chain)
            summary[i,5] = calc_pvalue_norm(summary[i,3] ./ summary[i,4])
            summary[i,6] = gewekediag(chain)[2]
        end
        writedlm(basedir .* "/summary." .* pheno .* ".spliceQTL.txt", summary, '\t')
    end
end

max_stars_repo_name: egmreynolds/MANA

max_issues_repo_name: egmreynolds/MANA

max_forks_repo_name: egmreynolds/MANA


text:
module IdealGas
include("Constants.jl")
include("Utils.jl")
using LightXML

abstract type ComponentDefinition end
abstract type ThermoData end

struct Gasphase <: ComponentDefinition
    species::Array{String,1}
end

struct NASAThermo <: ThermoData
    name::String
    phase::String
    composition::Dict{String,Int64}
    molWt::Float64
    ltl::Float64
    htl::Float64
    cmt::Float64
    ltp::Array{Float64,1}
    htp::Array{Float64,1}
end

#This array will be order as per ig.species
struct SpeciesThermoObj
    molwt::Array{Float64,1}
    thermo_all::Array{NASAThermo}
end
export SpeciesThermoObj

"""
    Tuple of element weights. This is used for the calculation of molecular weights
"""
elementWeight = (    H = 1.00794e-3,  He = 4.002602e-3,  Li = 6.941e-3,  Be = 9.012182e-3,  B = 1.811e-3,
                C = 12.011e-3,  N = 14.00674e-3,  O = 15.9994e-3,  F = 18.9984032e-3,  Ne = 20.1797e-3,
                Na = 22.98977e-3, Mg = 24.3050e-3, Al = 26.98154e-3, Si = 28.0855e-3, P = 30.97376e-3,
                S = 32.066e-3, Cl = 35.4527e-3, Ar = 39.948e-3, K = 39.0983e-3, Ca = 40.078e-3,
                Sc = 44.95591e-3, Ti = 47.88e-3, V = 50.9415e-3, Cr = 51.9961e-3, Mn = 54.9381e-3,
                Fe = 55.847e-3, Co = 58.9332e-3, Ni = 58.69e-3, Cu = 63.546e-3, Zn = 65.39e-3,
                Ga = 69.723e-3, Ge = 72.61e-3, As = 74.92159e-3, Se = 78.96e-3, Br = 79.904e-3,
                Kr = 83.80e-3,  Rb = 85.4678e-3, Sr = 87.62e-3, Y = 88.90585e-3, Zr = 91.224e-3,
                Nb = 92.90638e-3, Mo = 95.94e-3, Tc = 97.9072e-3, Ru = 101.07e-3, Rh = 102.9055e-3,
                Pd = 106.42e-3, Ag = 107.8682e-3,   Cd = 112.411e-3,    In = 114.82e-3, Sn = 118.710e-3,
                Sb = 121.75e-3, Te = 127.6e-3,  I = 126.90447e-3,   Xe = 131.29e-3,   Cs = 132.90543e-3,
                Ba = 137.327e-3, La = 138.9055e-3, Ce = 140.115e-3, Pr = 140.90765e-3, Nd = 144.24e-3,
                Pm = 144.9127e-3, Sm = 150.36e-3,   Eu = 151.965e-3,    Gd = 157.25e-3, Tb = 158.92534e-3,
                Dy = 162.50e-3, Ho = 164.93032e-3,  Er = 167.26e-3, Tm = 168.93421e-3,  Yb = 173.04e-3,
                Lu = 174.967e-3,    Hf = 178.49e-3, Ta = 180.9479e-3,   W = 183.85e-3, Re = 186.207e-3,
                Pt = 195.08e-3, Au = 196.96654e-3,  Hg = 200.59e-3, Tl = 204.3833e-3, Pb = 207.2e-3,
                Bi = 208.98037e-3,  Po = 208.9824e-3,  At = 209.9871e-3,  Rn = 222.0176e-3,
                Fr = 223.0197e-3, Ra = 226.0254e-3, Ac = 227.0279e-3, Th = 232.0381e-3,
                Pa = 231.03588e-3, U = 238.0508e-3, Np = 237.0482e-3, Pu = 244.0482e-3
            )

"""
get_element_weight(el)

function to return the element weights
#   Usage:
    get_element_weight(el)
-   'el::String' : Element name
"""
get_element_weight(el::String) = elementWeight[Symbol(titlecase(strip(el)))]

"""
Function to create thermo object. The function reads the therm.dat file and
parses the content based on the ideal gas object to create the thermo data object
The function returns SpeciesThermoObj
#   Usage:
    create_thermo() or
    create_thermo(file_path)
-   'ig::Gasphase' : Gasphase object obtained using the call create_gasphase
-   'file_path': relative path to the therm.dat file
"""
function create_thermo(species::Array{T,1}, thermo_file::AbstractString ) where T <: AbstractString
    #Create an array of thermo objects which stores the thermo data of all species
    species_thermo = Dict{String,NASAThermo}()
    molecular_weights_vector = Array{Float64,1}()

    open(thermo_file) do io
        lno = 0
        local  species_name, phase, mol_wt,LT,HT,CT
        comp = Dict{String,Int64}()
        local poly_coeff = Array{Float64,1}()
        while !eof(io)
            data_string = readline(io)
            if length(data_string) != 0
                if SubString(data_string,1,3) == "END"
                    break
                else
                    sp = uppercase(strip(SubString(data_string,1,18)))
                    if sp in species
                        lno = 1
                        species_name, phase,comp, mol_wt,LT,HT,CT = parase_thermo_species_data(data_string)
                        lno += 1
                        empty!(poly_coeff)
                    elseif lno > 1 && lno < 5
                        append!(poly_coeff,parse_thermo_polynomials(lno,data_string))
                        lno += 1
                        if lno > 4
                            nasa_thermo = NASAThermo(species_name,phase,comp ,mol_wt,LT,HT,CT,poly_coeff[8:14],poly_coeff[1:7])
                            species_thermo[nasa_thermo.name] = nasa_thermo
                        end
                    end
                end
            end
        end
    end
    species_thermo_array = Array{NASAThermo,1}()
    for sp in species
        push!(species_thermo_array,species_thermo[sp])
    end
    for td in species_thermo_array
        push!(molecular_weights_vector,td.molWt)
    end

    return SpeciesThermoObj(molecular_weights_vector, species_thermo_array)
end

"""
Function for parsing the first line of thermo data.
    It returns the species name, phase, molecular weight, high
    temperature, low temperature and Common temperature limits
    Not for external calls
#   Usage:
    parase_thermo_species_data(data_string)
"""
function parase_thermo_species_data(data_string::AbstractString)
    species_name = uppercase(strip(SubString(data_string,1,18)))
    a1 = split(strip(SubString(data_string,25,29)))
    a2 = split(strip(SubString(data_string,30,34)))
    a3 = split(strip(SubString(data_string,35,39)))
    a4 = split(strip(SubString(data_string,40,44)))
    elements_data = (a1,a2,a3,a4)
    mol_wt::Float64 = 0
    composition = Dict{String,Int64}()
    for el in elements_data
        if !isempty(el)
            mol_wt += parse(Float64,el[2])*get_element_weight(String(el[1]))
            composition[String(el[1])] = parse(Int64,el[2])
        end
    end
    phase = String(SubString(data_string,45,45))
    if length(phase)==0
        throw(error("Phase not specified for species $spname"))
    end
    LTs = String(SubString(data_string,46,57))
    if length(LTs) == 0
        throw(error("Low temperature limit not specified for species $spname"))
    end
    HTs = String(SubString(data_string,58,67))
    if length(HTs) == 0
        throw(error("High temperature limit not specified for species $spname"))
    end
    CTs = String(SubString(data_string,68,75))
    if length(CTs)==0
        throw(error("Common temperature limit not specified for species $spname"))
    end
    LT = parse(Float64, LTs)
    HT = parse(Float64, HTs)
    CT = parse(Float64, CTs)
    return (species_name,phase, composition ,mol_wt,LT,HT,CT)
end

"""
This function will parse the polynomials
    Not for external calls
"""
function parse_thermo_polynomials(lno::Int64, thdata::AbstractString)
    poly_data = Array{Float64,1}()
    if lno == 2
        push!(poly_data,parse(Float64,String(SubString(thdata,1,15))))
        push!(poly_data,parse(Float64,String(SubString(thdata,16,30))))
        push!(poly_data,parse(Float64,String(SubString(thdata,31,45))))
        push!(poly_data,parse(Float64,String(SubString(thdata,46,60))))
        push!(poly_data,parse(Float64,String(SubString(thdata,61,75))))
    elseif lno == 3
        push!(poly_data,parse(Float64,String(SubString(thdata,1,15))))
        push!(poly_data,parse(Float64,String(SubString(thdata,16,30))))
        push!(poly_data,parse(Float64,String(SubString(thdata,31,45))))
        push!(poly_data,parse(Float64,String(SubString(thdata,46,60))))
        push!(poly_data,parse(Float64,String(SubString(thdata,61,75))))
    elseif lno == 4
        push!(poly_data,parse(Float64,String(SubString(thdata,1,15))))
        push!(poly_data,parse(Float64,String(SubString(thdata,16,30))))
        push!(poly_data,parse(Float64,String(SubString(thdata,31,45))))
        push!(poly_data,parse(Float64,String(SubString(thdata,46,60))))
    end
    return poly_data
end

"""
Calculate the specific heat of pure species J/mol-K
#   Usage-1:
    cp(thermo,T)
-   'thermo::NASAThermo': NASAThermo of the species
-   'T::Float64': Temperature in K at which the property is required
# Usage-2:
    cp(sp,T,thermo,ig)
-   'sp::String' : species name
-   'T::Float64' : Temperature K
-   'thermo::SpeciesThermoObj' : Structure of SpeciesThermoObj
-   'ig::Gasphase ' : Gasphase object
"""
function cp(thermo::NASAThermo, T::Float64)
    TVec = [T^i for i in 0:4]
    T < thermo.cmt ? sum(thermo.ltp[1:5] .* TVec)*R : sum(thermo.htp[1:5] .* TVec)*R
end
cp(sp::String,T::Float64,thermoObj::SpeciesThermoObj,species_list::Array{String,1}) = cp(thermoObj.thermo_all[get_index(sp,species_list)],T)

"""
Calculates the enthalpy of pure species J/mol
#   Usage-1:
    H(thermo,T)
-   'thermo::NASAThermo': NASAThermo of the species
-   'T::Float64': Temperature in K at which the property is required
# Usage-2:
    H(sp,T,thermo,ig)
-   'sp::String' : species name
-   'T::Float64' : Temperature K
-   'thermo::SpeciesThermoObj' : Structure of SpeciesThermoObj
-   'ig::Gasphase ' : Gasphase object
"""
function H(thermo::NASAThermo, T::Float64)
    TVec = [1, T/2.0, T^2/3.0, T^3/4.0, T^4/5.0, 1.0/T]
    T < thermo.cmt ? sum(thermo.ltp[1:6] .* TVec)*R*T : sum(thermo.htp[1:6] .* TVec)*R*T
end
H(sp::String,T::Float64,thermoObj::SpeciesThermoObj,species_list::Array{String,1}) = H(thermoObj.thermo_all[get_index(sp,species_list)],T)

"""
Calculates the entropy of pure species J/mol-K
#   Usage-1:
    S(thermo,T)
-   'thermo::NASAThermo': NASAThermo of the species
-   'T::Float64': Temperature in K at which the property is required
# Usage-2:
    S(sp,T,thermo,ig)
-   'sp::String' : species name
-   'T::Float64' : Temperature K
-   'thermo::SpeciesThermoObj' : Structure of SpeciesThermoObj
-   'ig::Gasphase ' : Gasphase object
"""
function S(thermo::NASAThermo, T::Float64)
    TVec = [log(T), T, T^2/2.0, T^3/3.0, T^4/4.0]
    if T < thermo.cmt
        return (sum(thermo.ltp[1:5] .* TVec) + thermo.ltp[7])*R
    else
        return (sum(thermo.htp[1:5] .* TVec) + thermo.htp[7])*R
    end
end
S(sp::String,T::Float64,thermoObj::SpeciesThermoObj,species_list::Array{String,1}) = S(thermoObj.thermo_all[get_index(sp,species_list)],T)

"""
Calculates the specific heat of all species in J/mol-K
#   Usage
    cp_all(td,T)
-   'thermoObj::SpeciesThermoObj' : Structure of SpeciesThermoObj
-   'T::Float64' : Temperature in K at which the property is rquired
"""
function cp_all(thermoObj::SpeciesThermoObj,T::Float64)
    return map(x->cp(x,T),thermoObj.thermo_all)
end

"""
Calculates the enthalpy of all species in J/mol
#   Usage
    H_all(td,T)
-   'thermoObj::SpeciesThermoObj' : Structure of SpeciesThermoObj
-   'T::Float64' : Temperature in K at which the property is rquired
"""
function H_all(thermoObj::SpeciesThermoObj,T::Float64)
    return map(x->H(x,T),thermoObj.thermo_all)
end

"""
Calculates the entropy of all species in J/mol-K
#   Usage
    S_all(td,T)
-   'thermoObj::SpeciesThermoObj' : Structure of SpeciesThermoObj
-   'T::Float64' : Temperature in K at which the property is rquired
"""
function S_all(thermoObj::SpeciesThermoObj,T::Float64)
    return map(x->S(x,T),thermoObj.thermo_all)
end

"""
Calculates the enthalpy of a mixture J/mol
#   Usage
    Hmix(td,T,mlf)
-   'thermoObj::SpeciesThermoObj' : Structure of SpeciesThermoObj
-   'T::Float64' : Temperature in K
-   'mlf::Array{Float64,1}' : species mole fractions
"""
function Hmix(thermoObj::SpeciesThermoObj,T::Float64,mlf::Array{Float64,1})
    all_species_h = H_all(thermoObj,T)
    return sum(all_species_h .* mlf)
end

"""
Calculates the specific heat of a mixture in J/mol-K
#   Usage
    cpmix(td,T,mlf)
-   'thermoObj::SpeciesThermoObj' : Structure of SpeciesThermoObj
-   'T::Float64' : Temperature in K
-   'mlf::Array{Float64,1}' : species mole fractions
"""
function cpmix(thermoObj::SpeciesThermoObj,T::Float64,mlf::Array{Float64,1})
    all_species_cp = S_all(thermoObj,T)
    return sum(all_species_cp .* mlf)
end

"""
Calculates the entropy of a muxture in J/mol-K
#   Usage
    Smix(thermoObj,T,p,mlf)
-   'thermoObj::SpeciesThermoObj' : Structure of SpeciesThermoObj
-   'T::Float64' : Temperature in K
-   'p::Float64' : total pressure Pa
-   'mlf::Array{Float64,1} ': mole fractions
"""
function Smix(thermoObj::SpeciesThermoObj,T::Float64,p::Float64,mlf::Array{Float64,1})
    all_species_s = S_all(thermoObj,T)
    xp = mlf .* (p/p_std)
    #logterm = [x <= 0 ? 0 : log(x) for x in xp]
    logterm = collect(map(x->x <= 0 ? 0 : log(x) ,xp))
    return sum((all_species_s - (logterm .* R)) .* mlf)
end

"""
Calculates the Gibbs free energy of a muxture in J/mol
#   Usage
    Gmix(thermoObj,T,p,mlf)
-   'thermoObj::SpeciesThermoObj' : Structure of SpeciesThermoObj
-   'T::Float64' : Temperature in K
-   'p::Float64' : total pressure Pa
-   'mlf::Array{Float64,1} ': mole fractions
"""
function  Gmix(thermoObj::SpeciesThermoObj, T::Float64,p::Float64,mlf::Array{Float64,1})
    hmix = Hmix(thermoObj,T,mlf)
    smix = Smix(thermoObj,T,p,mlf)
    return hmix - T*smix
end

#end of module IdealGas
end


text:
"""
Main structure for holding attributes in e.g. plots!
"""
struct ObservableAttributes
    # name, for better error messages!
    name::String
    # We dont have one node per value anymore, but instead one node
    # that gets triggered on any setindex!, or whenever an input attribute node changes
    # This makes it easier to layer Observable{ObservableAttributes}()
    on_change::Observable{Pair{Symbol, Any}}
    # The supported fields, so we can throw an error, whenever fields are not supported
    supported_fields::Set{Symbol}
    # The attributes given at construction time, taking the highest priority
    from_user::Dict{Symbol, Any}
    # attributes filled in by e.g. a theme, or other processes in the pipeline
    from_theme::Dict{Symbol, Any}
    # The "global" default values for this specific attribute instance
    # Will be immutabe (maybe not a Dict then?) and shared between all similar objects
    default_values::Dict{Symbol, Any}

    function ObservableAttributes(name::String, default_values::Dict{Symbol, Any}, from_user::Dict{Symbol, Any})
        on_change = Observable{Pair{Symbol, Any}}()
        supported_fields = Set(keys(default_values))
        return new(name, on_change, supported_fields, from_user,
                   Dict{Symbol, Any}(), default_values)
    end
end

function ObservableAttributes(;kw...)
    return ObservableAttributes(kw)
end

value_convert2(x::Observables.AbstractObservable) = x[]
value_convert2(@nospecialize(x)) = x
function value_convert2(x::NTuple{N, Union{Any, Observables.AbstractObservable}}) where N
    return to_value.(x)
end
value_convert2(x::NamedTuple) = ObservableAttributes(x)

function ObservableAttributes(kw)
    defaults = Dict{Symbol, Any}()
    for (k, v) in pairs(kw)
        xx = value_convert2(v)
        defaults[k] = xx
        @assert !(xx isa Observable)
    end
    attributes = ObservableAttributes("", defaults, Dict{Symbol, Any}())
    for (k, v) in pairs(kw)
        if v isa Observables.AbstractObservable
            on(v) do value
                setproperty!(attributes, k, value)
            end
        end
    end
    return attributes
end

function Base.getindex(attributes::ObservableAttributes, field::Symbol)
    return getproperty(attributes, field)
end

function Base.getindex(attributes::ObservableAttributes, field::Symbol, b::Symbol)
    x = getproperty(attributes, b)
    return getproperty(x, field)
end

function Base.propertynames(attributes::ObservableAttributes)
    return getfield(attributes, :supported_fields)
end

function on_change(x::ObservableAttributes)
    return getfield(x, :on_change)
end

function on_change(f, x::ObservableAttributes)
    return on(f, on_change)
end

"""
attributes.attribute returns an observable.
Since we don't actually convert all values to Observables anymore,
we'll need to create new Observables on getproperty.
With OnFieldUpdate, we can do that lazily, store them in `listeners(onchange(obs))`,
and only create a new one for fields that aren't in listeners yet.
"""
struct OnFieldUpdate
    field::Symbol
    observable::Observable
end

function (of::OnFieldUpdate)(field_value::Pair{Symbol, <: Any})
    if field_value[1] === of.field
        of.observable[] = field_value[2]
    end
    return
end

function OnFieldUpdate(attributes::ObservableAttributes, field::Symbol)
    # We lazily store observables on field updates in the listeners
    # If we already have it in there, we just return the one we have.
    onchange = on_change(attributes)
    for listener in Observables.listeners(onchange)
        if listener isa OnFieldUpdate && listener.field === field
            return listener
        end
    end
    # we haven't found a listener, so we create a new one!
    result = Observable{Any}(get_value(attributes, field))
    of = OnFieldUpdate(field, result)
    on(of, onchange)
    return of
end

"""
    get_value(attributes::ObservableAttributes, field::Symbol)
Gets the value for `field`.
The values are looked up in the following order:
    1) user given at creation time
    2) theme given
    3) global defaults
"""
function get_value(attributes::ObservableAttributes, field::Symbol)
    name = getfield(attributes, :name)
    if field in propertynames(attributes)
        # The priority is:
        # User given
        from_user = getfield(attributes, :from_user)
        haskey(from_user, field) && return from_user[field]
        # ObservableAttributes given
        from_theme = getfield(attributes, :from_theme)
        haskey(from_theme, field) && return from_theme[field]
        # Construction defaults
        default_values = getfield(attributes, :default_values)
        haskey(default_values, field) && return default_values[field]
        error("Incorrectly constructed ObservableAttributes ($(name))! No value found for $(field)")
    else
        from_user = getfield(attributes, :from_user)
        default_values = getfield(attributes, :default_values)
        pn = propertynames(attributes)
        error("Field $(field) not in $(pn) $(from_user) $(default_values)!")
    end
end

function Base.getproperty(attributes::ObservableAttributes, field::Symbol)
    val = get_value(attributes, field)
    val isa ObservableAttributes && return val
    of = OnFieldUpdate(attributes, field)
    return of.observable
end

function Base.setproperty!(attributes::ObservableAttributes, field::Symbol, value)
    name = getfield(attributes, :name)
    # we always set the users data, since setting this is done by the user!
    from_user = getfield(attributes, :from_user)
    from_user[field] = value
    on_change = getfield(attributes, :on_change)
    # trigger change!
    on_change[] = field => value
    if !(field in propertynames(attributes))
        # TODO, dont let anyone change propertynames!
        push!(propertynames(attributes), field)
        default_values = getfield(attributes, :default_values)
        default_values[field] = value
    end
    return value
end

function Base.setproperty!(attributes::ObservableAttributes, field::Symbol, value::Observable)
    setproperty!(attributes, field, value[])
    on(value) do new_value
        setproperty!(attributes, field, value[])
    end
    return value
end

function Base.setindex!(attributes::ObservableAttributes, value, field::Symbol)
    setproperty!(attributes, field, value)
end

Base.broadcastable(x::ObservableAttributes) = Ref(x)

#######
## Dict interface
Base.keys(x::ObservableAttributes) = propertynames(x)
Base.values(x::ObservableAttributes) = (getproperty(x, field) for field in propertynames(x))
function Base.pop!(x::ObservableAttributes, field::Symbol)
    value = getproperty(x, field)
    delete!(x, field)
    return value
end

Base.haskey(x::ObservableAttributes, key::Symbol) = key in keys(x)

Base.filter(f, x::ObservableAttributes) = ObservableAttributes(filter(f, attributes(x)))
Base.empty!(x::ObservableAttributes) = (empty!(attributes(x)); x)
Base.length(x::ObservableAttributes) = length(attributes(x))

function Base.delete!(attributes::ObservableAttributes, field::Symbol)
    # The priority is:
    # User given
    from_user = getfield(attributes, :from_user)
    haskey(from_user, field) && delete!(from_user, field)
    # ObservableAttributes given
    from_theme = getfield(attributes, :from_theme)
    haskey(from_theme, field) && delete!(from_theme, field)
    # Construction defaults
    default_values = getfield(attributes, :default_values)
    haskey(default_values, field) && delete!(default_values, field)
    delete!(propertynames(attributes), field)
    return
end

function Base.iterate(x::ObservableAttributes, state...)
    s = iterate(keys(x), state...)
    s === nothing && return nothing
    return (s[1] => x[s[1]], s[2])
end

function Base.copy(attributes::ObservableAttributes)
    return ObservableAttributes(attributes)
end

Base.merge(target::ObservableAttributes, args::ObservableAttributes...) = merge!(copy(target), args...)

function merge_attributes!(input::ObservableAttributes, theme::ObservableAttributes)
    for (key, value) in theme
        if !haskey(input, key)
            input[key] = value
        else
            current_value = input[key]
            if value isa ObservableAttributes && current_value isa ObservableAttributes
                # if nested attribute, we merge recursively
                merge_attributes!(current_value, value)
            elseif value isa ObservableAttributes || current_value isa ObservableAttributes
                error("""
                Type missmatch while merging plot attributes with theme for key: $(key).
                Found $(to_value(value)) in theme, while attributes contains: $(current_value)
                """)
            else
                # we're good! input already has a value, can ignore theme
            end
        end
    end
    return input
end

function Base.merge!(target::ObservableAttributes, args::ObservableAttributes...)
    for elem in args
        merge_attributes!(target, elem)
    end
    return target
end

function Base.show(io::IO,::MIME"text/plain", attr::ObservableAttributes)
    d = Dict()
    for (k, v) in attr
        d[k] = to_value(v)
    end
    show(IOContext(io, :limit => false), MIME"text/plain"(), d)
end

Base.show(io::IO, attr::ObservableAttributes) = show(io, MIME"text/plain"(), attr)

"""
    get_attribute(dict::ObservableAttributes, key::Key)
Gets the attribute at `key`, converts it and extracts the value
"""
function get_attribute(dict::ObservableAttributes, key::Symbol)
    return convert_attribute(get_value(dict, key), Key{key}())
end

"""
    get_attribute(dict::ObservableAttributes, key::Key)
Gets the attribute at `key` as a converted signal
"""
function get_lifted_attribute(dict::ObservableAttributes, key::Symbol)
    return lift(x-> convert_attribute(x, Key{key}()), getproperty(dict, key))
end


text:
# # Control

#-

# A simple control problem on a system usually involves a variable $x(t)$
# that denotes the state of the system over time, and a variable $u(t)$ that
# denotes the input into the system over time. Linear constraints are used to
# capture the evolution of the system over time:
#
# $$
# x(t) = Ax(t - 1) + Bu(t), \ \text{for} \ t = 1,\ldots, T,
# $$
#
# where the numerical matrices $A$ and $B$ are called the dynamics and input matrices,
# respectively.
#
# The goal of the control problem is to find a sequence of inputs
# $u(t)$ that will allow the state $x(t)$ to achieve specified values
# at certain times. For example, we can specify initial and final states of the system:
#
# $$
# \begin{aligned}
#   x(0) &= x_i \\
#   x(T) &= x_f
# \end{aligned}
# $$
#
# Additional states between the initial and final states can also be specified. These
# are known as waypoint constraints. Often, the input and state of the system will
# have physical meaning, so we often want to find a sequence inputs that also
# minimizes a least squares objective like the following:
#
# $$
#   \sum_{t = 0}^T \|Fx(t)\|^2_2 + \sum_{t = 1}^T\|Gu(t)\|^2_2,
# $$
#
# where $F$ and $G$ are numerical matrices.
#
# We'll now apply the basic format of the control problem to an example of controlling
# the motion of an object in a fluid over $T$ intervals, each of $h$ seconds.
# The state of the system at time interval $t$ will be given by the position and the velocity of the
# object, denoted $p(t)$ and $v(t)$, while the input will be forces
# applied to the object, denoted by $f(t)$.
# By the basic laws of physics, the relationship between force, velocity, and position
# must satisfy:
#
# $$
#   \begin{aligned}
#     p(t+1) &= p(t) + h v(t) \\
#     v(t+1) &= v(t) + h a(t)
#   \end{aligned}.
# $$
#
# Here, $a(t)$ denotes the acceleration at time $t$, for which we we use
# $a(t) = f(t) / m + g - d v(t)$,
# where $m$, $d$, $g$ are constants for the mass of the object, the drag
# coefficient of the fluid, and the acceleration from gravity, respectively.
#
# Additionally, we have our initial/final position/velocity conditions:
#
# $$
#   \begin{aligned}
#     p(1) &= p_i\\
#     v(1) &= v_i\\
#     p(T+1) &= p_f\\
#     v(T+1) &= 0
#   \end{aligned}
# $$
#
# One reasonable objective to minimize would be
#
# $$
#   \text{objective} = \mu \sum_{t = 1}^{T+1} (v(t))^2 + \sum_{t = 1}^T (f(t))^2
# $$
#
# We would like to keep both the forces small to perhaps save fuel, and keep
# the velocities small for safety concerns.
# Here $\mu$ serves as a parameter to control which part of the objective we
# deem more important, keeping the velocity small or keeping the force small.
#
# The following code builds and solves our control example:

using Convex, SCS, Plots

## Some constraints on our motion
## The object should start from the origin, and end at rest
initial_velocity = [-20; 100]
final_position = [100; 100]

T = 100 # The number of timesteps
h = 0.1 # The time between time intervals
mass = 1 # Mass of object
drag = 0.1 # Drag on object
g = [0, -9.8] # Gravity on object

## Declare the variables we need
position = Variable(2, T)
velocity = Variable(2, T)
force = Variable(2, T - 1)

## Create a problem instance
mu = 1

## Add constraints on our variables
constraints = Constraint[
    position[:, i+1] == position[:, i] + h * velocity[:, i] for i in 1:T-1
]

for i in 1:T-1
    acceleration = force[:, i] / mass + g - drag * velocity[:, i]
    push!(constraints, velocity[:, i+1] == velocity[:, i] + h * acceleration)
end

## Add position constraints
push!(constraints, position[:, 1] == 0)
push!(constraints, position[:, T] == final_position)

## Add velocity constraints
push!(constraints, velocity[:, 1] == initial_velocity)
push!(constraints, velocity[:, T] == 0)

## Solve the problem
problem = minimize(sumsquares(force), constraints)
solve!(problem, SCS.Optimizer; silent_solver = true)

# We can plot the trajectory taken by the object.

pos = evaluate(position)
plot([pos[1, 1]], [pos[2, 1]], st = :scatter, label = "initial point")
plot!([pos[1, T]], [pos[2, T]], st = :scatter, label = "final point")
plot!(pos[1, :], pos[2, :], label = "trajectory")

# We can also see how the magnitude of the force changes over time.

plot(vec(sum(evaluate(force) .^ 2, dims = 1)), label = "force (magnitude)")


text:
@testset "De Bruijn Neighbors" begin
    @test collect(neighbors(DNAMer("ACG")))  == map(DNAMer, ["CGA",  "CGC",  "CGG",  "CGT"])
    @test collect(neighbors(DNAMer("GGGG"))) == map(DNAMer, ["GGGA", "GGGC", "GGGG", "GGGT"])
    @test collect(neighbors(RNAMer("ACG")))  == map(RNAMer, ["CGA",  "CGC",  "CGG",  "CGU"])
    @test collect(neighbors(RNAMer("GGGG"))) == map(RNAMer, ["GGGA", "GGGC", "GGGG", "GGGU"])

    @test collect(neighbors(BigDNAMer("ACG")))  == map(BigDNAMer, ["CGA",  "CGC",  "CGG",  "CGT"])
    @test collect(neighbors(BigDNAMer("GGGG"))) == map(BigDNAMer, ["GGGA", "GGGC", "GGGG", "GGGT"])
    @test collect(neighbors(BigRNAMer("ACG")))  == map(BigRNAMer, ["CGA",  "CGC",  "CGG",  "CGU"])
    @test collect(neighbors(BigRNAMer("GGGG"))) == map(BigRNAMer, ["GGGA", "GGGC", "GGGG", "GGGU"])
end