text:
import Base:
    ==,
    hash

##### state, action, and observation spaces and related functions #####

struct RockSampleState
    index::Int64
end
RockSampleState() = RockSampleState(-1)

struct RockSampleAction
    index::Int64
end
RockSampleAction() = RockSampleAction(-1)

struct RockSampleObs
    index::Int64
end
RockSampleObs() = RockSampleObs(-1)

## spaces
struct RockSampleStateSpace
    min_index::Int64
    max_index::Int64
end

struct RockSampleActionSpace
    min_index::Int64
    max_index::Int64
end

struct RockSampleObsSpace
    min_index::Int64
    max_index::Int64
end

## space iterator types
struct RockSampleStateIterator
    min_index::Int64
    max_index::Int64
end

struct RockSampleActionIterator
    min_index::Int64
    max_index::Int64
end

struct RockSampleObsIterator
    min_index::Int64
    max_index::Int64
end

# iterator enabling functions
Base.start(::RockSampleStateIterator)   = RockSampleState(0)
Base.start(::RockSampleActionIterator)  = RockSampleAction(0)
Base.start(::RockSampleObsIterator)     = RockSampleObs(0)

Base.next(::RockSampleStateIterator, state::RockSampleState)        =
    (state, RockSampleState(state.index+1))
Base.next(::RockSampleActionIterator, action::RockSampleAction)     =
    (action, RockSampleAction(action.index+1))
Base.next(::RockSampleObsIterator, obs::RockSampleObs)              =
    (obs, RockSampleObs(obs.index+1))

Base.done(iter::RockSampleStateIterator, state::RockSampleState)    =
    (state.index > iter.max_index)
Base.done(iter::RockSampleActionIterator, action::RockSampleAction) =
    (action.index > iter.max_index)
Base.done(iter::RockSampleObsIterator, obs::RockSampleObs)          =
    (obs.index > iter.max_index)

## iterator creation functions
POMDPs.iterator(space::RockSampleStateSpace)                        =
    RockSampleStateIterator(space.min_index, space.max_index)
POMDPs.iterator(space::RockSampleActionSpace)                       =
    RockSampleActionIterator(space.min_index, space.max_index)
POMDPs.iterator(space::RockSampleObsSpace)                          =
    RockSampleObsIterator(space.min_index, space.max_index)

==(x::RockSampleState, y::RockSampleState)      = (x.index == y.index)
==(x::RockSampleAction, y::RockSampleAction)    = (x.index == y.index)
==(x::RockSampleObs, y::RockSampleObs)          = (x.index == y.index)

Base.hash(x::RockSampleState,  h::Int64 = 0)    = x.index
Base.hash(x::RockSampleAction, h::Int64 = 0)    = x.index
Base.hash(x::RockSampleObs, h::Int64 = 0)       = x.index

######## RockSample type definition ########
mutable struct RockSample <: POMDPs.POMDP{RockSampleState, RockSampleAction, RockSampleObs}
    #problem parameters
    grid_size::Int64
    n_rocks::Int64
    rand_max::Int64
    model_seed ::UInt32 # random seed to construct arbitrary size scenarios
    belief_seed::UInt32 # random seed used in initial belief construction,
                        # if needed

    #problem properties
    n_cells::Int64
    n_actions::Int64
    n_states::Int64
    n_observations::Int64
    robot_start_cell::Int64
    half_eff_distance::Int64
    discount::Float64

    #internal variables and structures
    rock_set_start::Int64
    rock_at_cell::Array{Int64,1}
    cell_to_coords::Array{Vector{Int64},1}
    observation_effectiveness::Array{Float64,2}
    rocks::Array{Int64,1}
    T::Array{RockSampleState,2}
    R::Array{DESPOTReward,2}
    actions::Array{Int64,1} # needed for large problems

    #observation aliases
    #TODO: think how to best convert to const (or just wait for struct fields to be implemented...)
    BAD_OBS::Int64
    GOOD_OBS::Int64
    NONE_OBS::Int64
    TERMINAL_OBS::Int64

    # the const value for 'seed' is meant to provide compatibility with the C++ version of DESPOT
    function RockSample(grid_size::Int64 = 4,
                        n_rocks::Int64 = 4;
                        rand_max::Int64 = 2^31-1,
                        belief_seed::UInt32 = convert(UInt32, 479),
                        model_seed::UInt32  = convert(UInt32, 476),
                        discount::Float64 = 0.95)

          this = new()
          # problem parameters
          this.grid_size = grid_size
          this.n_rocks = n_rocks
          this.rand_max = rand_max
          this.belief_seed = belief_seed
          this.model_seed  = model_seed
          this.discount = discount

          # problem properties
          this.n_cells = grid_size*grid_size
          this.n_actions = n_rocks + 5
          this.n_states = (grid_size*grid_size+1)*(1 << n_rocks)
          this.n_observations = 4
          this.robot_start_cell = 1
          this.half_eff_distance = 20

          #internal variables and structures
          this.rock_set_start = 0
          this.rock_at_cell = Array{Int64}(this.n_cells)
          this.cell_to_coords = Array{Vector{Int64}}(this.n_cells)
          this.observation_effectiveness = Array{Float64}(this.n_cells, this.n_cells)
          this.rocks = Array{Int64}(this.n_rocks)                       # locations
          this.T = Array{RockSampleState}(this.n_states, this.n_actions)
          this.R = Array{Float64}(this.n_states, this.n_actions)
          this.actions = collect(1:n_rocks + 5) # default ordering
          this.BAD_OBS      = 0
          this.GOOD_OBS     = 1
          this.NONE_OBS     = 2
          this.TERMINAL_OBS = 3

          init_problem(this)
          return this
     end
end

POMDPs.state_index(pomdp::RockSample, state::RockSampleState)   = state.index
POMDPs.action_index(pomdp::RockSample, action::RockSampleAction) = action.index
POMDPs.obs_index(pomdp::RockSample, obs::RockSampleObs)       = obs.index

## distribution types
mutable struct RockSampleTransitionDistribution
    pomdp::RockSample
    state::RockSampleState
    action::RockSampleAction
end

mutable struct RockSampleObsDistribution
    pomdp::RockSample
    state::RockSampleState
    action::RockSampleAction
    next_state::RockSampleState
    debug::Int64 #TODO: consider removing

    #TODO: consider removing, not really needed except for debugging
    function RockSampleObsDistribution(pomdp::RockSample,
                                               state::RockSampleState,
                                               action::RockSampleAction,
                                               next_state::RockSampleState,
                                               debug::Int64 = 0)
        this = new()
        this.pomdp = pomdp
        this.state = state
        this.action = action
        this.next_state = next_state
        this.debug = debug

        return this
    end
end

start_state(pomdp::RockSample) =
    RockSampleState(make_state_index(pomdp, pomdp.robot_start_cell, pomdp.rock_set_start))

# Creates a default belief structure to store the problem's initial belief
DESPOT.create_belief(pomdp::RockSample)         =
    ParticleBelief(Array{Particle{RockSampleState}}(0))

create_transition_distribution(pomdp::RockSample)    =
    RockSampleTransitionDistribution(pomdp, RockSampleState(-1), RockSampleAction(-1))
create_observation_distribution(pomdp::RockSample)   =
    RockSampleObsDistribution(pomdp,
                              RockSampleState(-1),
                              RockSampleAction(-1),
                              RockSampleState(-1))

# default initial state distribution represented as a set of generic particles
function POMDPs.initial_state_distribution(pomdp::RockSample)
    n_states = 1 << pomdp.n_rocks
    pool = Array{POMDPToolbox.Particle{RockSampleState}}(0)

    p = 1.0/n_states
    for k = 0:n_states-1 #TODO: can make faster, potentially
        push!(pool, POMDPToolbox.Particle{RockSampleState}(RockSampleState(make_state_index(pomdp, pomdp.robot_start_cell, k)), p))
    end

    return POMDPToolbox.ParticleDistribution(pool)
end

## accessor functions
POMDPs.n_states(pomdp::RockSample)          = pomdp.n_states
POMDPs.n_actions(pomdp::RockSample)         = pomdp.n_actions
POMDPs.n_observations(pomdp::RockSample)    = pomdp.n_observations

# 0-based indexing in the following functions
POMDPs.states(pomdp::RockSample)    = RockSampleStateSpace(0, pomdp.n_states-1)
POMDPs.actions(pomdp::RockSample)   = RockSampleActionSpace(0, pomdp.n_actions-1)
POMDPs.observations(pomdp::RockSample) = RockSampleActionSpace(0, pomdp.n_observations-1)

POMDPs.discount(pomdp::RockSample) = pomdp.discount

## RockSample initialization
function init_4_4(pomdp::RockSample)
  pomdp.rocks[1] = cell_num(pomdp,0,2) # rocks is an array
  pomdp.rocks[2] = cell_num(pomdp,2,2)
  pomdp.rocks[3] = cell_num(pomdp,3,2)
  pomdp.rocks[4] = cell_num(pomdp,3,3)
  pomdp.robot_start_cell = cell_num(pomdp,2,0)
end

function init_7_8(pomdp::RockSample)
  pomdp.rocks[1] = cell_num(pomdp,0,1)
  pomdp.rocks[2] = cell_num(pomdp,1,5)
  pomdp.rocks[3] = cell_num(pomdp,2,2)
  pomdp.rocks[4] = cell_num(pomdp,2,3)
  pomdp.rocks[5] = cell_num(pomdp,3,6)
  pomdp.rocks[6] = cell_num(pomdp,5,0)
  pomdp.rocks[7] = cell_num(pomdp,5,3)
  pomdp.rocks[8] = cell_num(pomdp,6,2)
  pomdp.robot_start_cell = cell_num(pomdp,3,0)
end

function init_11_11(pomdp::RockSample)
  pomdp.rocks[1] = cell_num(pomdp,7,0)
  pomdp.rocks[2] = cell_num(pomdp,3,0)
  pomdp.rocks[3] = cell_num(pomdp,2,1)
  pomdp.rocks[4] = cell_num(pomdp,6,2)
  pomdp.rocks[5] = cell_num(pomdp,7,3)
  pomdp.rocks[6] = cell_num(pomdp,2,3)
  pomdp.rocks[7] = cell_num(pomdp,7,4)
  pomdp.rocks[8] = cell_num(pomdp,2,5)
  pomdp.rocks[9] = cell_num(pomdp,9,6)
  pomdp.rocks[10] = cell_num(pomdp,7,9)
  pomdp.rocks[11] = cell_num(pomdp,1,9)
  pomdp.robot_start_cell = cell_num(pomdp,5,0)
end

function init_15_15(pomdp::RockSample)
  pomdp.rocks[1] = cell_num(pomdp,7,0)
  pomdp.rocks[2] = cell_num(pomdp,3,0)
  pomdp.rocks[3] = cell_num(pomdp,2,1)
  pomdp.rocks[4] = cell_num(pomdp,6,2)
  pomdp.rocks[5] = cell_num(pomdp,7,3)
  pomdp.rocks[6] = cell_num(pomdp,2,3)
  pomdp.rocks[7] = cell_num(pomdp,7,4)
  pomdp.rocks[8] = cell_num(pomdp,2,5)
  pomdp.rocks[9] = cell_num(pomdp,9,6)
  pomdp.rocks[10] = cell_num(pomdp,7,9)
  pomdp.rocks[11] = cell_num(pomdp,1,9)
  pomdp.rocks[12] = cell_num(pomdp,8,11)
  pomdp.rocks[13] = cell_num(pomdp,10,13)
  pomdp.rocks[14] = cell_num(pomdp,9,14)
  pomdp.rocks[15] = cell_num(pomdp,2,12)
  pomdp.robot_start_cell = cell_num(pomdp,5,0)
end

function init_general(pomdp::RockSample, seed::Array{UInt32,1})

    rockIndex::Int64 = 1 # rocks is an array
    if !is_linux()
        srand(seed[1])
    end

    while rockIndex <= pomdp.n_rocks
        if is_linux()
            cell = ccall((:rand_r, "libc"), Int, (Ptr{Cuint},), seed) % pomdp.n_cells
        else
            cell = Base.rand(0:pomdp.rand_max) % pomdp.n_cells
        end

        if findfirst(pomdp.rocks, cell) == 0
            pomdp.rocks[rockIndex] = cell
            rockIndex += 1
        end
    end
    pomdp.robot_start_cell = cell_num(pomdp, round(Integer, pomdp.grid_size/2), 0)
end

function init_problem(pomdp::RockSample)

    pomdp.rocks = Array{Int64}(pomdp.n_rocks)
    seed = Cuint[convert(UInt32, pomdp.model_seed)]

    if pomdp.grid_size == 4 && pomdp.n_rocks == 4
        init_4_4(pomdp)
    elseif pomdp.grid_size == 7 && pomdp.n_rocks == 8
        init_7_8(pomdp)
    elseif pomdp.grid_size == 11 && pomdp.n_rocks == 11
        init_11_11(pomdp)
    elseif pomdp.grid_size == 15 && pomdp.n_rocks == 15
        init_15_15(pomdp)
    else
        init_general(pomdp, seed)
    end

    # Compute rock set start
    pomdp.rock_set_start = 0

    if !is_linux()
        srand(seed[1])
    end

    for i in 0 : pomdp.n_rocks-1
        if is_linux()
            rand_num = ccall((:rand_r, "libc"), Int, (Ptr{Cuint},), seed)
        else #Windows, etc
            rand_num = Base.rand(0:pomdp.rand_max)
        end

        if (rand_num & 1) == 1
            pomdp.rock_set_start |= (1 << i)
        end
    end

    # Fill in cellToCoord and init rock_at_cell mappings
    fill!(pomdp.rock_at_cell, -1)

    for i in 0 : pomdp.n_cells-1
        pomdp.cell_to_coords[i+1] = [trunc(Integer, i/pomdp.grid_size), i % pomdp.grid_size]
    end

    for i in 0 : pomdp.n_rocks-1
        pomdp.rock_at_cell[pomdp.rocks[i+1]+1] = i # rock_at_cell and rocks are arrays
    end

    # T and R - ALL INDICES BELOW ARE OFFSET BY +1 (for 1-based array indexing)
    for cell in 0 : pomdp.n_cells-1
        for rock_set = 0:(1 << pomdp.n_rocks)-1
        s_index = make_state_index(pomdp, cell, rock_set)

            #initialize transition and rewards with default values
            for a_index in 0:pomdp.n_actions-1
                pomdp.T[s_index+1,a_index+1] =
                    RockSampleState(s_index)
                pomdp.R[s_index+1,a_index+1] = 0.
            end

            row, col = pomdp.cell_to_coords[cell+1]
            # North
            if row == 0
                pomdp.T[s_index+1,1] =
                    RockSampleState(s_index)
                pomdp.R[s_index+1,1] = -100.
            else
                pomdp.T[s_index+1,1] =
                    RockSampleState(make_state_index(pomdp, cell_num(pomdp,row-1,col), rock_set))
                pomdp.R[s_index+1,1] = 0
            end

            # South
            if row == pomdp.grid_size-1
                pomdp.T[s_index+1,2] =
                    RockSampleState(s_index)
                pomdp.R[s_index+1,2] = -100.
            else
                pomdp.T[s_index+1,2] =
                    RockSampleState(make_state_index(pomdp, cell_num(pomdp,row+1,col), rock_set))
                pomdp.R[s_index+1,2] = 0.
            end

            # East
            if col == pomdp.grid_size-1
                pomdp.T[s_index+1,3] =
                    RockSampleState(make_state_index(pomdp, pomdp.n_cells, rock_set))
                pomdp.R[s_index+1,3] = 10.
            else
                pomdp.T[s_index+1,3] =
                    RockSampleState(make_state_index(pomdp, cell_num(pomdp,row,col+1), rock_set))
                pomdp.R[s_index+1,3] = 0.
            end

            # West
            if col == 0
                pomdp.T[s_index+1, 4] = RockSampleState(s_index)
                pomdp.R[s_index+1, 4] = -100.
            else
                pomdp.T[s_index+1, 4] =
                    RockSampleState(make_state_index(pomdp, cell_num(pomdp,row,col-1), rock_set))
                pomdp.R[s_index+1, 4] = 0.
            end

            # Sample
            rock = pomdp.rock_at_cell[cell+1] # array
            if rock != -1
                if rock_status(rock, rock_set)
                    pomdp.T[s_index+1, 5] =
                        RockSampleState(make_state_index(pomdp, cell, sample_rock_set(rock, rock_set)))
                    pomdp.R[s_index+1, 5] = +10.
                else
                    pomdp.T[s_index+1, 5] = RockSampleState(s_index)
                    pomdp.R[s_index+1, 5] = -10.
                end
            else
                pomdp.T[s_index+1, 5] = RockSampleState(s_index)
                pomdp.R[s_index+1, 5] = -100.
            end

            # Check
            for a_index in 5:pomdp.n_actions-1
                pomdp.T[s_index+1, a_index+1] = RockSampleState(s_index)
                pomdp.R[s_index+1, a_index+1] = 0.
            end
        end
    end

    # Terminal states
    for k = 0:(1 << pomdp.n_rocks)-1
        s_index = make_state_index(pomdp, pomdp.n_cells, k);
        for a_index in 0:pomdp.n_actions-1
            pomdp.T[s_index+1, a_index+1] = RockSampleState(s_index)
            pomdp.R[s_index+1, a_index+1] = 0.
        end
    end

    # precompute observation effectiveness table
    for i in 0 : pomdp.n_cells-1
        for j in 0 : pomdp.n_cells-1
        agent = pomdp.cell_to_coords[i+1]
        other = pomdp.cell_to_coords[j+1]
        dist = sqrt((agent[1] - other[1])^2 + (agent[2]-other[2])^2)
        pomdp.observation_effectiveness[i+1,j+1] =
            (1 + 2^(-dist / pomdp.half_eff_distance)) * 0.5 # Array indexing starts from 1.
                                                    # Remember to subtract one to go back
        end
    end
end

## utility functions

# True for good rock, false for bad rock, x can be a rock set or state index
rock_status(rock::Int64, x::Int64) = (((x >>> rock) & 1) == 1 ? true : false)

cell_num(pomdp::RockSample, row::Int64, col::Int64) = row * pomdp.grid_size + col

make_state_index(pomdp::RockSample, cell::Int64, rock_set::Int64) =
    convert(Int64, (cell << pomdp.n_rocks) + rock_set)

POMDPs.reward(pomdp::RockSample, s::RockSampleState, a::RockSampleAction) =
    pomdp.R[s.index+1, a.index+1]

function POMDPs.transition(
                    pomdp::RockSample,
                    state::RockSampleState,
                    action::RockSampleAction)

    distribution = create_transition_distribution(pomdp)
    distribution.pomdp = pomdp
    distribution.state = state
    distribution.action = action

    return distribution
end

function POMDPs.observation(
                    pomdp::RockSample,
                    state::RockSampleState,
                    action::RockSampleAction,
                    next_state::RockSampleState)

    distribution = create_observation_distribution(pomdp)
    distribution.pomdp = pomdp
    distribution.state = next_state
    distribution.action = action
    distribution.next_state = next_state

    return distribution
end

function POMDPs.rand(
                    rng::AbstractRNG,
                    distribution::RockSampleTransitionDistribution,
                    sample=nothing)

    return RockSampleState(
        distribution.pomdp.T[distribution.state.index+1, distribution.action.index+1].index)
end

function POMDPs.rand(
                    rng::AbstractRNG,
                    distribution::RockSampleObsDistribution,
                    sample=nothing)

    # generate a new random number regardless of whether it's used below or not
    rand_num::Array{Float64} = Array{Float64}(1)
    rand!(rng, rand_num)

    if (distribution.action.index < 5)
        obs = isterminal(distribution.pomdp, distribution.next_state) ?
                    RockSampleObs(distribution.pomdp.TERMINAL_OBS) :
                    RockSampleObs(distribution.pomdp.NONE_OBS) # rs.T is an array
    else
        rock_cell = distribution.pomdp.rocks[distribution.action.index - 4] # would be [action-5] with 0-based indexing
        agent_cell = cell_of(distribution.pomdp, distribution.state)
        eff = distribution.pomdp.observation_effectiveness[agent_cell+1, rock_cell+1]

        if (rand_num[1] <= eff) == rock_status(distribution.action.index - 5, distribution.state.index)
            obs = RockSampleObs(distribution.pomdp.GOOD_OBS)
        else
            obs = RockSampleObs(distribution.pomdp.BAD_OBS)
        end
    end

    return obs
end

function POMDPs.rand(
                    rng::AbstractRNG,
                    state_space::RockSampleStateSpace,
                    sample=nothing)

    if is_linux()
        random_number = ccall((:rand_r, "libc"), Int, (Ptr{Cuint},), rng.seed) / rng.rand_max
    else #Windows, etc
        srand(seed)
        random_number = Base.rand()
    end

    return RockSampleState(floor(
           random_number*(state_space.max_index-state_space.min_index)))
end

function POMDPs.pdf(distribution::RockSampleObsDistribution,
                    obs::RockSampleObs)
  # Terminal state should match terminal obs
  if isterminal(distribution.pomdp, distribution.next_state)
      if obs.index == distribution.pomdp.TERMINAL_OBS
          return 1.
      else
          return 0.
      end
  end

  if (distribution.action.index < 5)
      if obs.index == distribution.pomdp.NONE_OBS
          return 1.
      else
          return 0.
      end
  end

  if ((obs.index != distribution.pomdp.GOOD_OBS) && (obs.index != distribution.pomdp.BAD_OBS))
    return 0.
  end

  rock::Int64       = distribution.action.index - 5
  rockCell::Int64   = distribution.pomdp.rocks[rock+1]
  agentCell::Int64  = cell_of(distribution.pomdp, distribution.next_state)

  eff::Float64 = distribution.pomdp.observation_effectiveness[agentCell+1, rockCell+1]

  rstatus::Bool = rock_status(rock, distribution.next_state.index)
  if ((obs.index == distribution.pomdp.GOOD_OBS) && (rstatus == true)) ||
     ((obs.index == distribution.pomdp.BAD_OBS) && (rstatus == false))
    return eff
  else
    return 1. - eff
  end
end

POMDPs.isterminal(pomdp::RockSample, s::RockSampleState)        =
    cell_of(pomdp, s) == pomdp.n_cells

POMDPs.isterminal_obs(pomdp::RockSample, obs::RockSampleObs)    =
    obs.index == pomdp.TERMINAL_OBS

# Which cell the agent is in
cell_of(pomdp::RockSample, s::RockSampleState) = (s.index >>> pomdp.n_rocks)

# The rock set after sampling a rock from it
sample_rock_set(rock::Int64, rock_set::Int64) = (rock_set & ~(1 << rock))

# The set of rocks in the state
rock_set_of(pomdp::RockSample, s::RockSampleState) =
    s.index & ((1 << pomdp.n_rocks)-1)

function show_state(pomdp::RockSample, s::RockSampleState)
  ac = cell_of(pomdp, s)
  for i in 0:pomdp.grid_size-1
    for j in 0:pomdp.grid_size-1
      if ac == cell_num(pomdp,i,j)
        if pomdp.rock_at_cell[ac+1] == -1 # array
          print("R ")
        elseif rock_status(pomdp.rock_at_cell[ac+1], rock_set_of(pomdp,s))
          print("G ")
        else
          print("B ")
        end # if rock_at_cell[ac] == -1
        continue
      end # if ac == cell_num(i, j)
      if pomdp.rock_at_cell[cell_num(pomdp,i,j)+1] == -1
        print(". ")
      elseif (rock_status(pomdp.rock_at_cell[cell_num(pomdp,i,j)+1], rock_set_of(pomdp,s)))
        print("1 ")
      else
        print("0 ")
      end # if rock_at_cell[cell_num(i,j)] == -1

    end # for j in 1:grid_size
    println("")
  end # i in 1:grid_size
end

function show_obs(pomdp::RockSample, obs::RockSampleObs)
    if obs.index == pomdp.NONE_OBS
        println("NONE")
    elseif obs.index == pomdp.GOOD_OBS
        println("GOOD")
    elseif obs.index == pomdp.BAD_OBS
        println("BAD")
    elseif obs.index == pomdp.TERMINAL_OBS
        println("TERMINAL")
    else
        println("UNKNOWN")
    end
end


text:
import Pkg

function maketempdir()::String
    dir::String = mktempdir()
    atexit(() -> rm(dir; force = true, recursive = true,))
    return dir
end

mutable struct DummyOutputWrapperStruct{I, F, S, O}
    previous_time_seconds::I
    f::F
    interval_seconds::I
    dummy_output::S
    io::O
end

function DummyOutputWrapperStruct(
        ;
        interval_seconds::I = 60,
        initial_offset_seconds::I = interval_seconds,
        f::F,
        dummy_output::S = "This is a dummy line of output",
        io::O = Base.stdout,
        )::DummyOutputWrapperStruct{I, F, S, O} where
            I <: Integer where
            F <: Function where
            S <: AbstractString where
            O <: IO
    current_time_seconds::I = floor(I, time())
    initial_time_seconds::I = current_time_seconds + initial_offset_seconds
    wrapper_struct::DummyOutputWrapperStruct{I, F, S} =
        DummyOutputWrapperStruct(
            initial_time_seconds,
            f,
            interval_seconds,
            dummy_output,
            io,
            )
    return wrapper_struct
end

function (x::DummyOutputWrapperStruct{I, F, S, O})() where
        I <: Integer where
        F <: Function where
        S <: AbstractString where
        O <: IO
    current_time_seconds::I = floor(I, time())
    previous_time_seconds::I = x.previous_time_seconds
    f::F = x.f
    interval_seconds::I = x.interval_seconds
    dummy_output::S = x.dummy_output
    io::O = x.io
    elapsed_seconds::Int = current_time_seconds - previous_time_seconds
    print_dummy_output::Bool = elapsed_seconds > interval_seconds
    if print_dummy_output
        println(io, dummy_output)
        x.previous_time_seconds = current_time_seconds
    end
    f_result = f()
    return f_result
end

function dummy_output_wrapper(
        ;
        f::F,
        interval_seconds::I = 60,
        initial_offset_seconds::I = interval_seconds,
        dummy_output::S = "This is a dummy line of output",
        io::O = Base.stdout,
        ) where
            I <: Integer where
            F <: Function where
            S <: AbstractString where
            O <: IO
    wrapper_struct::DummyOutputWrapperStruct{I, F, S, O} =
        DummyOutputWrapperStruct(
            ;
            f = f,
            interval_seconds = interval_seconds,
            initial_offset_seconds = initial_offset_seconds,
            dummy_output = dummy_output,
            )
    function my_wrapper_function()
        result = wrapper_struct()
        return result
    end
    return my_wrapper_function
end

function command_ran_successfully(
        cmd::Base.AbstractCmd;
        max_attempts::Integer = 10,
        max_seconds_per_attempt::Real = 540,
        seconds_to_wait_between_attempts::Real = 30,
        error_on_failure::Bool = true,
        last_resort_run::Bool = true,
        before::Function = () -> (),
        )::Bool
    success_bool::Bool = false

    my_false = dummy_output_wrapper(
        ;
        f = () -> false,
        interval_seconds = 60,
        initial_offset_seconds = 60,
        dummy_output = "Still waiting between attempts...",
        io = Base.stdout,
        )

    for attempt = 1:max_attempts
        if success_bool
        else
            @debug(string("Attempt $(attempt)"))
            if attempt > 1
                timedwait(
                    () -> my_false(),
                    float(seconds_to_wait_between_attempts);
                    pollint = float(1.0),
                    )
            end
            before()
            p = run(cmd; wait = false,)
            my_process_exited = dummy_output_wrapper(
                ;
                f = () -> process_exited(p),
                interval_seconds = 60,
                initial_offset_seconds = 60,
                dummy_output = "The process is still running...",
                io = Base.stdout,
                )
            timedwait(
                () -> my_process_exited(),
                float(max_seconds_per_attempt),
                pollint = float(1.0),
                )
            if process_running(p)
                success_bool = false
                try
                    kill(p, Base.SIGTERM)
                catch exception
                    @warn("Ignoring exception: ", exception)
                end
                try
                    kill(p, Base.SIGKILL)
                catch exception
                    @warn("Ignoring exception: ", exception)
                end
            else
                success_bool = try
                    success(p)
                catch exception
                    @warn("Ignoring exception: ", exception)
                    false
                end
            end
        end
    end
    if !success_bool && last_resort_run
        @debug(string("Attempting the last resort run..."))
        run(cmd)
        success_bool = true
    end
    if success_bool
        @debug(string("Command ran successfully."),)
    else
        if error_on_failure
            error(string("Command did not run successfully."),)
        else
            @warn(string("Command did not run successfully."),)
        end
    end
    return success_bool
end

function retry_function_until_success(
        f::Function;
        max_attempts::Integer = 10,
        seconds_to_wait_between_attempts::Real = 30,
        )
    success_bool::Bool = false
    f_result = nothing

    my_false = dummy_output_wrapper(
        ;
        f = () -> false,
        interval_seconds = 60,
        initial_offset_seconds = 60,
        dummy_output = "Still waiting between attempts...",
        io = Base.stdout,
        )

    for attempt = 1:max_attempts
        if success_bool
        else
            @debug(string("Attempt $(attempt)"))
            if attempt > 1
                timedwait(
                    () -> my_false(),
                    float(seconds_to_wait_between_attempts);
                    pollint = float(1.0),
                    )
            end
            @debug(string("Running the provided function..."))
            success_bool = true
            f_result = try
                f()
            catch exception
                success_bool = false
                @warn("Ignoring exception: ", exception)
                nothing
            end
        end
    end

    if success_bool
        @debug(string("Function ran successfully."),)
        return f_result
    else
        error(string("Function did not run successfully."),)
    end
end

function git_clone_registry(url::AbstractString)::String
    tmp = maketempdir()
    Base.shred!(LibGit2.CachedCredentials()) do creds
        LibGit2.with(
            Pkg.GitTools.clone(
                url,
                tmp;
                header = "registry from $(repr(url))",
                credentials = creds,
                )
            ) do repo
        end
    end
    return tmp
end

function git_clone_repo(url::AbstractString)::String
    tmp = maketempdir()
    Base.shred!(LibGit2.CachedCredentials()) do creds
        LibGit2.with(
            Pkg.GitTools.clone(
                url,
                tmp;
                header = "git-repo from $(repr(url))",
                credentials = creds,
                )
            ) do repo
        end
    end
    return tmp
end

function Pkg_add_name_ignore_julia_version_error(
        name::AbstractString,
        )::Nothing
    try
        Pkg.add(name)
    catch e
        if occursin("Unsatisfiable requirements detected", repr(a)) &&
                occursin("restricted by julia compatibility", repr(a))
            @debug("ignoring error: ", e,)
        else
            @warn("rethrowing error: ", e,)
        end
    end
    return nothing
end

ext: jl

max_forks_repo_name: aluthge/OfflineRegistry


text:
module STMOZOO

# execute your source file and export the module you made
include("example.jl")
export Example

end # module


text:
module AndExport

##############################
# * Exporting
#----------------------------

export @xport

using MacroTools

# Convenient iteration flatten
# My version of flatten1, which doesn't remove blocks at the top level (this ruins structs)
function my_flatten1(ex)
  isexpr(ex, :block) || return ex
  #ex′ = :(;)
  ex′ = Expr(:block)
  for x in ex.args
    isexpr(x, :block) ? append!(ex′.args, x.args) : push!(ex′.args, x)
  end
  return ex′
end
iterflatten(ex) = MacroTools.postwalk(my_flatten1, block(ex)).args

macro xport(expr::Expr)
    expr = macroexpand(__module__, expr)

    exports = Expr[]

    for item in iterflatten(expr)
        if @capture (item) ( (function name_(args__) where {T__} content__ end) |
                                   (function name_(args__) content__ end) |
                                   (function name_ end) |
                                   (name_(args__) where {T__} = content__) |
                                   (name_(args__) = content__) |
                                   (mutable struct name_ content__ end) |
                                   (struct name_ content__ end) |
                                   (abstract type name_ end) |
                                   (const name_ = val__) |
                                   (name_ = val__) )
            if @capture (name) ( (name2_{T__} <: Z_ where T2_) |
                                 (name2_ <: Z_ where T2_) |
                                 (name2_{T__} where T2_) |
                                 (name2_ where T2_) |
                                 (name2_{T__} <: Z_) |
                                 (name2_ <: Z_) |
                                 (name2_{T__}) |
                                 name2_ )
                # Ignore (::Type{T}) functions for exports
                # if @capture (name2) (::Type{T__})
                #     @info "Ignoring" name
                if name2 isa Expr# && name2.head == Symbol("::")
                    # @info "Ignoring" name
                else
                    push!(exports, :(export $(esc(name2))))
                end
            end
        elseif @capture (item) (macro name_(args__) content__ end)
            name = Symbol("@" * string(name))
            push!(exports, :(export $(esc(name))))
        end
        item
    end

    if isempty(exports)
        buf = IOBuffer()
        dump(buf, expr, maxdepth=1)
        repr = String(take!(buf))
        error("Unknown object to export: $repr.")
    end

    expr = :(Base.@__doc__ $(esc(expr)))
    fullcode = Expr(:block, exports..., expr)
end

end # module


text:
using LightGraphs
using MetaGraphs
using LinearAlgebra
using SparseArrays
using SuiteSparse
using Combinatorics

include("./misc.jl")

function read_sparsity_pattern(instance_path::String)
    data = load_matpower(instance_path)
    ## Bus load and shunt information
    i_debut, i_fin = find_numarray(1, data)
    checkfor(data, i_debut-1, "mpc.bus")
    nb_bus = i_fin-i_debut+1
    index_bus = Dict( data[i+i_debut-1,1] => i for i=1:nb_bus)
    ## Bus generator information
    i_debut, i_fin = find_numarray(i_fin+1, data)
    checkfor(data, i_debut-1, "mpc.gen")
    #initialize network graph G
    sp = spzeros(nb_bus,nb_bus)
    ## Link information
    i_debut, i_fin = find_numarray(i_fin+1, data)
    checkfor(data, i_debut-1, "mpc.branch")
    for i=i_debut:i_fin
        if data[i, 11] == 0
        else
            orig = index_bus[data[i,1]]
            dest = index_bus[data[i,2]]
            sp[orig,dest] = 1
        end
    end
    sp_sym = sp + sp'
    diag = zeros(nb_bus)
    for i=1:nb_bus
        sum_col = sum(sp_sym[i,:])
        diag[i] = Int(sum_col + 1)
    end
    return sp_sym + Diagonal(diag)
end

function read_sparsity_pattern(graph::AbstractGraph)
    sparsity_pattern = adjacency_matrix(graph)
    for i in 1:nv(graph)
        sparsity_pattern[i, i] = length(neighbors(graph, i)) + 1
    end
    return sparsity_pattern
end

max_stars_repo_name: charlyalizadeh/OPFdecgen

max_issues_repo_name: charlyalizadeh/OPFdecgen

max_forks_repo_name: charlyalizadeh/OPFdecgen


text:
module GnssWrapper

using JavaCall

const AbsoluteDate = @jimport org.orekit.time.AbsoluteDate
const AbstractCycleSlipDetector = @jimport org.orekit.estimation.measurements.gnss.AbstractCycleSlipDetector
const AbstractDualFrequencyCombination = @jimport org.orekit.estimation.measurements.gnss.AbstractDualFrequencyCombination
const AbstractLambdaMethod = @jimport org.orekit.estimation.measurements.gnss.AbstractLambdaMethod
const AbstractMeasurement = @jimport org.orekit.estimation.measurements.AbstractMeasurement
const AbstractMeasurementBuilder = @jimport org.orekit.estimation.measurements.generation.AbstractMeasurementBuilder
const AbstractSingleFrequencyCombination = @jimport org.orekit.estimation.measurements.gnss.AbstractSingleFrequencyCombination
const AlternatingSampler = @jimport org.orekit.estimation.measurements.gnss.AlternatingSampler
const AmbiguityAcceptance = @jimport org.orekit.estimation.measurements.gnss.AmbiguityAcceptance
const AmbiguitySolver = @jimport org.orekit.estimation.measurements.gnss.AmbiguitySolver
const Class = @jimport java.lang.Class
const CombinationType = @jimport org.orekit.estimation.measurements.gnss.CombinationType
const CombinedObservationData = @jimport org.orekit.gnss.CombinedObservationData
const CombinedObservationDataSet = @jimport org.orekit.gnss.CombinedObservationDataSet
const ComparableMeasurement = @jimport org.orekit.estimation.measurements.ComparableMeasurement
const CorrelatedRandomVectorGenerator = @jimport org.hipparchus.random.CorrelatedRandomVectorGenerator
const CycleSlipDetectorResults = @jimport org.orekit.estimation.measurements.gnss.CycleSlipDetectorResults
const CycleSlipDetectors = @jimport org.orekit.estimation.measurements.gnss.CycleSlipDetectors
const Enum = @jimport java.lang.Enum
const EstimatedMeasurement = @jimport org.orekit.estimation.measurements.EstimatedMeasurement
const EstimationModifier = @jimport org.orekit.estimation.measurements.EstimationModifier
const FieldAbsoluteDate = @jimport org.orekit.time.FieldAbsoluteDate
const FieldVector3D = @jimport org.hipparchus.geometry.euclidean.threed.FieldVector3D
const Frequency = @jimport org.orekit.gnss.Frequency
const GRAPHICCombination = @jimport org.orekit.estimation.measurements.gnss.GRAPHICCombination
const GeometryFreeCombination = @jimport org.orekit.estimation.measurements.gnss.GeometryFreeCombination
const GeometryFreeCycleSlipDetector = @jimport org.orekit.estimation.measurements.gnss.GeometryFreeCycleSlipDetector
const GroundStation = @jimport org.orekit.estimation.measurements.GroundStation
const IntegerBootstrapping = @jimport org.orekit.estimation.measurements.gnss.IntegerBootstrapping
const IntegerLeastSquareSolution = @jimport org.orekit.estimation.measurements.gnss.IntegerLeastSquareSolution
const IntegerLeastSquareSolver = @jimport org.orekit.estimation.measurements.gnss.IntegerLeastSquareSolver
const IonosphereFreeCombination = @jimport org.orekit.estimation.measurements.gnss.IonosphereFreeCombination
const JString = @jimport java.lang.String
const LambdaMethod = @jimport org.orekit.estimation.measurements.gnss.LambdaMethod
const List = @jimport java.util.List
const Map = @jimport java.util.Map
const MeasurementCombination = @jimport org.orekit.estimation.measurements.gnss.MeasurementCombination
const MeasurementCombinationFactory = @jimport org.orekit.estimation.measurements.gnss.MeasurementCombinationFactory
const MelbourneWubbenaCombination = @jimport org.orekit.estimation.measurements.gnss.MelbourneWubbenaCombination
const ModifiedLambdaMethod = @jimport org.orekit.estimation.measurements.gnss.ModifiedLambdaMethod
const NarrowLaneCombination = @jimport org.orekit.estimation.measurements.gnss.NarrowLaneCombination
const Object = @jimport java.lang.Object
const ObservableSatellite = @jimport org.orekit.estimation.measurements.ObservableSatellite
const ObservationData = @jimport org.orekit.gnss.ObservationData
const ObservationDataSet = @jimport org.orekit.gnss.ObservationDataSet
const Optional = @jimport java.util.Optional
const ParameterDriver = @jimport org.orekit.utils.ParameterDriver
const Phase = @jimport org.orekit.estimation.measurements.gnss.Phase
const PhaseBuilder = @jimport org.orekit.estimation.measurements.gnss.PhaseBuilder
const PhaseMinusCodeCombination = @jimport org.orekit.estimation.measurements.gnss.PhaseMinusCodeCombination
const PhaseMinusCodeCycleSlipDetector = @jimport org.orekit.estimation.measurements.gnss.PhaseMinusCodeCycleSlipDetector
const RealFieldElement = @jimport org.hipparchus.RealFieldElement
const RealMatrix = @jimport org.hipparchus.linear.RealMatrix
const SatelliteSystem = @jimport org.orekit.gnss.SatelliteSystem
const SimpleRatioAmbiguityAcceptance = @jimport org.orekit.estimation.measurements.gnss.SimpleRatioAmbiguityAcceptance
const SpacecraftState = @jimport org.orekit.propagation.SpacecraftState
const TimeStampedFieldPVCoordinates = @jimport org.orekit.utils.TimeStampedFieldPVCoordinates
const TimeStampedPVCoordinates = @jimport org.orekit.utils.TimeStampedPVCoordinates
const Vector3D = @jimport org.hipparchus.geometry.euclidean.threed.Vector3D
const WideLaneCombination = @jimport org.orekit.estimation.measurements.gnss.WideLaneCombination
const WindUp = @jimport org.orekit.estimation.measurements.gnss.WindUp
const WindUpFactory = @jimport org.orekit.estimation.measurements.gnss.WindUpFactory

export AbstractCycleSlipDetector
export AbstractDualFrequencyCombination
export AbstractLambdaMethod
export AbstractSingleFrequencyCombination
export AlternatingSampler
export AmbiguityAcceptance
export AmbiguitySolver
export CombinationType
export CycleSlipDetectorResults
export CycleSlipDetectors
export GRAPHICCombination
export GeometryFreeCombination
export GeometryFreeCycleSlipDetector
export IntegerBootstrapping
export IntegerLeastSquareSolution
export IntegerLeastSquareSolver
export IonosphereFreeCombination
export LambdaMethod
export MeasurementCombination
export MeasurementCombinationFactory
export MelbourneWubbenaCombination
export ModifiedLambdaMethod
export NarrowLaneCombination
export Phase
export PhaseBuilder
export PhaseMinusCodeCombination
export PhaseMinusCodeCycleSlipDetector
export SimpleRatioAmbiguityAcceptance
export WideLaneCombination
export WindUp
export WindUpFactory
export accept
export add_modifier
export build
export combine
export compare_to
export describe_constable
export detect
export equals
export estimate
export fix_integer_ambiguities
export generate_next
export get_all_ambiguity_drivers
export get_base_weight
export get_begin_date
export get_class
export get_coordinates
export get_current
export get_cycle_slip_map
export get_date
export get_declaring_class
export get_dimension
export get_end_date
export get_geometry_free_combination
export get_graphic_combination
export get_ionosphere_free_combination
export get_melbourne_wubbena_combination
export get_mid_point
export get_modifiers
export get_name
export get_narrow_lane_combination
export get_observed_value
export get_parameters_drivers
export get_phase_minus_code_combination
export get_satellite_name
export get_satellites
export get_solution
export get_squared_distance
export get_station
export get_theoretical_standard_deviation
export get_wavelength
export get_wide_lane_combination
export get_wind_up
export hash_code
export in_range
export init
export is_enabled
export modify
export name
export notify
export notify_all
export number_of_candidates
export ordinal
export set_enabled
export set_radius
export signal_time_of_flight
export solve_ils
export to_string
export un_fix_ambiguity
export value_of
export values
export wait

include("abstract_cycle_slip_detector.jl")
include("abstract_dual_frequency_combination.jl")
include("abstract_lambda_method.jl")
include("abstract_single_frequency_combination.jl")
include("alternating_sampler.jl")
include("ambiguity_acceptance.jl")
include("ambiguity_solver.jl")
include("combination_type.jl")
include("cycle_slip_detector_results.jl")
include("cycle_slip_detectors.jl")
include("geometry_free_combination.jl")
include("geometry_free_cycle_slip_detector.jl")
include("graphic_combination.jl")
include("integer_bootstrapping.jl")
include("integer_least_square_solution.jl")
include("integer_least_square_solver.jl")
include("ionosphere_free_combination.jl")
include("lambda_method.jl")
include("measurement_combination.jl")
include("measurement_combination_factory.jl")
include("melbourne_wubbena_combination.jl")
include("modified_lambda_method.jl")
include("narrow_lane_combination.jl")
include("phase.jl")
include("phase_builder.jl")
include("phase_minus_code_combination.jl")
include("phase_minus_code_cycle_slip_detector.jl")
include("simple_ratio_ambiguity_acceptance.jl")
include("wide_lane_combination.jl")
include("wind_up.jl")
include("wind_up_factory.jl")

end


text:
module TestDoctest

import DataTools
using Documenter: doctest
using Test

@testset "doctest" begin
    doctest(DataTools; manual = false)
end

end  # module


text:
using Documenter
using OrbitalTrajectories

doctest(OrbitalTrajectories)


text:
"""
```
estimate(m, data; verbose=:low, proposal_covariance=Matrix(), method=:SMC)
```

Estimate the DSGE parameter posterior distribution.

### Arguments:
- `m::AbstractDSGEModel`: model object

### Optional Arguments:
- `data`: well-formed data as `Matrix` or `DataFrame`. If this is not provided, the `load_data` routine will be executed.

### Keyword Arguments:
- `verbose::Symbol`: The desired frequency of function progress messages printed to standard out.
   - `:none`: No status updates will be reported.
   - `:low`: Status updates will be provided in csminwel and at each block in
     Metropolis-Hastings.
   - `:high`: Status updates provided at each iteration in Metropolis-Hastings.
- `proposal_covariance::Matrix`: Used to test the metropolis_hastings algorithm with a precomputed
  covariance matrix for the proposal distribution. When the Hessian is singular,
  eigenvectors corresponding to zero eigenvectors are not well defined, so eigenvalue
  decomposition can cause problems. Passing a precomputed matrix allows us to ensure that
  the rest of the routine has not broken.
- `method::Symbol`: The method to use when sampling from the posterior distribution. Can
    be either `:MH` for standard Metropolis Hastings Markov Chain Monte Carlo, or `:SMC`
    for Sequential Monte Carlo.
- `mle`: Set to true if parameters should be estimated by maximum likelihood directly.
    If this is set to true, this function will return after estimating parameters.
- `sampling`: Set to false to disable sampling from the posterior.
"""
function estimate(m::AbstractDSGEModel, df::DataFrame;
                  verbose::Symbol = :low,
                  proposal_covariance::Matrix = Matrix(undef, 0, 0),
                  mle::Bool = false,
                  sampling::Bool = true,
                  filestring_addl::Vector{String} = Vector{String}(),
                  continue_intermediate::Bool = false,
                  intermediate_stage_start::Int = 0,
                  intermediate_stage_increment::Int = 10,
                  save_intermediate::Bool = false)
    data = df_to_matrix(m, df)
    estimate(m, data; verbose = verbose, proposal_covariance = proposal_covariance,
             mle = mle, sampling = sampling,
             intermediate_stage_increment = intermediate_stage_increment,
             save_intermediate = save_intermediate)
end

function estimate(m::AbstractDSGEModel;
                  verbose::Symbol = :low,
                  proposal_covariance::Matrix = Matrix(undef, 0, 0),
                  mle::Bool = false,
                  sampling::Bool = true,
                  filestring_addl::Vector{String} = Vector{String}(),
                  continue_intermediate::Bool = false,
                  intermediate_stage_start::Int = 0,
                  intermediate_stage_increment::Int = 10,
		          save_intermediate::Bool = false)
    # Load data
    df = load_data(m; verbose = verbose)
    estimate(m, df; verbose = verbose, proposal_covariance = proposal_covariance,
             mle = mle, sampling = sampling,
             intermediate_stage_increment = intermediate_stage_increment,
	         save_intermediate = save_intermediate)
end

function estimate(m::AbstractDSGEModel, data::AbstractArray;
                  verbose::Symbol = :low,
                  proposal_covariance::Matrix = Matrix(undef, 0,0),
                  mle::Bool = false,
                  sampling::Bool = true,
                  filestring_addl::Vector{String} = Vector{String}(),
                  continue_intermediate::Bool = false,
                  intermediate_stage_start::Int = 0,
                  intermediate_stage_increment::Int = 10,
		          save_intermediate::Bool = false)

    if !(get_setting(m, :sampling_method) in [:SMC, :MH])
        error("method must be :SMC or :MH")
    else
        method = get_setting(m, :sampling_method)
    end

    ########################################################################################
    ### Step 1: Find posterior/likelihood mode (if reoptimizing, run optimization routine)
    ########################################################################################

    # Specify starting mode

    vint = get_setting(m, :data_vintage)
    if reoptimize(m) && method == :MH
        println("Reoptimizing...")

        # Inputs to optimization algorithm
        n_iterations       = get_setting(m, :optimization_iterations)
        ftol               = get_setting(m, :optimization_ftol)
        xtol               = get_setting(m, :optimization_xtol)
        gtol               = get_setting(m, :optimization_gtol)
        step_size          = get_setting(m, :optimization_step_size)
        converged          = false

        # If the algorithm stops only because we have exceeded the maximum number of
        # iterations, continue improving guess of modal parameters
        total_iterations = 0
        optimization_time = 0
        max_attempts = get_setting(m, :optimization_attempts)
        attempts = 1

        while !converged
            begin_time = time_ns()
            out, H = optimize!(m, data;
                               method = get_setting(m, :optimization_method),
                               ftol = ftol, grtol = gtol, xtol = xtol,
                               iterations = n_iterations, show_trace = true, step_size = step_size,
                               verbose = verbose,
                               mle = mle)

            attempts += 1
            total_iterations += out.iterations
            converged = out.converged || attempts > max_attempts

            end_time = (time_ns() - begin_time)/1e9
            println(verbose, :low, @sprintf "Total iterations completed: %d\n" total_iterations)
            println(verbose, :low, @sprintf "Optimization time elapsed: %5.2f\n" optimization_time += end_time)

            # Write params to file after every `n_iterations` iterations
            params = map(θ->θ.value, m.parameters)
            h5open(rawpath(m, "estimate", "paramsmode.h5"),"w") do file
                file["params"] = params
            end
        end

        # write parameters to file one last time so we have the final mode
        h5open(rawpath(m, "estimate", "paramsmode.h5"),"w") do file
            file["params"] = params
        end
    end

    params = map(θ->θ.value, m.parameters)

    # Sampling does not make sense if mle=true
    if mle || !sampling
        return nothing
    end

    if get_setting(m,:sampling_method) == :MH
        ########################################################################################
        ### Step 2: Compute proposal distribution for Markov Chain Monte Carlo (MCMC)
        ###
        ### In Metropolis-Hastings, we draw sample parameter vectors from
        ### the proposal distribution, which is a degenerate multivariate
        ### normal centered at the mode. Its variance is the inverse of
        ### the hessian. We find the inverse via eigenvalue decomposition.
        ########################################################################################

        ## Calculate the Hessian at the posterior mode
        hessian = if calculate_hessian(m)
            println(verbose, :low, "Recalculating Hessian...")

            hessian, _ = hessian!(m, params, data; verbose=verbose)

            h5open(rawpath(m, "estimate","hessian.h5"),"w") do file
                file["hessian"] = hessian
            end

            hessian

        ## Read in a pre-calculated Hessian
        else
            fn = hessian_path(m)
            println(verbose, :low, "Using pre-calculated Hessian from $fn")

            hessian = h5open(fn,"r") do file
                read(file, "hessian")
            end

            hessian
	    end

        # Compute inverse hessian and create proposal distribution, or
        # just create it with the given cov matrix if we have it
        propdist = if isempty(proposal_covariance)
            # Make sure the mode and hessian have the same number of parameters
            n = length(params)
            @assert (n, n) == size(hessian)

            # Compute the inverse of the Hessian via eigenvalue decomposition
            #S_diag, U = eigen(hessian)
            F = svd(hessian)

            big_eig_vals = findall(x -> x > 1e-6, F.S)
            hessian_rank = length(big_eig_vals)

            S_inv = zeros(n, n)
            #for i = (n-hessian_rank+1):n
            for i = 1:hessian_rank
                S_inv[i, i] = 1/F.S[i]
            end

            #hessian_inv = U*sqrt.(S_inv) # this is the inverse of the hessian
            hessian_inv = F.V * S_inv * F.U'#sqrt.(S_inv) * F.U'
            DegenerateMvNormal(params, hessian_inv, hessian, diag(S_inv))
        else
            DegenerateMvNormal(params, proposal_covariance, pinv(proposal_covariance),
                               eigen(proposal_covariance).values)
        end

        if rank(propdist) != n_parameters_free(m)
            println("problem –    shutting down dimensions")
        end

        ########################################################################################
        ### Step 3: Sample from posterior using Metropolis-Hastings algorithm
        ########################################################################################

        # Set the jump size for sampling
        cc0 = get_setting(m, :mh_cc0)
        cc  = get_setting(m, :mh_cc)

        metropolis_hastings(propdist, m, data, cc0, cc; verbose = verbose,
                            filestring_addl = filestring_addl);

    elseif get_setting(m, :sampling_method) == :SMC
        ########################################################################################
        ### Step 3: Run Sequential Monte Carlo (SMC)
        ###
        ### In Sequential Monte Carlo, a large number of Markov Chains
        ### are simulated iteratively to create a particle approximation
        ### of the posterior. Portions of this method are executed in
        ### parallel.
        ########################################################################################
        smc2(m, data; verbose = verbose, filestring_addl = filestring_addl,
            continue_intermediate = continue_intermediate,
            intermediate_stage_start = intermediate_stage_start,
            save_intermediate = save_intermediate,
            intermediate_stage_increment = intermediate_stage_increment)
    end

    ########################################################################################
    ### Step 4: Calculate and save parameter covariance matrix
    ########################################################################################

    compute_parameter_covariance(m, filestring_addl = filestring_addl)

    return nothing
end

"""
```
compute_parameter_covariance(m::AbstractDSGEModel)
```

Calculates the parameter covariance matrix from saved parameter draws, and writes it to the
parameter_covariance.h5 file in the `workpath(m, "estimate")` directory.

### Arguments
* `m::AbstractDSGEModel`: the model object
"""
function compute_parameter_covariance(m::AbstractDSGEModel;
                                      filestring_addl::Vector{String} = Vector{String}(undef, 0))

    sampling_method = get_setting(m, :sampling_method)
    if sampling_method ∉ [:MH, :SMC]
        throw("Invalid sampling method specified in setting :sampling_method")
    end

    prefix           = sampling_method == :MH ? "mh" : "smc"
    param_draws_path = rawpath(m, "estimate", prefix * "save.h5", filestring_addl)
    savepath         = workpath(m, "estimate", "parameter_covariance.h5", filestring_addl)

    return compute_parameter_covariance(param_draws_path, sampling_method; savepath = savepath)
end

"""
```
compute_parameter_covariance(param_draws_path::String, sampling_method::Symbol;
                             savepath = "parameter_covariance.h5")
```

Generic function calculates the parameter covariance matrix from saved parameter draws,
writes it to the specified savepath.

### Arguments
* `param_draws_path::String`: Path to file where parameter draws are stored.
* `sampling_method::Symbol`: Sampling method used for estimation.
```
   - `:MH`: Metropolis-Hastings
   - `:SMC`: Sequential Monte Carlo
```
### Optional Arguments
* `savepath::String`: Where parameter covariance matrix is to be saved. Will default
    to "parameter_covariance.h5" if unspecified.
"""
function compute_parameter_covariance(param_draws_path::String, sampling_method::Symbol;
                                      savepath::String = "parameter_covariance.h5")
    if sampling_method ∉ [:MH, :SMC]
        throw("Invalid sampling method specified in setting :sampling_method")
    end
    prefix = sampling_method == :MH ? "mh" : "smc"

    if !isfile(param_draws_path)
        @printf stderr "Saved parameter draws not found.\n"
        return
    end

    param_draws = h5open(param_draws_path, "r") do f
        read(f, prefix * "params")
    end

    # Calculate covariance matrix
    param_covariance = cov(param_draws)

    # Write to file
    h5open(savepath, "w") do f
        f[prefix * "cov"] = param_covariance
    end
end

"""
```
get_estimation_output_files(m)
```

Returns a `Dict{Symbol, String}` with all files created by `estimate(m)`.
"""
function get_estimation_output_files(m::AbstractDSGEModel)
    output_files = Dict{Symbol, String}()

    for file in [:paramsmode, :hessian, :mhsave]
        output_files[file] = rawpath(m, "estimate", "$file.h5")
    end

    for file in [:paramsmean, :parameter_covariance]
        output_files[file] = workpath(m, "estimate", "$file.h5")
    end

    for file in [:priors, :prior_posterior_means, :moments]
        output_files[file] = tablespath(m, "estimate", "$file.tex")
    end

    return output_files
end


text:
## Matrix whole space
function sample_network_wholespace!(
  res::MCMCSREvaluationCache,
  prob::HermitianMatrixProblem,
  net, σ)

  lnψ, ∇lnψ = logψ_and_∇logψ!(res.∇lnψ, net, σ)
  E         = compute_Cloc(prob, net, σ, lnψ, res.σ)

  prob = exp(2*real(lnψ))
  res.Zave   += prob
  res.Eave   += prob  * E
  push!(res.Evalues, prob*E)

  for (i, _∇lnψ)= enumerate(∇lnψ.tuple_all_weights)
    res.Oave[i]  .+= prob .* _∇lnψ
    res.OOave[i] .+= prob .* conj.(_∇lnψ) .* transpose(_∇lnψ)
    res.EOave[i] .+= prob .* conj.(_∇lnψ) .* E
  end
  return res
end

## Matrix sampled
function sample_network!(
  res::MCMCSREvaluationCache,
  prob::HermitianMatrixProblem,
  net, σ)

  lnψ, ∇lnψ = logψ_and_∇logψ!(res.∇lnψ, net, σ)
  E         = compute_Cloc(prob, net, σ, lnψ, res.σ)

  res.Eave   += E
  res.Zave   += 1.0 #exp(2*real(lnψ))
  push!(res.Evalues, E)

  for (i,_∇lnψ)=enumerate(∇lnψ.tuple_all_weights)
    res.Oave[i]  .+= _∇lnψ
    res.OOave[i] .+= conj.(_∇lnψ) .* transpose(_∇lnψ)
    res.EOave[i] .+= conj.(_∇lnψ) .* E
  end
  return res
end

function sample_network!(res::MCMCSRLEvaluationCache,
                         prob::LRhoSquaredProblem,
                         net, σ, wholespace=false)
  CLO_i = res.LLO_i
  update_lookup!(σ, net)

  lnψ, ∇lnψ = logψ_and_∇logψ!(res.∇lnψ, net, σ)
  C_loc = compute_Cloc!(CLO_i, res.∇lnψ2, prob, net, σ, lnψ, res.σ)

  prob = wholespace ? exp(2*real(lnψ)) : 1.0
  E = abs(C_loc)^2

  res.Zave += prob
  res.Eave += prob * E
  push!(res.Evalues, prob*E)

  for (i, _∇lnψ)= enumerate(∇lnψ.tuple_all_weights)
    res.Oave[i]  .+= prob .* _∇lnψ
    res.OOave[i] .+= prob .* conj.(_∇lnψ) .* transpose(_∇lnψ)
    res.EOave[i] .+= prob .* conj.(_∇lnψ) .* E
    res.LLOave[i].+= prob .* conj.(C_loc)  .* CLO_i[i]
  end
  return res
end

#=
function sample_network!(res::MCMCSREvaluationCache,
                         prob::OpenTimeEvolutionProblem,
                         net, σ, wholespace=false)
  lnψ, ∇lnψ = logψ_and_∇logψ!(res.∇lnψ, net, σ)
  E         = compute_Cloc(prob, net, σ, lnψ)

  prob = exp(2*real(lnψ))
  res.Zave   += prob
  res.Eave   += prob  * E
  push!(res.Evalues, prob*E)

  for (i, _∇lnψ)= enumerate(∇lnψ.tuple_all_weights)
    res.Oave[i]  .+= prob .* _∇lnψ
    res.OOave[i] .+= prob .* conj.(_∇lnψ) .* transpose(_∇lnψ)
    res.EOave[i] .+= prob .* conj.(_∇lnψ) .* E
  end
  return res
end=#


text:
using PyCall
@pyimport numpy

x = numpy.load("test_x.npy")
y = numpy.load("test_y.npy")

num_tests = 2246
println("Loading model...")
using Keras
m = Keras.load("model-structure.json", "model-weights.h5")
println("Model loaded, Testing....")
count = 0
for i=1:num_tests
    if (findmax(y[i,:])[2] == findmax(m(x[i,:]))[2])
        count += 1
    end
end

println("Model tested: $(count*100/num_tests)% accuracy achieved.")


text:
include("nnet.jl")

nnet = NNet("../nnet/TestNetwork.nnet")
println("Num Inputs: "*string(num_inputs(nnet)))
println("Num Outputs: "*string(num_outputs(nnet)))
println("One evaluation:")
println(evaluate_network(nnet,[15299.0,-1.142,-1.142,600.0,500.0]))
println("\nMultiple evaluations at once:")
println(evaluate_network_multiple(nnet,[[15299.0, -1.142, -1.142, 600.0, 500.0] [15299.0,-1.142,-1.142,600.0,500.0]]))


text:
using DistributionParameters
using Test
using ForwardDiff, LinearAlgebra, PaddedMatrices, StructuredMatrices, Random, ProbabilityModels,
        ProbabilityModels, LoopVectorization, VectorizationBase, SLEEFPirates, SIMDPirates

using DistributionParameters: Bounds, RealScalarParameter, RealArrayParameter, LKJCorrCholeskyParameter, SimplexParameter, ParamDescriptionCore, LengthParamDescription

N = 200;
dv = [LengthParamDescription(rand(0:7), ParamDescriptionCore(Bounds(-Inf,Inf),i,RealArrayParameter)) for i ∈ 1:N];
dvc = copy(dv);
DistributionParameters.sort_by_rem!(dvc, 8); dvc;
po = DistributionParameters.parameter_offsets(dvc);
sum(iszero.(po .% 4)) / N

@testset "DistributionParameters.jl" begin
    # Write your own tests here.

K = 7
S = (@Constant randn(K,4K)) |> x -> x * x'
S *= (1/32)
pS = StructuredMatrices.SymmetricMatrixL(S)
L = StructuredMatrices.lower_chol(S);
    L * L'
T = 36; δₜ = (1/16) * reduce(+, (@Constant randexp(T-1)) for i ∈ 1:8);
    times = vcat(zero(ConstantFixedSizeVector{1,Float64}), cumsum(δₜ));
    ρtup = (0.35,0.4,0.45,0.5,0.55,0.60,0.65);
    ρ = ConstantFixedSizeVector{7}(ρtup);

KT = K*T;

    ARs = PaddedMatrices.MutableFixedSizeArray{Tuple{T,T,K},Float64}(undef);

    using PaddedMatrices: DynamicPaddedArray

workspace = (
    Sigfull = PaddedMatrices.DynamicPaddedArray{Float64}(undef, (K*T,K*T), (KT+15) & -8),
    ARs = ARs,
    ∂ARs = PaddedMatrices.MutableFixedSizeArray{Tuple{T,T,K},Float64}(undef)
);
fill!(workspace.Sigfull, 19.5);
Cov, pAR, pL = DistributionParameters.∂DynamicCovarianceMatrix(ρ, L, times, workspace, Val((true,true)));
    using BenchmarkTools
#pL |> typeof
#spL |> typeof
#spL
    mtimes = MutableFixedSizeVector(times);
    sp, (spC, spAR, spL) = DistributionParameters.∂CovarianceMatrix(ProbabilityModels.STACK_POINTER, ρ, L, mtimes, Val((true,true)));
#    sp0 = PaddedMatrices.StackPointer(Libc.malloc(8*1<<28));

#    sp, (spC, spAR, spL) = DistributionParameters.∂DynamicCovarianceMatrix(sp0, ρ, L, mtimes, Val((true,true)));
    @benchmark DistributionParameters.∂DynamicCovarianceMatrix($ρ, $L, $times, $workspace, Val((true,true)))
    @benchmark DistributionParameters.∂CovarianceMatrix(ProbabilityModels.STACK_POINTER, $ρ, $L, $mtimes, Val((true,true)))
#exit()
    @code_warntype DistributionParameters.∂DynamicCovarianceMatrix(sp0, ρ, L, mtimes, Val((true,true)),Val(true));
#typeof(L)
spL.ARs
pL.ARs

@benchmark PaddedMatrices.PtrArray{Tuple{36,36,7}}($sp)
@benchmark PaddedMatrices.PtrArray{Tuple{36,36,7},Float64,3,40}($sp)
@benchmark PaddedMatrices.PtrArray{Tuple{36,36,7},Float64,3,40}(Base.unsafe_convert(Ptr{Float64},$sp.ptr))
    @benchmark pointer($sp, Float64) + 121
@code_warntype pointer(sp,Float64)

@code_warntype PaddedMatrices.PtrArray{Tuple{36,36,7},Float64,3,40}(sp)
@code_warntype PaddedMatrices.PtrArray{Tuple{36,36,7},Float64,3,40}(Base.unsafe_convert(Ptr{Float64},sp.ptr))

@code_warntype PaddedMatrices.PtrArray{Tuple{36,36,7},Float64,3,40}(sp0)

    aspl = Array(spL.ARs);
    apl = Array(pL.ARs);
    @test all(aspl .≈ apl)
    asl = Array(spL.LKJ);
    al = Array(pL.LKJ);
    @test all(asl .≈ al)
    asC = Array(Symmetric(spC,:L));
    aC = Array(Symmetric(Cov,:L));
    @test all(asC .≈ aC)
#    L
#    L_K
    mtimes = MutableFixedSizeVector(times);
    Lm = MutableLowerTriangularMatrix(Array(L))
    ρ, ρs

    Cov, pAR, pL = DistributionParameters.∂DynamicCovarianceMatrix(ρ, L, times, workspace, Val((true,true)));
    sp, (spC, spAR, spL) = DistributionParameters.∂CovarianceMatrix(ProbabilityModels.STACK_POINTER, ρ, Lm, mtimes, Val((true,true)));
    asC = Array(Symmetric(spC,:L));
    aC = Array(Symmetric(Cov,:L));
    findall(Symmetric(Array(Cov),:L) .≉  Symmetric(Array(spC),:L) )
    cholesky(Symmetric(aC))

    times, time

    size(spC)
    Symmetric(spC, :L)
    Symmetric(Cov, :L)
spAR
    size(spC)
    size(Cov)
    function fill_block_diag(B::AbstractArray{T,3}) where {T}
        M, N, K = size(B)
        A = zeros(M*K,N*K)
        for k ∈ 0:K-1
            A[1+k*M:(k+1)*M,1+k*N:(k+1)*N] .= @view B[:,:,k+1]
        end
        A
    end
    Lfull = kron(L, Matrix{Float64}(I, T, T)); ARfull = fill_block_diag(pL.ARs);

    @test ( Lfull * ARfull * Lfull' .≈ aC ) |> all
    @test ( Lfull * ARfull * Lfull' .≈ asC ) |> all
    cholesky(Symmetric(aC))

    Cov = DistributionParameters.CovarianceMatrix(ρ, L, times, workspace); Cov
#    Cov1s = PaddedMatrices.PaddedArray{Float64}(undef, (K*T,K*T)); Cov1s .= 1;
using Statistics

    # some disagreements!
    #

    using ForwardDiff, StaticArrays

     function DistributionParameters.CovarianceMatrix(
         rhos::PaddedMatrices.AbstractFixedSizeVector{K,T1},
         L::StructuredMatrices.AbstractLowerTriangularMatrix{K,T2},
         times::ConstantFixedSizeVector{nT}
     ) where {K,T1,T2,nT}

         T = promote_type(T1, T2)
         workspace = (
             Sigfull = PaddedMatrices.DynamicPaddedArray{T}(undef, (K*nT,K*nT)),
             ARs = PaddedMatrices.MutableFixedSizeArray{Tuple{nT,nT,K},T}(undef),
         )
         DistributionParameters.DynamicCovarianceMatrix(rhos, L, times, workspace)

     end

    ForwardDiff.gradient(r -> sum(Symmetric(DistributionParameters.CovarianceMatrix( ConstantFixedSizeVector{7}(ntuple(i -> r[i], Val(7))), L, times))), SVector(ρtup))

    zd = ones(K*T,K*T) ; #.- LowerTriangular(ones(K*T,K*T));
    zr = randn(KT,KT); @. zr = zr + zr';

    @test all(
        ForwardDiff.gradient(r -> sum(LowerTriangular(zd .* DistributionParameters.CovarianceMatrix( ConstantFixedSizeVector{7}(ntuple(i -> r[i], Val(7))), L, times))), SVector(ρtup))'
        .≈
        Array(zd * spAR)
    )
    @test all(
        ForwardDiff.gradient(r -> sum(LowerTriangular(zr .* DistributionParameters.CovarianceMatrix( ConstantFixedSizeVector{7}(ntuple(i -> r[i], Val(7))), L, times))), SVector(ρtup))'
        .≈
        Array(zr * spAR)
    )
zr * pAR
zr * spAR

@generated function StructuredMatrices.LowerTriangularMatrix(A::SMatrix{K,K,T}) where {K,T}
        outtup = Expr(:tuple, )
        for k ∈ 1:K
            push!(outtup.args, :(A[$k,$k]))
        end
        for kc ∈ 1:K, kr ∈ kc+1:K
            push!(outtup.args, :(A[$kr,$kc]))
        end
        :(@inbounds LowerTriangularMatrix{$K,$T,$(length(outtup.args))}( $outtup ) )
    end

using StaticArrays
    sL = SMatrix{K,K}(Array(L))
    LowerTriangularMatrix(sL).data
    L
    L.data
    ForwardDiff.gradient(r -> sum(Symmetric(DistributionParameters.CovarianceMatrix( ρ, LowerTriangularMatrix(r), times ))), sL )
    ForwardDiff.gradient(r -> sum(UpperTriangular(DistributionParameters.CovarianceMatrix( ρ, LowerTriangularMatrix(r), times ))), sL )

sL
ForwardDiff.gradient(r -> (Lf = kron(LowerTriangular(r), Matrix{Float64}(I,T,T)); sum(Lf * ARfull * Lf')), sL)
    @test all(
        LowerTriangular(ForwardDiff.gradient(r -> (Lf = kron(LowerTriangular(r), Matrix{Float64}(I,T,T)); sum(Lf * ARfull * Lf')), sL))
        .≈
        Array( zd * pL )
    )

    ForwardDiff.gradient(r -> sum( Symmetric(zr .* DistributionParameters.CovarianceMatrix( ρ, LowerTriangularMatrix(r), times ))), sL )

ForwardDiff.gradient(r -> (Lf = kron(LowerTriangular(r), Matrix{Float64}(I,T,T)); sum(zr .* (Lf * ARfull * Lf'))), sL)
@test all(
        LowerTriangular(ForwardDiff.gradient(r -> (Lf = kron(LowerTriangular(r), Matrix{Float64}(I,T,T)); sum(zr .* (Lf * ARfull * Lf'))), sL))
        .≈
        Array( zr * pL )
    )
#zr
        Array( LowerTriangular(zr) * pL )
        Array( zr * pL )
        LowerTriangular(ForwardDiff.gradient(r -> (Lf = kron(LowerTriangular(r), Matrix{Float64}(I,T,T)); sum((zr .* (Lf * ARfull * Lf')))), sL))
        LowerTriangular(ForwardDiff.gradient(r -> (Lf = kron(LowerTriangular(r), Matrix{Float64}(I,T,T)); sum(LowerTriangular(zr .* (Lf * ARfull * Lf')))), sL))

    sum(ARs, dims = (1,2))
zrpL = zr * pL
zrspL = zr * spL
Lm
L_K
    Ctest = zeros(KT,KT); Ctest[1+(K-1)*T:KT, 1+(K-1)*T:KT] .= pAR.∂ARs[:,:,K];
    (Lfull * Ctest * Lfull') |> sum

    ForwardDiff.gradient(
        r -> sum(UpperTriangular(zd .* DistributionParameters.CovarianceMatrix( ConstantFixedSizeVector{4}(ntuple(i -> r[i], Val(4))), L, times))), SVector(0.5,0.65,0.7,0.75))

vunstable(x::AbstractVector{T}) where {T} = ConstantFixedSizeArray{Tuple{length(x)},T,1}(ntuple(i -> x[i], length(x)))
x = randn(16,16);
ltx = LowerTriangular(x);
ltxv = vec(ltx);

function trisub2(N)
    out = collect(1:(N+1):N^2)
    ind = 0
    for n ∈ 0:N-1
        ind += n+1
        r = N-n-1
        for i ∈ 1:r
            ind += 1
            push!(out, ind)
        end
    end
    out
end
ltxv2 = ltx[trisub2(16)];
trivec = vunstable(ltxv2);

r2 = @Constant randn(120);
@generated function load_lkj(a::PaddedMatrices.AbstractFixedSizeVector{L,T}) where {L,T}
    fp = quote end
    sp = quote end
    out = gensym(:L)
    M = (Int(sqrt(1 + 8L))-1)>>1
    DistributionParameters.load_parameter(fp.args, sp.args, out, DistributionParameters.LKJ_Correlation_Cholesky{M+1}, false)
    quote
        $(Symbol("##θparameter##")) = VectorizationBase.vectorizable(a);
        target = zero($T)
        $fp
        $sp
        $out, target
    end
end
@generated function ∂load_lkj(a::PaddedMatrices.AbstractFixedSizeVector{L,T}, seedout = (M = (Int(sqrt(1 + 8L))-1)>>1; fill(1.0, PaddedMatrices.Static{(((M+2)*(M+1))>>1)}()))) where {L,T}
    fp = quote end
    sp = quote end
    out = gensym(:L)
    M = (Int(sqrt(1 + 8L))-1)>>1
    DistributionParameters.load_parameter(fp.args, sp.args, out, DistributionParameters.LKJ_Correlation_Cholesky{M+1}, true)
    quote
        $(Symbol("##∂θparameter##m")) = PaddedMatrices.MutableFixedSizeVector{L,T}(undef);
        $(Symbol("##∂θparameter##")) = VectorizationBase.vectorizable($(Symbol("##∂θparameter##m")))
        $(Symbol("##θparameter##")) = VectorizationBase.vectorizable(a);
        target = zero($T)
        $fp
        $(Symbol("###seed###", out)) = seedout'
        $sp
        $out, target, $(Symbol("##∂θparameter##m"))
    end
end
Lm, t, partial = ∂load_lkj(r2, trivec); Lm

r2v = Array(r2);
nr2vl = @. 1 / (1 + exp(r2v));
r2vl = 1 .- nr2vl;
zvals = @. 1 - 2 * nr2vl;
r2vlpartial = @. r2vl * nr2vl;
r2l = vunstable(zvals);

lkj, lkjac, lkjldgrad, lkjgrad = DistributionParameters.∂lkj_constrain(r2l);
j = ForwardDiff.jacobian(x -> Array(DistributionParameters.lkj_constrain(vunstable(x))[1]), r2l);
ltxvj = ltxv' * j;
@test all(ltxvj .≈ trivec' * lkjgrad) # matches
g = ForwardDiff.gradient(x -> DistributionParameters.lkj_constrain(vunstable(x))[2], r2l); g'
@test all(lkjldgrad .≈ g)

end


text:
using Test, YaoBase, YaoArrayRegister, LinearAlgebra, LuxurySparse, SparseArrays
using YaoBase.Const

# NOTE: we don't have block here, feel safe to use
using YaoBase.Const

@testset "test general unitary instruction" begin
    U1 = randn(ComplexF64, 2, 2)
    ST = randn(ComplexF64, 1 << 4)
    REG = ArrayReg(ST)
    I2 = IMatrix(2)
    M = kron(I2, U1, I2, I2) * ST

    @test instruct!(copy(ST), U1, 3) ≈ M ≈ instruct!(reshape(copy(ST), :, 1), U1, 3)

    U2 = rand(ComplexF64, 4, 4)
    M = kron(I2, U2, I2) * ST
    @test instruct!(copy(ST), U2, (2, 3)) ≈ M

    @test instruct!(copy(ST), kron(U1, U1), (3, 1)) ≈ instruct!(instruct!(copy(ST), U1, 3), U1, 1)
    @test instruct!(copy(REG), kron(U1, U1), (3, 1)) ≈ instruct!(instruct!(copy(REG), U1, 3), U1, 1)
    @test instruct!(transpose_storage(REG), kron(U1, U1), (3, 1)) ≈ instruct!(
        instruct!(copy(REG), U1, 3),
        U1,
        1,
    )
    @test instruct!(transpose_storage(REG), kron(U1, U1), (3, 1)) ≈ instruct!(
        instruct!(transpose_storage(REG), U1, 3),
        U1,
        1,
    )

    @test instruct!(reshape(copy(ST), :, 1), kron(U1, U1), (3, 1)) ≈ instruct!(
        instruct!(reshape(copy(ST), :, 1), U1, 3),
        U1,
        1,
    )

    U2 = sprand(ComplexF64, 8, 8, 0.1)
    ST = randn(ComplexF64, 1 << 5)
    M = kron(I2, U2, I2) * ST
    @test instruct!(copy(ST), U2, (2, 3, 4)) ≈ M

    @test instruct!(copy(ST), I2, (1,)) ≈ ST
end

@testset "test auto conversion" begin
    v = rand(ComplexF32, 1 << 8)
    @test_logs (
        :warn,
        "Element Type Mismatch: register Complex{Float32}, operator Complex{Float64}. Converting operator to match, this may cause performance issue",
    ) instruct!(v, Const.CNOT, (1, 2))
end

@testset "test general control unitary operator" begin
    ST = randn(ComplexF64, 1 << 5)
    U1 = randn(ComplexF64, 2, 2)
    instruct!(copy(ST), U1, (3,), (1,), (1,))

    @test instruct!(copy(ST), U1, (3,), (1,), (1,)) ≈ general_controlled_gates(
        5,
        [P1],
        [1],
        [U1],
        [3],
    ) * ST
    @test instruct!(copy(ST), U1, (3,), (1,), (0,)) ≈ general_controlled_gates(
        5,
        [P0],
        [1],
        [U1],
        [3],
    ) * ST

    # control U2
    U2 = kron(U1, U1)
    @test instruct!(copy(ST), U2, (3, 4), (1,), (1,)) ≈ general_controlled_gates(
        5,
        [P1],
        [1],
        [U2],
        [3],
    ) * ST

    # multi-control U2
    @test instruct!(copy(ST), U2, (3, 4), (5, 1), (1, 0)) ≈ general_controlled_gates(
        5,
        [P1, P0],
        [5, 1],
        [U2],
        [3],
    ) * ST
end

@testset "test Pauli instructions" begin
    @testset "test $G instructions" for (G, M) in zip((:X, :Y, :Z), (X, Y, Z))
        @test linop2dense(s -> instruct!(s, Val(G), (1,)), 1) == M
        @test linop2dense(s -> instruct!(s, Val(G), (1, 2, 3)), 3) == kron(M, M, M)
    end

    @testset "test controlled $G instructions" for (G, M) in zip((:X, :Y, :Z), (X, Y, Z))
        @test linop2dense(s -> instruct!(s, Val(G), 4, (2, 1), (0, 1)), 4) ≈ general_controlled_gates(
            4,
            [P0, P1],
            [2, 1],
            [M],
            [4],
        )

        @test linop2dense(s -> instruct!(s, Val(G), 1, 2, 0), 2) ≈ general_controlled_gates(
            2,
            [P0],
            [2],
            [M],
            [1],
        )
    end
end

@testset "single qubit instruction" begin
    ST = randn(ComplexF64, 1 << 4)
    Pm = pmrand(ComplexF64, 2)
    Dv = Diagonal(randn(ComplexF64, 2))

    @test instruct!(copy(ST), Pm, 3) ≈ kron(I2, Pm, I2, I2) * ST ≈
          instruct!(reshape(copy(ST), :, 1), Pm, 3)
    @test instruct!(copy(ST), Dv, 3) ≈ kron(I2, Dv, I2, I2) * ST ≈
          instruct!(reshape(copy(ST), :, 1), Dv, 3)
end

@testset "swap instruction" begin
    ST = randn(ComplexF64, 1 << 2)
    @test instruct!(copy(ST), Val(:SWAP), (1, 2)) ≈ SWAP * ST
end

@testset "pswap instruction" begin
    ST = randn(ComplexF64, 1 << 2)
    θ = π / 3
    @test instruct!(copy(ST), Val(:PSWAP), (1, 2), θ) ≈ (cos(θ / 2) * IMatrix{4}() -
                                                         im * sin(θ / 2) * SWAP) * ST

    T = ComplexF64
    theta = 0.5
    for (R, G) in [(:Rx, X), (:Ry, Y), (:Rz, Z), (:PSWAP, SWAP)]
        @test rot_mat(T, Val(R), theta) ≈ rot_mat(T, G, theta)
    end
    @test rot_mat(T, Val(:CPHASE), theta) ≈ rot_mat(T, Diagonal([1, 1, 1, -1]), theta) *
                                            exp(im * theta / 2)
    for ST in [randn(ComplexF64, 1 << 5), randn(ComplexF64, 1 << 5, 10)]
        for R in [:Rx, :Ry, :Rz]
            @test instruct!(copy(ST), Val(R), (4,), θ) ≈ instruct!(
                copy(ST),
                Matrix(rot_mat(T, Val(R), θ)),
                (4,),
            )
            @test instruct!(copy(ST), Val(R), (4,), (1,), (0,), θ) ≈ instruct!(
                copy(ST),
                Matrix(rot_mat(T, Val(R), θ)),
                (4,),
                (1,),
                (0,),
            )
        end
        for R in [:CPHASE, :PSWAP]
            @test instruct!(copy(ST), Val(R), (4, 2), θ) ≈ instruct!(
                copy(ST),
                Matrix(rot_mat(T, Val(R), θ)),
                (4, 2),
            )
            instruct!(copy(ST), Val(R), (4, 2), (1,), (0,), θ)
            instruct!(copy(ST), Matrix(rot_mat(T, Val(R), θ)), (4, 2), (1,), (0,))
            @test instruct!(copy(ST), Val(R), (4, 2), (1,), (0,), θ) ≈ instruct!(
                copy(ST),
                Matrix(rot_mat(T, Val(R), θ)),
                (4, 2),
                (1,),
                (0,),
            )
        end
    end
end

@testset "Yao.jl/#189" begin
    st = rand(1 << 4)
    @test instruct!(st, IMatrix{2,Float64}(), 1) == st
end

@testset "test empty locs" begin
    st = rand(ComplexF64, 1 << 4)
    pm = pmrand(ComplexF64, 2)
    @test instruct!(copy(st), pm, ()) == st

    for G in [:Z, :S, :T, :Sdag, :Tdag]
        @test instruct!(copy(st), Val(G), ()) == st
    end
end


text:
"""
	PMA <: Factorization

The output of `pma` representing the Principal Moment Analysis factorization of a matrix `A`.
If `F::PMA` is the factorization object, `U`, `S`, `V` and `Vt` can be obtained via `F.U`, `F.S`, `F.V` and `F.Vt` such as A ≈ U * Diagonal(S) * Vt'.

See also `pma`.
"""
struct PMA{T} <: Factorization{T}
	U::Matrix{T}
	S::Vector{T}
	Vt::Matrix{T}
	VV::Matrix{T}
end

Base.size(F::PMA, dim::Integer) = dim == 1 ? size(F.U, dim) : size(F.Vt, dim)
Base.size(F::PMA) = (size(F, 1), size(F, 2))

function Base.getproperty(F::PMA,name::Symbol)
	name==:V && return F.Vt'
	getfield(F,name)
end
Base.propertynames(::PMA) = (:U,:S,:V,:Vt,:VV)

function simplices2kernelmatrix(sg::SimplexGraph)
	A = simplices2kernelmatrixroot(sg, simplify=false)
	A*A'
end

function simplices2kernelmatrixroot(sg::SimplexGraph; simplify=true)::SparseMatrixCSC
	s = vec(sum(sg.G,dims=1))

	L = sg.G * Diagonal(sqrt.(sg.w./(s.*(s.+1))))
	D = Diagonal(sqrt.(vec(sum(x->x^2,L,dims=2))))
	simplify || return [L D]

	QR = qr(vcat(sparse(L'), D))
	sparse(QR.R[:,invperm(QR.pcol)]')
end

function _pma(A::AbstractMatrix, S::AbstractMatrix; nsv::Integer=6)
	P,N = size(A)
	L = size(S,2)
	Y = A*S

	K = Symmetric(Y'Y)
	nsv = min(nsv, P, N)
	F = eigen(K, L-nsv+1:L)
	Σ = sqrt.(max.(0.,reverse(F.values)))
	VV = F.vectors[:,end:-1:1] # Coordinates of simplex equivalents, not interesting in practice.

	U = Y*VV ./ Σ'
	Vt = U'A ./ Σ # Coordinates of original sample points in low dimensional space.

	PMA(U,Σ,Vt,VV)
end

"""
	pma(A, G::SimplexGraph; nsv=6)

Computes the Principal Moment Analysis of the matrix `A` (variables × samples) using the sample SimplexGraph `G`.
Each column in `G` is a boolean vector representing one simplex. True means that a vertex is part of the simplex and false that it is not.
Set `nsv` to control the number of singular values and vectors returned.
Returns a `PMA` struct.

See also `PMA`, `SimplexGraph`.
"""
pma(A::AbstractMatrix, sg::SimplexGraph; kwargs...) = _pma(A, simplices2kernelmatrixroot(sg); kwargs...)

"""
	pma(A, G::AbstractMatrix{Bool}; nsv=6)

Computes the Principal Moment Analysis of the matrix `A` (variables × samples).
Each column in `G` is a boolean vector representing one simplex. True means that a vertex is part of the simplex and false that it is not.
Set `nsv` to control the number of singular values and vectors returned.
Returns a `PMA` struct.

See also `PMA`, `SimplexGraph`.
"""
pma(A::AbstractMatrix, G::AbstractMatrix{Bool}; kwargs...) = pma(A, SimplexGraph(G); kwargs...)


text:
Pkg.activate(dirname(@__DIR__))
using BenchmarkTools, VoxelSpace
using VoxelSpace: WIN_WIDTH, WIN_HEIGHT, render!

function bench_render()
    datac, datah = VoxelSpace.read_map("C1W")
    fbuffer = similar(datac,WIN_WIDTH,WIN_HEIGHT)
    hbuffer = Vector{Int}(undef,WIN_WIDTH)
    pos = (512f0,512f0,78f0)

    @benchmark render!($fbuffer,$hbuffer,$datac,$datah,$pos,θ) setup=θ=2*rand(Float32)*π
end

function profile_render()
    datac, datah = VoxelSpace.read_map("C1W")
    fbuffer = similar(datac,WIN_WIDTH,WIN_HEIGHT)
    hbuffer = Vector{Int}(undef,WIN_WIDTH)
    pos = (512f0,512f0,78f0)

    @profiler for i = 1:1000
        θ = 2*rand(Float32)*π
        render!(fbuffer,hbuffer,datac,datah,pos,θ)
    end
end

# VoxelSpace.run()

max_stars_repo_name: jebej/VoxelSpace

max_issues_repo_name: jebej/VoxelSpace

max_forks_repo_name: jebej/VoxelSpace


text:
using Pkg

cd("Julia")

Pkg.activate(".")
Pkg.instantiate()

using Revise

include("01_wrangling.jl")
include("../src/tweet_helpers_refac.jl")
include("../src/tweet_graphing.jl")

using Chain
using CairoMakie
using AlgebraOfGraphics
using DataFrames
using DataFramesMeta
using Dates
using Graphs
using SparseArrays
using DataStructures
using MetaGraphs
using GraphPlot
using Compose
import Cairo
import Compose

set_aog_theme!()

function graph2(df, date_func)
    fig = Figure()

    g1 = plot_graph(get_column_summary(df, :quote_count, date_func, sum),
        ["Day of year", "N quotes", "Daily quote volume"], false, fig[1,1])
    g2 = plot_graph(get_column_summary(df, :retweet_count, date_func, sum),
        ["Day of year", "N retweets", "Daily retweet volume"], false, fig[1,2])
    g3 = plot_graph(get_column_summary(df, :reply_count, date_func, sum),
        ["Day of year", "N replies", "Daily reply volume"], false, fig[2,1])

    return fig
end

function graph4(df, date_func)
    fig = Figure()

    g1 = plot_graph(@transform(get_column_summary(df, :quote_count, date_func, sum), :y = :y/1440),
        ["Day of year", "Avg n quotes per min", "Daily quote rate"], true, fig[1,1])
    g2 = plot_graph(@transform(get_column_summary(df, :retweet_count, date_func, sum), :y = :y/1440),
        ["Day of year", "Avg n retweets per min", "Daily retweet rate"], true, fig[1,2])
    g3 = plot_graph(@transform(get_column_summary(df, :reply_count, date_func, sum), :y = :y/1440),
        ["Day of year", "Avg n replies per min", "Daily reply rate"], true, fig[2,1])

    return fig
end

function graph5(df, date_func)
    fig = Figure()

    g1 = plot_graph(cumulative_sorting(df, :quote_count),
        ["Author", "Cumulative n quotes", "Cumulative quotes by author"], false, fig[1,1], true)
    g2 = plot_graph(cumulative_sorting(df, :retweet_count),
        ["Author", "Cumulative n retweets", "Cumulative retweets by author"], false, fig[1,2], true)
    g3 = plot_graph(cumulative_sorting(df, :reply_count),
        ["Author", "Cumulative n replies", "Cumulative replies by author"], false, fig[2,1], true)

    return fig
end

function _plot_stats(df, date_func, dir)
    save(dir * "/graph1.png", plot_graph(get_column_summary(df, :author_id, date_func, length ∘ unique), ["Day of year", "N unique users", "Daily unique users"], false), px_per_unit = 3)
    save(dir * "/graph2.png", graph2(df, date_func), px_per_unit = 3)
    save(dir * "/graph3.png", plot_graph(@transform(get_column_summary(df, :author_id, date_func, length ∘ unique), :y = :y/1440),  ["Day of year", "N unique users per min", "Daily unique user rate"], true), px_per_unit = 3)
    save(dir * "/graph4.png", graph4(df, date_func), px_per_unit = 3)
    save(dir * "/graph5.png", graph5(df, date_func), px_per_unit = 3)
end

function _plot_graphs(df, dir)
    g,userL = generate_author_retweet_source_graph(df)
    author_mg = generate_unipartite_projection(g, userL)
    Compose.draw(Compose.PNG(dir * "/graph6.png", 16cm, 16cm), generate_most_connected_user_author_plot(author_mg))
    Compose.draw(Compose.PNG(dir * "/graph7.png", 16cm, 16cm), generate_most_connected_user_author_plot(author_mg))
end

function gen_stat_plots(date_func, out_dir, dfs = read_files("Data"))
    scandals = read_names("Data")

    for s in scandals
        if !isdir(out_dir * "/" * s)
            mkdir(out_dir * "/" * s)
        end
    end

    for i = 1:length(dfs)
        _plot_stats(dfs[i], date_func, out_dir * "/" * scandals[i])
    end
end

function gen_graph_plots(out_dir)
    dfs = read_files("Data")
    scandals = read_names("Data")

    for s in scandals
        if !isdir(out_dir * "/" * s)
            mkdir(out_dir * "/" * s)
        end
    end

    for i = 1:length(dfs)
        _plot_graphs(dfs[i], out_dir * "/" * scandals[i])
    end
end

function _plot_net_met(df, met_funcs, per, cum, dir)
    g = generate_reply_quote_graph(df)
    mg = add_meta_reply_quote_graph(df, g)
    for met_func in met_funcs
        save(dir * "/" * string(Symbol(met_func)) * ".png", plot_graph_metric(mg, met_func, per, Date.(ZonedDateTime.(df.created_at))[1], Date.(ZonedDateTime.(df.created_at))[size(df,1)], cum), px_per_unit = 3)
    end
end

function gen_net_plots(met_funcs, per, cum, out_dir)
    dfs = read_files("Data")
    scandals = read_names("Data")

    for s in scandals
        if !isdir(out_dir * "/" * s)
            mkdir(out_dir * "/" * s)
        end
    end

    for i = 1:length(dfs)
        _plot_net_met(dfs[i], met_funcs, per, cum, out_dir * "/" * scandals[i])
    end
end

function graph_metric(x)
    return mean(degree(x))
end

functions = [graph_metric]
#gen_net_plots(functions, Day(1), false, "Plots2")

dfs = [unique(df,[:tweet_id,:author_id]) for df in read_files("Data")]

get_date(x) = Date(x)

gen_stat_plots(dayofyear, "PlotsDate", dfs)
#gen_graph_plots("Plots")


text:
##############################################################################
##
## We use SubDataFrame's to maintain a reference to a subset of a DataFrame
## without making copies.
##
##############################################################################

struct SubDataFrame{T <: AbstractVector{Int}} <: AbstractDataFrame
    parent::DataFrame
    rows::T # maps from subdf row indexes to parent row indexes

    function SubDataFrame{T}(parent::DataFrame, rows::T) where {T <: AbstractVector{Int}}
        if length(rows) > 0
            rmin, rmax = extrema(rows)
            if rmin < 1 || rmax > size(parent, 1)
                throw(BoundsError())
            end
        end
        new(parent, rows)
    end
end

"""
A view of row subsets of an AbstractDataFrame

A `SubDataFrame` is meant to be constructed with `view`.  A
SubDataFrame is used frequently in split/apply sorts of operations.

```julia
view(d::AbstractDataFrame, rows)
```

### Arguments

* `d` : an AbstractDataFrame
* `rows` : any indexing type for rows, typically an Int,
  AbstractVector{Int}, AbstractVector{Bool}, or a Range

### Notes

A `SubDataFrame` is an AbstractDataFrame, so expect that most
DataFrame functions should work. Such methods include `describe`,
`dump`, `nrow`, `size`, `by`, `stack`, and `join`. Indexing is just
like a DataFrame; copies are returned.

To subset along columns, use standard column indexing as that creates
a view to the columns by default. To subset along rows and columns,
use column-based indexing with `view`.

### Examples

```julia
df = DataFrame(a = repeat([1, 2, 3, 4], outer=[2]),
               b = repeat([2, 1], outer=[4]),
               c = randn(8))
sdf1 = view(df, 1:6)
sdf2 = view(df, df[:a] .> 1)
sdf3 = view(df[[1,3]], df[:a] .> 1)  # row and column subsetting
sdf4 = groupby(df, :a)[1]  # indexing a GroupedDataFrame returns a SubDataFrame
sdf5 = view(sdf1, 1:3)
sdf1[:,[:a,:b]]
```

"""
SubDataFrame

function SubDataFrame(parent::DataFrame, rows::T) where {T <: AbstractVector{Int}}
    return SubDataFrame{T}(parent, rows)
end

function SubDataFrame(parent::DataFrame, rows::Colon)
    return SubDataFrame(parent, 1:nrow(parent))
end

function SubDataFrame(parent::DataFrame, row::Integer)
    return SubDataFrame(parent, [Int(row)])
end

function SubDataFrame(parent::DataFrame, rows::AbstractVector{<:Integer})
    return SubDataFrame(parent, convert(Vector{Int}, rows))
end

function SubDataFrame(parent::DataFrame, rows::AbstractVector{Bool})
    return SubDataFrame(parent, findall(rows))
end

function SubDataFrame(sdf::SubDataFrame, rowinds::Union{T, AbstractVector{T}}) where {T <: Integer}
    return SubDataFrame(sdf.parent, sdf.rows[rowinds])
end

function SubDataFrame(sdf::SubDataFrame, rowinds::Colon)
    return sdf
end

function Base.view(adf::AbstractDataFrame, rowinds::AbstractVector{T}) where {T >: Missing}
    # Vector{>:Missing} need to be checked for missings
    any(ismissing, rowinds) && throw(MissingException("missing values are not allowed in indices"))
    return SubDataFrame(adf, convert(Vector{Missings.T(T)}, rowinds))
end

function Base.view(adf::AbstractDataFrame, rowinds::Any)
    return SubDataFrame(adf, rowinds)
end

function Base.view(adf::AbstractDataFrame, rowinds::Any, colinds::Union{Colon, AbstractVector})
    return SubDataFrame(adf[colinds], rowinds)
end

function Base.view(adf::AbstractDataFrame, rowinds::Any, colind::ColumnIndex)
    return SubDataFrame(adf[[colind]], rowinds)
end

##############################################################################
##
## AbstractDataFrame interface
##
##############################################################################

index(sdf::SubDataFrame) = index(sdf.parent)

# TODO: Remove these
nrow(sdf::SubDataFrame) = ncol(sdf) > 0 ? length(sdf.rows)::Int : 0
ncol(sdf::SubDataFrame) = length(index(sdf))

function Base.getindex(sdf::SubDataFrame, colinds::Any)
    return sdf.parent[sdf.rows, colinds]
end

function Base.getindex(sdf::SubDataFrame, rowinds::Any, colinds::Any)
    return sdf.parent[sdf.rows[rowinds], colinds]
end

function Base.setindex!(sdf::SubDataFrame, val::Any, colinds::Any)
    sdf.parent[sdf.rows, colinds] = val
    return sdf
end

function Base.setindex!(sdf::SubDataFrame, val::Any, rowinds::Any, colinds::Any)
    sdf.parent[sdf.rows[rowinds], colinds] = val
    return sdf
end

##############################################################################
##
## Miscellaneous
##
##############################################################################

Base.map(f::Function, sdf::SubDataFrame) = f(sdf) # TODO: deprecate

without(sdf::SubDataFrame, c) = view(without(sdf.parent, c), sdf.rows)