text:
#=
This code reads a FOO.jld generated by isingED.jl, generates N random
spin chain measurements in the computational basis, and saves them to
a text file FOO.csv.
=#

using ArgParse, HDF5, JLD

function parse_commandline()
    s = ArgParseSettings()
    @add_arg_table s begin
      "--infilename", "-i"
        help =  "Path to the .jld where the wavefunction is stored."
        arg_type = AbstractString
        required = true
      "-N"
        help = "Number of configurations to generate."
        arg_type = Int
        default = 20
      "--overwrite"
        help = "If infilename.out already exists and this flag is set it will
                be overwritten (by default it will be appended to)."
        action = :store_true

    end
    return parse_args(s)
end

"""
Generates a sample spin chain measurement from the wavefunction wf.
The idea is to simulate measuring
each site in sequence. After each measurement, the wavefunction is projected
onto the appropriate Hilbert subspace according to the Born rule.
"""
function measure(wf::AbstractArray)
  #rng = MersenneTwister()
  Nsites = convert(Int64, log2(length(wf)))
  measurement = zeros(Int, Nsites)

  #shape = tuple(collect(repeated(2,Nsites))...)
  for n = 1:Nsites
    shape=(2^(Nsites-n), 2)
    wf = reshape(wf, shape)
    #println(abs2(wf[:,1]))
    P0 = sum(abs2.(wf[:,1]))
    @assert P0>=0 && P0<(1.0+1e-12) "P0 was $P0 which is not a probability"
    roll = rand()
    if roll < P0
      outcome = 0
      wf = wf[:,1] / sqrt(P0) #this is the Born rule
    else
      outcome = 1
      wf = wf[:,2] / sqrt(1-P0)
    end
    measurement[Nsites-n+1] = outcome
  end
  return measurement
end

"""
***UNIMPLEMENTED***
Transforms the wavefunction so that computational-basis measurements upon
it are equivalent to measurements of the given string of Pauli operators
(0 -> I, 1->X, 2->Y, 3->Z).

For now this is just a null-op.
"""
function transformtopauli(wavefunction::AbstractArray, paulis::AbstractVector)
    return wavefunction
end

function writetofile(f, array)
    N = length(array)
    for i = 1:N
      write(f,  string(array[i]))
      if i!=N
        write(f, ",")
      else
        write(f, ";")
      end
    end
end

"""
Generates N measurements from the wavefunction and saves them to disk.
The overwrite flag determines whether the output file is overwritten or
appended to.
"""
function savemeasurements(wavefunction::AbstractArray,
                          outfilename::AbstractString, N::Int, overwrite::Bool)
  if overwrite
    mode = "w"
  else
    mode = "a"
  end
  open(outfilename, mode) do f
    for n = 1:N
      Nsites = convert(Int64, log2(length(wavefunction)))
      paulis = [rand(0:3) for j in 1:Nsites]
      pauliwf = transformtopauli(wavefunction, paulis)
      measurement = measure(pauliwf)
      println("Pauli ops: ", paulis)
      writetofile(f, paulis)
      println("Outcome: ", measurement)
      writetofile(f, measurement)
      write(f, "\n")
    end
  end
end

function main()
  parsed_args = parse_commandline()
  println("Hi!! Generating field configurations from:")
  for (arg, val) in parsed_args
    println("   $arg    =>    $val")
  end
  infilename = parsed_args["infilename"]
  wavefunction = load(infilename, "wavefunction")
  outfilename = infilename[1:end-3]*string("csv")
  savemeasurements(wavefunction, outfilename, parsed_args["N"],
                   parsed_args["overwrite"])
  println("Done.")

end
main()

max_stars_repo_name: alewis/isingED

max_issues_repo_name: alewis/isingED

max_forks_repo_name: alewis/isingED


text:
using Makie
using AbstractPlotting: textslider

function target_function(x, y, a, b, c)
        (a*tanh(x) + b*tanh(y) + c)^2
end

function plotsome(z::Array{Float32,2}, a, b, c)
    lspace = range(-1., stop=1., step=0.05)
    z = Float32[target_function(x, y, a, b, c) for x in lspace, y in lspace]
    z
end

function main()
    lspace = range(-1., stop=1., step=0.05)
    N = length(collect(lspace))
    s_a, obs_a = textslider(lspace, "a", start=0.1)
    s_b, obs_b = textslider(lspace, "b", start=0.1)
    s_c, obs_c = textslider(lspace, "c", start=0.1)
    args_n = (obs_a, obs_b, obs_c)
    update_rule = plotsome(zeros(Float32, (N, N)), to_value.(args_n)...)
    positions = lift(plotsome, Node(update_rule), args_n...)

    sur = surface(lspace, lspace, positions)
    axis = sur[Axis]
    axis[:names][:axisnames] = ("x", "y", "(a*tanh(x)+b*tanh(y)+c)^2")

    RecordEvents(
        hbox(
            sur,
            vbox(s_a, s_b, s_c),
            parent=Scene()
        ),
        "output"
    )

end

main()

ext: jl

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
using Revise
using HybridFunnelGraphs

graph = create_funnel_graph("test/pddl/domain.pddl", "test/pddl/problem.pddl"; max_levels=60)
graph.num_levels


text:
function square_of_sum(n)
    # https://docs.julialang.org/en/v1/manual/control-flow/#man-conditional-evaluation
    # https://docs.julialang.org/en/v1/base/collections/#Base.collect-Tuple{Any}
    # https://www.simonwenkel.com/2019/02/03/project-euler-problem-0006.html
    return n > 1 ? sum(1:n)^2 : n
end

function sum_of_squares(n)
    # https://docs.julialang.org/en/v1/base/collections/#Base.mapreduce-Tuple{Any,%20Any,%20Any}
    # https://github.com/programacaodinamica/exercism-julia/blob/main/difference-of-squares/difference-of-squares.jl
    return n > 1 ? mapreduce(x -> x^2, +, 1:n) : n
end

function difference(n)
    return square_of_sum(n) - sum_of_squares(n)
end

max_stars_repo_name: joaopalmeiro/learning-julia

max_issues_repo_name: joaopalmeiro/learning-julia

max_forks_repo_name: joaopalmeiro/learning-julia


text:
## A SortedDict is a wrapper around balancedTree with
## methods similiar to those of Julia container Dict.

type SortedDict{K, D, Ord <: Ordering} <: Associative{K,D}
    bt::BalancedTree23{K,D,Ord}
end

typealias SDSemiToken IntSemiToken

typealias SDToken{K,D, Ord <: Ordering} Token{SortedDict{K,D,Ord}, SDSemiToken}

## This constructor takes an ordering object which defaults
## to Forward

function SortedDict{K,D, Ord <: Ordering}(d::Associative{K,D}, o::Ord=Forward)
    bt1 = BalancedTree23{K,D,Ord}(o)
    for pr in d
        insert!(bt1, pr[1], pr[2], false)
    end
    SortedDict(bt1)
end

## This function implements m[k]; it returns the
## data item associated with key k.

function getindex{K,D, Ord <: Ordering}(m::SortedDict{K,D,Ord}, k_)
    i, exactfound = findkey(m.bt, convert(K,k_))
    !exactfound && throw(KeyError(k_))
    return m.bt.data[i].d
end

## This function implements m[k]=d; it sets the
## data item associated with key k equal to d.

function setindex!{K, D, Ord <: Ordering}(m::SortedDict{K,D,Ord}, d_, k_)
    insert!(m.bt, convert(K,k_), convert(D,d_), false)
    m
end

## Functions setindex! and getindex for semitokens:

function getindex{K, D, Ord <: Ordering}(m::SortedDict{K,D,Ord}, i::SDSemiToken)
    addr = i.address
    has_data(SDToken{K,D,Ord}(m,i))
    return m.bt.data[addr].d
end

function setindex!{K,D,Ord <: Ordering}(m::SortedDict{K,D,Ord},
                                        d_,
                                        i::SDSemiToken)
    addr = i.address
    has_data(SDToken{K,D,Ord}(m,i))
    m.bt.data[addr] = KDRec{K,D}(m.bt.data[addr].parent,
                                 m.bt.data[addr].k,
                                 convert(D,d_))
    m
end

sdtoken_construct{K, D, Ord <: Ordering}(m::SortedDict{K,D,Ord},int1::Int) =
    SDToken{K,D,Ord}(m, SDSemiToken(int1))

## This function looks up a key in the tree;
## if not found, then it returns a marker for the
## end of the tree.

function find{K,D,Ord <: Ordering}(m::SortedDict{K,D,Ord}, k_)
    ll, exactfound = findkey(m.bt, convert(K,k_))
    sdtoken_construct(m, exactfound? ll : 2)
end

## This function inserts an item into the tree.
## Unlike m[k]=d, it also returns a bool and a token.
## The bool is true if the inserted item is new.
## It is false if there was already an item
## with that key.
## The token points to the newly inserted item.

function insert!{K,D, Ord <: Ordering}(m::SortedDict{K,D,Ord}, k_, d_)
    b, i = insert!(m.bt, convert(K,k_), convert(D,d_), false)
    b, sdtoken_construct(m, i)
end

## delete! deletes an item given a token.

function delete!{K, D, Ord <: Ordering}(ii::SDToken{K,D,Ord})
    has_data(ii)
    delete!(ii.container.bt, ii.semitoken.address)
end

## Function startof returns the token that points
## to the first sorted order of the tree.  It returns
## the past-end token if the tree is empty.

startof(m::SortedDict) = sdtoken_construct(m, beginloc(m.bt))

## Function pastendtoken returns the otken past the end of the data.

pastendtoken(m::SortedDict) = sdtoken_construct(m,2)

## Function beforestarttoken returns the token before the start of the data.

beforestarttoken(m::SortedDict) = sdtoken_construct(m,1)

## Function advance takes a token and returns the
## next token in the sorted order.

function advance{K, D, Ord <: Ordering}(ii::SDToken{K,D,Ord})
    not_pastend(ii)
    sdtoken_construct(ii.container, nextloc0(ii.container.bt, ii.semitoken.address))
end

## Function regress takes a token and returns the
## previous token in the sorted order.

function regress{K, D, Ord <: Ordering}(ii::SDToken{K,D,Ord})
    not_beforestart(ii)
    sdtoken_construct(ii.container, prevloc0(ii.container.bt, ii.semitoken.address))
end

## Endof returns the token of the last item in the sorted order,
## or the before-start marker if the SortedDict is empty.

endof(m::SortedDict) = sdtoken_construct(m,endloc(m.bt))

## First and last return the first and last (key,data) pairs
## in the SortedDict.  It is an error to invoke them on an
## empty SortedDict.

function first(m::SortedDict)
    i = beginloc(m.bt)
    i == 2 && throw(BoundsError())
    return m.bt.data[i].k, m.bt.data[i].d
end

function last(m::SortedDict)
    i = endloc(m.bt)
    i == 1 && throw(BoundsError())
    return m.bt.data[i].k, m.bt.data[i].d
end

## Function deref(ii), where ii is a token, returns the
## (k,d) pair indexed by ii.

function deref{K, D, Ord <: Ordering}(ii::SDToken{K,D,Ord})
    has_data(ii)
    return ii.container.bt.data[ii.semitoken.address].k,
           ii.container.bt.data[ii.semitoken.address].d
end

## Function deref_key(ii), where ii is a token, returns the
## key indexed by ii.

function deref_key{K, D, Ord <: Ordering}(ii::SDToken{K,D,Ord})
    has_data(ii)
    return ii.container.bt.data[ii.semitoken.address].k
end

## Function deref_value(ii), where ii is a token, returns the
## value indexed by ii.

function deref_value{K, D, Ord <: Ordering}(ii::SDToken{K,D,Ord})
    has_data(ii)
    return ii.container.bt.data[ii.semitoken.address].d
end

## This function takes a key and returns the token
## of the first item in the tree that is >= the given
## key in the sorted order.  It returns the past-end marker
## if there is none.

function searchsortedfirst{K, D, Ord <: Ordering}(m::SortedDict{K,D,Ord}, k_)
    i, exactfound = findkey(m.bt, convert(K,k_))
    sdtoken_construct(m, exactfound? i : nextloc0(m.bt, i))
end

## This function takes a key and returns a token
## to the first item in the tree that is > the given
## key in the sorted order.  It returns the past-end marker
## if there is none.

function searchsortedafter{K, D, Ord <: Ordering}(m::SortedDict{K,D,Ord}, k_)
    i, exactfound = findkey(m.bt, convert(K,k_))
    sdtoken_construct(m, nextloc0(m.bt, i))
end

## This function takes a key and returns a token
## to the last item in the tree that is <= the given
## key in the sorted order.  It returns the before-start marker
## if there is none.

function searchsortedlast{K,D,Ord <: Ordering}(m::SortedDict{K,D,Ord}, k_)
    i, exactfound = findkey(m.bt, convert(K,k_))
    sdtoken_construct(m, i)
end

empty!(m::SortedDict) =  empty!(m.bt)
length(m::SortedDict) = length(m.bt.data) - length(m.bt.freedatainds) - 2
isempty(m::SortedDict) = length(m) == 0

## The next block of code handles iterations (for loops, i.e., start/next/done
## functions).  There are multiple ways to iterate over a loop; see the documentation.

## All the loops maintain a state which is an object of the
## following type.
immutable SDIterationState{K, D, Ord <: Ordering}
    m::SortedDict{K,D,Ord}
    next::Int
    final::Int
end

## This holds an object describing an exclude-last
## iteration.

immutable ExcludeLast{K, D, Ord <: Ordering}
    m::SortedDict{K, D, Ord}
    first::Int
    pastlast::Int
end

## This holds an object describing an include-last
## (i.e., colon operator) iteration.

immutable IncludeLast{K, D, Ord <: Ordering}
    m::SortedDict{K, D, Ord}
    first::Int
    last::Int
end

## The basic iterations are either over the whole sorted dict, an
## exclude-last object or include-last object.

typealias SDIterableTypesBase Union(SortedDict,ExcludeLast,IncludeLast)

## The compound iterations are obtained by applying keys(..) or tokens(..)
## to the basic iterations.  Furthermore, tokens(..) can be applied
## to either a basic iteration or a keys/values iteration.

immutable SDKeyIteration{T <: SDIterableTypesBase}
    base::T
end

immutable SDValIteration{T <: SDIterableTypesBase}
    base::T
end

immutable SDTokenIteration{T <: SDIterableTypesBase}
    base::T
end

immutable SDTokenKeyIteration{T <: SDIterableTypesBase}
    base::T
end

immutable SDTokenValIteration{T <: SDIterableTypesBase}
    base::T
end

typealias SDCompoundIterable Union(SDKeyIteration,
                                   SDValIteration, SDTokenIteration,
                                   SDTokenKeyIteration, SDTokenValIteration)

typealias SDAllIterable Union(SDIterableTypesBase, SDCompoundIterable)

done(::SDAllIterable, state::SDIterationState) = state.next == state.final

function nexthelper(state::SDIterationState)
    m = state.m
    sn = state.next
    (sn < 3 || !(sn in m.bt.useddatacells)) && throw(BoundsError())
    m.bt.data[sn], sdtoken_construct(m,sn), SDIterationState(m, nextloc0(m.bt, sn), state.final)
end

## The next function returns different objects depending on whether
## it is a basic iteration, a key iteration, a values iterations,
## a tokens/basic iteration, a tokens/key iteration, or tokens/values
## iteration.

function next(::SDIterableTypesBase, state::SDIterationState)
    dt, t, ni = nexthelper(state)
    (dt.k, dt.d), ni
end

function next(::SDKeyIteration, state::SDIterationState)
    dt, t, ni = nexthelper(state)
    dt.k, ni
end

function next(::SDValIteration, state::SDIterationState)
    dt, t, ni = nexthelper(state)
    dt.d, ni
end

function next(::SDTokenIteration, state::SDIterationState)
    dt, t, ni = nexthelper(state)
    (t, (dt.k, dt.d)), ni
end

function next(::SDTokenKeyIteration, state::SDIterationState)
    dt, t, ni = nexthelper(state)
    (t, dt.k), ni
end

function next(::SDTokenValIteration, state::SDIterationState)
    dt, t, ni = nexthelper(state)
    (t, dt.d), ni
end

function excludelast{K, D, Ord <: Ordering}(i1::SDToken{K,D,Ord}, i2::SDToken{K,D,Ord})
    !(i1.container === i2.container) &&
        throw(ArgumentError("SDToken range constructor requires tokens for the same container"))
    ExcludeLast(i1.container, i1.semitoken.address, i2.semitoken.address)
end

function colon{K, D, Ord <: Ordering}(i1::SDToken{K,D,Ord}, i2::SDToken{K,D,Ord})
    !(i1.container === i2.container) &&
        throw(ArgumentError("SDToken range constructor requires tokens for the same container"))
    IncludeLast(i1.container, i1.semitoken.address, i2.semitoken.address)
end

# need the next two to resolve ambiguity with keys(Associative)
keys(m::SortedDict) = SDKeyIteration(m)
values(m::SortedDict) = SDValIteration(m)
keys{T <: SDIterableTypesBase}(ba::T) = SDKeyIteration(ba)
values{T <: SDIterableTypesBase}(ba::T) = SDValIteration(ba)
tokens{T <: SDIterableTypesBase}(ba::T) = SDTokenIteration(ba)
tokens{T <: SDIterableTypesBase}(ki::SDKeyIteration{T}) = SDTokenKeyIteration(ki.base)
tokens{T <: SDIterableTypesBase}(vi::SDValIteration{T}) = SDTokenValIteration(vi.base)

start(m::SortedDict) = SDIterationState(m, nextloc0(m.bt,1), 2)

function start(e::ExcludeLast)
    (!(e.first in e.m.bt.useddatacells) || e.first == 1 ||
        !(e.pastlast in e.m.bt.useddatacells)) &&
        throw(BoundsError())
    if compareInd(e.m.bt, e.first, e.pastlast) < 0
        return SDIterationState(e.m, e.first, e.pastlast)
    else
        return SDIterationState(e.m, 2, 2)
    end
end

function start(e::IncludeLast)
    (!(e.first in e.m.bt.useddatacells) || e.first == 1 ||
        !(e.last in e.m.bt.useddatacells) || e.last == 2) &&
        throw(BoundsError())
    if compareInd(e.m.bt, e.first, e.last) <= 0
        return SDIterationState(e.m, e.first, nextloc0(e.m.bt, e.last))
    else
        return SDIterationState(e.m, 2, 2)
    end
end

start(e::SDCompoundIterable) = start(e.base)

function isless{K, D, Ord <: Ordering}(s::SDToken{K,D,Ord}, t::SDToken{K,D,Ord})
    !(s.container === t.container) &&
        throw(ArgumentError("SDToken isless requires tokens for the same container"))
    return compareInd(s.container.bt,
                      s.semitoken.address,
                      t.semitoken.address) < 0
end

function isequal{K, D, Ord <: Ordering}(s::SDToken{K,D,Ord}, t::SDToken{K,D,Ord})
    !(s.container === t.container) &&
        throw(ArgumentError("SDToken isequal requires tokens for the same container"))
    return s.semitoken.address == t.semitoken.address
end

function in{K,D,Ord <: Ordering}(pr::(@compat Tuple{Any,Any}), m::SortedDict{K,D,Ord})
    i, exactfound = findkey(m.bt,convert(K,pr[1]))
    return exactfound && isequal(m.bt.data[i].d,convert(D,pr[2]))
end

eltype{K,D,Ord <: Ordering}(m::SortedDict{K,D,Ord}) = (K,D)

orderobject(m::SortedDict) = m.bt.ord

function haskey{K,D,Ord <: Ordering}(m::SortedDict{K,D,Ord}, k_)
    i, exactfound = findkey(m.bt,convert(K,k_))
    exactfound
end

function get{K,D,Ord <: Ordering}(m::SortedDict{K,D,Ord}, k_, default_)
    i, exactfound = findkey(m.bt, convert(K,k_))
   return exactfound? m.bt.data[i].d : convert(D,default_)
end

function get!{K,D,Ord <: Ordering}(m::SortedDict{K,D,Ord}, k_, default_)
    k = convert(K,k_)
    i, exactfound = findkey(m.bt, k)
    if exactfound
        return m.bt.data[i].d
    else
        default = convert(D,default_)
        insert!(m.bt,k, default, false)
        return default
    end
end

function getkey{K,D,Ord <: Ordering}(m::SortedDict{K,D,Ord}, k_, default_)
    i, exactfound = findkey(m.bt, convert(K,k_))
    exactfound? m.bt.data[i].k : convert(K,default_)
end

## Function delete! deletes an item at a given
## key

function delete!{K,D,Ord <: Ordering}(m::SortedDict{K,D,Ord}, k_)
    i, exactfound = findkey(m.bt,convert(K,k_))
    !exactfound && throw(KeyError(k_))
    delete!(m.bt, i)
    m
end

function pop!{K,D,Ord <: Ordering}(m::SortedDict{K,D,Ord}, k_)
    i, exactfound = findkey(m.bt,convert(K,k_))
    !exactfound && throw(KeyError(k_))
    d = m.bt.data[i].d
    delete!(m.bt, i)
    d
end

## Check if two SortedDicts are equal in the sense of containing
## the same (K,D) pairs.  This sense of equality does not mean
## that semitokens valid for one are also valid for the other.

function isequal(m1::SortedDict, m2::SortedDict)
    ord = orderobject(m1)
    if !isequal(ord, orderobject(m2)) || !isequal(eltype(m1), eltype(m2))
        error("Cannot use isequal for two SortedDicts unless their element types and ordering objects are equal")
    end
    p1 = startof(m1)
    p2 = startof(m2)
    while true
        if p1 == pastendtoken(m1)
            return p2 == pastendtoken(m2)
        end
        if p2 == pastendtoken(m2)
            return false
        end
        k1,d1 = deref(p1)
        k2,d2 = deref(p2)
        if !eq(ord,k1,k2) || !isequal(d1,d2)
            return false
        end
        p1 = advance(p1)
        p2 = advance(p2)
    end
end

function mergetwo!{K,D,Ord <: Ordering}(m::SortedDict{K,D,Ord},
                                        m2::SortedDict{K,D,Ord})
    for (k,v) in m2
        m[k] = v
    end
end

function packcopy{K,D,Ord <: Ordering}(m::SortedDict{K,D,Ord})
    w = SortedDict(Dict{K,D}(),orderobject(m))
    mergetwo!(w,m)
    w
end

function packdeepcopy{K,D,Ord <: Ordering}(m::SortedDict{K,D,Ord})
    w = SortedDict(Dict{K,D}(),orderobject(m))
    for (k,v) in m
        newk = deepcopy(k)
        newv = deepcopy(v)
        w[newk] = newv
    end
    w
end

function merge!{K,D,Ord <: Ordering}(m::SortedDict{K,D,Ord},
                                     others::SortedDict{K,D,Ord}...)
    for o in others
        mergetwo!(m,o)
    end
end

function merge{K,D,Ord <: Ordering}(m::SortedDict{K,D,Ord},
                                    others::SortedDict{K,D,Ord}...)
    result = packcopy(m)
    merge!(result, others...)
    result
end

## status of a token is 0 if the token is invalid, 1 if it points to
## ordinary data, 2 if it points to the before-start location and 3 if
## it points to the past-end location.

status{K, D, Ord <: Ordering}(i::SDToken{K,D,Ord}) =
       !(i.semitoken.address in i.container.bt.useddatacells)? 0 :
                        (i.semitoken.address == 1? 2 : (i.semitoken.address == 2? 3 : 1))

## The next three are non-exported correctness checking routines.

not_beforestart{K, D, Ord <: Ordering}(i::SDToken{K,D,Ord}) =
    (!(i.semitoken.address in i.container.bt.useddatacells) ||
     i.semitoken.address == 1) && throw(BoundsError())

not_pastend{K, D, Ord <: Ordering}(i::SDToken{K,D,Ord}) =
    (!(i.semitoken.address in i.container.bt.useddatacells) ||
     i.semitoken.address == 2) &&
       throw(BoundsError())

has_data{K, D, Ord <: Ordering}(i::SDToken{K,D,Ord}) =
    (!(i.semitoken.address in i.container.bt.useddatacells) ||
     i.semitoken.address < 3) &&
       throw(BoundsError())

max_stars_repo_name: ambujone/DIDebug

max_issues_repo_name: ambujone/DIDebug

max_forks_repo_name: ambujone/DIDebug


text:
#
# Autogenerated by Thrift Compiler (1.0.0-dev)
#
# DO NOT EDIT UNLESS YOU ARE SURE THAT YOU KNOW WHAT YOU ARE DOING

module proxy

using Thrift
import Thrift.process, Thrift.meta, Thrift.distribute

export meta
export PartialKey # enum
export TablePermission # enum
export SystemPermission # enum
export ScanType # enum
export ScanState # enum
export ConditionalStatus # enum
export Durability # enum
export CompactionType # enum
export CompactionReason # enum
export IteratorScope # enum
export TimeType # enum
export Key # struct
export ColumnUpdate # struct
export DiskUsage # struct
export KeyValue # struct
export ScanResult # struct
export Range # struct
export ScanColumn # struct
export IteratorSetting # struct
export ScanOptions # struct
export BatchScanOptions # struct
export KeyValueAndPeek # struct
export KeyExtent # struct
export Column # struct
export Condition # struct
export ConditionalUpdates # struct
export ConditionalWriterOptions # struct
export ActiveScan # struct
export ActiveCompaction # struct
export WriterOptions # struct
export CompactionStrategyConfig # struct
export UnknownScanner # exception
export UnknownWriter # exception
export NoMoreEntriesException # exception
export AccumuloException # exception
export AccumuloSecurityException # exception
export TableNotFoundException # exception
export TableExistsException # exception
export MutationsRejectedException # exception
export AccumuloProxyProcessor, AccumuloProxyClient, AccumuloProxyClientBase, login, addConstraint, addSplits, attachIterator, checkIteratorConflicts, clearLocatorCache, cloneTable, compactTable, cancelCompaction, createTable, deleteTable, deleteRows, exportTable, flushTable, getDiskUsage, getLocalityGroups, getIteratorSetting, getMaxRow, getTableProperties, importDirectory, importTable, listSplits, listTables, listIterators, listConstraints, mergeTablets, offlineTable, onlineTable, removeConstraint, removeIterator, removeTableProperty, renameTable, setLocalityGroups, setTableProperty, splitRangeByTablets, tableExists, tableIdMap, testTableClassLoad, pingTabletServer, getActiveScans, getActiveCompactions, getSiteConfiguration, getSystemConfiguration, getTabletServers, removeProperty, setProperty, testClassLoad, authenticateUser, changeUserAuthorizations, changeLocalUserPassword, createLocalUser, dropLocalUser, getUserAuthorizations, grantSystemPermission, grantTablePermission, hasSystemPermission, hasTablePermission, listLocalUsers, revokeSystemPermission, revokeTablePermission, createBatchScanner, createScanner, hasNext, nextEntry, nextK, closeScanner, updateAndFlush, createWriter, update, flush, closeWriter, updateRowConditionally, createConditionalWriter, updateRowsConditionally, closeConditionalWriter, getRowRange, getFollowing # service AccumuloProxy

include("proxy_constants.jl")
include("proxy_types.jl")
include("proxy_impl.jl")  # server methods to be hand coded
include("AccumuloProxy.jl")

end # module proxy


text:
#=

  Euler #33 in Julia.

  Problem 33:
  """
  The fraction 49/98 is a curious fraction, as an inexperienced mathematician in
  attempting to simplify it may incorrectly believe that 49/98 = 4/8, which is correct,
  is obtained by cancelling the 9s.

  We shall consider fractions like, 30/50 = 3/5, to be trivial examples.

  There are exactly four non-trivial examples of this type of fraction, less than
  one in value, and containing two digits in the numerator and denominator.

  If the product of these four fractions is given in its lowest common terms, find
  the value of the denominator.
  """

  This Julia program was created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#

include("Euler.jl")

# 0.00000079s
function euler33a()
    s = 1;
    for y in 1:9
        for z in 1:9
            x = 9.0*y*z/(10.0*y-z)
            if floor(x)==x && y/z < 1.0 && x < 10.0
                s = (s*y)/z;
            end
        end
    end

    return 1/s;
end

run_euler(euler33a)

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
NOnes{N} = NTuple{N,StaticInt{1}}
@inline NOnes(::StaticInt{N}) where {N} = ntuple(Returns(static(1)), Val(N))

@generated function _flatten_indices(inds::I) where {I}
    out = Expr(:block, Expr(:meta, :inline))
    t = Expr(:tuple)
    any_flattened = false
    for i in known_length(I)
        if I.parameters[i] <: AbstractCartesianIndex
            isym = gensym(i)
            push!(out.args, :(@inbounds(Tuple(getfield(inds, $i)))))
            for j in known_length(I.parameters[i])
                push!(t.args, :(@inbounds(getfield($isym, $i))))
            end
            any_flattened = true
        elseif I.parameters[i] <:  CartesianIndices
            isym = gensym(i)
            push!(out.args, :(@inbounds(axes(getfield(inds, $i)))))
            for j in known_length(I.parameters[i])
                push!(t.args, :(@inbounds(getfield($isym, $i))))
            end
            any_flattened = true
        else
            push!(t.args, :(@inbounds(getfield(inds, $i))))
        end
    end
    push!(out.args, t)
    if any_flattened
        return out
    else
        return :inds
    end
end

const SymbolType = Union{Symbol,StaticSymbol}
const BoolType = Union{Bool,True,False}

const DynamicField = Union{Symbol,Int}
const StaticField{F} = Union{StaticSymbol{F},StaticInt{F}}

@inline get_field(x::X, f::DynamicField) where {X} = getfield(x, f)
@inline function get_field(x::X, f::StaticField{F}) where {X,F}
    out = known(field_type(X, f))
    if out === missing
        return getfield(x, F)
    else
        return static(known(field_type(X, f)))
    end
end
const Field{X} = Base.Fix2{typeof(get_field),X}
Field(f) = Base.Fix2(get_field, f)

_check_uplo(::StaticSymbol{S}) where {S} = _check_uplo(S)
function _check_uplo(uplo::Symbol)
    if uplo === :U || uplo === :L || uplo === :L
        return nothing
    else
        throw(ArgumentError("uplo argument must be :U (upper), :L (lower), or :D (diagonal)."))
    end
end


text:
lsv = reval("ls()")
@test length(lsv) == 0
@test isa(lsv, RObject{StrSxp})

lsd = reval("ls(\"package:datasets\")")
@test isa(lsv, RObject{StrSxp})
@test length(lsd) > 50
@test "airmiles" in rcopy(lsd)

langsexp = rlang(:solve, RObject([1 2; 0 4]))
@test length(langsexp) == 2
@test rcopy(reval(langsexp)) == [1 -0.5; 0 0.25]
@test rcopy(langsexp[1]) == :solve
langsexp[1] = RObject(:det)
langsexp[2] = RObject([1 2; 0 0])
@test rcopy(reval(langsexp))[1] == 0

globalEnv[:x] = RObject([1,2,3])
@test rcopy(globalEnv[:x]) == [1,2,3]
globalEnv[:y] = RObject([4,5,6])
@test rcopy(rcall(Symbol("+"),:x,:y)) == [5,7,9]

@test sprint(io -> rprint(io, RObject([1,2,3]))) == "[1] 1 2 3\n"
@test_logs (:warn, "RCall.jl: Warning: hello\n") reval("warning('hello')")
@test_throws RCall.REvalError reval("stop('hello')")

# @test rcopy(rcall(:besselI, 1.0, 2.0)) ≈ besseli(2.0,1.0)
# @test rcopy(rcall(:besselI, 1.0, 2.0, var"expon.scaled"=true)) ≈ besselix(2.0,1.0)

@test isna(R"list(a=1, b=NA)") == [false, true]
@test isna(R"list(a=1, b=NA)", 1) == false
@test isna(R"list(a=1, b=NA)", 2) == true

@test length(R"mtcars") == 11
@test size(R"mtcars") == (32, 11)

# setindex with missing and nothing
a = R"1:10"
a[2] = missing
@test isna(a, 2)

a = R"c('a', 'b', 'c')"
a[2] = missing
@test isna(a) == [false, true, false]

a = R"list(x=1, y=2)"
a[2] = nothing
@test isnull(a[2])
a[:x] = nothing
@test isnull(a[:x])

env = reval("new.env()")
env[:x] = 1
@test rcopy(env[:x]) == 1
env[:x] = nothing
@test isnull(env[:x])

# rparse
@test_throws RCall.RParseError rparse(raw"'\g'")
@test_throws RCall.RParseError rparse("``")

# callbacks
function testfn(x,y;a=3,b=4)
    [x;y;a;b]
end

r = rcall(testfn, 1, 2)
@test isa(r,RObject{IntSxp})
@test rcopy(r) == [1,2,3,4]

r = rcall(testfn, 1, 2,b=6)
@test isa(r,RObject{IntSxp})
@test rcopy(r) == [1,2,3,6]

r = rcall(:optimize,sin,[-2,0])
@test isapprox(r[:minimum][1], -pi/2, atol=eps()^0.25)
r = rcall(:optimize,sin,[0,2],maximum=true)
@test isapprox(r[:maximum][1], pi/2, atol=eps()^0.25)

nullfn() = nothing
@test isa(rcall(nullfn), RObject{NilSxp})

# S4 slots
t1 = R"""
track <- setClass("track", slots = c(x="numeric", y="numeric"))
track(x = 1:10, y = 1:10 + rnorm(10))
"""
@test rcopy(t1[:x]) == collect(1:10)
t1[:x] = 2:11
@test rcopy(t1[:x]) == collect(2:11)
@test_throws Exception t1[:x] = "a"

# graphics
RCall.rgui_init()
let f = tempname()
  rcall(:png,f)
  rcall(:plot,1:10)
  rcall(Symbol("dev.off"))
  @test isfile(f)
  rm(f)
  @test !RCall.rgui_start(true)
  @test_throws ErrorException RCall.rgui_start()
  @test RCall.rgui_stop()
end

# S4 rprint
@test occursin("An object of class",
  sprint(io -> rprint(io, reval("""
   setClass("Foo", representation(x = "numeric"))
   foo <- new("Foo", x = 20)
"""))))

# S3 rprint
@test occursin("hello",
  sprint(io -> rprint(io, reval("""
   print.Bar <- function(x) print("hello")
   bar <- 1
   class(bar) <- "Bar"
   bar
"""))))

@test_logs (:warn, "RCall.jl: Warning in print.Bar(x) : hello\n") rprint(reval("""
   print.Bar <- function(x) warning("hello")
   bar <- 1
   class(bar) <- "Bar"
   bar
"""))

@test_throws RCall.REvalError rprint(reval("""
   print.Bar <- function(x) stop("hello")
   bar <- 1
   class(bar) <- "Bar"
   bar
"""))

# operators
a = reval("a=c(1,2,3)")
b = reval("b=c(4,5,6)")
@test rcopy(a+b)==rcopy(R"a+b")
@test rcopy(a-b)==rcopy(R"a-b")
@test rcopy(a*b)==rcopy(R"a*b")
@test rcopy(a/b)==rcopy(R"a/b")
@test rcopy(a^b)==rcopy(R"a^b")

# misc
iris = rcopy(reval(:iris))
model =  R"lm(Sepal_Length ~ Sepal_Width,data=$iris)"
@test rcopy(RCall.getclass(model)) == "lm"
@test isapprox(rcopy(R"sum($iris$Sepal_Length)"), sum(iris[!, :Sepal_Length]), rtol=4*eps())
@test rcopy(R"factor(rep(1,10))") == fill("1",10)
@test rcopy(R"data.frame(a=rep('test',10), stringsAsFactors = TRUE)")[!, :a] == fill("test",10)

R"""
data(iris)
model <- lm(Sepal.Width ~ Petal.Length, iris)
"""
model = rcopy(R"model")
@test typeof(model[:call]) <: Expr

# getclass
@test rcopy(getclass(reval("1"))) == "numeric"
@test rcopy(getclass(reval("1L"))) == "integer"
@test rcopy(getclass(reval("complex(1,2)"))) == "complex"
if Rversion < v"4"
  @test rcopy(getclass(reval("matrix(1)"))) == "matrix"
else
  @test rcopy(getclass(reval("matrix(1)"))) == ["matrix", "array"]
end
@test rcopy(getclass(reval("function(x) x"))) == "function"
@test rcopy(getclass(reval("data.frame(x=1)"))) == "data.frame"
@test rcopy(getclass(reval("quote(zzz)"))) == "name"


text:
using BitInformation
using Test
import StatsBase.entropy
import Random

include("round_nearest.jl")
include("shave_set_groom.jl")
include("xor_transpose.jl")
include("information.jl")


text:
using GTrends
using Test

@testset "GTrends.jl" begin
    # Write your tests here.
end


text:
module PredatorPreyExample

export run_predatorprey_model!

using Agents
using Random
using InteractiveDynamics
using CairoMakie

mutable struct SheepWolf <: AbstractAgent
    id::Int
    pos::Dims{2}
    type::Symbol
    energy::Float64
    reproduction_prob::Float64
    Δenergy::Float64
end

Sheep(id, pos, energy, reprob, Δe) = SheepWolf(id, pos, :sheep, energy, reprob, Δe)
Wolf(id, pos, energy, reprob, Δe) = SheepWolf(id, pos, :wolf, energy, reprob, Δe)

function initialize_sheep!(model, n_sheep; start_id=0, Δenergy_sheep=4, sheep_reproduce=0.04)
    id = start_id
    for _ in 1:n_sheep
        id += 1
        energy = rand(model.rng, 1:(Δenergy_sheep*2)) - 1
        sheep = Sheep(id, (0, 0), energy, sheep_reproduce, Δenergy_sheep)
        add_agent!(sheep, model)
    end
end

function initialize_wolves!(model, n_wolves; start_id=0, Δenergy_wolf=20, wolf_reproduce=0.05)
    id = start_id
    for _ in 1:n_wolves
        id += 1
        energy = rand(model.rng, 1:(Δenergy_wolf*2)) - 1
        wolf = Wolf(id, (0, 0), energy, wolf_reproduce, Δenergy_wolf)
        add_agent!(wolf, model)
    end
end

function initialize_grass!(model; regrowth_time=30)
    for p in positions(model) # random grass initial growth
        fully_grown = rand(model.rng, Bool)
        countdown = fully_grown ? regrowth_time : rand(model.rng, 1:regrowth_time) - 1
        model.countdown[p...] = countdown
        model.fully_grown[p...] = fully_grown
    end
end

function initialize_model(;
    n_sheep=100,
    n_wolves=50,
    dims=(20, 20),
    regrowth_time=30,
    Δenergy_sheep=4,
    Δenergy_wolf=20,
    sheep_reproduce=0.04,
    wolf_reproduce=0.05,
    seed=23182
)

    rng = MersenneTwister(seed)
    space = GridSpace(dims, periodic=false)
    # Model properties contain the grass as two arrays: whether it is fully grown
    # and the time to regrow. Also have static parameter `regrowth_time`.
    # Notice how the properties are a `NamedTuple` to ensure type stability.
    properties = (
        fully_grown=falses(dims),
        countdown=zeros(Int, dims),
        regrowth_time=regrowth_time,
    )
    model = ABM(SheepWolf, space; properties, rng, scheduler=Schedulers.randomly)
    initialize_sheep!(model, n_sheep;
        Δenergy_sheep=Δenergy_sheep,
        sheep_reproduce=sheep_reproduce
    )
    initialize_wolves!(model, n_wolves;
        start_id=n_sheep,
        Δenergy_wolf=Δenergy_wolf,
        wolf_reproduce=wolf_reproduce
    )
    initialize_grass!(model; regrowth_time=regrowth_time)
    return model
end

function sheepwolf_step!(agent::SheepWolf, model)
    if agent.type == :sheep
        sheep_step!(agent, model)
    else # then `agent.type == :wolf`
        wolf_step!(agent, model)
    end
end

function sheep_step!(sheep, model)
    walk!(sheep, rand, model)
    sheep.energy -= 1
    sheep_eat!(sheep, model)
    if sheep.energy < 0
        kill_agent!(sheep, model)
        return
    end
    if rand(model.rng) <= sheep.reproduction_prob
        reproduce!(sheep, model)
    end
end

function nearest_sheep_in_range(wolf, model; range=3)
    agents = collect(nearby_agents(wolf,model,range))
    sheeps = filter!(x -> x.type == :sheep, agents)
    nearest_sheep = nothing
    nearest_distance = Inf
    for sheep in sheeps
        edist = edistance(wolf, sheep, model)
        if edist < nearest_distance
            nearest_distance = edist
            nearest_sheep = sheep
        end
    end
    return nearest_sheep
end

# A minor improvement on wolf_step! :
# if there are no sheep in the current position,
# the wolf will walk towards a sheep in the vision range
# instead of walking in a random direction
function wolf_step!(wolf, model)
    wolf.energy -= 1
    agents = collect(agents_in_position(wolf.pos, model))
    dinner = filter!(x -> x.type == :sheep, agents)
    if isempty(dinner)
        walk!(wolf, rand, model)
    else
        ns = nearest_sheep_in_range(wolf,model)
        if ns === nothing
            walk!(wolf, rand, model)
        else
            nsd = get_direction(wolf.pos,ns.pos,model)
            walk!(wolf, nsd, model)
        end
    end
    wolf_eat!(wolf, dinner, model)
    if wolf.energy < 0
        kill_agent!(wolf, model)
        return
    end
    if rand(model.rng) <= wolf.reproduction_prob
        reproduce!(wolf, model)
    end
end

function sheep_eat!(sheep, model)
    if model.fully_grown[sheep.pos...]
        sheep.energy += sheep.Δenergy
        model.fully_grown[sheep.pos...] = false
    end
end

function wolf_eat!(wolf, sheep, model)
    if !isempty(sheep)
        dinner = rand(model.rng, sheep)
        kill_agent!(dinner, model)
        wolf.energy += wolf.Δenergy
    end
end

function reproduce!(agent, model)
    agent.energy /= 2
    id = nextid(model)
    offspring = SheepWolf(
        id,
        agent.pos,
        agent.type,
        agent.energy,
        agent.reproduction_prob,
        agent.Δenergy,
    )
    add_agent_pos!(offspring, model)
    return
end

function grass_step!(model)
    @inbounds for p in positions(model) # we don't have to enable bound checking
        if !(model.fully_grown[p...])
            if model.countdown[p...] ≤ 0
                model.fully_grown[p...] = true
                model.countdown[p...] = model.regrowth_time
            else
                model.countdown[p...] -= 1
            end
        end
    end
end

offset(a) = a.type == :sheep ? (-0.7, -0.5) : (-0.3, -0.5)
ashape(a) = a.type == :sheep ? :circle : :utriangle
acolor(a) = a.type == :sheep ? :white : :black
grasscolor(model) = model.countdown ./ model.regrowth_time
heatkwargs = (colormap = [:brown, :green], colorrange = (0, 1))
sheep(a) = a.type == :sheep
wolves(a) = a.type == :wolf
count_grass(model) = count(model.fully_grown)

# data to collect
adata = [(sheep, count), (wolves, count)]
mdata = [count_grass]

function run_predatorprey_model!(;file="sheepwolf.mp4",frames=500,spf=1)
    model = initialize_model(
    n_wolves = 20,
    dims = (25, 25),
    Δenergy_sheep = 5,
    sheep_reproduce = 0.2,
    wolf_reproduce = 0.08,
    seed = 7758
    )
    adf = init_agent_dataframe(model, adata)
    mdf = init_model_dataframe(model, mdata)
    # Create 'Step N' as the Observable title
    s = Observable(0) # counter of current step
    t = lift(x -> "step = " * string(x), s)
    axiskwargs = (title = t, titlealign = :left) #title and position
    fig, abmstepper = abm_plot(
        model;
        ac = acolor,
        as = 15,
        am = ashape,
        offset = offset,
        heatarray = grasscolor,
        heatkwargs = heatkwargs,
        axiskwargs = axiskwargs
    )
    record(fig, file; framerate = 8) do io
        for j in 1:frames-1
            recordframe!(io)
            # MUST USE THE SAME ABMSTEPPER! otherwise frame won't update
            step!(abmstepper, model, sheepwolf_step!, grass_step!, spf)
            collect_model_data!(mdf,model,mdata,j)
            collect_agent_data!(adf,model,adata,j)
            s[] += spf
        end
        recordframe!(io) # last frame
    end
    isfile(file) || error("file not created")
    return (file, adf, mdf)
end

function plot_population_timeseries(adf, mdf; figfile="predatorprey.png")
    figure = Figure(resolution = (600, 400))
    ax = figure[1, 1] = Axis(figure; xlabel = "Step", ylabel = "Population")
    sheepl = lines!(ax, adf.step, adf.count_sheep, color = :blue)
    wolfl = lines!(ax, adf.step, adf.count_wolves, color = :orange)
    grassl = lines!(ax, mdf.step, mdf.count_grass, color = :green)
    figure[1, 2] = Legend(figure, [sheepl, wolfl, grassl], ["Sheep", "Wolves", "Grass"])
    save(figfile,figure)
    isfile(figfile) || error("plot file not created")
    return figure
end

# _, adf, mdf = run_predatorprey_model!()
# plot_population_timeseries(adf, mdf)

end


text:
"""
# IAG: The Identity and Access Governance Julia module
"""
module IAG

include("AbstractThing.jl");
export AbstractThing;
include("AbstractObject.jl");
export AbstractObject;
include("AbstractSubject.jl");
export AbstractSubject;
include("AbstractRight.jl");
export AbstractRight;

include("Object.jl");
export Object;
include("Subject.jl");
export Subject;
include("Right.jl");
export Right;
include("AccessControlSystem.jl");
export AccessControlSystem;

include("get_label.jl");
export get_label;
export get_label!;
include("set_label.jl");
export set_label;
export set_label!;

end; # module

using .IAG;

max_stars_repo_name: daviddoret/julia-cybersecurity

max_issues_repo_name: daviddoret/julia-cybersecurity

max_forks_repo_name: daviddoret/julia-cybersecurity


text:
import NMRDataSetup

using FFTW
import PyPlot
import BSON

import Base.@kwdef

PyPlot.close("all")
fig_num = 1

PyPlot.matplotlib["rcParams"][:update](["font.size" => 22, "font.family" => "serif"])

@kwdef mutable struct StoreConfigType
    save_dir::String
    project_name::String
    experiment_full_path::String
    solvent_ppm_guess::Float64 = 4.7
    solvent_window_ppm::Float64 = 0.1
end
# function StoreConfigType(x,y,z)
#     return StoreConfigType(x,y,z,4.7,0.1)
# end

metadata_set = Vector{StoreConfigType}(undef, 0)

### user inputs.
save_dir = "/home/roy/MEGAsync/outputs/NMR/calibrate/final"

push!(metadata_set, StoreConfigType( save_dir = save_dir,
project_name = "L-Alanine-700",
experiment_full_path = "/home/roy/Documents/repo/NMRData/combination/BMRB-700/L-Alanine"))

push!(metadata_set, StoreConfigType( save_dir = save_dir,
project_name = "L-Arginine-500-2mM",
experiment_full_path = "/home/roy/Documents/repo/NMRData/combination/BMRB-500-2mM/L-Arginine"))

# missing aspartic.

push!(metadata_set, StoreConfigType( save_dir = save_dir,
project_name = "L-Cysteine-500-2mM",
experiment_full_path = "/home/roy/Documents/repo/NMRData/combination/BMRB-500-2mM/L-Cysteine"))

push!(metadata_set, StoreConfigType( save_dir = save_dir,
project_name = "L-Glutamic acid-700",
experiment_full_path = "/home/roy/Documents/repo/NMRData/combination/BMRB-700/L-Glutamic acid"))

push!(metadata_set, StoreConfigType( save_dir = save_dir,
project_name = "L-Glutamine-700",
experiment_full_path = "/home/roy/Documents/repo/NMRData/combination/BMRB-700/L-Glutamine"))

push!(metadata_set, StoreConfigType( save_dir = save_dir,
project_name = "Glycine-500-2mM",
experiment_full_path = "/home/roy/Documents/repo/NMRData/combination/BMRB-500-2mM/Glycine"))

push!(metadata_set, StoreConfigType( save_dir = save_dir,
project_name = "L-Histidine-700",
experiment_full_path = "/home/roy/Documents/repo/NMRData/combination/BMRB-700/L-Histidine"))

push!(metadata_set, StoreConfigType( save_dir = save_dir,
project_name = "L-Isoleucine-700",
experiment_full_path = "/home/roy/Documents/repo/NMRData/combination/BMRB-700/L-Isoleucine"))

push!(metadata_set, StoreConfigType( save_dir = save_dir,
project_name = "L-Leucine-500-2mM",
experiment_full_path = "/home/roy/Documents/repo/NMRData/combination/BMRB-500-2mM/L-Leucine"))

push!(metadata_set, StoreConfigType( save_dir = save_dir,
project_name = "L-Lysine-500-2mM",
experiment_full_path = "/home/roy/Documents/repo/NMRData/combination/BMRB-500-2mM/L-Lysine"))

push!(metadata_set, StoreConfigType( save_dir = save_dir,
project_name = "L-Methionine-2mM",
experiment_full_path = "/home/roy/Documents/repo/NMRData/combination/BMRB-500-2mM/L-Methionine"))

push!(metadata_set, StoreConfigType( save_dir = save_dir,
project_name = "L-Phenylalanine-700",
experiment_full_path = "/home/roy/Documents/repo/NMRData/combination/BMRB-700/L-Phenylalanine"))

push!(metadata_set, StoreConfigType( save_dir = save_dir,
project_name = "L-Proline-500-2mM",
experiment_full_path = "/home/roy/Documents/repo/NMRData/combination/BMRB-500-2mM/L-Proline"))

push!(metadata_set, StoreConfigType( save_dir = save_dir,
project_name = "L-Serine-700",
experiment_full_path = "/home/roy/Documents/repo/NMRData/combination/BMRB-700/L-Serine"))

push!(metadata_set, StoreConfigType( save_dir = save_dir,
project_name = "L-Threonine-700",
experiment_full_path = "/home/roy/Documents/repo/NMRData/combination/BMRB-700/L-Threonine"))

push!(metadata_set, StoreConfigType( save_dir = save_dir,
project_name = "L-Tryptophan-700",
experiment_full_path = "/home/roy/Documents/repo/NMRData/combination/BMRB-700/L-Tryptophan"))

push!(metadata_set, StoreConfigType( save_dir = save_dir,
project_name = "L-Tyrosine-500-2mM",
experiment_full_path = "/home/roy/Documents/repo/NMRData/combination/BMRB-500-2mM/L-Tyrosine"))

push!(metadata_set, StoreConfigType( save_dir = save_dir,
project_name = "L-Valine-700",
experiment_full_path = "/home/roy/Documents/repo/NMRData/combination/BMRB-700/L-Valine"))

push!(metadata_set, StoreConfigType( save_dir = save_dir,
project_name = "D-(+)-Glucose-700",
experiment_full_path = "/home/roy/Documents/repo/NMRData/combination/BMRB-700/D-(+)-Glucose"))

push!(metadata_set, StoreConfigType( save_dir = save_dir,
project_name = "D-(+)-Glucose-NRC-600",
experiment_full_path = "/home/roy/Documents/repo/NMRData/src/experiments/NRC/misc/glucose/Sep-25-2018"))

### end inputs.

## load.
function batchstoreexp(metadata_set)

    for i = 1:length(metadata_set)
        save_dir = metadata_set[i].save_dir
        project_name = metadata_set[i].project_name
        experiment_full_path = metadata_set[i].experiment_full_path
        solvent_ppm_guess = metadata_set[i].solvent_ppm_guess
        solvent_window_ppm = metadata_set[i].solvent_window_ppm

        s_t, S, hz2ppmfunc, ppm2hzfunc, ν_0ppm, fs, SW, α_0ppm, β_0ppm, λ_0ppm, Ω_0ppm,
            results_0ppm, dic, α_solvent, β_solvent, λ_solvent, Ω_solvent,
            results_solvent = NMRDataSetup.loadspectrum(experiment_full_path;
            solvent_ppm = solvent_ppm_guess,
            solvent_window_ppm = solvent_window_ppm)

        ## store.
        save_folder_path = joinpath(save_dir, project_name)
        isdir(save_folder_path) || mkpath(save_folder_path); # make save folder if it doesn't exist.

        save_path = joinpath(save_folder_path, "$(project_name).bson")
        BSON.bson(save_path,
        s_t = s_t,
        fs = fs,
        SW = SW,
        ν_0ppm = ν_0ppm)
    end
end

### batch.
### end batch.

### singular.
resize!(metadata_set, 0)

push!(metadata_set, StoreConfigType( save_dir = save_dir,
project_name = "Nam2022_Serine",
experiment_full_path = "/home/roy/MEGAsync/data/NMR/NRC/Nam_4_amino_acid_mixture_Jan_2022/NHK220113_4AA_glu_PBS_d2o/1"))
### end singular.

# store.
batchstoreexp(metadata_set)


text:
"""
    RDA (R data archive) reading context.

    * Stores flags that define how R objects are read and converted
      into Julia objects.
    * Maintains the list of R objects that could be referenced later in
      the RDA stream.
"""
type RDAContext{T <: RDAIO}
    io::T                      # RDA input stream

    # RDA format properties
    fmtver::UInt32             # format version
    Rver::VersionNumber        # R version used to write the file
    Rmin::VersionNumber        # minimal R version to read the file

    kwdict::Dict{Symbol,Any}   # options defining RDA deserialization behaviour

    # intermediate data
    ref_tab::Vector{RSEXPREC}  # SEXP array for references

    function RDAContext(io::T, kwoptions::Vector{Any})
        fmtver = readint32(io)
        rver = readint32(io)
        rminver = readint32(io)
        kwdict = Dict{Symbol,Any}(kwoptions)
        new(io,
            fmtver,
            VersionNumber(div(rver,65536), div(rver%65536, 256), rver%256),
            VersionNumber(div(rminver,65536), div(rminver%65536, 256), rminver%256),
            kwdict,
            RSEXPREC[])
    end
end

RDAContext{T <: RDAIO}(io::T, kwoptions::Vector{Any}) = RDAContext{T}(io, kwoptions)

"""
    Registers R object, so that it could be referenced later
    (by its index in the reference table).
"""
function registerref(ctx::RDAContext, obj::RSEXPREC)
    push!(ctx.ref_tab, obj)
    obj
end


text:
#
#    formats = (
#   '{{first_name}} {{last_,"{last_name}}'",
#   '{{first_name}} {{first_,"{last_name}}'",
#   '{{first_name}} {{first_,"{last_name}} {{last_name}},"'{{first_name}} {{last_name}}'",
#   '{{prefix}} {{first_,"{last_name}}'",
#    )

    first_names = ("Abel","Abelardo","Abigail","Abraham","Abril","Adalberto","Adán","Adela","Adriana","Aida","Alejandra","Agustín","Alberto","Aldonza","Alicia","Alta  Gracia","Alonso","Aldo","Alejandro","Alfonso","Alfredo","Alma","Alvaro","Amalia","Amanda","Amador","Amelia","Ana","Anabel","Ana Luisa","Ana María","Anel","Andrea","Andrés","Ángel","Ángela","Angélica","Antonia","Antonio","Araceli","Arcelia","Ariadna","Armando","Arturo","Asunción","Augusto","Aurora","Aurelio","Barbara","Beatriz'","Berta","Benito","Benjamín","Bernardo","Bernabé","Bianca","Blanca","Bruno","Camila","Camilo","Caridad","Carla","Carlos","Carlota","Carmen","Carolina","Catalina","César","Cecilia","Celia","Citlali","Clara","Claudia","Claudio","Clemente","Concepción","Conchita","Cornelio","Cristian","Cristal","Cristina","Cristobal","Cynthia","Dalia","Daniel","Daniela","Darío","David","Débora","Delia","Diana","Diego","Dolores","Dulce","Dulce María","Eduardo","Elena","Elias","Elisa","Eloisa","Elsa","Elvia","Elvira","Eloy","Emilia","Emiliano","Emilio","Enrique","Eric","Ernesto","Esmeralda","Esteban","Estefanía","Estela","Esparta","Espartaco","Esperanza","Estela","Esther","Eugenia","Eugenio","Eva","Evelio","Fabiola","Federico","Felipe","Fernando","Felix","Fidel","Flavio","Florencia","Francisco","Francisco Javier","Francisca","Frida","Gabino","Gabriela","Gabriel","Genaro","Georgina","Gerardo","Gerónimo","Germán","Gilberto","Guillermina","Gloria","Gonzalo","Graciela","Gregorio","Guillermo","Guadalupe","Gustavo","Héctor","Helena","Hermelinda' 'Hernán","Hilda","Homero","Horacio","Hugo","Humberto","Ignacio","Ilse","Indira","Inés","Irene","Irma","Itzel","Isaac","Isabel","Isabela","Israel","Iván","Ivonne","Jacinto' 'Jacobo","Jaime","Javier","Jaqueline","Jerónimo","Jesús","Joaquín","Jonás","Jorge","Jorge Luis","Jos","José","Josefina","José Carlos","José Eduardo","José Emilio","José Luis","José Manuél","José María","Juan","Juana","Juan Carlos","Judith","Julia","Julio","Julio César","Laura","Leonardo","Leonel","Leonor","Karla","Karina' 'Leticia","Lorenzo","Lucas","Lilia","Liliana","Linda","Lorena","Lourdes","Lucía","Luisa","Luz","Luis","Luis Miguel","Luis Manuel","Magdalena","Manuel","Marco Antonio","Marcela","Marcos","Margarita","María","Marisela","Marisol","María del Carmen","María Cristina","María Elena","María Eugenia","María José","María Luisa","María Teresa","Marisol","Martha","Mayte","Mariano","Mariana","Mario","Martín","Mateo","Mauro","Mauricio","Maximiliano","Mercedes","Micaela","Minerva","Mitzy","Miguel","Miguel Ángel","Miriam","Modesto","Mónica","Nadia","Natalia","Natividad","Nancy","Nayeli","Nelly","Noelia","Noemí","Norma","Nicolás","Octavio","Ofelia","Olivia","Óliver","Olga","Óscar","Oswaldo","Omar","Pablo","Paola","Patricia","Pamela","Patricio","Pascual","Paulina","Pedro","Perla","Pilar","Porfirio","Rafaél","Ramiro","Ramón","Raúl","Raquel","Rebeca","Reina","Renato","René","Reynaldo","Ricardo","Roberto","Rodolfo","Rocío","Rodrigo","Rolando","Rosa","Rosalia","Rosario","Rubén","Rufino","Ruby","Salvador","Salma","Samuel","Sandra","Santiago","Sara","Sessa","Sergio","Serafín","Silvano","Silvia","Sofía","Socorro","Soledad","Sonia","Susana","Tania","Teresa","Teodoro","Timoteo","Tomás","Trinidad","Verónica","Vicente","Violeta","Uriel","Úrsula","Vanesa","Víctor","Victoria","Virginia","Wilfrido","Wendolin","Yeni'","Yolanda","Yuridia","Zacarías","Zeferino","Zoé")

    last_names = ( "Abrego", "Abreu", "Acevedo", "Acosta", "Acuña","Adame", "Aguayo", "Aguilar", "Aguilera", "Aguirre", "Alarcón", "Alba","Alcala", "Alcántar", "Alcaraz", "Alejandro", "Alemán", "Alfaro","Almanza", "Almaraz", "Almonte","Alonso", "Alonzo", "Altamirano","Alva", "Alvarado", "Alvarez", "Amador", "Amaya", "Anaya", "Anguiano","Angulo","Aparicio", "Apodaca", "Aponte", "Aragón", "Aranda", "Arce","Archuleta", "Arellano", "Arenas", "Arevalo", "Arguello", "Arias","Armas", "Armendáriz", "Armenta", "Arredondo", "Arreola", "Arriaga","Arroyo", "Arteaga", "Ávalos", "Ávila", "Avilés", "Ayala", "Baca","Badillo", "Báez", "Baeza", "Bahena", "Balderas", "Ballesteros","Bañuelos", "Barajas", "Barela", "Barragán", "Barraza", "Barrera","Barreto", "Barrientos", "Barrios", "Batista", "Becerra", "Beltrán","Benavides", "Benavídez", "Benítez", "Bermúdez", "Bernal", "Berríos","Bétancourt", "Blanco", "Bonilla", "Borrego", "Botello", "Bravo","Briones", "Briseño", "Brito", "Bueno", "Burgos", "Bustamante","Bustos", "Caballero", "Cabán", "Cabrera", "Cadena", "Caldera","Calderón", "Calvillo", "Camacho", "Camarillo", "Campos", "Canales","Candelaria", "Cano", "Cantú", "Caraballo", "Carbajal", "Cardenas","Cardona", "Carmona", "Carranza", "Carrasco", "Carreón", "Carrera","Carrero", "Carrillo", "Carrión", "Carvajal", "Casanova", "Casares","Casárez", "Casas", "Casillas", "Castañeda", "Castellanos", "Castillo","Castro", "Cavazos", "Cazares", "Ceballos", "Cedillo", "Ceja","Centeno", "Cepeda", "Cervantes", "Cervántez", "Chacón", "Chapa","Chavarría", "Chávez", "Cintrón", "Cisneros", "Collado", "Collazo","Colón", "Colunga", "Concepción", "Contreras", "Cordero", "Córdova","Cornejo", "Corona", "Coronado", "Corral", "Corrales", "Correa","Cortés", "Cortez", "Cotto", "Covarrubias", "Crespo", "Cruz", "Cuellar","Curiel", "Dávila", "de Anda", "de Jesús", "de la Crúz", "de la Fuente","de la Garza", "de la O", "de la Rosa", "de la Torre", "de León","Delgadillo", "Delgado", "del Río", "del Valle", "Díaz", "Domínguez","Duarte", "Dueñas", "Durán", "Echeverría", "Elizondo", "Enríquez","Escalante", "Escamilla", "Escobar", "Escobedo", "Esparza", "Espinal","Espino", "Espinosa", "Espinoza", "Esquibel", "Esquivel", "Estévez","Estrada", "Fajardo", "Farías", "Feliciano", "Fernández", "Ferrer","Fierro", "Figueroa", "Flores", "Flórez", "Fonseca", "Franco", "Frías","Fuentes", "Gaitán", "Galarza", "Galindo", "Gallardo", "Gallegos","Galván", "Gálvez", "Gamboa", "Gamez", "Gaona", "Garay", "García","Garibay", "Garica", "Garrido", "Garza", "Gastélum", "Gaytán", "Gil","Girón", "Godínez", "Godoy", "Gómez", "Gonzales", "González", "Gollum","Gracia", "Granado", "Granados", "Griego", "Grijalva", "Guajardo","Guardado", "Guerra", "Guerrero", "Guevara", "Guillen", "Gurule","Gutiérrez", "Guzmán", "Haro", "Henríquez", "Heredia", "Hernádez","Hernandes", "Hernández", "Herrera", "Hidalgo", "Hinojosa", "Holguín","Huerta", "Hurtado", "Ibarra", "Iglesias", "Irizarry", "Jaime","Jaimes", "Jáquez", "Jaramillo", "Jasso", "Jiménez", "Jimínez","Juárez", "Jurado", "Laboy", "Lara", "Laureano", "Leal", "Lebrón","Ledesma", "Leiva", "Lemus", "León", "Lerma", "Leyva", "Limón","Linares", "Lira", "Llamas", "Loera", "Lomeli", "Longoria", "López","Lovato", "Loya", "Lozada", "Lozano", "Lucero", "Lucio", "Luevano","Lugo", "Luna", "Macías", "Madera", "Madrid", "Madrigal", "Maestas","Magaña", "Malave", "Maldonado", "Manzanares", "Mares", "Marín","Márquez", "Marrero", "Marroquín", "Martínez", "Mascareñas", "Mata","Mateo", "Matías", "Matos", "Maya", "Mayorga", "Medina", "Medrano","Mejía", "Meléndez", "Melgar", "Mena", "Menchaca", "Méndez", "Mendoza","Menéndez", "Meraz", "Mercado", "Merino", "Mesa", "Meza", "Miramontes","Miranda", "Mireles", "Mojica", "Molina", "Mondragón", "Monroy","Montalvo", "Montañez", "Montaño", "Montemayor", "Montenegro","Montero", "Montes", "Montez", "Montoya", "Mora", "Morales", "Moreno","Mota", "Moya", "Munguía", "Muñiz", "Muñoz", "Murillo", "Muro","Nájera", "Naranjo", "Narváez", "Nava", "Navarrete", "Navarro","Nazario", "Negrete", "Negrón", "Nevárez", "Nieto", "Nieves", "Niño","Noriega", "Núñez", "Ocampo", "Ocasio", "Ochoa", "Ojeda", "Olivares","Olivárez", "Olivas", "Olivera", "Olivo", "Olmos", "Olvera","Ontiveros", "Oquendo", "Ordóñez", "Orellana", "Ornelas", "Orosco","Orozco", "Orta", "Ortega", "Ortiz", "Osorio", "Otero", "Ozuna","Pabón", "Pacheco", "Padilla", "Padrón", "Páez", "Palacios", "Palomino","Palomo", "Pantoja", "Paredes", "Parra", "Partida", "Patiño", "Paz","Pedraza", "Pedroza", "Pelayo", "Peña", "Perales", "Peralta", "Perea","Peres", "Pérez", "Pichardo", "Piña", "Pineda", "Pizarro", "Polanco","Ponce", "Porras", "Portillo", "Posada", "Prado", "Preciado", "Prieto","Puente", "Puga", "Pulido", "Quesada", "Quezada", "Quiñones","Quiñónez", "Quintana", "Quintanilla", "Quintero", "Quiroz", "Rael","Ramírez", "Ramón", "Ramos", "Rangel", "Rascón", "Raya", "Razo","Regalado", "Rendón", "Rentería", "Reséndez","Reyes", "Reyna","Reynoso", "Rico", "Rincón", "Riojas", "Ríos", "Rivas", "Rivera","Rivero", "Robledo", "Robles", "Rocha", "Rodarte", "Rodrígez","Rodríguez", "Rodríquez", "Rojas", "Rojo", "Roldán", "Rolón", "Romero","Romo", "Roque", "Rosado", "Rosales", "Rosario", "Rosas", "Roybal","Rubio", "Ruelas", "Ruiz", "Saavedra", "Sáenz", "Saiz", "Salas","Salazar", "Salcedo", "Salcido", "Saldaña", "Saldivar", "Salgado","Salinas", "Samaniego", "Sanabria", "Sanches", "Sánchez", "Sandoval","Santacruz", "Santana", "Santiago", "Santillán", "Sarabia", "Sauceda","Saucedo", "Segovia","Segura", "Sepúlveda", "Serna", "Serrano","Serrato", "Sevilla", "Sierra", "Sisneros", "Solano", "Solís", "Soliz","Solorio", "Solorzano", "Soria", "Sosa", "Sotelo", "Soto", "Suárez","Tafoya", "Tamayo", "Tamez", "Tapia", "Tejada", "Tejeda", "Téllez","Tello", "Terán", "Terrazas", "Tijerina", "Tirado", "Toledo", "Toro","Torres", "Tórrez", "Tovar", "Trejo", "Treviño", "Trujillo", "Ulibarri","Ulloa", "Urbina", "Ureña", "Urías", "Uribe", "Urrutia", "Vaca","Valadez", "Valdés", "Valdez", "Valdivia", "Valencia", "Valentín","Valenzuela", "Valladares", "Valle", "Vallejo", "Valles", "Valverde","Vanegas", "Varela", "Vargas", "Vásquez", "Vázquez", "Vega", "Vela","Velasco", "Velásquez", "Velázquez", "Vélez", "Véliz", "Venegas","Vera", "Verdugo", "Verduzco", "Vergara", "Viera", "Vigil", "Villa","Villagómez", "Villalobos", "Villalpando", "Villanueva", "Villareal","Villarreal", "Villaseñor", "Villegas", "Yáñez", "Ybarra", "Zambrano","Zamora", "Zamudio", "Zapata", "Zaragoza", "Zarate", "Zavala", "Zayas","Zedillo", "Zelaya", "Zepeda", "Zúñiga" )

prefixes = ("Sr(a).","Lic.","Ing.","Dr(a).")

first_name() = first_names[rand(1:length(first_names))]

last_name() = last_names[rand(1:length(last_names))]

prefixe() = prefixes[rand(1:length(prefixes))]


text:
import Base: *, size, ∈, isempty

export SparseMatrixExp,
       ExponentialMap,
       size,
       ProjectionSparseMatrixExp,
       ExponentialProjectionMap,
       get_row,
       get_rows,
       get_column,
       get_columns

# --- SparseMatrixExp & ExponentialMap ---

"""
    SparseMatrixExp{N}

Type that represents the matrix exponential, ``\\exp(M)``, of a sparse matrix.

### Fields

- `M` -- sparse matrix

### Examples

Take for exammple a random sparse matrix:

```jldoctest SparseMatrixExp_constructor
julia> A = sprandn(100, 100, 0.1);

julia> E = SparseMatrixExp(A);

julia> size(E)
(100, 100)
```

Now, `E` is a lazy representation of ``\\exp(A)``. To compute with `E`, use
`get_row` and `get_column` (or `get_rows` and `get_columns`;
they return row and column vectors (or matrices). For example:

```jldoctest SparseMatrixExp_constructor
julia> get_row(E, 10); # compute E[10, :]

julia> get_column(E, 10); # compute E[:, 10]

julia> get_rows(E, [10]); # same as get_row(E, 10) but a 1x100 matrix is returned

julia> get_columns(E, [10]); # same as get_column(E, 10) but a 100x1 matrix is returned
```

### Notes

This type is provided for use with very large and very sparse matrices.
The evaluation of the exponential matrix action over vectors relies on the
[Expokit](https://github.com/acroy/Expokit.jl) package.
"""
struct SparseMatrixExp{N, MN<:AbstractSparseMatrix{N}} <: AbstractMatrix{N}
    M::MN
end

SparseMatrixExp(M::Matrix) =
        error("only sparse matrices can be used to create a `SparseMatrixExp`")

Base.IndexStyle(::Type{<:SparseMatrixExp}) = IndexCartesian()
Base.getindex(spmexp::SparseMatrixExp, I::Vararg{Int, 2}) = get_column(spmexp, I[2])[I[1]]

function size(spmexp::SparseMatrixExp)::Tuple{Int, Int}
    return size(spmexp.M)
end

function size(spmexp::SparseMatrixExp, ax::Int)::Int
    return size(spmexp.M, ax)
end

function load_expokit_sparsematrixexp()
return quote

using .Expokit: expmv

function get_column(spmexp::SparseMatrixExp{N}, j::Int)::Vector{N} where {N}
    n = size(spmexp, 1)
    aux = zeros(N, n)
    aux[j] = one(N)
    return expmv(one(N), spmexp.M, aux)
end

function get_columns(spmexp::SparseMatrixExp{N},
                     J::AbstractArray)::Matrix{N} where {N}

    n = size(spmexp, 1)
    aux = zeros(N, n)
    ans = zeros(N, n, length(J))
    count = 1
    one_N = one(N)
    zero_N = zero(N)
    @inbounds for j in J
        aux[j] = one_N
        ans[:, count] = expmv(one_N, spmexp.M, aux)
        aux[j] = zero_N
        count += 1
    end
    return ans
end

"""
    get_row(spmexp::SparseMatrixExp{N}, i::Int) where {N}

Return a single row of a sparse matrix exponential.

### Input

- `spmexp` -- sparse matrix exponential
- `i`      -- row index

### Output

A row vector corresponding to the `i`th row of the matrix exponential.

### Notes

This function uses Julia's `transpose` function to create the result.
The result is of type `Transpose`; in Julia versions older than v0.7, the result
was of type `RowVector`.
"""
function get_row(spmexp::SparseMatrixExp{N}, i::Int) where {N}
    n = size(spmexp, 1)
    aux = zeros(N, n)
    aux[i] = one(N)
    return transpose(expmv(one(N), transpose(spmexp.M), aux))
end

function get_rows(spmexp::SparseMatrixExp{N},
                  I::AbstractArray{Int})::Matrix{N} where {N}
    n = size(spmexp, 1)
    aux = zeros(N, n)
    ans = zeros(N, length(I), n)
    Mtranspose = transpose(spmexp.M)
    count = 1
    one_N = one(N)
    zero_N = zero(N)
    @inbounds for i in I
        aux[i] = one_N
        ans[count, :] = expmv(one_N, Mtranspose, aux)
        aux[i] = zero_N
        count += 1
    end
    return ans
end

end end  # quote / load_expokit_sparsematrixexp

"""
    ExponentialMap{N<:Real, S<:LazySet{N}} <: LazySet{N}

Type that represents the action of an exponential map on a convex set.

### Fields

- `spmexp` -- sparse matrix exponential
- `X`      -- convex set

### Examples

The `ExponentialMap` type is overloaded to the usual times `*` operator when the
linear map is a lazy matrix exponential. For instance,

```jldoctest
julia> A = sprandn(100, 100, 0.1);

julia> E = SparseMatrixExp(A);

julia> B = BallInf(zeros(100), 1.);

julia> M = E * B; # represents the image set: exp(A) * B

julia> M isa ExponentialMap
true

julia> dim(M)
100
```
"""
struct ExponentialMap{N, S<:LazySet{N}} <: LazySet{N}
    spmexp::SparseMatrixExp{N}
    X::S
end

"""
```
    *(spmexp::SparseMatrixExp{N},
      X::LazySet{N})::ExponentialMap{N} where {N<:Real}
```

Return the exponential map of a convex set from a sparse matrix exponential.

### Input

- `spmexp` -- sparse matrix exponential
- `X`      -- convex set

### Output

The exponential map of the convex set.
"""
function *(spmexp::SparseMatrixExp{N},
           X::LazySet{N})::ExponentialMap{N} where {N<:Real}
    return ExponentialMap(spmexp, X)
end

"""
    dim(em::ExponentialMap)::Int

Return the dimension of an exponential map.

### Input

- `em` -- an ExponentialMap

### Output

The ambient dimension of the exponential map.
"""
function dim(em::ExponentialMap)::Int
    return size(em.spmexp.M, 1)
end

function load_expokit_exponentialmap()
return quote

using .Expokit: expmv

"""
    σ(d::AbstractVector{N}, em::ExponentialMap{N}) where {N<:Real}

Return the support vector of the exponential map.

### Input

- `d`  -- direction
- `em` -- exponential map

### Output

The support vector in the given direction.
If the direction has norm zero, the result depends on the wrapped set.

### Notes

If ``E = \\exp(M)⋅S``, where ``M`` is a matrix and ``S`` is a convex set, it
follows that ``σ(d, E) = \\exp(M)⋅σ(\\exp(M)^T d, S)`` for any direction ``d``.

We allow sparse direction vectors, but will convert them to dense vectors to be
able to use `expmv`.
"""
function σ(d::AbstractVector{N}, em::ExponentialMap{N}) where {N<:Real}
    d_dense = d isa Vector ? d : Vector(d)
    v = expmv(one(N), transpose(em.spmexp.M), d_dense) # v   <- exp(M') * d
    return expmv(one(N), em.spmexp.M, σ(v, em.X)) # res <- exp(M) * σ(v, S)
end

"""
    ρ(d::AbstractVector{N}, em::ExponentialMap{N}) where {N<:Real}

Return the support function of the exponential map.

### Input

- `d`  -- direction
- `em` -- exponential map

### Output

The support function in the given direction.

### Notes

If ``E = \\exp(M)⋅S``, where ``M`` is a matrix and ``S`` is a convex set, it
follows that ``ρ(d, E) = ρ(\\exp(M)^T d, S)`` for any direction ``d``.

We allow sparse direction vectors, but will convert them to dense vectors to be
able to use `expmv`.
"""
function ρ(d::AbstractVector{N}, em::ExponentialMap{N}) where {N<:Real}
    d_dense = d isa Vector ? d : Vector(d)
    v = expmv(one(N), transpose(em.spmexp.M), d_dense) # v <- exp(M^T) * d
    return ρ(v, em.X)
end

"""
    ∈(x::AbstractVector{N}, em::ExponentialMap{N})::Bool where {N<:Real}

Check whether a given point is contained in an exponential map of a convex set.

### Input

- `x`  -- point/vector
- `em` -- exponential map of a convex set

### Output

`true` iff ``x ∈ em``.

### Algorithm

This implementation exploits that ``x ∈ \\exp(M)⋅S`` iff ``\\exp(-M)⋅x ∈ S``.
This follows from ``\\exp(-M)⋅\\exp(M) = I`` for any ``M``.

### Examples

```jldoctest
julia> using SparseArrays

julia> em = ExponentialMap(
        SparseMatrixExp(sparse([1, 2], [1, 2], [2.0, 1.0], 2, 2)),
        BallInf([1., 1.], 1.));

julia> ∈([-1.0, 1.0], em)
false
julia> ∈([1.0, 1.0], em)
true
```
"""
function ∈(x::AbstractVector{N}, em::ExponentialMap{N})::Bool where {N<:Real}
    @assert length(x) == dim(em)
    return ∈(expmv(-one(N), em.spmexp.M, x), em.X)
end

"""
    vertices_list(em::ExponentialMap{N})::Vector{Vector{N}} where {N<:Real}

Return the list of vertices of a (polytopic) exponential map.

### Input

- `em` -- exponential map

### Output

A list of vertices.

### Algorithm

We assume that the underlying set `X` is polytopic.
Then the result is just the exponential map applied to the vertices of `X`.
"""
function vertices_list(em::ExponentialMap{N})::Vector{Vector{N}} where {N<:Real}
    # collect low-dimensional vertices lists
    vlist_X = vertices_list(em.X)

    # create resulting vertices list
    vlist = Vector{Vector{N}}()
    sizehint!(vlist, length(vlist_X))
    for v in vlist_X
        push!(vlist, expmv(one(N), em.spmexp.M, v))
    end

    return vlist
end

end end  # quote / load_expokit_exponentialmap

"""
    isbounded(em::ExponentialMap)::Bool

Determine whether an exponential map is bounded.

### Input

- `em` -- exponential map

### Output

`true` iff the exponential map is bounded.
"""
function isbounded(em::ExponentialMap)::Bool
    return isbounded(em.X)
end

"""
    isempty(em::ExponentialMap)::Bool

Return if an exponential map is empty or not.

### Input

- `em` -- exponential map

### Output

`true` iff the wrapped set is empty.
"""
function isempty(em::ExponentialMap)::Bool
    return isempty(em.X)
end

# --- ProjectionSparseMatrixExp & ExponentialProjectionMap ---

"""
    ProjectionSparseMatrixExp{N<:Real}

Type that represents the projection of a sparse matrix exponential, i.e.,
``L⋅\\exp(M)⋅R`` for a given sparse matrix ``M``.

### Fields

- `L` -- left multiplication matrix
- `E` -- sparse matrix exponential
- `R` -- right multiplication matrix
"""
struct ProjectionSparseMatrixExp{N<:Real, MN1<:AbstractSparseMatrix{N},
                                 MN2<:AbstractSparseMatrix{N},
                                 MN3<:AbstractSparseMatrix{N}}
    L::MN1
    spmexp::SparseMatrixExp{N, MN2}
    R::MN3
end

"""
    ExponentialProjectionMap{N<:Real, S<:LazySet{N}} <: LazySet{N}

Type that represents the application of a projection of a sparse matrix
exponential to a convex set.

### Fields

- `spmexp` -- projection of a sparse matrix exponential
- `X`      -- convex set
"""
struct ExponentialProjectionMap{N<:Real, S<:LazySet{N}} <: LazySet{N}
    projspmexp::ProjectionSparseMatrixExp
    X::S
end

"""
```
    *(projspmexp::ProjectionSparseMatrixExp,
      X::LazySet)::ExponentialProjectionMap
```

Return the application of a projection of a sparse matrix exponential to a
convex set.

### Input

- `projspmexp` -- projection of a sparse matrix exponential
- `X`          -- convex set

### Output

The application of the projection of a sparse matrix exponential to the convex
set.
"""
function *(projspmexp::ProjectionSparseMatrixExp,
           X::LazySet)::ExponentialProjectionMap
    return ExponentialProjectionMap(projspmexp, X)
end

"""
    dim(eprojmap::ExponentialProjectionMap)::Int

Return the dimension of a projection of an exponential map.

### Input

- `eprojmap` -- projection of an exponential map

### Output

The ambient dimension of the projection of an exponential map.
"""
function dim(eprojmap::ExponentialProjectionMap)::Int
    return size(eprojmap.projspmexp.L, 1)
end

function load_expokit_exponentialprojectionmap()
return quote

using .Expokit: expmv

"""
    σ(d::AbstractVector{N},
      eprojmap::ExponentialProjectionMap{N}) where {N<:Real}

Return the support vector of a projection of an exponential map.

### Input

- `d`        -- direction
- `eprojmap` -- projection of an exponential map

### Output

The support vector in the given direction.
If the direction has norm zero, the result depends on the wrapped set.

### Notes

If ``S = (L⋅M⋅R)⋅X``, where ``L`` and ``R`` are matrices, ``M`` is a matrix
exponential, and ``X`` is a set, it follows that
``σ(d, S) = L⋅M⋅R⋅σ(R^T⋅M^T⋅L^T⋅d, X)`` for any direction ``d``.

We allow sparse direction vectors, but will convert them to dense vectors to be
able to use `expmv`.
"""
function σ(d::AbstractVector{N},
           eprojmap::ExponentialProjectionMap{N}) where {N<:Real}
    d_dense = d isa Vector ? d : Vector(d)
    daux = transpose(eprojmap.projspmexp.L) * d_dense
    aux1 = expmv(one(N), transpose(eprojmap.projspmexp.spmexp.M), daux)
    daux = _At_mul_B(eprojmap.projspmexp.R, aux1)
    svec = σ(daux, eprojmap.X)

    aux2 = eprojmap.projspmexp.R * svec
    daux = expmv(one(N), eprojmap.projspmexp.spmexp.M, aux2)
    return eprojmap.projspmexp.L * daux
end

end end  # quote / load_expokit_exponentialprojectionmap

"""
    isbounded(eprojmap::ExponentialProjectionMap)::Bool

Determine whether an exponential projection map is bounded.

### Input

- `eprojmap` -- exponential projection map

### Output

`true` iff the exponential projection map is bounded.

### Algorithm

We first check if the left or right projection matrix is zero or the wrapped set
is bounded.
Otherwise, we check boundedness via [`isbounded_unit_dimensions`](@ref).
"""
function isbounded(eprojmap::ExponentialProjectionMap)::Bool
    if iszero(eprojmap.projspmexp.L) || iszero(eprojmap.projspmexp.R) ||
            isbounded(eprojmap.X)
        return true
    end
    return isbounded_unit_dimensions(eprojmap)
end

"""
    isempty(eprojmap::ExponentialProjectionMap)::Bool

Return if an exponential projection map is empty or not.

### Input

- `eprojmap` -- exponential projection map

### Output

`true` iff the wrapped set is empty.
"""
function isempty(eprojmap::ExponentialProjectionMap)::Bool
    return isempty(eprojmap.X)
end


text:
function Base.convert(::Type{StateSpace}, t::TransferFunction)
    if !isproper(t)
        error("System is improper, a state-space representation is impossible")
    end
    ny, nu = size(t)
    mat = t.matrix
    # TODO : These are added due to scoped for blocks, but is a hack. This
    # could be much cleaner.
    Ac = Bc = Cc = Dc = A = B = C = D = Array(Float64, 0, 0)
    for i=1:nu
        for j=1:ny
            a, b, c, d = siso_tf_to_ss(mat[j, i])
            if j > 1
                # vcat
                Ac = blkdiag(Ac, a)
                Bc = vcat(Bc, b)
                Cc = blkdiag(Cc, c)
                Dc = vcat(Dc, d)
            else
                Ac, Bc, Cc, Dc = (a, b, c, d)
            end
        end
        if i > 1
            # hcat
            A = blkdiag(A, Ac)
            B = blkdiag(B, Bc)
            C = hcat(C, Cc)
            D = hcat(D, Dc)
        else
            A, B, C, D = (Ac, Bc, Cc, Dc)
        end
    end
    A, B, C = balance_statespace(A, B, C)[1:3]
    return ss(A, B, C, D, t.Ts, inputnames=t.inputnames, outputnames=t.outputnames)
end

Base.convert(::Type{StateSpace}, t::Real) = ss(t)

Base.promote_rule{T<:SisoTf}(::Type{StateSpace}, ::Type{TransferFunction{T}}) = StateSpace
Base.promote_rule{T<:Real}(::Type{StateSpace}, ::Type{T}) = StateSpace

siso_tf_to_ss(t::SisoTf) = siso_tf_to_ss(convert(SisoRational, t))

function siso_tf_to_ss(t::SisoRational)
    t = normalize_tf(t)
    tnum = num(t)
    tden = den(t)
    len = length(tden)
    d = Array(Float64, 1, 1)
    d[1] = tnum[1]

    if len==1 || tnum == zero(Poly{Float64})
        a = zeros(0, 0)
        b = zeros(0, 1)
        c = zeros(1, 0)
    else
        tden = tden[2:end]
        a = [-tden' ; eye(len - 2, len - 1)]
        b = eye(len - 1, 1)
        c = tnum[2:len]' - d * tden[:]'
    end
    return float64mat(a), float64mat(b), float64mat(c), d
end

function normalize_tf(t::SisoRational)
    d = t.den[1]
    return SisoTf(t.num/d, t.den/d)
end

"""
`A, B, C, T = balance_statespace{S}(A::Matrix{S}, B::Matrix{S}, C::Matrix{S}, perm::Bool=false)`

`sys, T = balance_statespace(sys::StateSpace, perm::Bool=false)`

Computes a balancing transformation `T` that attempts to scale the system so
that the row and column norms of [T*A/T T*B; C/T 0] are approximately equal.
If `perm=true`, the states in `A` are allowed to be reordered.

This is not the same as finding a balanced realization with equal and diagonal observability and reachability gramians, see `balreal`
"""
function balance_statespace{S}(A::Matrix{S}, B::Matrix{S},
        C::Matrix{S}, perm::Bool=false)
    nx = size(A, 1)
    nu = size(B,2)
    ny = size(C,1)

    # Compute the transformation matrix
    mag_A = abs(A)
    mag_B = maximum([abs(B)  zeros(S, nx, 1)], 2)
    mag_C = maximum([abs(C); zeros(S, 1, nx)], 1)
    T = balance_transform(mag_A, mag_B, mag_C, perm)

    # Perform the transformation
    A = T*A/T
    B = T*B
    C = C/T

    return A, B, C, T
end

function balance_statespace(sys::StateSpace, perm::Bool=false)
    A, B, C, T = balance_statespace(sys.A,sys.B,sys.C, perm)
    return ss(A,B,C,sys.D), T
end

"""
`T = balance_transform{R}(A::Matrix{R}, B::Matrix{R}, C::Matrix{R}, perm::Bool=false)`

`T = balance_transform(sys::StateSpace, perm::Bool=false) = balance_transform(A,B,C,perm)`

Computes a balancing transformation `T` that attempts to scale the system so
that the row and column norms of [T*A/T T*B; C/T 0] are approximately equal.
If `perm=true`, the states in `A` are allowed to be reordered.

This is not the same as finding a balanced realization with equal and diagonal observability and reachability gramians, see `balreal`
See also `balance_statespace`, `balance`
"""
function balance_transform{R}(A::Matrix{R}, B::Matrix{R}, C::Matrix{R}, perm::Bool=false)
    nx = size(A, 1)
    # Compute a scaling of the system matrix M
    S = diag(balance([A B; C zeros(size(C*B))], false)[1])
    Sx = S[1:nx]
    Sio = S[nx+1]
    # Compute permutation of x (if requested)
    pvec = perm ? balance(A, true)[2] * [1:nx;] : [1:nx;]
    # Compute the transformation matrix
    T = zeros(R, nx, nx)
    T[pvec, :] = Sio * diagm(1./Sx)
    return T
end

balance_transform(sys::StateSpace, perm::Bool=false) = balance_transform(sys.A,sys.B,sys.C,perm)

@doc """`sys = ss2tf(s::StateSpace)`, ` sys = ss2tf(A, B, C, Ts = 0; inputnames = "", outputnames = "")`

Convert a `StateSpace` realization to a `TransferFunction`""" ->
function ss2tf(s::StateSpace)
    return ss2tf(s.A, s.B, s.C, s.D, s.Ts, inputnames=s.inputnames, outputnames=s.outputnames)
end

function ss2tf(A, B, C, D, Ts = 0; inputnames = "", outputnames = "")
    nu,ny = size(B,2),size(C,1)
    ubernum = Matrix{Vector}(ny,nu)
    uberden = Matrix{Vector}(ny,nu)
    for i = 1:nu, j=1:ny
        ubernum[j,i],uberden[j,i] = sisoss2tf(A, B[:,i], C[j,:]', D[j,i])
    end
    tf(ubernum,uberden, Ts, inputnames=inputnames, outputnames=outputnames)
end

function sisoss2tf(A, B, C, D)
    charpolA = charpoly(A)
    numP = charpoly(A-B*C) - charpolA + D*charpolA
    denP = charpolA
    return numP[1:length(numP)], denP[1:length(denP)]
end

tf(sys::StateSpace) = ss2tf(sys)
zpk(sys::StateSpace) = zpk(ss2tf(sys))

function charpoly(A)
    λ = eigvals(A);
    p = reduce(*,ControlSystems.Poly([1.]), ControlSystems.Poly[ControlSystems.Poly([1, -λᵢ]) for λᵢ in λ]);
    if maximum(imag(p[:])./(1+abs(real(p[:])))) < sqrt(eps())
        for i = 1:length(p)
            p[i] = real(p[i])
        end
    else
        error("Characteristic polynomial should be real")
    end
    p
end


text:
using Revise
using ADCME
using ADCMEKit
using NNFEM
using PyPlot
using ProgressMeter
close("all")
include("load_domain_function.jl")

NT = 100
Δt = 20/NT

## solve static elastic response
domain = load_crack_domain()
H = domain.elements[1].mat[1].H
EBC_func = nothing
FBC_func = nothing
body_func = nothing
globaldata = GlobalData(missing, missing, missing, missing, domain.neqs, EBC_func, FBC_func, body_func)
Fext = compute_external_force(globaldata, domain)
d = LinearStaticSolver(globaldata, domain, domain.state, H, Fext)
sess = Session(); init(sess)
domain.state = run(sess, d)

figure()
visualize_boundary(domain)
figure()
subplot(211)
visualize_scalar_on_scoped_body(domain.state[1:domain.nnodes], zeros(size(domain.state)...), domain)
subplot(212)
visualize_scalar_on_scoped_body(domain.state[domain.nnodes+1:end], zeros(size(domain.state)...), domain)

## solve viscoelastic relaxation
strain = getStrain(domain)
stress = getStress(domain)
state = domain.state

domain = load_crack_domain(option="mixed")
ts = GeneralizedAlphaSolverTime(Δt, NT)
ubd, abd = compute_boundary_info(domain, globaldata, ts)
Fext = compute_external_force(domain, globaldata, ts)

d0 = state
v0 = zeros(2domain.nnodes)
a0  = zeros(2domain.nnodes)
# σ0 = zeros(getNGauss(domain),3)
# ϵ0 = zeros(getNGauss(domain),3)
σ0 = stress
ϵ0 = strain

μ = zeros(getNGauss(domain))
λ = zeros(getNGauss(domain))
η = zeros(getNGauss(domain))
k = 0
for i = 1:domain.neles
  e = domain.elements[i]
  for mat in e.mat
    global k += 1
    μ[k] = mat.μ
    λ[k] = mat.λ
    η[k] = mat.η
  end
end

d, v, a, σ, ϵ = ViscoelasticitySolver(
  globaldata, domain, d0, v0, a0, σ0, ϵ0, Δt, NT, μ, λ, η, Fext, ubd, abd
)

sess = Session(); init(sess)

d_, σ_ = run(sess, [d, σ])
domain.history["state"] = []
domain.history["stress"] = []
for i = 1:NT+1
  push!(domain.history["state"], d_[i,:])
  push!(domain.history["stress"], σ_[i,:,:])
end

### visulalize results

# figure()
# p = visualize_scalar_on_scoped_body(σ_[:, 1:domain.nnodes, 2], d_, domain, scale_factor=1.0)
# # p = visualize_von_mises_stress_on_scoped_body(d_, domain, scale_factor=10.0)
# # p = visualize_total_deformation_on_scoped_body(d_, domain, scale_factor=10.0)
# saveanim(p, "ad_solver.gif")
# close("all")

###################### check surface ########################
# idx = []
# k = 0
# pts = getGaussPoints(domain)
# xcoord = Float64[]
# for e in domain.elements
#     global k
#     if sum(map(x->x≈0.0, e.coords[:,2]))≠2
#         k += length(e.weights)
#         continue
#     end
#     l = argmin(pts[k+1:k+length(e.weights), 2])
#     push!(idx, k + l)
#     push!(xcoord, pts[k+l, 1])
#     k += length(e.weights)
# end
# ii = sortperm(xcoord)
# xcoord = xcoord[ii]
# idx = idx[ii]

# σ_surface = σ_[:, idx, :]
# s = zeros(size(σ_surface, 1), size(σ_surface, 2))
# for i = 1:size(σ_surface, 1)
#     for j = 1:size(σ_surface, 2)
#         s[i,j] = NNFEM.postprocess_stress(σ_surface[i,j,:], "vonMises")
#     end
# end

# figure(figsize=(8,3))
# semilogy(s[1, :], "-", label="x0")
# semilogy(s[NT+1, :], "--", label="x")
# legend()
# title("Stress")
# savefig("Stress_test.png")

# figure(figsize=(8,3))
# semilogy(s[:, 13], "-")
# title("Stress")

slip_idx = findall(domain.nodes[:,2] .≈ 0.0)
coords = domain.nodes[slip_idx,1]
ii = sortperm(coords)
coords = coords[ii]
x_ = d_[:, slip_idx[ii]]
y_ = d_[:, slip_idx[ii] .+ domain.nnodes]

x0 = d0[slip_idx[ii]]
y0 = d0[slip_idx[ii] .+ domain.nnodes]

figure(figsize=(8,3))
subplot(211)
plot((x_' .- x0)[:,1:10:end])
subplot(212)
plot(x0, label="x")
legend()

figure(figsize=(8,3))
subplot(211)
S = 1:20:NT
for i = 1:length(S)
  plot(-(y_' .- y0)[:,S[i]],"C$i",  label="$i")
end
legend()

subplot(212)
S = 1:20:NT
for i = 1:length(S)
  plot(-y_'[:,S[i]],"C$i",  label="$i")
end
plot(-y0, label="y")
# plot(y_[1, :] .- y0, "-", label="x0")
# plot(y_[NT+1, :] .- y0, "--", label="x")
# plot(y0, ":", label="static")
legend()
savefig("disp.png")

figure()
plot(y_[:, 10])
# plot(y_[1, :] .- y0, "-", label="x0")
# plot(y_[NT+1, :] .- y0, "--", label="x")
# plot(y0, ":", label="static")
legend()

max_stars_repo_name: ADCMEMarket/ADCMEImages

max_issues_repo_name: ADCMEMarket/ADCMEImages

max_forks_repo_name: ADCMEMarket/ADCMEImages


text:
module HybridFunnelGraphs

using LinearAlgebra
using Combinatorics
using HPD
using Julog

include("types.jl")
# include("graph.jl")
include("hpd_graph.jl")

export create_funnel_graph,
        get_preconditions,
        get_effects,
        Graph
#        Region,
#        intersects
end


text:
function print(lp::StandardLpData)
    println("Number of rows   : $(lp.nrows)")
    println("Number of columns: $(lp.ncols)")
    println("Constraint matrix: $(Matrix(lp.A))")
    println("Row lower bounds: $(lp.bl)")
    println("Row upper bounds: $(lp.bu)")
    println("Objective coefficients: $(lp.c)")
    println("Column lower bound: $(lp.xl)")
    println("Column upper bound: $(lp.xu)")
end

function print(lp::CanonicalLpData)
    println("Number of rows   : $(lp.nrows)")
    println("Number of columns: $(lp.ncols)")
    println("Constraint matrix: $(Matrix(lp.A))")
    println("Row bounds: $(lp.b)")
    println("Objective coefficients: $(lp.c)")
    println("Column lower bound: $(lp.xl)")
    println("Column upper bound: $(lp.xu)")
end

function summary(lp::StandardLpData)
    println("Number of rows   : $(lp.nrows)")
    println("Number of columns: $(lp.ncols)")
    # println("A: Min $(minimum(nonzeros(lp.A))), Max $(maximum(nonzeros(lp.A)))")
    println("bl: Min $(minimum(lp.bl)), Max $(maximum(lp.bl))")
    println("bu: Min $(minimum(lp.bu)), Max $(maximum(lp.bu))")
    println("xl: Min $(minimum(lp.xl)), Max $(maximum(lp.xl))")
    println("xu: Min $(minimum(lp.xu)), Max $(maximum(lp.xu))")
    println("c: Min $(minimum(lp.c)), Max $(maximum(lp.c))")

    cons = Dict("EQ"=>0, "GT"=>0, "LT"=>0, "RNG"=>0, "FR"=>0)
    vars = Dict("EQ"=>0, "GT"=>0, "LT"=>0, "RNG"=>0, "FR"=>0)
    for i = 1:lp.nrows
        if lp.bl[i] == lp.bu[i]
            cons["EQ"] += 1
        elseif lp.bl[i] > -Inf
            if lp.bu[i] < Inf
                cons["RNG"] += 1
            else
                cons["GT"] += 1
            end
        elseif lp.bu[i] <= Inf
            cons["LT"] += 1
        else
            cons["FR"] += 1
        end
    end

    for j = 1:lp.ncols
        if lp.xl[j] == lp.xu[j]
            vars["EQ"] += 1
        elseif lp.xl[j] > -Inf
            if lp.xu[j] < Inf
                vars["RNG"] += 1
            else
                vars["GT"] += 1
            end
        elseif lp.xu[j] <= Inf
            vars["LT"] += 1
        else
            vars["FR"] += 1
        end
    end

    println("Constraints:")
    for (k,v) = cons
        println(" $k = $v")
    end
    println("Variables:")
    for (k,v) = vars
        println(" $k = $v")
    end

    # for j = 1:lp.ncols
    #     if lp.xl[j]
    # end
    # println("Free: ")
end


text:
function __init__()
    set_num_threads(Threads.nthreads())
end


text:
module RegisterMismatch

import Base: copy, eltype, isnan, ndims

using ImageCore
using RFFT, FFTW
using RegisterCore, PaddedViews, MappedArrays
using Printf
using Reexport
@reexport using RegisterMismatchCommon
import RegisterMismatchCommon: mismatch0, mismatch, mismatch_apertures

export
    CMStorage,
    fillfixed!,
    mismatch0,
    mismatch,
    mismatch!,
    mismatch_apertures,
    mismatch_apertures!

"""
The major types and functions exported are:

- `mismatch` and `mismatch!`: compute the mismatch between two images
- `mismatch_apertures` and `mismatch_apertures!`: compute apertured mismatch between two images
- `mismatch0`: simple direct mismatch calculation with no shift
- `nanpad`: pad the smaller image with NaNs
- `highpass`: highpass filter an image
- `correctbias!`: replace corrupted mismatch data (due to camera bias inhomogeneity) with imputed data
- `truncatenoise!`: threshold mismatch computation to prevent problems from roundoff
- `aperture_grid`: create a regular grid of apertures
- `allocate_mmarrays`: create storage for output of `mismatch_apertures!`
- `CMStorage`: a type that facilitates re-use of intermediate storage during registration computations
"""
RegisterMismatch

FFTW.set_num_threads(min(Sys.CPU_THREADS, 8))
set_FFTPROD([2,3])

mutable struct NanCorrFFTs{T<:AbstractFloat,N,RCType<:RCpair{T,N}}
    I0::RCType
    I1::RCType
    I2::RCType
end

copy(x::NanCorrFFTs) = NanCorrFFTs(copy(x.I0), copy(x.I1), copy(x.I2))

"""
    CMStorage{T}(undef, aperture_width, maxshift; flags=FFTW.ESTIMATE, timelimit=Inf, display=true)

Prepare for FFT-based mismatch computations over domains of size `aperture_width`, computing the
mismatch up to shifts of size `maxshift`.  The keyword arguments allow you to control the planning
process for the FFTs.
"""
mutable struct CMStorage{T<:AbstractFloat,N,RCType<:RCpair{T,N},FFT<:Function,IFFT<:Function}
    aperture_width::Vector{Float64}
    maxshift::Vector{Int}
    getindices::Vector{UnitRange{Int}}   # indices for pulling padded data, in source-coordinates
    padded::Array{T,N}
    fixed::NanCorrFFTs{T,N,RCType}
    moving::NanCorrFFTs{T,N,RCType}
    buf1::RCType
    buf2::RCType
    # the next two store the result of calling plan_fft! and plan_ifft!
    fftfunc!::FFT
    ifftfunc!::IFFT
    shiftindices::Vector{Vector{Int}} # indices for performing fftshift & snipping from -maxshift:maxshift
end

function CMStorage{T,N}(::UndefInitializer, aperture_width::NTuple{N,<:Real}, maxshift::Dims{N}; flags=FFTW.ESTIMATE, timelimit=Inf, display=true) where {T,N}
    blocksize = map(x->ceil(Int,x), aperture_width)
    padsz = padsize(blocksize, maxshift)
    padded = Array{T}(undef, padsz)
    getindices = padranges(blocksize, maxshift)
    maxshiftv = [maxshift...]
    region = findall(maxshiftv .> 0)
    fixed  = NanCorrFFTs(RCpair{T}(undef, padsz, region), RCpair{T}(undef, padsz, region), RCpair{T}(undef, padsz, region))
    moving = NanCorrFFTs(RCpair{T}(undef, padsz, region), RCpair{T}(undef, padsz, region), RCpair{T}(undef, padsz, region))
    buf1 = RCpair{T}(undef, padsz, region)
    buf2 = RCpair{T}(undef, padsz, region)
    tcalib = 0
    if display && flags != FFTW.ESTIMATE
        print("Planning FFTs (maximum $timelimit seconds)...")
        flush(stdout)
        tcalib = time()
    end
    fftfunc = plan_rfft!(fixed.I0, flags=flags, timelimit=timelimit/2)
    ifftfunc = plan_irfft!(fixed.I0, flags=flags, timelimit=timelimit/2)
    if display && flags != FFTW.ESTIMATE
        dt = time()-tcalib
        @printf("done (%.2f seconds)\n", dt)
    end
    shiftindices = Vector{Int}[ [size(padded,i).+(-maxshift[i]+1:0); 1:maxshift[i]+1] for i = 1:length(maxshift) ]
    CMStorage{T,N,typeof(buf1),typeof(fftfunc),typeof(ifftfunc)}(Float64[aperture_width...], maxshiftv, getindices, padded, fixed, moving, buf1, buf2, fftfunc, ifftfunc, shiftindices)
end

CMStorage{T}(::UndefInitializer, aperture_width::NTuple{N,<:Real}, maxshift::Dims{N}; kwargs...) where {T<:Real,N} =
    CMStorage{T,N}(undef, aperture_width, maxshift; kwargs...)

eltype(cms::CMStorage{T,N}) where {T,N} = T
 ndims(cms::CMStorage{T,N}) where {T,N} = N

"""
    mm = mismatch([T], fixed, moving, maxshift; normalization=:intensity)

Compute the mismatch between `fixed` and
`moving` as a function of translations (shifts) up to size `maxshift`.
Optionally specify the element-type of the mismatch arrays (default
`Float32` for Integer- or FixedPoint-valued images) and the
normalization scheme (`:intensity` or `:pixels`).

`fixed` and `moving` must have the same size; you can pad with
`NaN`s as needed. See `nanpad`.
"""
function mismatch(::Type{T}, fixed::AbstractArray, moving::AbstractArray, maxshift::DimsLike; normalization = :intensity) where T<:Real
    msz = 2 .* maxshift .+ 1
    mm = MismatchArray(T, msz...)
    cms = CMStorage{T}(undef, size(fixed), maxshift)
    fillfixed!(cms, fixed)
    erng = shiftrange.((cms.getindices...,), first.(axes(fixed)) .- 1)  # expanded rng
    mpad = PaddedView(convert(T, NaN), of_eltype(T, moving), erng)
    mismatch!(mm, cms, mpad, normalization=normalization)
    return mm
end

"""
`mismatch!(mm, cms, moving; [normalization=:intensity])`
computes the mismatch as a function of shift, storing the result in
`mm`. The `fixed` image has been prepared in `cms`, a `CMStorage` object.
"""
function mismatch!(mm::MismatchArray, cms::CMStorage, moving::AbstractArray; normalization = :intensity)
    # Pad the moving snippet using any available data, including
    # regions that might be in the parent Array but are not present
    # within the boundaries of the SubArray. Use NaN only for pixels
    # truly lacking data.
    checksize_maxshift(mm, cms.maxshift)
    copyto!(cms.padded, CartesianIndices(cms.padded), moving, CartesianIndices(moving))
    fftnan!(cms.moving, cms.padded, cms.fftfunc!)
    # Compute the mismatch
    f0 = complex(cms.fixed.I0)
    f1 = complex(cms.fixed.I1)
    f2 = complex(cms.fixed.I2)
    m0 = complex(cms.moving.I0)
    m1 = complex(cms.moving.I1)
    m2 = complex(cms.moving.I2)
    tnum   = complex(cms.buf1)
    tdenom = complex(cms.buf2)
    if normalization == :intensity
        for i = 1:length(tnum)
            c = 2*conj(f1[i])*m1[i]
            q = conj(f2[i])*m0[i] + conj(f0[i])*m2[i]
            tdenom[i] = q
            tnum[i] = q - c
        end
    elseif normalization == :pixels
        for i = 1:length(tnum)
            tdenom[i] = conj(f0[i])*m0[i]
            tnum[i] = conj(f2[i])*m0[i] - 2*conj(f1[i])*m1[i] + conj(f0[i])*m2[i]
        end
    else
        error("normalization $normalization not recognized")
    end
    cms.ifftfunc!(cms.buf1)
    cms.ifftfunc!(cms.buf2)
    copyto!(mm, (view(real(cms.buf1), cms.shiftindices...), view(real(cms.buf2), cms.shiftindices...)))
    mm
end

"""
`mms = mismatch_apertures([T], fixed, moving, gridsize, maxshift;
[normalization=:pixels], [flags=FFTW.MEASURE], kwargs...)` computes
the mismatch between `fixed` and `moving` over a regularly-spaced grid
of aperture centers, effectively breaking the images up into
chunks. The maximum-allowed shift in any aperture is `maxshift`.

`mms = mismatch_apertures([T], fixed, moving, aperture_centers,
aperture_width, maxshift; kwargs...)` computes the mismatch between
`fixed` and `moving` over a list of apertures of size `aperture_width`
at positions defined by `aperture_centers`.

`fixed` and `moving` must have the same size; you can pad with `NaN`s
as needed to ensure this.  You can optionally specify the real-valued
element type mm; it defaults to the element type of `fixed` and
`moving` or, for Integer- or FixedPoint-valued images, `Float32`.

On output, `mms` will be an Array-of-MismatchArrays, with the outer
array having the same "grid" shape as `aperture_centers`.  The centers
can in general be provided as an vector-of-tuples, vector-of-vectors,
or a matrix with each point in a column.  If your centers are arranged
in a rectangular grid, you can use an `N`-dimensional array-of-tuples
(or array-of-vectors) or an `N+1`-dimensional array with the center
positions specified along the first dimension. See `aperture_grid`.
"""
function mismatch_apertures(::Type{T},
                            fixed::AbstractArray,
                            moving::AbstractArray,
                            aperture_centers::AbstractArray,
                            aperture_width::WidthLike,
                            maxshift::DimsLike;
                            normalization = :pixels,
                            flags = FFTW.MEASURE,
                            kwargs...) where T
    nd = sdims(fixed)
    (length(aperture_width) == nd && length(maxshift) == nd) || error("Dimensionality mismatch")
    mms = allocate_mmarrays(T, aperture_centers, maxshift)
    cms = CMStorage{T}(undef, aperture_width, maxshift; flags=flags, kwargs...)
    mismatch_apertures!(mms, fixed, moving, aperture_centers, cms, normalization=normalization)
end

"""
`mismatch_apertures!(mms, fixed, moving, aperture_centers, cms;
[normalization=:pixels])` computes the mismatch between `fixed` and
`moving` over a list of apertures at positions defined by
`aperture_centers`.  The parameters and working storage are contained
in `cms`, a `CMStorage` object. The results are stored in `mms`, an
Array-of-MismatchArrays which must have length equal to the number of
aperture centers.
"""
function mismatch_apertures!(mms, fixed, moving, aperture_centers, cms::CMStorage{T}; normalization=:pixels) where T
    N = ndims(cms)
    fillvalue = convert(T, NaN)
    getinds = (cms.getindices...,)::NTuple{ndims(fixed),UnitRange{Int}}
    fixedT, movingT = of_eltype(T, fixed), of_eltype(T, moving)
    for (mm,center) in zip(mms, each_point(aperture_centers))
        rng = aperture_range(center, cms.aperture_width)
        fsnip = PaddedView(fillvalue, fixedT, rng)
        erng = shiftrange.(getinds, first.(rng) .- 1)  # expanded rng
        msnip = PaddedView(fillvalue, movingT, erng)
        # Perform the calculation
        fillfixed!(cms, fsnip)
        mismatch!(mm, cms, msnip; normalization=normalization)
    end
    mms
end

# Calculate the components needed to "nancorrelate"
function fftnan!(out::NanCorrFFTs{T}, A::AbstractArray{T}, fftfunc!::Function) where T<:Real
    I0 = real(out.I0)
    I1 = real(out.I1)
    I2 = real(out.I2)
    _fftnan!(parent(I0), parent(I1), parent(I2), A)
    fftfunc!(out.I0)
    fftfunc!(out.I1)
    fftfunc!(out.I2)
    out
end

function _fftnan!(I0, I1, I2, A::AbstractArray{T}) where T<:Real
    @inbounds for i in CartesianIndices(size(A))
        a = A[i]
        f = !isnan(a)
        I0[i] = f
        af = f ? a : zero(T)
        I1[i] = af
        I2[i] = af*af
    end
end

function fillfixed!(cms::CMStorage{T}, fixed::AbstractArray) where T
    fill!(cms.padded, NaN)
    pinds = CartesianIndices(ntuple(d->(1:size(fixed,d)).+cms.maxshift[d], ndims(fixed)))
    copyto!(cms.padded, pinds, fixed, CartesianIndices(fixed))
    fftnan!(cms.fixed, cms.padded, cms.fftfunc!)
end

#### Utilities

Base.isnan(A::Array{Complex{T}}) where {T} = isnan(real(A)) | isnan(imag(A))
function sumsq_finite(A)
    s = 0.0
    for a in A
        if isfinite(a)
            s += a*a
        end
    end
    if s == 0
        error("No finite values available")
    end
    s
end

### Deprecations

function CMStorage{T}(::UndefInitializer, aperture_width::WidthLike, maxshift::WidthLike; kwargs...) where {T<:Real}
    Base.depwarn("CMStorage with aperture_width::$(typeof(aperture_width)) and maxshift::$(typeof(maxshift)) is deprecated, use tuples instead", :CMStorage)
    (N = length(aperture_width)) == length(maxshift) || error("Dimensionality mismatch")
    return CMStorage{T,N}(undef, (aperture_width...,), (maxshift...,); kwargs...)
end

@deprecate CMStorage(::Type{T}, aperture_width, maxshift; kwargs...) where T   CMStorage{T}(undef, aperture_width, maxshift; kwargs...)

end


text:
abstract type IntegratorRK{dType, tType} <: DeterministicIntegrator{dType, tType} end
abstract type IntegratorPRK{dType, tType} <: IntegratorRK{dType, tType} end

equation(integrator::IntegratorRK) = integrator.params.equ
timestep(integrator::IntegratorRK) = integrator.params.Δt
tableau(integrator::IntegratorRK)  = integrator.params.tab
nstages(integrator::IntegratorRK)  = integrator.params.tab.s


text:
#
# Correctness Tests
#

my_tests = ["test/bayes.jl"]

@printf "Running tests:\n"

for my_test in my_tests
    @printf " * %s\n" my_test
    include(my_test)
end