text:
# simple adjacency list, which represents both nodes and timestamps as integers.

"""
    IntAdjacencyList(nv, nt)

Construct a graph represented by an adjacency list with `nv` nodes and `nt` timestamps, where both nodes and timestamps are represented by integers.

# Example

```jldoctest
julia> using EvolvingGraphs

julia> g = IntAdjacencyList(4,3)
Directed IntAdjacencyList (4 nodes, 0 static edges, 3 timestamps)

julia> add_edge!(g, 1, 2, 1)
Directed IntAdjacencyList (4 nodes, 1 static edges, 3 timestamps)

julia> add_edge!(g, 2, 3, 2)
Directed IntAdjacencyList (4 nodes, 2 static edges, 3 timestamps)

julia> num_edges(g)
2
```
"""
mutable struct IntAdjacencyList <: AbstractEvolvingGraph{Int, Tuple{Int,Int}, Int}
    is_directed::Bool
    nodes::UnitRange{Int}
    timestamps::Vector{Int}
    nnodes::Int      # number of nodes
    nedges::Int      # number of static edges
    forward_adjlist::Vector{Vector{Int}}
    backward_adjlist::Vector{Vector{Int}}
end
function IntAdjacencyList(nv::Int, nt::Int; is_directed::Bool = true)
    ts = Array{Int}(nv*nt)
    f_adj = Vector{Int}[]
    b_adj = Vector{Int}[]
    for i = 1:nv*nt
        push!(f_adj, Int[])
        push!(b_adj, Int[])
    end
    for t = 1:nt
        for v = 1:nv
            ts[v + nv*(t-1)] = t
        end
    end
    IntAdjacencyList(is_directed, 1:nv*nt, ts, nv, 0, f_adj, b_adj)
end

is_directed(g::IntAdjacencyList) = g.is_directed

"""
    evolving_graph_to_adj(g)

Convert an evolving graph `g` to an adjacency list.

# Example

```jldoctest
julia> using EvolvingGraphs

julia> g = EvolvingGraph()
Directed EvolvingGraph 0 nodes, 0 static edges, 0 timestamps

julia> add_bunch_of_edges!(g, [(1,2,2001),(2,3,2002), (2,4,2002), (3,1,2003)])
Directed EvolvingGraph 4 nodes, 4 static edges, 3 timestamps

julia> evolving_graph_to_adj(g)
Directed IntAdjacencyList (4 nodes, 4 static edges, 3 timestamps)
```
"""
function evolving_graph_to_adj(g::AbstractEvolvingGraph)
    g1 = IntAdjacencyList(num_nodes(g), num_timestamps(g), is_directed = is_directed(g))

    tmap = Dict(t => i for (i, t) in enumerate(unique_timestamps(g)))
    for e in edges(g)
        v1 = node_index(source(e))
        v2 = node_index(target(e))
        add_edge!(g1, v1, v2, tmap[e.timestamp])
    end
    g1
end

nodes(g::IntAdjacencyList) = collect(1:g.nnodes)
num_nodes(g::IntAdjacencyList) = g.nnodes
num_edges(g::IntAdjacencyList) = g.nedges
timestamps(g::IntAdjacencyList) = g.timestamps
unique_timestamps(g::IntAdjacencyList) = unique(g.timestamps)
num_timestamps(g::IntAdjacencyList) = length(unique_timestamps(g))
node_index(g::IntAdjacencyList, node::Int) = node

deepcopy(g::IntAdjacencyList) = IntAdjacencyList(is_directed(g),
                                             deepcopy(g.nodes),
                                             deepcopy(g.timestamps),
                                             g.nnodes,
                                             g.nedges,
                                             deepcopy(g.forward_adjlist),
                                             deepcopy(g.backward_adjlist))

function forward_neighbors(g::IntAdjacencyList, v::Int, t::Int)
    ns = g.nnodes
    n = v + ns*(t-1)
    nn = Tuple{Int, Int}[]
    for node in g.forward_adjlist[n]
        t = g.timestamps[node]
        v1 = node - ns*(t-1)
        push!(nn, (v1, t))
    end
    nn
end
forward_neighbors(g::IntAdjacencyList, vt::Tuple{Int,Int}) = forward_neighbors(g, vt[1], vt[2])

function backward_neighbors(g::IntAdjacencyList, v::Int, t::Int)
    ns = g.nnodes
    n = v + ns*(t-1)
    nn = Tuple{Int, Int}[]
    for node in g.backward_adjlist[n]
        t = g.timestamps[node]
        v1 = node - ns*(t-1)
        push!(nn, (v1, t))
    end
    nn
end
backward_neighbors(g::IntAdjacencyList, vt::Tuple{Int,Int}) = backward_neighbors(g, vt[1], vt[2])

function add_edge!(g::IntAdjacencyList, v1::Int, v2::Int, t::Int)
    ns = g.nnodes
    n1 = v1 + ns*(t-1)
    n2 = v2 + ns*(t-1)
    if _insertnode!(g.forward_adjlist[n1], n2)
        g.nedges += 1

        for (v, n) in ((v1, n1), (v2, n2))
            for i in 1:t-1
                len = length(g.forward_adjlist[v + ns*(i-1)]) +
                length(g.backward_adjlist[v + ns*(i-1)])
                if len > 0
                    _insertnode!(g.forward_adjlist[v + ns*(i-1)], n)
                    _insertnode!(g.backward_adjlist[n], (v+ns*(i-1)))
                end
            end
            for i in t+1:num_timestamps(g)
                #   println("node: $(v + ns*(i-1))")
                len = length(g.forward_adjlist[v + ns*(i-1)]) +
                length(g.backward_adjlist[v + ns*(i-1)])
                if len > 0
                    _insertnode!(g.forward_adjlist[n], (v+ns*(i-1)))
                    _insertnode!(g.backward_adjlist[v+ns*(i-1)], n)
                end
            end
        end
    end
    _insertnode!(g.backward_adjlist[n2], n1)
    g
end

# The function _insertnode! is from LightGraphs.jl
# The LightGraphs.jl package is licensed under the Simplified "2-clause" BSD License:

#Copyright (c) 2015: Seth Bromberger and other contributors.

#Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:

#Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
#Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
#THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

# returns true if insert succeeded, false if it was a duplicate
_insertnode!(v::Vector{Int}, x::Int) = isempty(splice!(v, searchsorted(v,x), x))

# has_edge on active nodes
function _has_edge(g::IntAdjacencyList, n1::Int, n2::Int)
    nn = length(g.nodes)
    if n1 > nn || n2 > nn
        return false
    end
    return n2 in g.forward_adjlist[n1]
end

function has_edge(g::IntAdjacencyList, v1::Int, v2::Int, t::Int)
    ns = g.nnodes
    n1 = v1 + ns*(t-1)
    n2 = v2 + ns*(t-1)
    return _has_edge(g, n1, n2)
end

function adjacency_matrix(g::IntAdjacencyList, t::Int)
    ns = nodes(g)
    num_ns = num_nodes(g)
    A = zeros(Float64, num_ns, num_ns)
    for n in ns
        fns = forward_neighbors(g, n, t)
        for (fn, ft) in fns
            if ft == t
                 A[n, fn] = 1.0
             end
         end
     end
     return A
end

function sparse_adjacency_matrix(g::IntAdjacencyList, t::Int)
    is = Int[]
    js = Int[]
    ns = nodes(g)
    num_ns = num_nodes(g)
    for n in ns
        fns = forward_neighbors(g, n, t)
        for (fn, ft) in fns
            if ft == t
                push!(is, n)
                push!(js, fn)
            end
        end
    end
    return sparse(is, js, ones(is, Float64), num_ns, num_ns)
end

function block_adjacency_matrix(g::IntAdjacencyList)
    ns = nodes(g)
    num_ns = num_nodes(g)
    uts = unique_timestamps(g)
    num_ts = num_timestamps(g)
    is = Int[]
    js = Int[]
    vs = Float64[]
    for n in ns
        for t in uts
            fns = forward_neighbors(g, n, t)
            for (fn, ft) in fns
                push!(is, n + num_ns*(t-1))
                push!(js, fn + num_ns*(ft-1))
                abs(ft - t) == 0 ? push!(vs, 1.0) : push!(vs, 1./abs(ft-t))
            end
        end
    end
    return sparse(is, js, vs, num_ns*num_ts, num_ns*num_ts)
end


text:
function f(x::Vector)
    (x[1] - 5.0)^4
end

function g!(x::Vector, storage::Vector)
    storage[1] = 4.0 * (x[1] - 5.0)^3
end

function h!(x::Vector, storage::Matrix)
    storage[1, 1] = 12.0 * (x[1] - 5.0)^2
end

d = TwiceDifferentiableFunction(f, g!, h!)

results = Optim.newton(d, [0.0])
@assert length(results.trace.states) == 0
@assert results.converged
@assert norm(results.minimum - [5.0]) < 0.01

eta = 0.9

function f(x::Vector)
  (1.0 / 2.0) * (x[1]^2 + eta * x[2]^2)
end

function g!(x::Vector, storage::Vector)
  storage[1] = x[1]
  storage[2] = eta * x[2]
end

function h!(x::Vector, storage::Matrix)
  storage[1, 1] = 1.0
  storage[1, 2] = 0.0
  storage[2, 1] = 0.0
  storage[2, 2] = eta
end

d = TwiceDifferentiableFunction(f, g!, h!)
results = Optim.newton(d, [127.0, 921.0])
@assert length(results.trace.states) == 0
@assert results.converged
@assert norm(results.minimum - [0.0, 0.0]) < 0.01

ext: jl


text:
mutable struct CProps
  ai::Vector{Int}
  bi::Vector{Int}
  ci::Vector{Int}
  nactiveA::Int
  nactiveB::Int
  nactiveC::Int
  AtoB::Vector{Int}
  AtoC::Vector{Int}
  BtoC::Vector{Int}
  permuteA::Bool
  permuteB::Bool
  permuteC::Bool
  dleft::Int
  dmid::Int
  dright::Int
  ncont::Int
  Acstart::Int
  Bcstart::Int
  Austart::Int
  Bustart::Int
  PA::Vector{Int}
  PB::Vector{Int}
  PC::Vector{Int}
  ctrans::Bool
  newArange::Vector{Int}
  newBrange::Vector{Int}
  newCrange::Vector{Int}
  function CProps(ai::Vector{Int},
                  bi::Vector{Int},
                  ci::Vector{Int})
    new(ai,bi,ci,0,0,0,Vector{Int}(),Vector{Int}(),Vector{Int}(),false,false,false,1,1,1,0,
        length(ai),length(bi),length(ai),length(bi),Vector{Int}(),Vector{Int}(),Vector{Int}(),
        false,Vector{Int}(),Vector{Int}(),Vector{Int}())
  end
end

function compute_perms!(props::CProps)::Nothing
  #Use !AtoB.empty() as a check to see if we've already run this
  length(props.AtoB)!=0 && return

  na = length(props.ai)
  nb = length(props.bi)
  nc = length(props.ci)

  props.AtoB = fill(0,na)
  props.AtoC = fill(0,na)
  props.BtoC = fill(0,nb)
  for i = 1:na
    for j = 1:nb
      if props.ai[i]==props.bi[j]
        props.ncont += 1
        #TODO: check this if this should be i,j or i-1,j-1 (0-index or 1-index)
        i<=props.Acstart && (props.Acstart = i)
        j<=props.Bcstart && (props.Bcstart = j)
        props.AtoB[i] = j
        break
      end
    end
  end

  for i = 1:na
    for k = 1:nc
      if props.ai[i]==props.ci[k]
        #TODO: check this if this should be i,j or i-1,j-1 (0-index or 1-index)
        i<=props.Austart && (props.Austart = i)
        props.AtoC[i] = k
        break
      end
    end
  end

  for j = 1:nb
    for k = 1:nc
      if props.bi[j]==props.ci[k]
        #TODO: check this if this should be i,j or i-1,j-1 (0-index or 1-index)
        j<=props.Bustart && (props.Bustart = j)
        props.BtoC[j] = k
        break
      end
    end
  end

end

function is_trivial_permutation(P::Vector{Int})::Bool
  for n = 1:length(P)
    P[n]!=n && return false
  end
  return true
end

function checkACsameord(props::CProps)::Bool
  props.Austart>=length(props.ai) && return true
  aCind = props.AtoC[props.Austart]
  for i = 1:length(props.ai)
    if !contractedA(props,i)
      props.AtoC[i]!=aCind && return false
      aCind += 1
    end
  end
  return true
end

function checkBCsameord(props::CProps)::Bool
  props.Bustart>=length(props.bi) && return true
  bCind = props.BtoC[props.Bustart]
  for i = 1:length(props.bi)
    if !contractedB(props,i)
      props.BtoC[i]!=bCind && return false
      bCind += 1
    end
  end
  return true
end

contractedA(props::CProps,i::Int) = (props.AtoC[i]<1)
contractedB(props::CProps,i::Int) = (props.BtoC[i]<1)
Atrans(props::CProps) = contractedA(props,1)
Btrans(props::CProps) = !contractedB(props,1)
Ctrans(props::CProps) = props.ctrans

function find_index(v::Vector{Int},t)::Int
  for i = 1:length(v)
    v[i]==t && return i
  end
  return -1
end

function permute_extents(R::Vector{Int},P::Vector{Int})::Vector{Int}
  Rb = fill(0,length(R))
  n = 1
  for pn in P
    Rb[pn] = R[n]
    n += 1
  end
  return Rb
end

function compute!(props::CProps,
                  A::Array,
                  B::Array,
                  C::Array)
  compute_perms!(props)

  #Use props.PC.size() as a check to see if we've already run this
  length(props.PC)!=0 && return

  ra = length(props.ai)
  rb = length(props.bi)
  rc = length(props.ci)

  props.PC = fill(0,rc)

  props.dleft = 1
  props.dmid = 1
  props.dright = 1
  c = 1
  for i = 1:ra
    if !contractedA(props,i)
      props.dleft *= size(A,i)
      props.PC[props.AtoC[i]] = c
      c += 1
    else
      props.dmid *= size(A,i)
    end
  end
  for j = 1:rb
    if !contractedB(props,j)
      props.dright *= size(B,j)
      props.PC[props.BtoC[j]] = c
      c += 1
    end
  end

  if !is_trivial_permutation(props.PC)
    props.permuteC = true
    if checkBCsameord(props) && checkACsameord(props)
      #Can avoid permuting C by
      #computing Bt*At = Ct
      props.ctrans = true
      props.permuteC = false
    end
  end

  #Check if A can be treated as a matrix without permuting
  props.permuteA = false
  if !(contractedA(props,1) || contractedA(props,ra))
    #If contracted indices are not all at front or back,
    #will have to permute A
    props.permuteA = true
  else
    #Contracted ind start at front or back, check if contiguous
    #TODO: check that the limits are correct (1-indexed vs. 0-indexed)
    for i = 1:props.ncont
      if !contractedA(props,props.Acstart+i-1)
        #Contracted indices not contiguous, must permute
        props.permuteA = true
        break
      end
    end
  end

  #Check if B is matrix-like
  props.permuteB = false
  if !(contractedB(props,1) || contractedB(props,rb))
    #If contracted indices are not all at front or back,
    #will have to permute B
    props.permuteB = true
  else
    #TODO: check that the limits are correct (1-indexed vs. 0-indexed)
    for i = 1:props.ncont
      if !contractedB(props,props.Bcstart+i-1)
        #Contracted inds not contiguous, permute
        props.permuteB = true
        break
      end
    end
  end

  if !props.permuteA && !props.permuteB
    #Check if contracted inds. in same order
    #TODO: check these limits are correct
    for i = 1:props.ncont
      if props.AtoB[props.Acstart+i-1]!=(props.Bcstart+i-1)
        #If not in same order,
        #must permute one of A or B
        #so permute the smaller one
        props.dleft<props.dright ? (props.permuteA = true) : (props.permuteB = true)
        break
      end
    end
  end

  if props.permuteC && !(props.permuteA && props.permuteB)
    PCost(d::Real) = d*d
    #Could avoid permuting C if
    #permute both A and B, worth it?
    pCcost = PCost(props.dleft*props.dright)
    extra_pABcost = 0
    !props.permuteA && (extra_pABcost += PCost(props.dleft*props.dmid))
    !props.permuteB && (extra_pABcost += PCost(props.dmid*props.dright))
    if extra_pABcost<pCcost
      props.permuteA = true
      props.permuteB = true
      props.permuteC = false
    end
  end

  if props.permuteA
    props.PA = fill(0,ra)
    #Permute contracted indices to the front,
    #in the same order as on B
    newi = 0
    #TODO: check this is correct for 1-indexing
    bind = props.Bcstart
    for i = 1:props.ncont
      while !contractedB(props,bind) bind += 1 end
      j = find_index(props.ai,props.bi[bind])
      newi += 1
      props.PA[newi] = j
      bind += 1
    end
    #Reset p.AtoC:
    fill!(props.AtoC,0)
    #Permute uncontracted indices to
    #appear in same order as on C
    #TODO: check this is correct for 1-indexing
    for k = 1:rc
      j = find_index(props.ai,props.ci[k])
      if j!=-1
        props.AtoC[newi] = k
        props.PA[newi+1] = j
        newi += 1
      end
      newi==ra && break
    end
  end

  ##Also update props.Austart,props.Acstart
  props.Acstart = ra+1
  props.Austart = ra+1
  #TODO: check this is correct for 1-indexing
  for i = 1:ra
    if contractedA(props,i)
      props.Acstart = min(i,props.Acstart)
    else
      props.Austart = min(i,props.Austart)
    end
    props.newArange = permute_extents([size(A)...],props.PA)
  end

  if(props.permuteB)
    props.PB = fill(0,rb)
    #TODO: check this is correct for 1-indexing
    newi = 0 #1
    if(props.permuteA)
      #A's contracted indices already set to
      #be in same order as B above, so just
      #permute contracted indices to the front
      #keeping relative order
      #TODO: how to translate this for loop?
      #for(int i = props.Bcstart; newi < props.ncont; ++newi)
      i = props.Bcstart
      #TODO: check this is correct for 1-indexing
      while newi < props.ncont
        while !contractedB(props,i) i += 1 end
        props.PB[newi+1] = i
        i += 1
        newi += 1
      end
    else
      #Permute contracted indices to the
      #front and in same order as on A
      aind = props.Acstart
      for i = 0:(props.ncont-1)
        while !contractedA(props,aind) aind += 1 end
        j = find_index(props.bi,props.ai[aind])
        newi += 1
        props.PB[newi] = j
        aind += 1
      end
    end
    #Reset p.BtoC:
    fill!(props.BtoC,0)
    #Permute uncontracted indices to
    #appear in same order as on C
    for k = 1:rc
      j = find_index(props.bi,props.ci[k])
      if j!=-1
        props.BtoC[newi] = k
        props.PB[newi+1] = j
        newi += 1
      end
      newi==rb && break
    end
    props.Bcstart = rb
    props.Bustart = rb
    for i = 1:rb
      if(contractedB(props,i))
          props.Bcstart = min(i,props.Bcstart)
      else
          props.Bustart = min(i,props.Bustart)
      end
    end
    props.newBrange = permute_extents([size(B)...],props.PB)
  end

  if props.permuteA || props.permuteB
    #Recompute props.PC
    c = 1
    #TODO: check this is correct for 1-indexing
    for i = 1:ra
      if !contractedA(props,i)
        props.PC[props.AtoC[i]] = c
        c += 1
      end
    end
    #TODO: check this is correct for 1-indexing
    for j = 1:rb
      if !contractedB(props,j)
        props.PC[props.BtoC[j]] = c
        c += 1
      end
    end
    props.ctrans = false
    if(is_trivial_permutation(props.PC))
      props.permuteC = false
    else
      props.permuteC = true
      #Here we already know since pc_triv = false that
      #at best indices from B precede those from A (on result C)
      #so if both sets remain in same order on C
      #just need to transpose C, not permute it
      if  checkBCsameord(props) && checkACsameord(props)
        props.ctrans = true
        props.permuteC = false
      end
    end
  end

  if props.permuteC
    Rb = Int[]
    if !props.permuteA
      #TODO: check this is correct for 1-indexing
      for i = 1:ra
        if !contractedA(props,i)
          push!(Rb,size(A,i))
        end
      end
    else
      #TODO: check this is correct for 1-indexing
      for i = 1:ra
        if !contractedA(props,i)
          push!(Rb,size(props.newArange,i))
        end
      end
    end
    if !props.permuteB
      #TODO: check this is correct for 1-indexing
      for j = 1:rb
        if !contractedB(props,j)
          push!(Rb,size(B,j))
        end
      end
    else
      #TODO: check this is correct for 1-indexing
      for j = 1:rb
        if !contractedB(props,j)
          push!(Rb,size(props.newBrange,j))
        end
      end
    end
    props.newCrange = Rb
  end

end

function contract!(C::Array{T},
                   p::CProps,
                   A::Array{T},
                   B::Array{T},
                   α::T=one(T),
                   β::T=zero(T)) where {T}

  tA = 'N'
  if p.permuteA
    aref = reshape(permutedims(A,p.PA),p.dmid,p.dleft)
    tA = 'T'
  else
    #A doesn't have to be permuted
    if Atrans(p)
      aref = reshape(A,p.dmid,p.dleft)
      tA = 'T'
    else
      aref = reshape(A,p.dleft,p.dmid)
    end
  end

  tB = 'N'
  if p.permuteB
    bref = reshape(permutedims(B,p.PB),p.dmid,p.dright)
  else
    if Btrans(p)
      bref = reshape(B,p.dright,p.dmid)
      tB = 'T'
    else
      bref = reshape(B,p.dmid,p.dright)
    end
  end

  if p.permuteC
    cref = reshape(copy(C),p.dleft,p.dright)
  else
    if Ctrans(p)
      cref = reshape(C,p.dleft,p.dright)
      if tA=='N' && tB=='N'
        (aref,bref) = (bref,aref)
        tA = tB = 'T'
      elseif tA=='T' && tB=='T'
        (aref,bref) = (bref,aref)
        tA = tB = 'N'
      end
    else
      cref = reshape(C,p.dleft,p.dright)
    end
  end

  #BLAS.gemm!(tA,tB,promote_type(T,Tα)(α),aref,bref,promote_type(T,Tβ)(β),cref)
  BLAS.gemm!(tA,tB,α,aref,bref,β,cref)

  if p.permuteC
    permutedims!(C,reshape(cref,p.newCrange...),p.PC)
  end
  return
end

#TODO: this should be optimized
function contract_scalar!(Cdata::Array,Clabels::Vector{Int},
                          Bdata::Array,Blabels::Vector{Int},α,β)
  p = calculate_permutation(Blabels,Clabels)
  if β==0
    if is_trivial_permutation(p)
      Cdata .= α.*Bdata
    else
      #TODO: make an optimized permutedims!() that also scales the data
      permutedims!(Cdata,α*Bdata)
    end
  else
    if is_trivial_permutation(p)
      Cdata .= α.*Bdata .+ β.*Cdata
    else
      #TODO: make an optimized permutedims!() that also adds and scales the data
      permBdata = permutedims(Bdata,p)
      Cdata .= α.*permBdata .+ β.*Cdata
    end
  end
  return
end

function contract!(Cdata::Array{T},Clabels::Vector{Int},
                   Adata::Array{T},Alabels::Vector{Int},
                   Bdata::Array{T},Blabels::Vector{Int},
                   α::T=one(T),β::T=zero(T)) where {T}
  if(length(Alabels)==0)
    contract_scalar!(Cdata,Clabels,Bdata,Blabels,α*Adata[1],β)
  elseif(length(Blabels)==0)
    contract_scalar!(Cdata,Clabels,Adata,Alabels,α*Bdata[1],β)
  else
    props = CProps(Alabels,Blabels,Clabels)
    compute!(props,Adata,Bdata,Cdata)
    contract!(Cdata,props,Adata,Bdata,α,β)
  end
  return
end

function contract(Cinds::IndexSet,
                  Clabels::Vector{Int},
                  Astore::Dense{SA},
                  Ainds::IndexSet,
                  Alabels::Vector{Int},
                  Bstore::Dense{SB},
                  Binds::IndexSet,
                  Blabels::Vector{Int}) where {SA<:Number,SB<:Number}
  SC = promote_type(SA,SB)

  # Convert the arrays to a common type
  # since we will call BLAS
  Astore = convert(Dense{SC},Astore)
  Bstore = convert(Dense{SC},Bstore)

  Adims = dims(Ainds)
  Bdims = dims(Binds)
  Cdims = dims(Cinds)

  # Create storage for output tensor
  Cstore = Dense{SC}(prod(Cdims))

  Adata = reshape(data(Astore),Adims)
  Bdata = reshape(data(Bstore),Bdims)
  Cdata = reshape(data(Cstore),Cdims)

  contract!(Cdata,Clabels,Adata,Alabels,Bdata,Blabels)
  return Cstore
end


text:
struct TopKBatchQs <: ExtendingBatchQs
    sequentialQs::SequentialPQs
    k::Int

    function TopKBatchQs(sequentialQs::SequentialPQs; k::Int)::TopKBatchQs
        (k < 1) && throw(ArgumentError("Invalid batch size k=$(k)."))
        return new(sequentialQs, k)
    end
end

"""
Use sequential strategy, select batch_size observations with highest usefulness.
"""
function select_batch(qs::TopKBatchQs, x::Array{T, 2}, labels::Dict{Symbol, Vector{Int}}, candidate_indices::Vector{Int})::Vector{Int} where T <: Real
    num_observations = length(candidate_indices)
    if num_observations <= qs.k
        return candidate_indices
    end

    sequential_strategy = qs.sequentialQs
    candidate_scores = qs_score(sequential_strategy, x, labels)[candidate_indices]
    descending_indices = sortperm(candidate_scores; rev=true)
    best_k = candidate_indices[descending_indices[1:qs.k]]
    return best_k
end


text:
module Axelrodnet

    using LightGraphs
    using GraphIO
    using DataFrames
    using Random
    using Feather
    using ParserCombinator
    using REPL
    using StatsBase

    export run_experiment
    export run
    export export_experiment

    include("agent.jl")
    include("assert.jl")
    include("io.jl")
    include("tick.jl")
    include("network.jl")
    include("utilities.jl")
    include("simulation.jl")

end

max_stars_repo_name: JohannesNakayama/axelrodnet

max_issues_repo_name: JohannesNakayama/axelrodnet

max_forks_repo_name: JohannesNakayama/axelrodnet


text:
# check out the following for active learning:
# https://ocw.mit.edu/courses/electrical-engineering-and-computer-science/6-867-machine-learning-fall-2006/lecture-notes/lec6.pdf
# says we need to pick the largest eigenvector of A=(X'X)^-1 where X has row vectors for all previous examples.
# note the closed form solution for linreg is w=AX'y
# all this assumes X'X is invertible.  In our case X is fat and it is probably not.

# https://en.wikipedia.org/wiki/Active_learning_(machine_learning)
# says svm active learning picks points closest to the hyperplane or furthest or alternate.
# closest suggests picking vecs orthogonal to wpred - in 10K dimensions any random vector will do.
# it may not be possible to beat random vector.

using Knet
isdefined(:MNIST) || include(Knet.dir("examples/mnist.jl"))
using MNIST: minibatch,xtst,ytst,xtrn,ytrn
(x0,y0) = minibatch(xtst, ytst, 10000; atype=KnetArray{Float32})[1]
loss(w) = -sum(y0 .* logp(w*x0,1))/size(y0,2)
gloss = grad(loss)
vcos(a,b)=dot(a,b)/(vecnorm(a)*vecnorm(b))

# This version gets stuck around loss=0.4 in 100 epochs with batch=100.

function gfo(w=nothing; batch=100, epochs=100, atype=KnetArray{Float64}, vstd=1e-5, lr_g=1.0, lr_h=1.0, hpred=0.01, pr=0, lr=0.05, decay=0.8)
    Knet.rng(true)
    dtrn = minibatch(xtrn,ytrn,batch;atype=atype)
    dtst = minibatch(xtst,ytst,length(ytst);atype=atype)
    f(w,x,y) = -sum(y .* logp(w*x,1))/size(y,2); gradf = grad(f)
    w==nothing && (w = convert(atype, zeros(10,784)))
    vrnd = zeros(w); gpred = zeros(w); hpred0=hpred; lr0=lr; loss0=loss1=Inf
    report(epoch)=println(map((epoch,
                               :lr, lr,
                               :loss, mean(f(w,x,y) for (x,y) in dtst),
                               :vcos, mean(vcos(gpred,gradf(w,x,y)) for (x,y) in dtst),
                               :ggold, mean(vecnorm(gradf(w,x,y)) for (x,y) in dtst),
                               :gpred, vecnorm(gpred), :hpred, hpred
                               )) do x; isa(x,AbstractFloat) ? parse(@sprintf("%g",x)) : x; end)
    report(0)
    approx(a,b)=isapprox(a,b;rtol=1e-3,atol=1e-12)
    for epoch=1:epochs
        for (x,y) in dtrn
            f0 = f(w,x,y)
            g0 = gradf(w,x,y)

            # adjust direction using small random step
            v = randn!(vrnd,0,vstd)
            vv = sumabs2(v)
            w1 = w+v
            f1 = f(w1,x,y)
            dfgold = f1-f0
            dfpred = dot(gpred,v)
            delta = (dfpred - dfgold)
            vcos1 = vcos(g0,gpred)
            gpred -= (lr_g * delta / vv) * v
            vcos2 = vcos(g0,gpred)
            approx(dfgold, dot(gpred,v)) || error("dfgold1=$dfgold dfpred1=$(dot(gpred,v))")

            # adjust size using small g step
            v = (-vecnorm(v)/vecnorm(gpred))*gpred
            approx(vv, sumabs2(v)) || error("vv=$vv sumabs2(v)=$(sumabs2(v))")
            w2 = w+v
            f2 = f(w2,x,y)
            # f2 < f0 || warn("f0=$f0 f2=$f2") # f may not decrease when w changes
            dfgold = f2-f0
            dfpred = dot(gpred,v)
            delta = (dfpred - dfgold)
            size1 = vecnorm(gpred)
            size0 = abs(dot(g0,gpred)/size1)
            gpred -= (lr_g * delta / vv) * v
            approx(dfgold, dot(gpred,v)) || error("dfgold2=$dfgold dfpred2=$(dot(gpred,v))")
            size2 = vecnorm(gpred)
            isapprox(size2, size0, rtol=.1, atol=1e-5) || error("size0=$size0 size2=$size2")

            # adjust curvature using newton step
            v = (-1/hpred)*gpred
            vv = sumabs2(v)
            w3 = w+v
            f3 = f(w3,x,y)
            dfgold = f3-f0
            dfpred = dot(gpred,v) + 0.5 * hpred * vv
            dfpred < 0 || error("dfpred3=$dfpred dot=$(dot(gpred,v)) hpred=$hpred vv=$vv")
            delta = (dfpred - dfgold) # =0.5*(hpred-hgold)*vv
            hpred1 = hpred
            hpred2 = hpred - (2 * delta / vv)
            approx(dfgold, dot(gpred,v)+0.5*hpred2*vv) || error("dfgold3=$dfgold dfpred3=$(dot(gpred,v)+0.5*hpred2*vv)")
            hpred3 = hpred - (lr_h * 2 * delta / vv)
            if hpred2 >= 0 # 0.001
                hpred = hpred2 # clamp(hpred3, 0.01, 10.0)
            elseif hpred2 < 0
                warn("hpred=$hpred hpred2=$hpred2 delta=$delta vv=$vv dfpred=$dfpred dfgold=$dfgold")
            end

            # finally take that step
            v = (-lr/hpred)*gpred
            w4 = w+v
            f4 = f(w4,x,y)
            if f4 <= f0
                w = w4
                gpred += hpred * v
            else
                warn("f0=$f0 f4=$f4 hpred=$hpred gpred=$(vecnorm(gpred))")
            end

            if epoch in pr
                @printf("f=%g->%g vcos=%g->%g size=%g->%g curv=%g->%g delta=%g\n",
                        f0,f4,vcos1,vcos2,size1,size2,hpred1,hpred2,delta)
            end

            # gpred -= (lr_g * delta / (2*vv)) * v
            # hpred -= (lr_h * delta / vv)

            # if 0.01 < hpred2 < 100.0    # do not let hpred get too small or big
            #     hpred = hpred2
            # end
            # if f2 <= f0
            #     w = w2
            #     gpred += (hpred - lr_g * delta / (2*vv)) * v
            #     # gpred -= (lr_g * delta / (2*vv)) * v
            #     # gpred += hpred * v    # we assume hpred stays constant, but gpred should be updated after the step
            # else
            #     fill!(gpred,0)
            #     delta < 0 || error("delta=$delta")
            #     # (0.01<hpred<100) && (hpred -= (lr_h * delta / vv))
            #     # hpred = hpred0
            #     # gpred *= -1
            #     # lr *= 0.99f0
            #     # lr < (1e-6) && (lr=lr0; fill!(gpred,0))
            # end
        end
        report(epoch)
        loss1 = mean(f(w,x,y) for (x,y) in dtst)
        loss1 < loss0 || (lr *= decay)
        loss0 = loss1
    end
    return w
end

# SGD that takes a newton step in gradient direction every iteration:
# gets stuck around 0.39 with lr=1.0.
# gets stuck around 0.31 with lr=0.5.
# gets stuck around 0.2775 with lr=0.1.
# reaches 0.2721 with lr=0.1, decay=0.5.
# vcos usually negative.

function supersgd(w=nothing; batch=100, epochs=100, atype=KnetArray{Float64}, hpred=0.4, pr=0, lr=1.0, decay=1.0)
    Knet.rng(true)
    dtrn = minibatch(xtrn,ytrn,batch;atype=atype)
    dtst = minibatch(xtst,ytst,length(ytst);atype=atype)
    f(w,x,y) = -sum(y .* logp(w*x,1))/size(y,2); gradf = grad(f)
    w==nothing && (w = convert(atype, zeros(10,784)))
    gpred = zeros(w); loss0=loss1=Inf
    report(epoch)=println(map((epoch,
                               :lr, lr,
                               :loss, mean(f(w,x,y) for (x,y) in dtst),
                               :vcos, mean(vcos(gpred,gradf(w,x,y)) for (x,y) in dtst),
                               :ggold, mean(vecnorm(gradf(w,x,y)) for (x,y) in dtst),
                               :gpred, vecnorm(gpred), :hpred, hpred
                               )) do x; isa(x,AbstractFloat) ? parse(@sprintf("%g",x)) : x; end)
    report(0)
    approx(a,b)=isapprox(a,b;rtol=1e-3,atol=1e-12)
    for epoch=1:epochs
        for (x,y) in dtrn
            f0 = f(w,x,y)
            gpred = gradf(w,x,y)   # perfect direction and size, only need curvature

            # adjust curvature using newton step
            v = (-1/hpred)*gpred
            vv = sumabs2(v)
            w3 = w+v
            f3 = f(w3,x,y)
            dfgold = f3-f0
            dfpred = dot(gpred,v) + 0.5 * hpred * vv
            dfpred < 0 || error("dfpred3=$dfpred dot=$(dot(gpred,v)) hpred=$hpred vv=$vv")
            delta = (dfpred - dfgold) # =0.5*(hpred-hgold)*vv
            hpred1 = hpred
            hpred2 = hpred - (2 * delta / vv)
            approx(dfgold, dot(gpred,v)+0.5*hpred2*vv) || error("dfgold3=$dfgold dfpred3=$(dot(gpred,v)+0.5*hpred2*vv)")
            if hpred2 >= 0 # 0.001
                hpred = hpred2 # clamp(hpred3, 0.01, 10.0)
            elseif hpred2 < 0
                warn("hpred=$hpred hpred2=$hpred2 delta=$delta vv=$vv dfpred=$dfpred dfgold=$dfgold")
            end

            # finally take that step
            v = (-lr/hpred)*gpred
            w4 = w+v
            f4 = f(w4,x,y)
            if f4 <= f0
                w = w4
                # gpred += hpred * v
            else
                warn("f0=$f0 f4=$f4 hpred=$hpred gpred=$(vecnorm(gpred))")
            end

            if epoch in pr
                @printf("f=%g->%g vcos=%g->%g size=%g->%g curv=%g->%g delta=%g\n",
                        f0,f4,vcos1,vcos2,size1,size2,hpred1,hpred2,delta)
            end

        end
        report(epoch)
        loss1 = mean(f(w,x,y) for (x,y) in dtst)
        loss1 < loss0 || (lr *= decay)
        loss0 = loss1
    end
    return w
end

# gets to loss=0.2730 in 100 epochs with batch=100 lr=0.1
# gets to loss=0.2715 with lr=0.1, decay=0.5.
# vcos between test gradient and training minibatches is around 0.1.
# vcos changes when lr changes, looks like a bug!

function sgd(; batch=100, lr=0.1, epochs=100, atype=KnetArray{Float32}, decay=0.5)
    dtrn = minibatch(xtrn,ytrn,batch;atype=atype)
    dtst = minibatch(xtst,ytst,length(ytst);atype=atype)
    f(w,x,y) = -sum(y .* logp(w*x,1))/size(y,2)
    g = grad(f)
    w = convert(atype, zeros(10,784))
    println((0, mean([f(w,x,y) for (x,y) in dtst])))
    loss1 = loss0 = Inf
    for epoch=1:epochs
        vcos1 = Any[]
        for (x,y) in dtrn
            dw = g(w,x,y)
            push!(vcos1, mean(vcos(dw,g(w,x,y)) for (x,y) in dtst))
            axpy!(-lr,dw,w)
        end
        loss1 = mean([f(w,x,y) for (x,y) in dtst])
        println((epoch, :lr, lr, :loss, loss1, :vcos, mean(vcos1)))
        loss1 > loss0 && (lr *= decay)
        loss0 = loss1
    end
    return w
end

# measure gradient variance
function gradvar(n=100)
    f(w,x,y) = -sum(y .* logp(w*x,1))/size(y,2)
    g = grad(f)
    data = minibatch(xtrn,ytrn,n; atype=Array{Float64})
    grads = Any[]
    w = oftype(data[1][1], zeros(10,784))
    for (x,y) in data
        push!(grads, g(w,x,y))
    end
    g0 = mean(grads)
    cosine = Any[]
    sqdiff = Any[]
    for g in grads
        push!(cosine, vcos(g,g0))
        push!(sqdiff, sumabs2(g-g0))
    end
    (:cos, mean(cosine), :std, sqrt(mean(sqdiff)), :nrm, mean(map(vecnorm,grads)))
end

function flinreg(w,x,y;l1=0,l2=0,l3=0)
    lss = zero(eltype(w))
    if l1 != 0; (lss += l1*sumabs(w)); end
    if l2 != 0; (lss += l2*sumabs2(w)); end
    if l3 != 0 # prevent l1 reg from changing the sign
        lss += sum(min(0.5*abs2(w),l3*abs(w)))
    end
    lss += abs2(y - sum(w.*x))
    return lss
end

glinreg = grad(flinreg)

# with no l1/l2: 10.18 10000 (:cos,0.74111754f0)
# best l1=1e-14 11.24 10000 (:cos,0.7515351f0)
# no meaningful improvement with l2.

function linreg(epochs=10000;l1=1e-14,l2=0,l3=0,lr=0.5)
    Knet.rng(true)
    w0 = oftype(x0, zeros(10,784))
    f0 = loss(w0)
    g0 = gloss(w0)
    gpred = oftype(x0, zeros(10,784))
    dw = similar(w0)
    progress(epoch)=(:cos,dot(gpred,g0)/(vecnorm(gpred)*vecnorm(g0)))
    report()
    for epoch=1:epochs
        randn!(dw, 0, 1e-6)
        w1 = w0 + dw
        f1 = loss(w1)
        dgold = f1 - f0
        ggrad = glinreg(gpred, dw, dgold; l1=l1, l2=l2, l3=l3)
        # axpy!(-lr, ggrad, gpred)
        # See the changelog 2016-11-03 for the sumabs2(dw) factor.
        axpy!(-lr/sumabs2(dw), ggrad, gpred)
        report(progress, epoch)
    end
    report(progress, epochs, final=true)
    return gpred
end

function flogreg(w,x,y,l1,l2)
    lss = zero(eltype(w))
    if l1 != 0; (lss += l1*sumabs(w)); end
    if l2 != 0; (lss += l2*sumabs2(w)); end
    prob = sigm(sum(w.*x)) #sigm(dot(w,x))
    if sign(y) < 0; (prob = 1-prob); end
    lss -= log(prob)
    return lss
end

glogreg = grad(flogreg)

# l1=3e-10 works best in epochs=1000,lr=1 -- 1.09 1000 (:cos,0.26973265f0)
# compared to no l1/l2: 1.15 1000 (:cos,0.26972535f0)
# compared to perceptron: 1.07 1000 (:cos,0.21155934f0,:acc,0.5769524461878847)
# l1=3e-9 works best in epochs=10000,lr=1 -- 10.69 10000 (:cos,0.6744472f0)
# l2=1e-6 works best in epochs=10000,lr=1 -- 10.68 10000 (:cos,0.66643745f0)
# compared to no l1/l2: 9.97 10000 (:cos,0.66643345f0)
# compared to perceptron: 20.90 10000 (:cos,0.5908562f0,:acc,0.7445132059813894)
function logreg(epochs=10000;l1=3e-9,l2=0,lr=1.0)
    Knet.rng(true)
    w0 = oftype(x0, zeros(10,784))
    f0 = loss(w0)
    g0 = gloss(w0)
    gpred = oftype(x0, zeros(10,784))
    dw = similar(w0)
    progress(epoch)=(:cos,dot(gpred,g0)/(vecnorm(gpred)*vecnorm(g0)))
    report()
    for epoch=1:epochs
        randn!(dw, 0, 1e-6)
        w1 = w0 + dw
        f1 = loss(w1)
        dgold = f1 - f0
        ggrad = glogreg(gpred, dw, dgold, l1, l2)
        axpy!(-lr, ggrad, gpred)
        report(progress, epoch)
    end
    report(progress, epochs, final=true)
    return gpred
end

function perceptron(epochs=100000)
    w0 = oftype(x0, zeros(10,784))
    f0 = loss(w0)
    g0 = gloss(w0)
    gpred = oftype(x0, randn(10,784)*1e-6)
    dw = similar(w0)
    acc = 0.5
    progress(epoch)=(:cos,dot(gpred,g0)/(vecnorm(gpred)*vecnorm(g0)),:acc,acc)
    report()
    for epoch=1:epochs
        randn!(dw, 0, 1e-6)
        w1 = w0 + dw
        f1 = loss(w1)
        dgold = f1 - f0
        dpred = dot(gpred,dw)
        acc = 0.999 * acc
        if dpred * dgold < 0
            gpred += sign(dgold) * dw
        else
            acc += 0.001
        end
        report(progress, epoch)
    end
    report(progress, epochs, final=true)
    return gpred
end

function spsa(epochs=100000)
    w0 = oftype(x0, zeros(10,784))
    f0 = loss(w0)
    g0 = gloss(w0)
    gpred = oftype(x0, randn(10,784)*1e-6)
    dw = similar(w0)
    acc = 0.5
    progress(epoch)=(:cos,dot(gpred,g0)/(vecnorm(gpred)*vecnorm(g0)),:acc,acc)
    report()
    for epoch=1:epochs
        dw = oftype(x0, rand(-1e6:2e6:1e6, size(w0)))
        w1 = w0 + dw
        f1 = loss(w1)
        dgold = f1 - f0
        dpred = dot(gpred,dw)
        acc = 0.999 * acc + 0.001 * (dpred * dgold >= 0)
        gpred += dgold ./ dw
        report(progress, epoch)
    end
    report(progress, epochs, final=true)
    return gpred
end

let time0=tnext=nnext=1
    global report
    report()=(time0=time();tnext=time0;nnext=1)
    function report(f,n; dt=10, dn=2, final=false)
        tnext == 1 && report()
        if final || n >= nnext || time() >= tnext
            t = time()
            @printf("%.2f %d %s\n", t-time0, n, f(n))
            t >= tnext && (tnext = t + dt)
            n >= nnext && (nnext *= dn)
        end
    end
end

macro cuda1(lib,fun,x...)
    if Libdl.find_library(["lib$lib"], []) == ""
        msg = "Cannot find lib$lib, please install it and rerun Pkg.build(\"Knet\")."
        :(error($msg))
    else
        f2 = ("$fun","lib$lib")
        fx = Expr(:ccall, f2, :UInt32, x...)
        err = "$lib.$fun error "
        quote
            local _r = $fx
            if _r != 0
                warn($err, _r)
                Base.show_backtrace(STDOUT, backtrace())
            end
        end
    end
end

using Knet: Cptr, rng, cublashandle

import Base: randn!
randn!(a::KnetArray{Float32},mean,stddev)=(@cuda1(curand,curandGenerateNormal,(Cptr,Ptr{Float32},Csize_t,Float32,Float32),rng(),a,length(a),mean,stddev); a)
randn!(a::KnetArray{Float64},mean,stddev)=(@cuda1(curand,curandGenerateNormalDouble,(Cptr,Ptr{Float64},Csize_t,Float64,Float64),rng(),a,length(a),mean,stddev); a)
randn!(a::Array,mean,stddev)=(randn!(a);stddev!=1 && scale!(a,stddev);mean!=0 && (a[:]+=mean);a)

import Base.LinAlg: dot
dot(x::KnetArray{Float32},y::KnetArray{Float32})=(c=Float32[0];@cuda1(cublas,cublasSdot_v2, (Cptr, Cint, Ptr{Float32}, Cint, Ptr{Float32}, Cint, Ptr{Float32}), cublashandle, length(x), x, 1, y, 1, c);c[1])
dot(x::KnetArray{Float64},y::KnetArray{Float64})=(c=Float64[0];@cuda1(cublas,cublasDdot_v2, (Cptr, Cint, Ptr{Float64}, Cint, Ptr{Float64}, Cint, Ptr{Float64}), cublashandle, length(x), x, 1, y, 1, c);c[1])

import Knet: sigmback
sigmback(a::Number,b::Number)=sigmback(promote(a,b)...)

max_stars_repo_name: JuliaTagBot/gfo

max_issues_repo_name: JuliaTagBot/gfo

max_forks_repo_name: JuliaTagBot/gfo


text:
module TestFileSystemCache

using Base.Test
using SimpleTasksTests.Utils.TestServices

import SimpleTasks.Services.Cache
import SimpleTasks.Services.FileSystemCache

function test_blank_directory()
    @test_throws ArgumentError FileSystemCache.FileSystemCacheService(" ")
end

function test_good_cache()
    @test FileSystemCache.FileSystemCacheService(TEST_BASE_DIRECTORY) !=
        nothing
end

function test_sanitize_key()
    @test FileSystemCache.sanitize("../../.ssh/id_rsa") == "//.ssh/id_rsa"
end

function test_to_filename()
    cache = make_valid_file_system_cache()
    @test FileSystemCache.to_filename(cache, "file") ==
        "$TEST_BASE_DIRECTORY/file"
end

function test_to_filename_sanitized()
    cache = make_valid_file_system_cache()
    @test FileSystemCache.to_filename(cache, "../../.ssh/id_rsa") ==
        "$TEST_BASE_DIRECTORY///.ssh/id_rsa"
end

function test_has_key()
    filename = "testfile"
    full_path_file_name = "$TEST_BASE_DIRECTORY/$filename"
    file = open(full_path_file_name, "w")
    close(file)
    cache = make_valid_file_system_cache()
    @test haskey(cache, filename)
    rm(full_path_file_name)
end

function test_not_exist()
    filename = "testfileDOES NOT EXIST"
    full_path_file_name = "$TEST_BASE_DIRECTORY/$filename"
    cache = make_valid_file_system_cache()
    @test !haskey(cache, filename)
end

function test_put()
    filename = "key"
    test_string = "TEST"
    buffer = IOBuffer(test_string)
    cache = make_valid_file_system_cache()
    cache[filename] = buffer
    file = open("$TEST_BASE_DIRECTORY/$filename", "r")
    @test isreadable(file)
    @test readall(file) == test_string
    close(file)
end

function test_put_path()
    filename = "a/b/c/key"
    test_string = "TEST"
    buffer = IOBuffer(test_string)
    cache = make_valid_file_system_cache()
    cache[filename] = buffer
    file = open("$TEST_BASE_DIRECTORY/$filename", "r")
    @test isreadable(file)
    @test readall(file) == test_string
    @test ispath("$TEST_BASE_DIRECTORY/a/b/c/")
    rm("$TEST_BASE_DIRECTORY/$filename")
    rm("$TEST_BASE_DIRECTORY/a/b/c/")
    rm("$TEST_BASE_DIRECTORY/a/b/")
    rm("$TEST_BASE_DIRECTORY/a/")
    close(file)
end

function test_get()
    filename = "a/b/c/key"
    mkpath("$TEST_BASE_DIRECTORY/a/b/c/")
    file = open("$TEST_BASE_DIRECTORY/$filename", "w")
    test_string = "TEST"
    write(file, "TEST")
    close(file)

    cache = make_valid_file_system_cache()
    buffer = Cache.get(cache, filename)
    @test readall(buffer) == test_string
    close(buffer)
end

function test_get_not_has_key()
    filename = "testfileDOES NOT EXIST"
    cache = make_valid_file_system_cache()
    @test_throws KeyError Cache.get(cache, filename)
end

function test_delete()
    filename = "key"
    test_string = "TEST"
    buffer = IOBuffer(test_string)
    cache = make_valid_file_system_cache()
    cache[filename] = buffer
    delete!(cache, filename)
    @test !isfile("$TEST_BASE_DIRECTORY/$filename")
end

function test_delete_not_has_key()
    filename = "keyDOES NOT EXIST"
    cache = make_valid_file_system_cache()
    # should not throw
    Cache.delete!(cache, filename)
end

function test_delete_path()
    filename = "a/b/c/key"
    test_string = "TEST"
    buffer = IOBuffer(test_string)
    cache = make_valid_file_system_cache()
    cache[filename] = buffer
    delete!(cache, filename)
    @test !isfile("$TEST_BASE_DIRECTORY/$filename")
end

function test_clear()
    cache = make_valid_file_system_cache()

    filename1 = "a/b/c/key1"
    test_string1 = "TEST1"
    buffer1 = IOBuffer(test_string1)

    filename2 = "a/b/c/key2"
    test_string2 = "TEST2"
    buffer2 = IOBuffer(test_string2)

    cache[filename1] = buffer1
    cache[filename2] = buffer2
    Cache.clear!(cache)

    @test isempty(readdir(cache.base_directory))
end

function __init__()
    test_blank_directory()
    test_good_cache()

    test_sanitize_key()
    test_to_filename()
    test_to_filename_sanitized()

    test_has_key()
    test_not_exist()

    test_put()
    test_put_path()

    test_get()
    test_get_not_has_key()

    test_delete()
    test_delete_not_has_key()
    test_delete_path()

    test_clear()
end

end # module TestFileSystemCache

max_forks_repo_name: seung-lab/SimpleTasks


text:
abstract type EsatMethod end
struct Sonntag1990 <: EsatMethod end
struct Alduchov1996 <: EsatMethod end
struct Allen1998 <: EsatMethod end

"""
    Esat_slope(Tair; Esat_formula, constants)
    Esat_from_Tair(Tair; Esat_formula, constants)
    Esat_from_Tair_deriv(Tair; Esat_formula, constants)

Saturation Vapor Pressure (Esat) and Slope of the Esat Curve

# Arguemtns
- `Tair`:      Air temperature (deg C)
- `Esat_formula=Sonntag1990()`:   Esat_formula to be used. Either
   `Sonntag1990()` (Default), `Alduchov1996()`, or `Allen1998()`.
- `constants=BigleafConstants()`: Dictionary with entry  :Pa2kPa

# Details
Esat (kPa) is calculated using the Magnus equation:

```Esat = a * exp((b * Tair) / (c + Tair)) / 1000```
where the coefficients a, b, c take different values depending on the Esat_formula use
The default values are from Sonntag 1990 (a=611.2, b=17.62, c=243.12). This versi
of the Magnus equation is recommended by the WMO (WMO 2008; p1.4-29). Alternativel
parameter values determined by Alduchov & Eskridge 1996 or Allen et al. 1998 can b
used (see references).
The slope of the Esat curve (``\\Delta``) is calculated as the first derivative of the function:

``\\Delta = {dEsat \\over dTair}``

# Value
- `Esat_from_Tair`: Saturation vapor pressure (kPa)
- `Esat_from_Tair_deriv`: Slope of the saturation vapor pressure curve (kPa K-1)
- `Esat_slope`: A tuple of the two above values
# References
Sonntag D. 1990: Important new values of the physical constants of 1986, vapor
pressure formulations based on the ITS-90 and psychrometric formulae.
Zeitschrift fuer Meteorologie 70, 340-344.

World Meteorological Organization 2008: Guide to Meteorological Instruments
and Methods of Observation (WMO-No.8). World Meteorological Organization,
Geneva. 7th Edition,

Alduchov, O. A. & Eskridge, R. E., 1996: Improved Magnus form approximation of
saturation vapor pressure. Journal of Applied Meteorology, 35, 601-609

Allen, R.G., Pereira, L.S., Raes, D., Smith, M., 1998: Crop evapotranspiration -
Guidelines for computing crop water requirements - FAO irrigation and drainage
paper 56, FAO, Rome.

```@example
Esat_from_Tair(20.0)          # Esat in kPa
Esat_from_Tair_deriv(20.0)    # its derivative to temperature in kPa K-1
Esat_slope(20.0)              # both as a tuple
```
"""
function Esat_slope(Tair::Number; Esat_formula=Sonntag1990(), constants=BigleafConstants())
  Esat = Esat_from_Tair(Tair; Esat_formula, constants)
  Delta = Esat_from_Tair_deriv(Tair; Esat_formula, constants)
  Esat, Delta
end,
function Esat_from_Tair(Tair; Esat_formula=Sonntag1990(), constants=BigleafConstants())
  a,b,c = get_EsatCoef(Esat_formula)
  Esat = a * exp((b * Tair) / (c + Tair)) * constants.Pa2kPa
end,
function Esat_from_Tair_deriv(Tair; Esat_formula=Sonntag1990(), constants=BigleafConstants())
  # slope of the saturation vapor pressure curve
  #Delta = eval(D(expression(a * exp((b * Tair) / (c + Tair))),name="Tair"))
  a,b,c = get_EsatCoef(Esat_formula)
  #Delta_Pa = @. a*(b / (Tair + c) + (-Tair*b) / ((Tair + c)^2))*exp((Tair*b) / (Tair + c))
  Delta_Pa = @. a * (exp((b * Tair)/(c + Tair)) * (b/(c + Tair) - (b * Tair)/(c + Tair)^2))
  Delta = Delta_Pa .* constants.Pa2kPa
end

get_EsatCoef(::Sonntag1990) = (a=611.2,b=17.62,c=243.12)
get_EsatCoef(::Alduchov1996) = (a=610.94,b=17.625,c=243.04)
get_EsatCoef(::Allen1998) = (a=610.8,b=17.27,c=237.3)

"""
    LE_to_ET(LE,Tair)
    ET_to_LE(ET,Tair)

Convert evaporative water flux from mass (ET=evapotranspiration)
             to energy (LE=latent heat flux) units, or vice versa.

# Arguments
- LE   Latent heat flux (W m-2)
- ET   Evapotranspiration (kg m-2 s-1)
- Tair Air temperature (deg C)

# Details
The conversions are given by:
- ``ET = LE/\\lambda``
- ``LE = \\lambda ET``

where ``\\lambda`` is the latent heat of vaporization (J kg-1) as calculated by
[`latent_heat_vaporization`](@ref).

# Examples
```jldoctest; output = false
# LE of 200 Wm-2 and air temperature of 25degC
ET = LE_to_ET(200,25)
≈(ET, 8.19e-5, atol =1e-7)
# output
true
```
"""
function LE_to_ET(LE,Tair)
  lambda = latent_heat_vaporization(Tair)
  ET     = LE/lambda
end,
function ET_to_LE(ET,Tair)
  lambda = latent_heat_vaporization(Tair)
  LE     = ET*lambda
end

"""
    ms_to_mol(G_ms,Tair,pressure; constants=BigleafConstants())
    mol_to_ms(G_mol,Tair,pressure; constants=BigleafConstants())

Converts conductances from mass (m s-1) to molar units (mol m-2 s-1), or vice versa

# Details
The conversions are given by
- ``G_{mol} = G_{ms}  \\, pressure / (Rgas  Tair)``
- ``G_{ms} = G_{mol}  \\, (Rgas  Tair) / pressure``

where Tair is in Kelvin and pressure in Pa (converted from kPa internally).

# References
Jones, H_G_ 1992_ Plants and microclimate: a quantitative approach to
  environmental plant physiology_
  2nd Edition_, Cambridge University Press, Cambridge_ 428

# Examples
```jldoctest; output = false
G_ms,Tair,pressure = 0.005,25.0,100.0
rmol = ms_to_mol(G_ms,Tair,pressure)
≈(rmol, 0.2017, atol =1e-4)
# output
true
```
"""
function ms_to_mol(G_ms,Tair,pressure; constants=BigleafConstants())
  Tair     = Tair + oftype(Tair,constants.Kelvin)
  pressure = pressure * constants.kPa2Pa
  G_mol  = G_ms * pressure / (oftype(Tair, constants.Rgas)  * Tair)
end,
function mol_to_ms(G_mol,Tair,pressure; constants=BigleafConstants())
  Tair     = Tair + oftype(Tair,constants.Kelvin)
  pressure = pressure * constants.kPa2Pa
  G_ms  = G_mol * (oftype(Tair, constants.Rgas) * Tair) / (pressure)
end

"""
    VPD_to_rH(VPD,Tair; ...)
    H_to_VPD(rH,Tair; ...)
    e_to_rH(e,Tair; ...)
    VPD_to_e(VPD,Tair; ...)
    e_to_VPD(e,Tair; ...)
    e_to_q(e,pressure; ...)
    q_to_e(q,pressure; ...)
    q_to_VPD(q,Tair,pressure; ...)
    VPD_to_q(VPD,Tair,pressure; ...)

Conversion between vapor pressure (e), vapor pressure deficit (VPD),
             specific humidity (q), and relative humidity (rH).

# Arguments
- Tair:      Air temperature (deg C)
- pressure:  Atmospheric pressure (kPa)
- e:         Vapor pressure (kPa)
- q:        Specific humidity (kg kg-1)
- VPD:       Vapor pressure deficit (kPa)
- rH:        Relative humidity (-)

All functions accept the optional arguemtns:
- `Esat_formula`: formula used in [`Esat_from_Tair`](@ref)
- `constants`: dictionary from [`BigleafConstants`](@ref) with entries
  eps and Pa2kPa

# Rreferences
Foken, T, 2008: Micrometeorology_ Springer, Berlin, Germany.
"""
function VPD_to_rH(VPD,Tair; Esat_formula=Sonntag1990(),
                      constants=BigleafConstants())
  esat = Esat_from_Tair(Tair; Esat_formula = Esat_formula,constants)
  rH   = 1 - VPD/esat
end,
function rH_to_VPD(rH,Tair; Esat_formula=Sonntag1990(),
                      constants=BigleafConstants())
  if !ismissing(rH) && rH > 1
    @warn("Expected relative humidity (rH) between 0 and 1, but was ", rH)
  end
  esat = Esat_from_Tair(Tair; Esat_formula = Esat_formula,constants)
  VPD  = esat - rH*esat
end,
function e_to_rH(e,Tair; Esat_formula=Sonntag1990(),
                    constants=BigleafConstants())
  esat = Esat_from_Tair(Tair; Esat_formula = Esat_formula,constants)
  if !ismissing(e) && (e > esat + sqrt(eps()))
    @warn("Provided vapour pressure that was higher than saturation_
             Returning rH=1 for those cases.")
  end
  rH  = min(1, e/esat)
end,
function VPD_to_e(VPD,Tair; Esat_formula=Sonntag1990(),
                     constants=BigleafConstants())
  esat = Esat_from_Tair(Tair; Esat_formula = Esat_formula,constants)
  e    = esat - VPD
end,
function e_to_VPD(e,Tair; Esat_formula=Sonntag1990(),
                     constants=BigleafConstants())
  esat = Esat_from_Tair(Tair; Esat_formula = Esat_formula,constants)
  VPD  = esat - e
end,
function e_to_q(e::FT,pressure; constants=BigleafConstants()) where FT
  q = (FT(constants.eps) * e / (pressure - (FT(1)-FT(constants.eps)) * e))
end,
function q_to_e(q::FT,pressure; constants=BigleafConstants()) where FT
  e = (q * pressure / ((FT(1)-FT(constants.eps)) * q + FT(constants.eps)))
end,
function q_to_VPD(q,Tair,pressure; Esat_formula=Sonntag1990(),
                     constants=BigleafConstants())
  esat = Esat_from_Tair(Tair; Esat_formula = Esat_formula,constants)
  e    = q_to_e(q,pressure; constants)
  VPD  = esat - e
end,
function VPD_to_q(VPD,Tair,pressure; Esat_formula=Sonntag1990(),
                     constants=BigleafConstants())
  esat = Esat_from_Tair(Tair; Esat_formula = Esat_formula,constants)
  e    = esat - VPD
  q    = e_to_q(e,pressure;constants)
end

"""
    Rg_to_PPFD(Rg,J_to_mol=4.6,frac_PAR=0.5)
    PPFD_to_Rg(PPFD,J_to_mol=4.6,frac_PAR=0.5)

Conversions between Global Radiation (W m-2) and Photosynthetic Photon Flux
Density (umol m-2 s-1)

# Arguments
- Rg:       Global radiation = incoming short-wave radiation at the surface (W m-2)
- PPFD:     Photosynthetic photon flux density (umol m-2 s-1)
- J_to_mol: Conversion factor from J m-2 s-1 (= W m-2) to umol (quanta) m-2 s-1
- frac_PAR: Fraction of incoming solar irradiance that is photosynthetical-
  active radiation (PAR); defaults to 0.5

# Details
The conversion is given by:

 ``PPFD = Rg * frac_PAR * J_to_mol``

by default, the combined conversion factor (`frac_PAR * J_to_mol`) is 2.3

# Examples
```@example
# convert a measured incoming short-wave radiation of 500 Wm-2 to
# PPFD in umol m-2 s-1 and backwards
Rg_to_PPFD(500)
PPFD_to_Rg(1150)
```
"""
function Rg_to_PPFD(Rg,J_to_mol=4.6,frac_PAR=0.5)
  PPFD = Rg * frac_PAR * J_to_mol
end,
function PPFD_to_Rg(PPFD,J_to_mol=4.6,frac_PAR=0.5)
  Rg = PPFD / frac_PAR / J_to_mol
end

"""
    kg_to_mol(mass, molarMass=BigleafConstants()[:H2Omol])

Conversion between Mass (kg) and Molar Units (mol).
"""
function kg_to_mol(mass, molarMass=BigleafConstants()[:H2Omol])
  moles = mass / molarMass
end

"""
    umolCO2_to_gC(CO2_flux; constants=BigleafConstants())
    gC_to_umolCO2(C_flux; constants=BigleafConstants())

Convert CO2 quantities from (umol CO2 m-2 s-1) to (g C m-2 d-1) and vice versa.

# Arguments
- CO2_flux  CO2 flux (umol CO2 m-2 s-1)
- C_flux    Carbon (C) flux (gC m-2 d-1)
- `constants`: dictionary from [`BigleafConstants`](@ref) with entries:
  Cmol, umol2mol, mol2umol, kg2g, g2kg, says2seconds

# Examples
```@example
umolCO2_to_gC(20)  # gC m-2 d-1
```
"""
function umolCO2_to_gC(CO2_flux::FT; constants=BigleafConstants()) where FT
  C_flux = CO2_flux * FT(constants.umol2mol) * FT(constants.Cmol) *
  FT(constants.kg2g) * FT(constants.days2seconds)
end,
function gC_to_umolCO2(C_flux::FT; constants=BigleafConstants()) where FT
  CO2_flux = (C_flux * FT(constants.g2kg)  / FT(constants.days2seconds)) /
  FT(constants.Cmol)  * FT(constants.mol2umol)
end


text:
# Butcher Tableau
const B = [[1/5],
     [3/40, 9/40],
     [44/45, -56/15, 32/9],
     [19372/6561, -25360/2187, 64448/6561, -212/729],
     [9017/3168, -355/33, 46732/5247, 49/176, -5103/18656],
     [35/384, 0, 500/1113, 125/192, -2187/6784, 11/84]
     ]

# Step size fractions
const nodes = [1/5, 3/10, 4/5, 8/9, 1, 1]

# Weights for 5th order method
const b5 = [5179/57600, 0, 7571/16695, 393/640, -92097/339200, 187/2100, 1/40]
# Weights for 4th order method
const b4 = [35/384, 0, 500/1113, 125/192, -2187/6784, 11/84, 0]
# Error estimate
const errest = b5 .- b4

#Interpolation coefficients
const interpC = hcat(
[1,       -183/64,      37/12,       -145/128],
[0,       0,             0,             0],
[0,        1500/371,     -1000/159,   1000/371],
[0,       -125/32,      125/12,      -375/64],
[0,       9477/3392,    -729/106,    25515/6784],
[0,       -11/7,        11/3,        -55/28],
[0,       3/2,          -4,           5/2],
)


text:
export numderiv

# --------------------- #
# Numerical Derivatives #
# --------------------- #

"""
    ngradient(f, xs::AbstractArray...)

Computes the numerical gradients of `f` w.r.t `xs`. The function `f` must return
a scalar value for this to work. This function is not meant for general usage as
the value of the parameter `δ` has been tuned for this package specifically.

Also, it should be noted that these gradients are highly unstable and should be
used only for confirming the values obtained through other methods. For meaningful
results be sure to use `Float64` as `Float32` is too numerically unstable.
"""
function ngradient(f, xs::AbstractArray...)
    grads = zero.(xs)
    for (x, Δ) in zip(xs, grads), i in 1:length(x)
        δ = 1.0e-11
        tmp = x[i]
        x[i] = tmp - δ/2
        y1 = f(xs...)
        x[i] = tmp + δ/2
        y2 = f(xs...)
        x[i] = tmp
        Δ[i] = (y2 - y1)/δ
    end
    return grads
end

# ---------- #
# Parameters #
# ---------- #

"""
    get_params(x)

Get the parameters from a struct that can be tuned. The output is in
the form of an array.

### Example:

```jldoctest
julia> using RayTracer;

julia> scene = Triangle(Vec3(-1.9f0, 1.6f0, 1.0f0), Vec3(1.0f0, 1.0f0, 0.0f0), Vec3(-0.5f0, -1.0f0, 0.0f0),
                        Material())
Triangle Object:
    Vertex 1 - x = -1.9, y = 1.6, z = 1.0
    Vertex 2 - x = 1.0, y = 1.0, z = 0.0
    Vertex 3 - x = -0.5, y = -1.0, z = 0.0
    Material{Array{Float32,1},Array{Float32,1},Nothing,Nothing,Nothing,Nothing}(x = 1.0, y = 1.0, z = 1.0, x = 1.0, y = 1.0, z = 1.0, x = 1.0, y = 1.0, z = 1.0, Float32[50.0], Float32[0.5], nothing, nothing, nothing, nothing)

julia> RayTracer.get_params(scene)
20-element Array{Float32,1}:
 -1.9
  1.6
  1.0
  1.0
  1.0
  0.0
 -0.5
 -1.0
  0.0
  1.0
  1.0
  1.0
  1.0
  1.0
  1.0
  1.0
  1.0
  1.0
 50.0
  0.5
```
"""
get_params(x::T, typehelp = Float32) where {T<:AbstractArray} = x

get_params(x::T, typehelp = Float32) where {T<:Real} = [x]

get_params(::Nothing, typehelp = Float32) = typehelp[]

get_params(x::T, typehelp = Float32) where {T} =
    foldl((a, b) -> [a; b], [map(i -> get_params(getfield(x, i), typehelp),
                                 fieldnames(T))...])

"""
    set_params!(x, y::AbstractArray)

Sets the tunable parameters of the struct `x`. The index of the last element
set into the struct is returned as output. This may be used to confirm that
the size of the input array was as expected.

### Example:

```jldoctest
julia> using RayTracer;

julia> scene = Triangle(Vec3(-1.9f0, 1.6f0, 1.0f0), Vec3(1.0f0, 1.0f0, 0.0f0), Vec3(-0.5f0, -1.0f0, 0.0f0),
                        Material())
Triangle Object:
    Vertex 1 - x = -1.9, y = 1.6, z = 1.0
    Vertex 2 - x = 1.0, y = 1.0, z = 0.0
    Vertex 3 - x = -0.5, y = -1.0, z = 0.0
    Material{Array{Float32,1},Array{Float32,1},Nothing,Nothing,Nothing,Nothing}(x = 1.0, y = 1.0, z = 1.0, x = 1.0, y = 1.0, z = 1.0, x = 1.0, y = 1.0, z = 1.0, Float32[50.0], Float32[0.5], nothing, nothing, nothing, nothing)

julia> new_params = collect(1.0f0:20.0f0)
20-element Array{Float32,1}:
  1.0
  2.0
  3.0
  4.0
  5.0
  6.0
  7.0
  8.0
  9.0
 10.0
 11.0
 12.0
 13.0
 14.0
 15.0
 16.0
 17.0
 18.0
 19.0
 20.0

julia> RayTracer.set_params!(scene, new_params);

julia> scene
Triangle Object:
    Vertex 1 - x = 1.0, y = 2.0, z = 3.0
    Vertex 2 - x = 4.0, y = 5.0, z = 6.0
    Vertex 3 - x = 7.0, y = 8.0, z = 9.0
    Material{Array{Float32,1},Array{Float32,1},Nothing,Nothing,Nothing,Nothing}(x = 10.0, y = 11.0, z = 12.0, x = 13.0, y = 14.0, z = 15.0, x = 16.0, y = 17.0, z = 18.0, Float32[19.0], Float32[20.0], nothing, nothing, nothing, nothing)
```
"""
function set_params!(x::AbstractArray, y::AbstractArray)
    l = length(x)
    x .= y[1:l]
    return l
end

# FIXME: This update strategy fails for real parameters. Will need to deal with
#        this later. But currently I am not much concered about these parameters.
set_params!(x::T, y::AbstractArray) where {T<:Real} = set_params!([x], y)

function set_params!(x, y::AbstractArray)
    start = 1
    for i in 1:nfields(x)
        isnothing(getfield(x, i)) && continue
        start += set_params!(getfield(x, i), y[start:end])
    end
    return start - 1
end

# ------------------------- #
# Numerical Differentiation #
# ------------------------- #
"""
    numderiv(f, θ)
    numderiv(f, θ::AbstractArray)
    numderiv(f, θ::Real)

Compute the numerical derivates wrt one of the scene parameters.
The parameter passed cannot be enclosed in an Array, thus making
this not a general method for differentiation.

!!! note
    This is not a generalized method for getting the gradients.
    For that please use Zygote. However, this can be used to
    debug you model, incase you feel the gradients obtained by
    other methods is sketchy.
"""
function numderiv(f, θ)
    arr = get_params(θ)
    function modified_f(x::AbstractArray)
        θ₂ = deepcopy(θ)
        set_params!(θ₂, x)
        return f(θ₂)
    end
    grads = ngradient(modified_f, arr)[1]
    grad_θ = deepcopy(θ)
    set_params!(grad_θ, grads)
    return grad_θ
end

numderiv(f, θ::AbstractArray) = ngradient(f, θ)

numderiv(f, θ::Real) = ngradient(f, [θ])


text:
LCR6xxx_Triggering_Dict = Dict()
# Set trigger type
LCR6xxx_Triggering_Dict["edge"]  = 0
LCR6xxx_Triggering_Dict["width"] = 1
# Set trigger source
LCR6xxx_Triggering_Dict["1"]  = 0
LCR6xxx_Triggering_Dict["2"]  = 1
LCR6xxx_Triggering_Dict["3"]  = 2
LCR6xxx_Triggering_Dict["4"]  = 3
LCR6xxx_Triggering_Dict["ext"]  = 4
# Set trigger mode
LCR6xxx_Triggering_Dict["auto"]   = 0
LCR6xxx_Triggering_Dict["normal"] = 1
LCR6xxx_Triggering_Dict["single"] = 2
LCR6xxx_Triggering_Dict["stop"] = 3
# Set trigger slope
LCR6xxx_Triggering_Dict["pos"]    = 0
LCR6xxx_Triggering_Dict["neg"]    = 1
LCR6xxx_Triggering_Dict["either"] = 2
# Trigger coupling
LCR6xxx_Triggering_Dict["dc"]    = 0
LCR6xxx_Triggering_Dict["ac"]    = 1
#LCR6xxx_Triggering_Dict["gnd"]   = 2
LCR6xxx_Triggering_Dict["lfrej"] = 2
LCR6xxx_Triggering_Dict["hfrej"] = 3

# Set trigger type
function set_trig_type(obj::T;tp::String="edge") where {T<:LCR6xxx}
	tp = LCR6xxx_Triggering_Dict[tp]
	write(obj,"VBS app.Acquisition.Trigger.Type=$tp")
end
# Set trigger source
function set_trig_src(obj::T;src::String="1") where {T<:LCR6xxx}
	src = LCR6xxx_Triggering_Dict[src]
	write(obj,"VBS app.Acquisition.Trigger.Edge.Source=$src")
end
# Set trigger mode
function set_trig_mode(obj::T;mode::String="normal") where {T<:LCR6xxx}
	mode = LCR6xxx_Triggering_Dict[mode]
	write(obj,"VBS app.Acquisition.TriggerMode=$mode")
end
# Set trigger slope
function set_trig_slope(obj::T;slp::String="pos") where {T<:LCR6xxx}
	slp = LCR6xxx_Triggering_Dict[slp]
	write(obj,"VBS app.Acquisition.Trigger.Edge.Slope=$slp")
end
# Set trigger level
function set_trig_lev(obj::T;lev=1) where {T<:LCR6xxx}
	write(obj,"VBS app.Acquisition.Trigger.Edge.Level=$lev")
end
# Set trigger holdoff
function set_trig_holdoff(obj::T;hf=1e-9) where {T<:LCR6xxx}
	write(obj,"VBS app.Acquisition.Trigger.Edge.HoldoffTime=$hf")
end
# Trigger delay time
function set_hoffs_t(obj::T;href=1e-3) where {T<:LCR6xxx}
	write(obj,"VBS app.Acquisition.Horizontal.HorOffset=$href")
end
# Set trigger coupling
function set_trig_cpl(obj::T;src::String="1",cpl::String="dc") where {T<:LCR6xxx}
	cpl = LCR6xxx_Triggering_Dict[cpl]
	write(obj,"VBS app.Acquisition.Trigger.C$(src)Coupling=$cpl")
end


text:
function largez(k::Integer)
#
#  test polynomial suggested by Z Zeng
#
    p = [
         1.0
        -0.7
        -0.19
         0.177
        -0.7364
        -0.43780
        -0.952494
        -0.2998258
        -0.00322203
        -0.328903811
        -0.4959527435
        -0.9616679762
         0.4410459281
         0.1090273141
         0.6868094008
         0.0391923826
         0.0302248540
         0.6603775863
        -0.1425784968
        -0.3437618593
         0.4357949015
        ];
    z = [
         0.5 + im
         0.5 - im
        -1 + 0.2im
        -1 - 0.2im
        -0.1 + im
        -0.1 - im
         0.8 + 0.6im
         0.8 - 0.6im
        -0.7 + 0.7im
        -0.7 - 0.7im
         1.4
        -0.4 + 0.9im
        -0.4 - 0.9im
         0.9
        -0.8 + 0.3im
        -0.8 - 0.3im
         0.3 + 0.8im
         0.3 - 0.8im
         0.6 + 0.4im
         0.6 - 0.4im
        ]

    mult = 1
    for i = 1:k
        p = conv(p, p)
        mult *= 2
    end

    z = [z mult*ones(20)]
    p, PolyZeros(z)
end

largez01() = largez(0)
largez02() = largez(1)
largez03() = largez(2)
largez04() = largez(3)
largez05() = largez(4)
largez06() = largez(5)


text:
module ValueNoise
using Base.Iterators: countfrom

include("value_noise.jl")
export value_noise

end # module


text:
include("./ParserCombinator.jl")
using ParserCombinator

@rule(E) do
  A()
end

@rule(A) do
  chainl(
    M(),
    (p_string("+") > (op) -> (lhs, rhs) -> lhs + rhs) /
    (p_string("-") > (op) -> (lhs, rhs) -> lhs - rhs)
  )
end

@rule(M) do
  chainl(
    P(),
    (p_string("*") > (op) -> (lhs, rhs) -> lhs * rhs) /
    (p_string("/") > (op) -> (lhs, rhs) -> lhs / rhs)
  )
end

@rule(P) do
  (
    (p_string("(") + E() + p_string(")")) >
    values -> values[1][2]
  ) / N()
end

@rule(N) do
  regex(r"[0-9]+") > (n) -> parse(Int64, n)
end

const parser = E()
println(parser.parse("(1+2)*(3+4)"))

max_stars_repo_name: kmizu/mlcombinators


text:
# Practical Machine learning
# Artificial Neural Network
# Chapter 11

function computeNumericalGradient(J, theta)

  # Notes: The following code implements numerical gradient checking, and
  #        returns the numerical gradient.It sets numgrad(i) to (a numerical
  #        approximation of) the partial derivative of J with respect to the
  #        i-th input argument, evaluated at theta. (i.e., numgrad(i) should
  #        be the (approximately) the partial derivative of J with respect
  #        to theta(i).)

  numgrad = COMPUTENUMERICALGRADIENT(J, theta)
  perturb = zeros(size(theta))
  e = 1e-4
  for p in 1:length(theta)
    # Set perturbation vector
    perturb[p] = e
    loss1, _ = J(theta - perturb)
    loss2, _ = J(theta + perturb)
    # Compute Numerical Gradient
    numgrad[p] = (loss2 - loss1) / 2e
    perturb[p] = 0
  end
  return numgrad
end


text:
function e004(N)
    return largest_palindrome(N)
end

function is_palindrome(i::Integer) :: Bool
    return parse(Int, reverse(repr(i))) == i
end

function largest_palindrome(factor_length::Integer)
    largest::Integer = 0
    for i::Integer = 10^(factor_length-1):10^(factor_length)-1
        for j::Integer = 10^(factor_length-1):i
            product = i*j
            if is_palindrome(product) && product > largest
                largest = product
            end
        end
    end
    return largest
end

function solution004()
    N = 3
    answer = e004(N)
    println("The answer to problem 004 is $answer")
end

max_stars_repo_name: mithodin/ttpe

max_issues_repo_name: mithodin/ttpe

max_forks_repo_name: mithodin/ttpe


text:
using Revise
using ADCME
using NNFEM
using JLD2
using PyPlot
reset_default_graph()

include("nnutil.jl")

loss = constant(0.0)
for tid = [1,2,4,5]
    strain, stress = read_strain_stress("Data/$(tid).dat")
    X, Y = prepare_strain_stress_data1D(strain, stress )
    y = squeeze(nn(constant(X[:,1:1]), constant(X[:,2:2]), constant(X[:,3:3])))
    global loss += sum((y - Y)^2)
end
sess = Session(); init(sess)
@show run(sess,loss)

if !isdir("Data/$(nntype)")
    mkdir("Data/$(nntype)")
end

loss_all = []
lossfile = matopen("Data/$(nntype)/learned_nn$(idx)_loss_$(nn_init_id).txt", "w")
for i = 1:10
    println("************************** Outer Iteration = $i ************************** ")
    loss_ = BFGS!(sess, loss, 3000)
    push!(loss_all, loss_)
    ADCME.save(sess, "Data/$(nntype)/learned_nn$(idx)_ite$(i)_$(nn_init_id).mat")
end

write(lossfile, "loss", vcat(loss_all...))
close(lossfile)

# error()

# tid = 3
# strain, stress = read_strain_stress("Data/$(tid).dat")
# X, Y = prepare_strain_stress_data1D(strain, stress )
# y = squeeze(nn(constant(X[:,1:1]), constant(X[:,2:2]), constant(X[:,3:3])))
# out = run(sess, y)
# plot(X[:,1], out,"+", label="NN")
# plot(X[:,1], Y, ".", label="Exact")
# #legend()
# savefig("nnlearn_$(nntype)_nn$(idx)_truss1d_stress$tid.png")

# @load "Data/domain.jld2" domain
# X, Y = prepare_strain_stress_data1D(domain)
# x = constant(X)
# y = squeeze(nn(constant(X[:,1]), constant(X[:,2]), constant(X[:,3])))
# sess = Session(); init(sess)
# close("all")
# ADCME.load(sess, "Data/learned_nn.mat")
# out = run(sess, y)
# plot(X[:,1], out,"+", label="NN")
# plot(X[:,1], Y, ".", label="Exact")
# ADCME.load(sess, "Data/trained_nn_fem.mat")
# out = run(sess, y)
# plot(X[:,1], out,">", label="End-to-end")
# legend()


text:
export IntersectionArray,
       array

# ================================
# Intersection of an array of sets
# ================================

"""
    IntersectionArray{N, S<:LazySet{N}} <: LazySet{N}

Type that represents the intersection of a finite number of sets.

### Fields

- `array` -- array of sets

### Notes

This type assumes that the dimensions of all elements match.

The `EmptySet` is the absorbing element for `IntersectionArray`.

The intersection preserves convexity: if the set arguments are convex, then
their intersection is convex as well.

Constructors:

- `IntersectionArray(array::Vector{<:LazySet})` -- default constructor

- `IntersectionArray([n]::Int=0, [N]::Type=Float64)`
 -- constructor for an empty sum with optional size hint and numeric type
"""
struct IntersectionArray{N, S<:LazySet{N}} <: LazySet{N}
   array::Vector{S}
end

isoperationtype(::Type{<:IntersectionArray}) = true
isconvextype(::Type{IntersectionArray{N, S}}) where {N, S} = isconvextype(S)

# constructor for an empty sum with optional size hint and numeric type
function IntersectionArray(n::Int=0, N::Type=Float64)
   arr = Vector{LazySet{N}}()
   sizehint!(arr, n)
   return IntersectionArray(arr)
end

# Universe is the neutral element for IntersectionArray
@neutral(IntersectionArray, Universe)

# EmptySet is the absorbing element for IntersectionArray
@absorbing(IntersectionArray, EmptySet)

# add functions connecting Intersection and IntersectionArray
@declare_array_version(Intersection, IntersectionArray)

"""
    array(ia::IntersectionArray)

Return the array of an intersection of a finite number of sets.

### Input

- `ia` -- intersection of a finite number of sets

### Output

The array of an intersection of a finite number of sets.
"""
function array(ia::IntersectionArray)
   return ia.array
end

# --- LazySet interface functions ---

"""
    dim(ia::IntersectionArray)

Return the dimension of an intersection of a finite number of sets.

### Input

- `ia` -- intersection of a finite number of sets

### Output

The ambient dimension of the intersection of a finite number of sets.
"""
function dim(ia::IntersectionArray)
   return length(ia.array) == 0 ? 0 : dim(ia.array[1])
end

"""
    σ(d::AbstractVector, ia::IntersectionArray)

Return the support vector of an intersection of a finite number of sets in a
given direction.

### Input

- `d`  -- direction
- `ia` -- intersection of a finite number of sets

### Output

The support vector in the given direction.
If the direction has norm zero, the result depends on the individual sets.
"""
function σ(d::AbstractVector, ia::IntersectionArray)
    X = concretize(ia)
    return σ(d, X)
end

"""
    isbounded(ia::IntersectionArray)

Determine whether an intersection of a finite number of sets is bounded.

### Input

- `ia` -- intersection of a finite number of sets

### Output

`true` iff the intersection is bounded.

### Algorithm

We first check if any of the wrapped sets is bounded.
Otherwise, we check boundedness via [`LazySets._isbounded_unit_dimensions`](@ref).
"""
function isbounded(ia::IntersectionArray)
   if any(isbounded, ia.array)
       return true
   end
   return _isbounded_unit_dimensions(ia)
end

function isboundedtype(::Type{<:IntersectionArray{N, S}}) where {N, S}
    return isboundedtype(S)
end

"""
    ∈(x::AbstractVector, ia::IntersectionArray)

Check whether a given point is contained in an intersection of a finite number
of sets.

### Input

- `x`  -- point/vector
- `ia` -- intersection of a finite number of sets

### Output

`true` iff ``x ∈ ia``.
"""
function ∈(x::AbstractVector, ia::IntersectionArray)
   for S in ia.array
       if x ∉ S
           return false
       end
   end
   return true
end

"""
    constraints_list(ia::IntersectionArray{N}) where {N}

Return the list of constraints of an intersection of a finite number of
(polyhedral) sets.

### Input

- `ia` -- intersection of a finite number of (polyhedral) sets

### Output

The list of constraints of the intersection.

### Notes

We assume that the underlying sets are polyhedral, i.e., offer a method
`constraints_list`.

### Algorithm

We create the polyhedron from the `constraints_list`s of the sets and remove
redundant constraints.
"""
function constraints_list(ia::IntersectionArray{N}) where {N}
   constraints = Vector{LinearConstraint{N, Vector{N}}}() # TODO: use vector type of ia
   for X in array(ia)
       clist_X = _normal_Vector(X)
       append!(constraints, clist_X)
   end
   remove_redundant_constraints!(constraints)
   return constraints
end

function concretize(ia::IntersectionArray)
    a = array(ia)
    @assert !isempty(a) "an empty intersection is not allowed"
    X = ia
    @inbounds for (i, Y) in enumerate(a)
        if i == 1
            X = concretize(Y)
        else
            X = intersection(X, concretize(Y))
        end
    end
    return X
end


text:
@testset "DMD" begin
    # Create some linear data
    A = [0.9 -0.2; 0.0 0.2]
    y = [[10.; -10.]]
    for i in 1:10
        push!(y, A*y[end])
    end
    X = hcat(y...)

    estimator = DMD(X)
    # TODO export
    #@test isstable(estimator)
    @test operator(estimator) ≈ A
    @test eigvals(estimator) ≈ eigvals(A)
    @test eigvecs(estimator) ≈ eigvecs(A)
    @test is_discrete(estimator)
    # TODO add Array{T,2}
    #@test estimator(X[:, 1]) ≈ X[:, 2]
    @test_nowarn update!(estimator, X[:, end-1], X[:,end])
    @test_throws AssertionError outputmap(estimator)

    @variables u[1:2]
    b = Basis(u, u)
    sys = EDMD(X[:, 1:3], b)
    @test_nowarn update!(sys, X[:, 3:end-1], X[:, 4:end], threshold = 1e-17)
    @test operator(sys) ≈ A

    # Test for a large system with reduced svd
    dd(u, p, t) = -0.9*ones(length(u))*u[1] - 0.05*u
    u0 = randn(100)*100.0
    u0[1] = 100.0
    prob = ODEProblem(dd, u0, (0.0, 100.0))
    sol = solve(prob, Tsit5(), saveat = 0.1)

    X = sol[:,:]
    DX = sol(sol.t, Val{1})[:,:]

    t = DMDSVD(0.01)
    d = gDMD(X, DX, alg = t)

    test = ODEProblem(d, u0, (0.0, 100.0))
    sol_ = solve(test, Tsit5(), saveat = 0.1)
    @test norm(sol-sol_, Inf) < 2.0

    d2 = gDMD(sol.t, X)
    d3 = gDMD(sol.t, X, dt = 0.2)
    @test generator(d) ≈ generator(d2) atol = 1e-1
    @test generator(d) ≈ generator(d3) atol = 1e-1

    K = -0.5*I + [0 0 -0.2; 0.1 0 -0.1; 0. -0.2 0]
    F = qr(randn(20, 3))
    Q = F.Q[:, 1:3]
    dudt(u, p, t) = K*u
    prob = ODEProblem(dudt, [10.0; 0.3; -5.0], (0.0, 10.0))
    sol_ = solve(prob, Tsit5(), saveat = 0.01)
    alg = TOTALDMD(3, DMDPINV())
    X = Q*sol_[:,:] + 1e-3*randn(20, 1001)
    DX = Q*sol_(sol_.t, Val{1})[:,:] + 1e-3*randn(20, 1001)
    approx = gDMD(X, DX, alg = alg)
    @test Q'*generator(approx)*Q ≈ K atol = 1e-1
    @test Q*K*Q' ≈ generator(approx) atol = 1e-1
    alg = TOTALDMD(0.01, DMDSVD())
    approx = gDMD(X, DX, alg = alg)
    @test Q'*generator(approx)*Q ≈ K atol = 1e-1
    @test Q*K*Q' ≈ generator(approx) atol = 1e-1

end

@testset "EDMD" begin
    # Test for linear system
    function linear_sys(u, p, t)
        x = -0.9*u[1]
        y = -0.3*u[2]
        return [x;y]
    end

    u0 = [π; 1.0]
    tspan = (0.0, 50.0)
    prob = DiscreteProblem(linear_sys, u0, tspan)
    sol = solve(prob,FunctionMap())

    @variables u[1:2]
    h = [u[1]; u[2]; sin(u[1]); cos(u[1]); u[1]*u[2]; u[2]^2]
    basis = Basis(h, u)

    estimator = EDMD(sol[:,:], basis, alg = DMDSVD())
    @test basis == estimator.basis
    basis_2 = reduce_basis(estimator, threshold = 1e-5)
    @test size(basis_2)[1] < size(basis)[1]

    estimator_2 = EDMD(sol[:,:], basis_2)
    p1 = DiscreteProblem(estimator, u0, tspan)
    s1 = solve(p1,FunctionMap())
    p2 = DiscreteProblem(estimator_2, u0, tspan)
    s2 = solve(p2,FunctionMap())
    # TODO add linear dynamics ?
    #p3 = DiscreteProblem(linear_dynamics(estimator_2), estimator_2(u0), tspan, [])
    #s3 = solve(p3,FunctionMap())

    @test sol[:,:] ≈ s1[:,:]
    @test sol[:,:] ≈ s2[:,:]
    @test eigvals(estimator_2) ≈ [-0.9; -0.3]
    @test isstable(estimator_2)

    # Test for nonlinear system
    function nonlinear_sys(du, u, p, t)
        du[1] = 0.9u[1] + 0.1u[2]^2
        du[2] = sin(u[1]) - 0.1u[1]
    end

    prob = DiscreteProblem(nonlinear_sys, u0, tspan)
    sol = solve(prob,FunctionMap())
    estimator = EDMD(sol[:,:], basis, alg = DMDPINV())
    sys = ODESystem(estimator)
    #@test isa(sys, ODESystem)
    dudt = ODEFunction(sys)
    p4 = DiscreteProblem(dudt, u0, tspan)

    s4 = solve(p4,FunctionMap())
    @test sol[:,:] ≈ s4[:,:]

    function slow_manifold(du, u, p, t)
      du[1] = p[1]*u[1]
      du[2] = p[2]*(u[2]-u[1]^2)
    end

    u0 = [3.0; -2.0]
    tspan = (0.0, 10.0)
    p = [-0.05, -0.9]

    problem = ODEProblem(slow_manifold, u0, tspan, p)
    sol_slow = solve(problem, Tsit5(), saveat = 0.2)

    X = Array(sol_slow)
    # This enforces more accurate results
    DX = similar(X)
    for (i,dx) in enumerate(eachcol(DX))
        slow_manifold(dx, X[:, i], p, sol_slow.t[i])
    end
    basis = Basis([u; u[1]^2], u)
    estimator2 = gEDMD(X, DX, basis)
    A_analytical = [p[1] 0 0; 0 p[2] -p[2]; 0 0 2*p[1]]
    outputmap(estimator2)
    generator(estimator2)
    @test generator(estimator2) ≈ [p[1] 0 0; 0 p[2] -p[2]; 0 0 2*p[1]] atol =  1e-3
    @test outputmap(estimator2) ≈ [1 0 0 ; 0 1 0] atol = 1e-3
    @test abs.(modes(estimator2)) ≈ abs.(eigvecs(A_analytical)) atol = 1e-3
    @test frequencies(estimator2) ≈ eigvals(A_analytical) atol = 1e-3

    function nonlinear_sys2(du, u, p, t)
        du[1] = u[2]
        du[2] = -0.9u[1]
    end

    prob_nl = ODEProblem(nonlinear_sys2, u0, (0.0, 10.0))
    sol_nl = solve(prob_nl, Tsit5())

    DX = sol_nl(sol_nl.t, Val{1})[:,:]
    X = Array(sol_nl)

    basis = Basis(u, u)
    estimator3 = gEDMD(X, DX, basis)
    @test generator(estimator3) ≈ [0  1.0; -0.9 0]
    @test outputmap(estimator3) ≈ I(2)

    d2 = gEDMD(sol_nl.t, X, basis)
    d3 = gEDMD(sol_nl.t, X, basis, dt = 0.1)
    @test generator(estimator3) ≈ generator(d2) atol = 1e-1
    @test generator(estimator3) ≈ generator(d3) atol = 1e-1

end

@testset "DMDc" begin
    # Define measurements from unstable system with known control input
    X = [4 2 1 0.5 0.25; 7 0.7 0.07 0.007 0.0007]
    U = [-4 -2 -1 -0.5]
    B = Float32[1; 0]

    # But with a little more knowledge
    sys = DMDc(X, U, B = B)
    @test operator(sys) ≈[1.5 0; 0 0.1]
    @test inputmap(sys) ≈ [1.0; 0.0]
    @test !isstable(sys)
    @test_nowarn eigen(sys)
    sys2 = DMDc(X, U, B = B, alg = DMDSVD())
    @test operator(sys2) ≈ operator(sys)
    sys3 = gDMDc(X[:, 1:end-1], X[:, 2:end], U, B = B)
    @test generator(sys3) ≈ [1.5 0; 0 0.1]

    A = [0.9 0.3; -0.2 0.9]
    B = [1.0 0;0 0.5]
    X = zeros(Float64, 2, 11)
    U = zeros(Float64, 2, 10)
    X[:, 1] = [1.0; -3.0]
    for i in 1:10
        U[:, i] .= [sin(2.0/5.0*i); cos(5.0/7.0*i)]
        X[:, i+1] .= A*X[:,i]+B*[sin(2.0/5.0*i); cos(5.0/7.0*i)]
    end

    sys4 = DMDc(X, U)
    @test operator(sys4) ≈ A
    @test inputmap(sys4) ≈ B
    sys5 = gDMDc(X[:, 1:end-1], X[:, 2:end], U[:, 1:end])
    @test_throws AssertionError operator(sys5)
    @test generator(sys5) ≈ A
    @test inputmap(sys5) ≈ B

    sys6 = gDMDc(collect(0.0:9.0), X[:, 1:end-1], U, B = B)
    @test exp(generator(sys6)) ≈ A atol = 1e-1
    @test inputmap(sys6) ≈ B

end


text:
using Plots, UnPack, RollingFunctions

"""
    p1,p2,p3=plot_final(outputs,parameters)
"""
function plot_final(outputs, parameters)

    @unpack U0, UA, THETA, SST, TO, ustar, psimh, psixh = outputs
    @unpack xi, rd, rh, re, SW, LW, SH, LH = outputs
    @unpack dt, n, ndays, cp, dz, rhow0 = parameters

    k1 = plot(
        rollmean(U0[1, :], Int(24 * 60 * 60 / dt)),
        xticks = (0:24*60*60/dt:n, string.(0:ndays)),
        ylabel = "U0",
    )
    k2 = plot(
        rollmean(UA[:, end], Int(24 * 60 * 60 / dt)),
        xticks = (0:24*60*60/dt:n, string.(0:ndays)),
        ylabel = "Utop",
    )
    k3 = plot(
        rollmean(UA[:, end] - U0[1, :], Int(24 * 60 * 60 / dt)),
        xticks = (0:24*60*60/dt:n, string.(0:ndays)),
        ylabel = "Utop-U0",
    )
    k4 = plot(
        rollmean(UA[:, 1], Int(24 * 60 * 60 / dt)),
        xticks = (0:24*60*60/dt:n, string.(0:ndays)),
        ylabel = "U",
    )
    k5 = plot(
        rollmean(THETA[:, 2], Int(24 * 60 * 60 / dt)),
        xticks = (0:24*60*60/dt:n, string.(0:ndays)),
        ylabel = "THETA",
    )
    k6 = plot(
        rollmean(SST, Int(24 * 60 * 60 / dt)),
        xticks = (0:24*60*60/dt:n, string.(0:ndays)),
        ylabel = "SST",
    )
    k7 = plot(
        rollmean(TO[:, 10], Int(24 * 60 * 60 / dt)),
        xticks = (0:24*60*60/dt:n, string.(0:ndays)),
        ylabel = "T",
    )
    #tt=string("DZ=",dz)
    plt1 = plot(k1, k2, k3, k4, k5, k6, k7, layout = (7, 1), legend = false)

    k1 = plot(
        rollmean(xi[1, :], Int(24 * 60 * 60 / dt)),
        xticks = (0:24*60*60/dt:n, string.(0:ndays)),
        title = "xi",
    )
    k2 = plot(
        rollmean(ustar[1, :], Int(24 * 60 * 60 / dt)),
        xticks = (0:24*60*60/dt:n, string.(0:ndays)),
        title = "ustar",
    )
    k3 = plot(
        rollmean(psimh[1, :], Int(24 * 60 * 60 / dt)),
        xticks = (0:24*60*60/dt:n, string.(0:ndays)),
        title = "psimh",
    )
    k4 = plot(
        rollmean(psixh[1, :], Int(24 * 60 * 60 / dt)),
        xticks = (0:24*60*60/dt:n, string.(0:ndays)),
        title = "psixh",
    )
    k5 = plot(
        rollmean((rd[1, :] .^ 2), Int(24 * 60 * 60 / dt)),
        xticks = (0:24*60*60/dt:n, string.(0:ndays)),
        title = "CD",
    )
    k6 = plot(
        rollmean((rh[1, :] .* rd[1, :]), Int(24 * 60 * 60 / dt)),
        xticks = (0:24*60*60/dt:n, string.(0:ndays)),
        title = "CH",
    )
    k7 = plot(
        rollmean((re[1, :] .* rd[1, :]), Int(24 * 60 * 60 / dt)),
        xticks = (0:24*60*60/dt:n, string.(0:ndays)),
        title = "CE",
    )
    #tt=string("DZ=",dz)
    plt2 = plot(k1, k2, k3, k4, k5, k6, k7, layout = (7, 1), legend = false)

    k1 = plot(
        rollmean(SW[:, 1] / cp / (dz * rhow0) * dt, Int(24 * 60 * 60 / dt)),
        xticks = (0:24*60*60/dt:n, string.(0:ndays)),
        title = "SW",
    )
    k2 = plot(
        rollmean(LW[1, :] / cp / (dz * rhow0) * dt, Int(24 * 60 * 60 / dt)),
        xticks = (0:24*60*60/dt:n, string.(0:ndays)),
        title = "LW",
    )
    k3 = plot(
        rollmean(SH[1, :] / cp / (dz * rhow0) * dt, Int(24 * 60 * 60 / dt)),
        xticks = (0:24*60*60/dt:n, string.(0:ndays)),
        title = "SH",
    )
    k4 = plot(
        rollmean(LH[1, :] / cp / (dz * rhow0) * dt, Int(24 * 60 * 60 / dt)),
        xticks = (0:24*60*60/dt:n, string.(0:ndays)),
        title = "LH",
    )
    k5 = plot(
        rollmean(
            (SW[:, 1] - LW[1, :] - SH[1, :] - LH[1, :]) / cp / (dz * rhow0) * dt,
            Int(24 * 60 * 60 / dt),
        ),
        xticks = (0:24*60*60/dt:n, string.(0:ndays)),
        title = "QNET",
    )
    plt3 = plot(k1, k2, k3, k4, k5, layout = (5, 1), legend = false)
    return plt1, plt2, plt3
end

function plot_state()
    f1 = plot(-ZO[1:end-1], TO[i, :], ylabel = ("TO"))#,ylims=(19,21)
    f2 = plot(-ZO[1:end-1], UO[i, :], ylabel = ("UO"))
    f3 = plot(-ZO[1:end-1], kom, ylabel = ("KOm"))
    f4 = plot(ZA, UA[i, :], ylabel = ("UA"))
    f5 = plot(ZA, Tv .- 273.15, ylabel = ("Tv")) #ylim([19 21]);...
    f6 = plot(ZA, qA[i, :], ylabel = ("q"))
    f7 = plot(ZA, kam, ylabel = ("KAm"))
    plt =
        plot(f1, f2, f3, f4, f5, f6, f7, layout = (7, 1), legend = false, titlefontsize = 6)
    display(plt)
end


text:
abstract type MHV2RuleModel <: MixingRule end

struct MHV2Rule{γ} <: MHV2RuleModel
    components::Array{String,1}
    activity::γ
    references::Array{String,1}
end

@registermodel MHV2Rule
export MHV2Rule
function MHV2Rule(components::Vector{String}; activity = Wilson, userlocations::Vector{String}=String[],activity_userlocations::Vector{String}=String[], verbose::Bool=false)
    init_activity = activity(components;userlocations = activity_userlocations,verbose)

    references = ["10.1016/0378-3812(90)85053-D"]
    model = MHV2Rule(components, init_activity,references)
    return model
end

MHV2q(::PRModel) = (-0.4347,-0.003654)
MHV2q(::RKModel) = (-0.4783,-0.0047)
function mixing_rule(model::Union{PRModel,RKModel},V,T,z,mixing_model::MHV2RuleModel,α,a,b,c)
    n = sum(z)
    invn = (one(n)/n)
    invn2 = invn^2
    g_E = excess_gibbs_free_energy(mixing_model.activity,1e5,T,z)*invn
    b̄ = dot(z,Symmetric(b),z) * invn2
    c̄ = dot(z,c)*invn
    #ᾱ = a.*sqrt.(α.*α')./(b*R̄*T)
    q1,q2 = MHV2q(model)
    Σlogb = zero(first(z))
    Σab = zero(T+first(z))
    Σab2 = Σab
    for i ∈ @comps
        ᾱ  = a[i,i]*α[i]/(b[i,i]*R̄*T)
        zia = z[i]*ᾱ
        Σab += zia
        Σab2 += zia*ᾱ
        Σlogb += z[i]*log(b̄/b[i,i])
    end
    Σab *= invn
    Σab2 *= invn
    Σlogb *= invn
    c  = -q1*Σab-q2*Σab2-g_E/(R̄*T)-Σlogb
    ā = b̄*R̄*T*(-q1-sqrt(q1^2-4*q2*c))/(2*q2)
    return ā,b̄,c̄
end
#=
function mixing_rule(model::PRModel,V,T,z,mixing_model::MHV2RuleModel,α,a,b,c)
    n = sum(z)
    x = z./n
    invn2 = (one(n)/n)^2
    g_E = excess_gibbs_free_energy(mixing_model.activity,1e5,T,z) / n
    b̄ = dot(z,Symmetric(b),z) * invn2
    c̄ = dot(z,c)/n

    ᾱ = a.*sqrt.(α.*α')./(b*R̄*T)

    q1 = -0.4347
    q2 = -0.003654
    c  = -q1*sum(x[i]*ᾱ[i,i] for i ∈ @comps)-q2*sum(x[i]*ᾱ[i,i]^2 for i ∈ @comps)-g_E/(R̄*T)-sum(x[i]*log(b̄/b[i,i]) for i ∈ @comps)

    ā = b̄*R̄*T*(-q1-sqrt(q1^2-4*q2*c))/(2*q2)
    return ā,b̄,c̄
end=#


text:
using DataFrames

"""
Holds information about a single write
"""
struct WriteOperation
  data::AbstractDataFrame
  format::String
  key::String
  filename::String
  function WriteOperation(data, format, key="", filename="data.$format")
    return new(data, format, key, filename)
  end
end

"""
Write a DataFrame to a path using a user-provided function
  - writer(data::DataFrame, path::String)
"""
function write(op::WriteOperation, path::String, writer::Function)
  dir = "$path/$(op.key)"
  !isdir(dir) && mkpath(dir)
  filepath = "$dir/$(op.filename)"
  try
    writer(op.data, filepath)
    return (true, nothing)
  catch err
    rm(filepath; force=true)
    return (false, err)
  end
end

"""
Container of multiple `WriteOpration`s.
"""
struct WriteTransaction
  ops::Vector{WriteOperation}
  format::String
  dest::String
  tmp::String
  bkp::String
  function WriteTransaction(ops, format, dest)
    tmp_dir = "$(ensure_path(dest))/.tmp"
    bkp_dir = "$(ensure_path(dest))/.bkp"
    return new(ops, format, dest, tmp_dir, bkp_dir)
  end
end

"""
Main function that runs a transaction.
"""
function execute(tr::WriteTransaction, writer::Function)
  has_backup, bkps = backup(tr)
  success, err = write_tmp(tr, writer)
  if !success
    cleanup(tr)
    return (
      false, ErrorException("Failed to execute transaction. Write to temp directory failed")
    )
  end
  try
    commit(tr)
    #@debug "Transaction at '$(tr.dest)' with $(length(tr.ops)) atomic operations executed sucesfully."
    return (true, nothing)
  catch err
    rollback(bkps)
    return (false, err)
  finally
    cleanup(tr)
  end
end

"""
Write data to a temporary location using a user-provided function.
"""
function write_tmp(tr::WriteTransaction, writer::Function)
  isdir(tr.tmp) && rm(tr.tmp; recursive=true, force=true)
  mkpath(tr.tmp)
  tmp_writes = Channel(length(tr.ops))
  @sync begin
    for op in tr.ops
      Threads.@spawn begin
        ok, err = write(op, tr.tmp, writer)
        push!(tmp_writes, (success=ok, error=err))
      end
    end
  end
  close(tmp_writes)
  if all(map(w -> w.success, tmp_writes))
    return (true, nothing)
  else
    msg = "At least one write failed. Can't commit transaction."
    @warn msg
    return (false, ErrorException(msg))
    # for write in tmp_writes
    #   !isnothing(write.error) && return (false, ErrorException(msg))
    # end
  end
end

"""
Create a pre-transaction backup by copying everything from src to destination.
"""
function backup(tr::WriteTransaction)
  !isdir(tr.dest) && return (false, nothing)
  try
    isdir(tr.bkp) && rm(tr.bkp; recursive=true, force=true)
    mkpath(tr.bkp)
    bkps = @chain begin
      tr.ops
      map(
        op -> (
          key=!isempty(op.key) ? "$(tr.bkp)/$(op.key)" : tr.bkp,
          src=if !isempty(op.key)
            "$(tr.dest)/$(op.key)/$(op.filename)"
          else
            "$(tr.dest)/$(op.filename)"
          end,
          dest=if !isempty(op.key)
            "$(tr.bkp)/$(op.key)/$(op.filename)"
          else
            "$(tr.bkp)/$(op.filename)"
          end,
        ),
        _,
      )
    end
    @sync begin
      for bkp in bkps
        Threads.@spawn begin
          mkpath(bkp.key)
          cp(bkp.src, bkp.dest)
        end
      end
    end
    return (true, bkps)
  catch err
    return (false, err)
  end
end

"""
Copies the file from temporary to permanent location.
"""
function commit(tr::WriteTransaction)
  dir_mapping = map(
    key -> (
      src=!isempty(key) ? "$(tr.tmp)/$key" : tr.tmp,
      dest=!isempty(key) ? "$(tr.dest)/$key" : tr.dest,
    ),
    readdir(tr.tmp),
  )
  @sync begin
    for x in dir_mapping
      Threads.@spawn begin
        if isfile(x.dest)
          cp(x.src, x.dest; force=true)
        else
          !isdir(x.dest) && mkpath(x.dest)
          cp(x.src, x.dest; force=true)
        end
      end
    end
  end
end

"""
Uses the backup to revert to the original state.
"""
function rollback(backups)
  isnothing(backups) && return nothing
  @chain begin
    backups
    foreach(x -> begin
      !isdir(x.src) && mkpath(x.src)
      cp(x.dest, x.src; force=true)
    end, _)
  end
end

"""
Removes all temporary and backup files.
"""
function cleanup(tr::WriteTransaction)
  @chain begin
    [tr.bkp, tr.tmp]
    foreach(dir -> isdir(dir) && rm(dir; recursive=true, force=true), _)
  end
end