text:
# POMDP model functions
"""
    POMDP{S,A,O}

Abstract base type for a partially observable Markov decision process.

    S: state type
    A: action type
    O: observation type
"""
abstract type POMDP{S,A,O} end

"""
    MDP{S,A}

Abstract base type for a fully observable Markov decision process.

    S: state type
    A: action type
"""
abstract type MDP{S,A} end

"""
    n_states(problem::POMDP)
    n_states(problem::MDP)

Return the number of states in `problem`. Used for discrete models only.
"""
function n_states end

"""
    n_actions(problem::POMDP)
    n_actions(problem::MDP)

Return the number of actions in `problem`. Used for discrete models only.
"""
function n_actions end

"""
    n_observations(problem::POMDP)

Return the number of observations in `problem`. Used for discrete models only.
"""
function n_observations end

"""
    discount(problem::POMDP)
    discount(problem::MDP)

Return the discount factor for the problem.
"""
function discount end

"""
    transition{S,A,O}(problem::POMDP{S,A,O}, state::S, action::A)
    transition{S,A}(problem::MDP{S,A}, state::S, action::A)

Return the transition distribution from the current state-action pair
"""
function transition end

"""
    observation{S,A,O}(problem::POMDP{S,A,O}, statep::S)
    observation{S,A,O}(problem::POMDP{S,A,O}, action::A, statep::S)
    observation{S,A,O}(problem::POMDP{S,A,O}, state::S, action::A, statep::S)

Return the observation distribution. You need only define the method with the fewest arguments needed to determine the observation distribution.

# Example
```julia
using POMDPToolbox # for SparseCat

struct MyPOMDP <: POMDP{Int, Int, Int} end

observation(p::MyPOMDP, sp::Int) = SparseCat([sp-1, sp, sp+1], [0.1, 0.8, 0.1])
```
"""
function observation end

"""
    observation{S,A,O}(problem::POMDP{S,A,O}, action::A, statep::S)

Return the observation distribution for the a-s' tuple (action and next state)
"""
observation(problem::POMDP, a, sp) = observation(problem, sp)
# @impl_dep observation(::P,::A,::S) where {P<:POMDP,S,A} observation(::P,::S)
@impl_dep observation(::P,::A,::S) where {P<:POMDP,S,A} observation(::P,::S)

"""
    observation{S,A,O}(problem::POMDP{S,A,O}, state::S, action::A, statep::S)

Return the observation distribution for the s-a-s' tuple (state, action, and next state)
"""
observation(problem::POMDP, s, a, sp) = observation(problem, a, sp)
@impl_dep observation(::P,::S,::A,::S) where {P<:POMDP,S,A} observation(::P,::A,::S)

"""
    reward{S,A,O}(problem::POMDP{S,A,O}, state::S, action::A)
    reward{S,A}(problem::MDP{S,A}, state::S, action::A)

Return the immediate reward for the s-a pair
"""
function reward end

"""
    reward{S,A,O}(problem::POMDP{S,A,O}, state::S, action::A, statep::S)
    reward{S,A}(problem::MDP{S,A}, state::S, action::A, statep::S)

Return the immediate reward for the s-a-s' triple
"""
reward(problem::Union{POMDP,MDP}, s, a, sp) = reward(problem, s, a)
@impl_dep reward(::P,::S,::A,::S) where {P<:Union{POMDP,MDP},S,A} reward(::P,::S,::A)

"""
    isterminal_obs{S,A,O}(problem::POMDP{S,A,O}, observation::O)

Check if an observation is terminal.
"""
isterminal_obs(problem::POMDP, observation) = false

"""
    isterminal_act{S,A,O}(problem::POMDP{S,A,O}, action::A)
    isterminal_act{S,A}(problem::MDP{S,A}, action::A)

Check if an action is terminal.
"""
isterminal_act(problem::Union{POMDP,MDP}, action) = false

"""
    isterminal{S,A,O}(problem::POMDP{S,A,O}, state::S)
    isterminal{S,A}(problem::MDP{S,A}, state::S)

Check if state s is terminal
"""
isterminal(problem::Union{POMDP,MDP}, state) = false

"""
    isterminal{S,A,O}(problem::POMDP{S,A,O}, state::S, action::A, observation::O)

Check if state s, action a or observation o is terminal
"""
isterminal(problem::POMDP, state, action, observation) = isterminal(problem, state, action) || isterminal_obs(problem, observation)

"""
    isterminal{S,A,O}(problem::POMDP{S,A,O}, state::S, action::A)
    isterminal{S,A}(problem::MDP{S,A}, state::S, action::A)

Check if state s or action a is terminal
"""
isterminal(problem::Union{POMDP,MDP}, state, action) = isterminal(problem, state) || isterminal_act(problem, action)

"""
    initialstate_distribution(pomdp::POMDP)
    initialstate_distribution(mdp::MDP)

Return a distribution of the initial state of the pomdp or mdp.
"""
function initialstate_distribution end
@deprecate initial_state_distribution initialstate_distribution

"""
    stateindex{S,A,O}(problem::POMDP{S,A,O}, s::S)
    stateindex{S,A}(problem::MDP{S,A}, s::S)

Return the integer index of state `s`. Used for discrete models only.
"""
function stateindex end
@deprecate state_index stateindex

"""
    actionindex{S,A,O}(problem::POMDP{S,A,O}, a::A)
    actionindex{S,A}(problem::MDP{S,A}, a::A)

Return the integer index of action `a`. Used for discrete models only.
"""
function actionindex end
@deprecate action_index actionindex

"""
    obsindex{S,A,O}(problem::POMDP{S,A,O}, o::O)

Return the integer index of observation `o`. Used for discrete models only.
"""
function obsindex end
@deprecate obs_index obsindex

"""
    convert_s(::Type{V}, s, problem::Union{MDP,POMDP}) where V<:AbstractArray
    convert_s(::Type{S}, vec::V, problem::Union{MDP,POMDP}) where {S,V<:AbstractArray}

Convert a state to vectorized form or vice versa.
"""
function convert_s end

convert_s(T::Type{A1}, s::A2, problem::Union{MDP, POMDP}) where {A1<:AbstractArray, A2<:AbstractArray} = convert(T, s)

convert_s(::Type{A}, s::Number, problem::Union{MDP,POMDP}) where A<:AbstractArray{F} where F<:Number = F[s]
convert_s(::Type{N}, v::AbstractArray{F}, problem::Union{MDP, POMDP}) where {N<:Number, F<:Number} = convert(N, first(v))

"""
    convert_a(::Type{V}, a, problem::Union{MDP,POMDP}) where V<:AbstractArray
    convert_a(::Type{A}, vec::V, problem::Union{MDP,POMDP}) where {A,V<:AbstractArray}

Convert an action to vectorized form or vice versa.
"""
function convert_a end

convert_a(T::Type{A1}, s::A2, problem::Union{MDP, POMDP}) where {A1<:AbstractArray, A2<:AbstractArray} = convert(T, s)

convert_a(::Type{A}, s::Number, problem::Union{MDP,POMDP}) where A<:AbstractArray{F} where F<:Number = F[s]
convert_a(::Type{N}, v::AbstractArray{F}, problem::Union{MDP, POMDP}) where {N<:Number, F<:Number} = convert(N, first(v))

"""
    convert_o(::Type{V}, o, problem::Union{MDP,POMDP}) where V<:AbstractArray
    convert_o(::Type{O}, vec::V, problem::Union{MDP,POMDP}) where {O,V<:AbstractArray}

Convert an observation to vectorized form or vice versa.
"""
function convert_o end

convert_o(T::Type{A1}, s::A2, problem::Union{MDP, POMDP}) where {A1<:AbstractArray, A2<:AbstractArray} = convert(T, s)

convert_o(::Type{A}, s::Number, problem::Union{MDP,POMDP}) where A<:AbstractArray{F} where F<:Number = F[s]
convert_o(::Type{N}, v::AbstractArray{F}, problem::Union{MDP, POMDP}) where {N<:Number, F<:Number} = convert(N, first(v))

# XXX DEPRECATED - remove in version 0.7
function Base.convert(T::Type, x::X, problem::Union{MDP,POMDP}) where X
    Base.depwarn("POMDPs.convert is deprecated. Please use convert_s, convert_a, or convert_o.", :convert)
    return convert_s(T, x, problem)
end

function max_possible_weight end
function new_particle end
function reset_distribution end


text:
using Plots
push!(LOAD_PATH, "./src")
using Basis2DTri
using Basis1D
theme(:wong)
gr(aspect_ratio=1,border=:none,legend=:bottomright,axis=nothing,ticks=nothing,legendfontsize=14,xlims=[-1,4],
   markerstrokewidth=0,markersize=1)
r1D, w1D = gauss_quad(0,0,2)
Nfp = length(r1D) # number of points per face
e = ones(Nfp) # vector of all ones
z = zeros(Nfp) # vector of all zeros
rf = [r1D; -r1D; -e];
sf = [-e; r1D; -r1D];

r_tri,s_tri,_ = Basis2DTri.quad_nodes_tri(4)
plot([-1;-1],[-1;1],seriescolor=:darkorange1,lw=5,label = nothing)
plot!([-1;1],[1;-1],seriescolor=:darkorange1,lw=5,label = nothing)
plot!([-1;1],[-1;-1],seriescolor=:darkorange1,lw=5,label = nothing)
plot!([-1;-1].+3,[-1;1].+2,seriescolor=:darkorange1,lw=5,label = nothing)
plot!([-1;1].+3,[1;-1].+2,seriescolor=:darkorange1,lw=5,label = nothing)
plot!([-1;1].+3,[-1;-1].+2,seriescolor=:darkorange1,lw=5,label = nothing)
plot!([1;4],[-1,1],seriescolor=:darkorange1,lw=5,label = nothing)
plot!([-1;2],[-1,1],seriescolor=:darkorange1,lw=5,ls=:dash,label = nothing)
plot!([-1;2],[1,3],seriescolor=:darkorange1,lw=5,label = nothing)
savefig("wedge.png")

ext: jl

max_stars_repo_name: yiminllin/ESDG-Fourier

max_issues_repo_name: yiminllin/ESDG-Fourier

max_forks_repo_name: yiminllin/ESDG-Fourier


text:
include("gp_brute_force_kernels.jl")

function make_kernel_data(x::AA1, kernel::Function; sigmasq_obs::AA2 = zeros(length(x)), sigmasq_cor::Real=1.0, rho::Real=1.0)  where { T1<:Real, AA1<:AbstractArray{T1,1},  T2<:Real, AA2<:AbstractArray{T2,1} }
  @assert length(x) == length(sigmasq_obs)
  K = diagm(sigmasq_obs)
  for i in 1:size(K,1)
      for j in 1:size(K,2)
	      K[i,j] += kernel(x[i]-x[j], sigmasq=sigmasq_cor, rho=rho)
	  end
  end
  return PDMat(K)
end

function make_kernel_obs_pred(xobs::AA1, xpred::AA2,
			kernel::Function; sigmasq_cor::Real=1.0, rho::Real=1.0)  where { T1<:Real, AA1<:AbstractArray{T1,1}, T2<:Real, AA2<:AbstractArray{T2,1}, T<:Real }
  K = zeros(length(xobs),length(xpred))
  for i in 1:length(xobs)
      for j in 1:length(xpred)
	      K[i,j] += kernel(xobs[i]-xpred[j], sigmasq=sigmasq_cor, rho=rho)
	  end
  end
  return K
end

ncalls = 0
function reset_ncalls()
	global ncalls = 0
end

function predict_mean(xobs::AA1, yobs::AA2, xpred::AA3; kernel::Function = matern52_sparse_kernel,
			sigmasq_obs::AA4 = 1e-16*ones(length(xobs)), sigmasq_cor::Real=1.0, rho::Real=1.0)  where {
				T1<:Real, AA1<:AbstractArray{T1,1}, T2<:Real, AA2<:AbstractArray{T2,1}, T3<:Real, AA3<:AbstractArray{T3,1}, T4<:Real, AA4<:AbstractArray{T4,1} }
	# global ncalls += 1
  	println("# size(xobs) = ",size(xobs), "  size(xpred) = ", size(xpred))
	kobs = make_kernel_data(xobs, kernel, sigmasq_obs=sigmasq_obs, sigmasq_cor=sigmasq_cor, rho=rho)
  	kobs_pred = make_kernel_obs_pred(xobs, xpred, kernel, sigmasq_cor=sigmasq_cor, rho=rho)
  	alpha = kobs \ yobs
	pred_mean = kobs_pred' * alpha
end

function predict_deriv(xobs::AA, yobs::AA, xpred::AA;
			kernel::Function = matern52_sparse_kernel, dkerneldx::Function = dkerneldx_matern52_sparse,
			sigmasq_obs::AA = 1e-16*ones(length(xobs)),	sigmasq_cor::Real=1.0, rho::Real=1.0)  where { T<:Real, AA<:AbstractArray{T,1} }
  kobs = make_kernel_data(xobs, kernel=kernel, sigmasq_obs=sigmasq_obs, sigmasq_cor=sigmasq_cor, rho=rho)
  kobs_pred_deriv = make_kernel_obs_pred(xobs,xpred, kernel=dkerneldx, sigmasq_cor=sigmasq_cor, rho=rho)
  alpha = kobs \ yobs
  pred_deriv = kobs_pred_deriv' * alpha
end

function predict_deriv2(xobs::AA, yobs::AA, xpred::AA;
			kernel::Function = matern52_sparse_kernel, d2kerneldx2::Function = d2kerneldx2_matern52_sparse,
			sigmasq_obs::AA = 1e-16*ones(length(xobs)),	sigmasq_cor::Real=1.0, rho::Real=1.0)  where { T<:Real, AA<:AbstractArray{T,1} }
  kobs = make_kernel_data(xobs, kernel=kernel, sigmasq_obs=sigmasq_obs, sigmasq_cor=sigmasq_cor, rho=rho)
  kobs_pred_deriv2 = make_kernel_obs_pred(xobs,xpred, kernel=d2kerneldx2, sigmasq_cor=sigmasq_cor, rho=rho)
  alpha = kobs \ yobs
  pred_deriv = kobs_pred_deriv2' * alpha
end

function predict_mean_and_deriv(xobs::AA, yobs::AA, xpred::AA;
			kernel::Function = matern52_sparse_kernel, dkerneldx::Function = dkerneldx_matern52_sparse,
			sigmasq_obs::AA = 1e-16*ones(length(xobs)),	sigmasq_cor::Real=1.0, rho::Real=1.0)  where { T<:Real, AA<:AbstractArray{T,1} }
  kobs = make_kernel_data(xobs, kernel=kernel, sigmasq_obs=sigmasq_obs, sigmasq_cor=sigmasq_cor, rho=rho)
  alpha = kobs \ yobs
  kobs_pred = make_kernel_obs_pred(xobs,xpred, kernel=kernel, sigmasq_cor=sigmasq_cor, rho=rho)
  pred_mean = kobs_pred' * alpha
  kobs_pred_deriv = make_kernel_obs_pred(xobs,xpred, kernel=dkerneldx, sigmasq_cor=sigmasq_cor, rho=rho)
  pred_deriv = kobs_pred_deriv' * alpha
  return (mean=pred_mean, deriv=pred_deriv)
end

function predict_mean_and_derivs(xobs::AA, yobs::AA, xpred::AA;
			kernel::Function = matern52_sparse_kernel, dkerneldx::Function = dkerneldx_matern52_sparse, d2kerneldx2::Function = d2kerneldx2_matern52_sparse,
			sigmasq_obs::AA = 1e-16*ones(length(xobs)),	sigmasq_cor::Real=1.0, rho::Real=1.0)  where { T<:Real, AA<:AbstractArray{T,1} }
  kobs = make_kernel_data(xobs, kernel=kernel, sigmasq_obs=sigmasq_obs, sigmasq_cor=sigmasq_cor, rho=rho)
  alpha = kobs \ yobs
  kobs_pred = make_kernel_obs_pred(xobs,xpred, kernel=kernel, sigmasq_cor=sigmasq_cor, rho=rho)
  pred_mean = kobs_pred' * alpha
  kobs_pred_deriv = make_kernel_obs_pred(xobs,xpred, kernel=dkerneldx, sigmasq_cor=sigmasq_cor, rho=rho)
  pred_deriv = kobs_pred_deriv' * alpha
  kobs_pred_deriv2 = make_kernel_obs_pred(xobs,xpred, kernel=d2kerneldx2, sigmasq_cor=sigmasq_cor, rho=rho)
  pred_deriv2 = kobs_pred_deriv2' * alpha
  return (mean=pred_mean, deriv=pred_deriv, deriv2=pred_deriv2)
end

function gp_marginal(xobs::AA, yobs::AA, kernel::Function;
			sigmasq_obs::AA = 1e-16*ones(length(xobs)),	sigmasq_cor::Real=1.0, rho::Real=1.0)  where { T<:Real, AA<:AbstractArray{T,1} }
  kobs = make_kernel_data(xobs, kernel=kernel, sigmasq_obs=sigmasq_obs, sigmasq_cor=sigmasq_cor, rho=rho)
  -0.5*( invquad(kobs, yobs) + logdet(kobs) + length(xobs)*log(2pi) )
end

#=
function calc_gp_marginal_on_segments(lambda::AA, flux::AA;
                                sigmasq_obs::AA = 1e-16*ones(length(lambda)),	sigmasq_cor::Real=1.0, rho::Real=1.0,
                                half_chunk_size::Integer = 100)  where { T<:Real, AA<:AbstractArray{T,1} }
  @assert length(lambda) == length(flux) == length(sigmasq_obs)
  println("# sigmasq_obs[1,1] = ", sigmasq_obs[1,1], " sigmasq_cor= ", sigmasq_cor, " rho= ", rho)
  output = 0.0
  num_seg = convert(Int64,ceil(length(lambda)/half_chunk_size)-1)
  for i in 1:num_seg
    idx_begin = 1+half_chunk_size*(i-1)
    idx_end = min(half_chunk_size*(i+1), length(lambda))
    write_idx_begin = idx_begin + div(half_chunk_size,2)
    write_idx_end = idx_end - div(half_chunk_size,2)
    if i==1 write_idx_begin=1 end
    if i==num_seg
      idx_begin = max(1,idx_end-2*half_chunk_size)
      write_idx_end=length(lambda)
    end
    #println("# i= ",i,": ", idx_begin, " - ", idx_end, " -> ", write_idx_begin, " - ", write_idx_end)
    #output[write_idx_begin:write_idx_end] = predict_gp(view(lambda,idx_begin:idx_end), view(flux,idx_begin:idx_end), view(lambda,write_idx_begin:write_idx_end), sigmasq_obs=view(sigmasq_obs,idx_begin:idx_end), sigmasq_cor=sigmasq_cor, rho=rho)
    output += gp_marginal(lambda[idx_begin:idx_end], flux[idx_begin:idx_end], sigmasq_obs=sigmasq_obs[idx_begin:idx_end], sigmasq_cor=sigmasq_cor, rho=rho)
  end
  output
end

function gp_marginal_wrapper(param::Vector)
   @assert length(param) == 2
   idx_min = 40000
   idx_max = 45000
   println("# wrapper: ", param)
   calc_gp_marginal_on_segments(lambda[idx_min:idx_max],vec(mean(obs[idx_min:idx_max],2)), sigmasq_obs = ones(length(lambda[idx_min:idx_max]))/150000,
                    sigmasq_cor = exp(param[1]), rho = exp(param[2]) )
end

function calc_doppler_component_gp(lambda::AA, flux::AA;
                                sigmasq_obs::AA = 1e-16*ones(length(lambda)),	sigmasq_cor::Real=1.0, rho::Real=1.0,
                                half_chunk_size::Integer = 100)  where { T<:Real, AA<:AbstractArray{T,1} }
   lambda.*calc_gp_on_segments(predict_deriv,lambda, flux, sigmasq_obs=sigmasq_obs,	sigmasq_cor=sigmasq_cor, rho=rho, half_chunk_size=half_chunk_size)
end

function calc_gp_on_segments(predict_gp::Function, lambda::AA, flux::AA;
                                sigmasq_obs::AA = 1e-16*ones(length(lambda)),	sigmasq_cor::Real=1.0, rho::Real=1.0,
                                half_chunk_size::Integer = 100)  where { T<:Real, AA<:AbstractArray{T,1} }
  @assert length(lambda) == length(flux) == length(sigmasq_obs)
  output = Array{eltype(flux)}(undef,length(lambda))
  num_seg = convert(Int64,ceil(length(lambda)/half_chunk_size)-1)
  for i in 1:num_seg
    idx_begin = 1+half_chunk_size*(i-1)
    idx_end = min(half_chunk_size*(i+1), length(lambda))
    write_idx_begin = idx_begin + div(half_chunk_size,2)
    write_idx_end = idx_end - div(half_chunk_size,2)
    if i==1 write_idx_begin=1 end
    if i==num_seg
      idx_begin = max(1,idx_end-2*half_chunk_size)
      write_idx_end=length(lambda)
    end
    #println("# i= ",i,": ", idx_begin, " - ", idx_end, " -> ", write_idx_begin, " - ", write_idx_end)
    #output[write_idx_begin:write_idx_end] = predict_gp(view(lambda,idx_begin:idx_end), view(flux,idx_begin:idx_end), view(lambda,write_idx_begin:write_idx_end), sigmasq_obs=view(sigmasq_obs,idx_begin:idx_end), sigmasq_cor=sigmasq_cor, rho=rho)
    output[write_idx_begin:write_idx_end] = predict_gp(lambda[idx_begin:idx_end], flux[idx_begin:idx_end], lambda[write_idx_begin:write_idx_end], sigmasq_obs=sigmasq_obs[idx_begin:idx_end], sigmasq_cor=sigmasq_cor, rho=rho)
  end
  output
end

=#

#=
function calc_doppler_component_gp(lambda::AA, flux::AbstractArray{T,2};
                                sigmasq_obs::AA = 1e-16*ones(length(lambda)),	sigmasq_cor::Real=1.0, rho::Real=1.0)  where { T<:Real, AA<:AbstractArray{T,1} }
  doppler_basis = calc_doppler_component_gp(lambda,vec(mean(flux,2)),sigmasq_obs=sigmasq_obs,sigmasq_cor=sigmasq_cor,rho=rho)
end

function calc_doppler_quadratic_term_gp(lambda::AA, flux::AA;
                                sigmasq_obs::AA = 1e-16*ones(length(lambda)),	sigmasq_cor::Real=1.0, rho::Real=1.0,
                                half_chunk_size::Integer = 100)  where { T<:Real, AA<:AbstractArray{T,1} }
   0.5*lambda.^2.*calc_gp_on_segments(predict_deriv2,lambda, flux, sigmasq_obs=sigmasq_obs,	sigmasq_cor=sigmasq_cor, rho=rho, half_chunk_size=half_chunk_size)
end

function calc_doppler_quadratic_term_gp(lambda::AA, flux::AbstractArray{T,2};
                                     sigmasq_obs::AA = 1e-16*ones(length(lambda)),	sigmasq_cor::Real=1.0, rho::Real=1.0)  where rho*bandwidth{ T<:Real, AA<:AbstractArray{T,1} }
  doppler_quad_term = calc_doppler_quadratic_term_gp(lambda,vec(mean(flux,2)),sigmasq_obs=sigmasq_obs,sigmasq_cor=sigmasq_cor,rho=rho)
end
=#

#=
function test_gp_deriv()
  #lambda_small = collect(linspace(1000.0,1100.0,101))
  #yobs_small = sin(2pi*lambda/40)+ 0.3*randn(length(lambda_small))
  #lambda_pred = collect(linspace(1000.0,1100.0,201))
  lambda_small = lambda[10000:10100]
  yobs_small = obs[10000:10100,1]
  lambda_pred = lambda_small
  (f,df,df2) = predict_mean_and_derivs(lambda_small,yobs_small,lambda_pred, sigmasq_obs=1.0/150000.0*ones(length(lambda_small)), sigmasq_cor=0.1,rho = 0.01)
  plot(lambda_small,(yobs_small),"r.")
  plot(lambda_pred,f,"b-")
  plot(lambda_pred,df,"g-")
  #plot(lambda_pred,df2,"m-")
end
=#


text:
function initialize_dynamic_device!(
    dynamic_device::DynamicWrapper{DynG},
    static::PSY.StaticInjection,
    initial_inner_vars::AbstractVector,
) where {DynG <: PSY.DynamicGenerator}
    #Obtain States
    device_states = zeros(PSY.get_n_states(dynamic_device))

    #Initialize Machine and Shaft: δ and ω
    initialize_mach_shaft!(device_states, static, dynamic_device, initial_inner_vars)
    #Initialize extra Shaft states
    initialize_shaft!(device_states, static, dynamic_device, initial_inner_vars)
    #Initialize AVR
    initialize_avr!(device_states, static, dynamic_device, initial_inner_vars)
    #Initialize TG
    initialize_tg!(device_states, static, dynamic_device, initial_inner_vars)
    #Initialize PSS
    initialize_pss!(device_states, static, dynamic_device, initial_inner_vars)

    return device_states
end

function initialize_dynamic_device!(
    dynamic_device::DynamicWrapper{DynI},
    static::PSY.StaticInjection,
    initial_inner_vars::AbstractVector,
) where {DynI <: PSY.DynamicInverter}
    #Obtain States
    device_states = zeros(PSY.get_n_states(dynamic_device))

    #Initialize Machine and Shaft: V and I
    initialize_filter!(device_states, static, dynamic_device, initial_inner_vars)
    #Initialize Converter
    initialize_converter!(device_states, static, dynamic_device, initial_inner_vars)
    #Initialize freq estimator
    initialize_frequency_estimator!(
        device_states,
        static,
        dynamic_device,
        initial_inner_vars,
    )
    #Initialize OuterLoop
    initialize_outer!(device_states, static, dynamic_device, initial_inner_vars)
    #Initialize DCside
    initialize_DCside!(device_states, static, dynamic_device, initial_inner_vars)
    #Initialize Converter
    initialize_converter!(device_states, static, dynamic_device, initial_inner_vars)
    #Initialize InnerLoop
    initialize_inner!(device_states, static, dynamic_device, initial_inner_vars)
    return device_states
end

function initialize_static_device!(
    ::StaticWrapper{PSY.PowerLoad, T},
) where {T <: LoadCategory}
    return
end

function initialize_static_device!(
    ::StaticWrapper{PSY.FixedAdmittance, T},
) where {T <: LoadCategory}
    return
end

function initialize_static_device!(
    device::StaticWrapper{PSY.Source, T},
) where {T <: BusCategory}
    #PowerFlow Data
    P0 = PSY.get_active_power(device)
    Q0 = PSY.get_reactive_power(device)
    Vm = PSY.get_magnitude(PSY.get_bus(device))
    θ = PSY.get_angle(PSY.get_bus(device))
    S0 = P0 + Q0 * 1im
    V_R = Vm * cos(θ)
    V_I = Vm * sin(θ)
    V = V_R + V_I * 1im
    I = conj(S0 / V)
    I_R = real(I)
    I_I = imag(I)
    R_th = PSY.get_R_th(device.device)
    X_th = PSY.get_X_th(device.device)
    Zmag = R_th^2 + X_th^2

    function f!(out, x)
        V_R_internal = x[1]
        V_I_internal = x[2]

        out[1] =
            R_th * (V_R_internal - V_R) / Zmag + X_th * (V_I_internal - V_I) / Zmag - I_R
        out[2] =
            R_th * (V_I_internal - V_I) / Zmag - X_th * (V_R_internal - V_R) / Zmag - I_I
    end
    x0 = [V_R, V_I]
    sol = NLsolve.nlsolve(f!, x0)
    if !NLsolve.converged(sol)
        @warn("Initialization in Source failed")
    else
        sol_x0 = sol.zero
        #Update terminal voltages
        V_internal = sqrt(sol_x0[1]^2 + sol_x0[2]^2)
        θ_internal = atan(sol_x0[2] / sol_x0[1])
        PSY.set_internal_voltage!(device.device, V_internal)
        PSY.set_internal_angle!(device.device, θ_internal)
        set_V_ref(device, PSY.get_internal_voltage(device.device))
        set_θ_ref(device, PSY.get_internal_angle(device.device))
    end
    return
end

function initialize_dynamic_device!(
    dynamic_device::DynamicWrapper{PSY.PeriodicVariableSource},
    source::PSY.Source,
    ::AbstractVector,
)
    device_states = zeros(PSY.get_n_states(dynamic_device))

    #PowerFlow Data
    P0 = PSY.get_active_power(source)
    Q0 = PSY.get_reactive_power(source)
    Vm = PSY.get_magnitude(PSY.get_bus(source))
    θ = PSY.get_angle(PSY.get_bus(source))
    S0 = P0 + Q0 * 1im
    V_R = Vm * cos(θ)
    V_I = Vm * sin(θ)
    V = V_R + V_I * 1im
    I = conj(S0 / V)
    I_R = real(I)
    I_I = imag(I)
    R_th = PSY.get_R_th(source)
    X_th = PSY.get_X_th(source)
    Zmag = R_th^2 + X_th^2
    function f!(out, x)
        V_R_internal = x[1]
        V_I_internal = x[2]

        out[1] =
            R_th * (V_R_internal - V_R) / Zmag + X_th * (V_I_internal - V_I) / Zmag - I_R
        out[2] =
            R_th * (V_I_internal - V_I) / Zmag - X_th * (V_R_internal - V_R) / Zmag - I_I
    end
    x0 = [V_R, V_I]
    sol = NLsolve.nlsolve(f!, x0)
    if !NLsolve.converged(sol)
        @warn("Initialization in Periodic Variable Source failed")
    else
        sol_x0 = sol.zero
        #Update terminal voltages
        V_internal = sqrt(sol_x0[1]^2 + sol_x0[2]^2)
        θ_internal = atan(sol_x0[2] / sol_x0[1])

        V_internal_freqs = 0.0
        V_freqs = PSY.get_internal_voltage_frequencies(get_device(dynamic_device))
        V_coeff = PSY.get_internal_voltage_coefficients(get_device(dynamic_device))
        for (ix, ω) in enumerate(V_freqs)
            V_internal_freqs += V_coeff[ix][2]     #sin(0) = 0; cos(0)=1
        end

        θ_internal_freqs = 0.0
        θ_freqs = PSY.get_internal_angle_frequencies(get_device(dynamic_device))
        θ_coeff = PSY.get_internal_angle_coefficients(get_device(dynamic_device))
        for (ix, ω) in enumerate(θ_freqs)
            θ_internal_freqs += θ_coeff[ix][2]     #sin(0) = 0; cos(0)=1
        end
        V_internal_bias = V_internal - V_internal_freqs
        θ_internal_bias = θ_internal - θ_internal_freqs

        device_states[1] = V_internal
        device_states[2] = θ_internal
        PSY.set_internal_voltage_bias!(get_device(dynamic_device), V_internal_bias)
        PSY.set_internal_angle_bias!(get_device(dynamic_device), θ_internal_bias)
    end
    return device_states
end

function initialize_dynamic_device!(branch::BranchWrapper)
    device_states = zeros(PSY.get_n_states(branch))
    #PowerFlow Data
    arc = PSY.get_arc(branch)
    Vm_from = PSY.get_magnitude(PSY.get_from(arc))
    θ_from = PSY.get_angle(PSY.get_from(arc))
    Vm_to = PSY.get_magnitude(PSY.get_to(arc))
    θ_to = PSY.get_angle(PSY.get_to(arc))
    V_R_from = Vm_from * cos(θ_from)
    V_I_from = Vm_from * sin(θ_from)
    V_R_to = Vm_to * cos(θ_to)
    V_I_to = Vm_to * sin(θ_to)
    R = PSY.get_r(branch)
    X = PSY.get_x(branch)
    Zmag = R^2 + X^2
    #Compute Current
    I_R = R * (V_R_from - V_R_to) / Zmag + X * (V_I_from - V_I_to) / Zmag
    I_I = R * (V_I_from - V_I_to) / Zmag - X * (V_R_from - V_R_to) / Zmag
    #Update Current
    device_states[1] = I_R
    device_states[2] = I_I
    return device_states
end


text:
using Documenter, ContMechTensors

makedocs(
    modules = [ContMechTensors],
    format = :html,
    sitename = "ContMechTensors.jl",
    doctest = true,
    strict = true,
    pages = Any[
        "Home" => "index.md",
        "Manual" => [
            "man/constructing_tensors.md",
            "man/indexing.md",
            "man/binary_operators.md",
            "man/other_operators.md",
            "man/storing_tensors.md",
            "man/automatic_differentiation.md",
            ],
        "Demos" => "demos.md"
    ]
)

deploydocs(
    repo = "github.com/KristofferC/ContMechTensors.jl.git",
    target = "build",
    julia = "0.5",
    deps = nothing,
    make = nothing
)


text:
# # Airplane
#
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/models/Airplane.ipynb)
#
# This example consists of a dynamical system that is a simple model of a flying
# airplane. There are 12 state variables:
#
# ```math
# \begin{aligned}
# [x, y, z, u, v, w, ϕ, θ, ψ, r, p, q]
# \end{aligned}
# ```
#
# where $(x, y, z)$ is the position of the center of gravity (C.G.),
# $(u, v, w)$ are the components of velocity in $(x, y, z)$ directions,
# $(p, q, r)$ are body rotation rates, and $(ϕ, θ, ψ)$ are the Euler angles.
#
# See also the file `airplane.png` for a visualization of the meaning of the
# coordinates with respect to the airplane.

module Airplane  #jl

using NeuralNetworkAnalysis
using NeuralNetworkAnalysis: SingleEntryVector

# The following option determines whether the falsification settings should be
# used or not. The falsification settings are sufficient to show that the safety
# property is violated. Concretely we start from an initial point and use a
# smaller time step.
const falsification = true;

# ## Model
#
# The equations of motion are reduced to:
#
# ```math
# \left\{ \begin{array}{lcl}
# u̇ &=& − g \sin(\theta) + \frac{F_x}{m} - qw + rv \\
# v̇ &=& g \cos(\theta) \sin(\phi) + \frac{F_y}{m} - ru + pw \\
# ẇ &=& g \cos(\theta) \cos(\phi) + \frac{F_z}{m} - pv + qu \\
# I_x \dot p + I_{xz} \dot r &=& M_x - (I_z - I_y) qr - I_{xz} pq \\
# I_y \dot q &=& M_y - I_{xz}(r^2 - p^2) - (I_x - I_z) pr \\
# I_{xz} \dot p + I_z \dot r &=& M_z - (I_y - I_x) qp - I_{xz} rq
# \end{array} \right.
# ```
#
# The mass of the airplane is denoted with m and $I_x$ , $I_y$ , $I_z$
# and $I_{xz}$ are the moment of inertia with respect to the indicated axis.
# The control parameters include three force components $F_x$ , $F_y$ and
# $F_z$ and three moment components $M_x$ , $M_y$ , $M_z$ . Notice that
# for simplicity we have assumed the aerodynamic forces are absorbed in the
# $F$'s. In addition to these six equations, we have six additional kinematic
# equations:
#
# ```math
# \begin{bmatrix}
# \dot x \\ \dot y \\ \dot z
# \end{bmatrix}
# =
# \begin{bmatrix}
# \cos(\psi) & -\sin(\psi) & 0 \\
# \sin(\psi) & \cos(\psi) & 0 \\
# 0 & 0 & 1
# \end{bmatrix}
# \begin{bmatrix}
# \cos(\theta) & 0 & \sin(\theta) \\
# 0 & 1 & 0 \\
# -\sin(\theta) & 0 & \cos(\theta)
# \end{bmatrix}
# \begin{bmatrix}
# 1 & 0 & 0 \\
# 0 & \cos(\phi) & -\sin(\phi) \\
# 0 & \sin(\phi) & \cos(\phi)
# \end{bmatrix}
# ```
#
# and
#
# ```math
# \begin{bmatrix}
# \phi \\ \theta \\ \psi
# \end{bmatrix}
# =
# \begin{bmatrix}
# 1 & \tan(\theta) \sin(\phi) & \tan(\theta) \cos(\phi) \\
# 0 & \cos(\phi) & -\sin(\phi) \\
# 0 & \sec(\theta) sin(\phi) & \sec(\theta) \cos(\phi)
# \end{bmatrix}
# \begin{bmatrix}
# p \\ q \\ r
# \end{bmatrix}
# ```
#
# For the simplicity of control design, the parameters have been chosen to
# have some nominal dimensionless values:
# $m = 1$, $I_x = I_y = I_z = 1$, $I_{xz} = 0$ and $g = 1$.

Tψ = ψ -> [cos(ψ)  -sin(ψ)  0;
           sin(ψ)   cos(ψ)  0;
                0        0  1]

Tθ = θ -> [ cos(θ)  0  sin(θ);
                 0  1       0;
           -sin(θ)  0  cos(θ)]

Tϕ = ϕ -> [1       0        0;
           0  cos(ϕ)  -sin(ϕ);
           0  sin(ϕ)   cos(ϕ)]

Rϕθ = (ϕ, θ) -> [1  tan(θ) * sin(ϕ)  tan(θ) * cos(ϕ);
                 0           cos(θ)          -sin(ϕ);
                 0  sec(θ) * sin(ϕ)  sec(θ) * cos(ϕ)]

## alternative matrix with only sin/cos but with postprocessing
## Rϕθ_ = (ϕ, θ) -> [cos(θ)  sin(θ) * sin(ϕ)   sin(θ) * cos(ϕ);
##                        0  cos(θ) * cos(ϕ)  -cos(θ) * sin(ϕ);
##                        0           sin(ϕ)            cos(ϕ)]

## model constants
const m = 1.0
const g = 1.0

## unused constants (terms are simplified instead)
const Ix = 1.0
const Iy = 1.0
const Iz = 1.0
const Ixz = 0.0

@taylorize function airplane!(dx, x, p, t)
    _x, y, z, u, v, w, ϕ, θ, ψ, r, _p, q, Fx, Fy, Fz, Mx, My, Mz = x

    T_ψ = Tψ(ψ)
    T_θ = Tθ(θ)
    T_ϕ = Tϕ(ϕ)
    mat_1 = T_ψ * T_θ * T_ϕ
    xyz = mat_1 * vcat(u, v, w)

    mat_2 = Rϕθ(ϕ, θ)
    ## mat_2 = 1 / cos(θ) * Rϕθ_(ϕ, θ)  # alternative matrix with postprocessing
    ϕθψ = mat_2 * vcat(_p, q, r)

    dx[1] = xyz[1]
    dx[2] = xyz[2]
    dx[3] = xyz[3]
    dx[4] = -g * sin(θ) + Fx / m - q * w + r * v
    dx[5] = g * cos(θ) * sin(ϕ) + Fy / m - r * u + _p * w
    dx[6] = g * cos(θ) * cos(ϕ) + Fz / m - _p * v + q * u
    dx[7] = ϕθψ[1]
    dx[8] = ϕθψ[2]
    dx[9] = ϕθψ[3]
    dx[10] = Mx  # simplified term
    dx[11] = My  # simplified term
    dx[12] = Mz  # simplified term
    dx[13] = zero(Fx)
    dx[14] = zero(Fy)
    dx[15] = zero(Fz)
    dx[16] = zero(Mx)
    dx[17] = zero(My)
    dx[18] = zero(Mz)
end;

# ## Specifications

controller = read_nnet(@modelpath("Airplane", "controller_airplane.nnet"))

X₀ = Hyperrectangle(low=[0.0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
                    high=[0.0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0])
if falsification
    ## choose a single point in the initial states (here: the top-most one)
    X₀ = Singleton(high(X₀))
end
U₀ = ZeroSet(6)

vars_idx = Dict(:state_vars=>1:12, :control_vars=>13:18)
ivp = @ivp(x' = airplane!(x), dim: 18, x(0) ∈ X₀ × U₀)

period = 0.1  # control period

prob = ControlledPlant(ivp, controller, vars_idx, period);

# Safety specification: $x_2$ should be in $[−0.5, 0.5]$ and $x_7, x_8, x_9$
# should be in $[-1, 1]$.
if falsification
    k = 4  # falsification can run for a shorter time horizon
else
    k = 20
end
T = k * period  # time horizon

safe_states = HPolyhedron([HalfSpace(SingleEntryVector(2, 18, 1.0), 0.5),
                           HalfSpace(SingleEntryVector(2, 18, -1.0), 0.5),
                           HalfSpace(SingleEntryVector(7, 18, 1.0), 1.0),
                           HalfSpace(SingleEntryVector(7, 18, -1.0), 1.0),
                           HalfSpace(SingleEntryVector(8, 18, 1.0), 1.0),
                           HalfSpace(SingleEntryVector(8, 18, -1.0), 1.0),
                           HalfSpace(SingleEntryVector(9, 18, 1.0), 1.0),
                           HalfSpace(SingleEntryVector(9, 18, -1.0), 1.0)])

## property for guaranteed violation
predicate = X -> isdisjoint(overapproximate(X, Hyperrectangle), safe_states)
predicate_sol = sol -> any(predicate(R) for F in sol for R in F);

# ## Results

import DifferentialEquations

alg = TMJets(abstol=1e-10, orderT=7, orderQ=1)
alg_nn = Ai2()

function benchmark(; silent::Bool=false)
    ## We solve the controlled system:
    silent || println("flowpipe construction")
    res_sol = @timed solve(prob, T=T, alg_nn=alg_nn, alg=alg)
    sol = res_sol.value
    silent || print_timed(res_sol)

    ## Next we check the property for an overapproximated flowpipe:
    silent || println("property checking")
    solz = overapproximate(sol, Zonotope)
    res_pred = @timed predicate_sol(solz)
    silent || print_timed(res_pred)
    if res_pred.value
        silent || println("The property is violated.")
    else
        silent || println("The property may be satisfied.")
    end

    ## We also compute some simulations:
    silent || println("simulation")
    trajectories = falsification ? 1 : 50
    res_sim = @timed simulate(prob, T=T, trajectories=trajectories,
                              include_vertices=!falsification)
    sim = res_sim.value
    silent || print_timed(res_sim)

    return solz, sim
end

benchmark(silent=true)  # warm-up
res = @timed benchmark()  # benchmark
sol, sim = res.value
println("total analysis time")
print_timed(res);
print(Main.io, "JuliaReach, Airplane, -, falsified, $(res.time * 1e9)\n")

# Finally we plot the results:

using Plots

## set more precise tolerance for plotting small sets correctly
LazySets.set_ztol(Float64, 1e-9)

function plot_helper(fig, vars)
    if vars[1] == 0
        safe_states_projected = project(safe_states, [vars[2]])
        time = Interval(0, T)
        safe_states_projected = cartesian_product(time, safe_states_projected)
    else
        safe_states_projected = project(safe_states, vars)
    end
    plot!(fig, safe_states_projected, color=:lightgreen, lab="safe states")
    if !falsification && 0 ∉ vars
        plot!(fig, project(initial_state(prob), vars), lab="X₀")
    end
    plot!(fig, sol, vars=vars, color=:yellow, lab="")
    lab_sim = falsification ? "simulation" : ""
    plot_simulation!(fig, sim; vars=vars, color=:black, lab=lab_sim)
end

vars = (2, 7)
fig = plot(xlab="y", ylab="ϕ", leg=:bottomleft)
plot_helper(fig, vars)
if falsification
    xlims!(-0.01, 0.65)
    ylims!(0.9, 1.01)
else
    xlims!(-1.8, 22.5)
    ylims!(-1.05, 1.05)
end
savefig("Airplane-x2-x7.png")

#-

vars = (8, 9)
fig = plot(xlab="θ", ylab="ψ", leg=:bottom)
plot_helper(fig, vars)
if falsification
    xlims!(0.999, 1.03)
    ylims!(0.99, 1.001)
else
    xlims!(-1.05, 1.2)
    ylims!(-1.05, 1.2)
end
savefig("Airplane-x8-x9.png")

#-

end  #jl
nothing  #jl


text:
"""
DioOption represents a set of options that is used in DIO, a fundamental
frequency analysis.

**Fields**

$(FIELDS)
"""
immutable DioOption
    f0floor::Cdouble
    f0ceil::Cdouble
    channels_in_octave::Cdouble
    "frame period in ms"
    period::Cdouble
    speed::Cint
    "added in v0.2.1-2 (WORLD 0.2.0_2)"
    allowed_range::Cdouble

    function DioOption(f0floor, f0ceil, channels_in_octave, period, speed,
                       allowed_range = 0.1)
        f0floor > f0ceil && throw(ArgumentError("F0floor must be larger than F0ceil"))
        f0floor < 0 && throw(ArgumentError("f0floor must be positive"))
        if channels_in_octave < 0
            throw(ArgumentError("channels_in_octave must be positive"))
        end
        period <= 0 && throw(ArgumentError("period must be positive"))
        (1 <= speed <= 12) || throw(ArgumentError("1 ≤ speed ≤ 12 is supprted"))
        allowed_range >= 0 || throw(ArgumentError("allowed_range >= 0 is supported"))
        new(f0floor, f0ceil, channels_in_octave, period, speed, allowed_range)
    end
end

"""
CheapTrick options

**Fields**

$(FIELDS)
"""
immutable CheapTrickOption
    q1::Cdouble
    CheapTrickOption(q1=-0.09) = new(q1)
end

"""
D4C options (nothing for now, but for future changes)

**Fields**

$(FIELDS)
"""
immutable D4COption
    dummy::Cdouble
    D4COption(dummy=0.0) = new(dummy)
end

# Note that the default options assume that the sampling frequency of a input
# speech signal is 44.1 kHz.
function DioOption(;
                   f0floor::Float64=71.0,
                   f0ceil::Float64=800.0,
                   channels_in_octave::Float64=2.0,
                   period::Float64=5.0,
                   speed::Integer=1,
                   allowed_range::Float64=0.1
                   )
    DioOption(f0floor, f0ceil, channels_in_octave, period, speed, allowed_range)
end

function get_samples_for_dio(fs::Real, len::Integer, period::Real)
    ccall((:GetSamplesForDIO, libworld), Cint,
          (Cint, Cint, Cdouble), fs, len, period)
end

"""
$(SIGNATURES)

Dio estimates F0 trajectory given a monoral input signal.

**Paremters**

- `x`  : Input signal
- `fs` : Sampling frequency
- `opt` : DioOption

**Returns**

- `time_axis`  : Temporal positions.
- `f0`         : F0 contour.

"""
function dio(x::StridedVector{Cdouble}, fs::Real, opt::DioOption=DioOption())
    expectedlen = get_samples_for_dio(fs, length(x), opt.period)
    f0 = Array(Cdouble, expectedlen)
    timeaxis = Array(Cdouble, expectedlen)
    # Note that value passinig of julia-type to C-struct doesn't work.
    ccall((:DioByOptPtr, libworld),  Void,
          (Ptr{Cdouble}, Cint, Cint, Ptr{DioOption}, Ptr{Cdouble}, Ptr{Cdouble}),
          x, length(x), fs, &opt, timeaxis, f0)
    f0, timeaxis
end

"""
$(SIGNATURES)

StoneMask refines the estimated F0 by Dio,

**Parameters**

- `x` : Input signal
- `fs` : Sampling frequency
- `time_axis` : Temporal information
- `f0` : f0 contour

**Returns**

- `refined_f0` : Refined F0
"""
function stonemask(x::StridedVector{Cdouble}, fs::Integer,
                   timeaxis::StridedVector{Cdouble},
                   f0::StridedVector{Cdouble})
    refinedF0 = Array(Cdouble, length(f0))
    ccall((:StoneMask, libworld),  Void,
          (Ptr{Cdouble}, Cint, Cint, Ptr{Cdouble}, Ptr{Cdouble}, Cint,
           Ptr{Cdouble}),
          x, length(x), fs, timeaxis, f0, length(f0), refinedF0)
    refinedF0
end

# Array{T,2} -> Array{Ptr{T}}
function ptrarray2d!{T<:Real}(dst::Array{Ptr{T},1}, src::Array{T,2})
    for i=1:size(src, 2)
        @inbounds dst[i] = pointer(view(src, 1:size(src, 1), i), 1)
    end
end

"""
$(SIGNATURES)

GetFFTSizeForCheapTrick calculates the FFT size based on the sampling
frequency and the lower limit of f0 (It is defined in world.h).

**Parameters**

- `fs`: Sampling frequency

**Returns**

- `fftsize` : FFT size
"""
function get_fftsize_for_cheaptrick(fs::Integer)
    fftsize = ccall((:GetFFTSizeForCheapTrick, libworld), Cint, (Cint,), fs)
    convert(Int, fftsize)
end

"""
$(SIGNATURES)

CheapTrick calculates the spectrogram that consists of spectral envelopes
estimated by CheapTrick.

**Parameters**

- `x` : Input signal
- `fs` : Sampling frequency
- `time_axis` : Time axis
- `f0` : F0 contour
- `opt` : CheapTrick option

**Returns**

- `spectrogram`  : Spectrogram estimated by CheapTrick.
"""
function cheaptrick(x::StridedVector{Cdouble}, fs::Integer,
                    timeaxis::StridedVector{Cdouble},
                    f0::StridedVector{Cdouble};
                    opt::CheapTrickOption=CheapTrickOption()
    )
    freqbins = get_fftsize_for_cheaptrick(fs)>>1 + 1
    spectrogram = Array(Cdouble, freqbins, length(f0))

    # Array{Cdouble,2} -> Array{Ptr{Cdouble}}
    cspectrogram = Array(Ptr{Cdouble}, size(spectrogram, 2))
    ptrarray2d!(cspectrogram, spectrogram)

    ccall((:CheapTrick, libworld), Void,
          (Ptr{Cdouble}, Cint, Cint, Ptr{Cdouble}, Ptr{Cdouble}, Cint,
           Ptr{CheapTrickOption}, Ptr{Ptr{Cdouble}}),
          x, length(x), fs, timeaxis, f0, length(f0), &opt, cspectrogram)

    # Array{Cdouble,2} <- Array{Ptr{Cdouble}}
    for i=1:length(f0), j=1:freqbins
        @inbounds spectrogram[j,i] = unsafe_load(cspectrogram[i], j)
    end

    spectrogram
end

"""
$(SIGNATURES)

D4C calculates the aperiodicity estimated by D4C.

**Parameters**

- `x` : Input signal
- `fs` : Sampling frequency
- `time_axis` : Time axis
- `f0` : F0 contour

**Returns**

- `aperiodicity` : Aperiodicity estimated by D4C.
"""
function d4c(x::StridedVector{Cdouble}, fs::Integer,
             timeaxis::StridedVector{Cdouble},
             f0::StridedVector{Cdouble};
             opt::D4COption=D4COption())
    fftsize = get_fftsize_for_cheaptrick(fs)
    freqbins = fftsize>>1 + 1
    aperiodicity = zeros(Cdouble, freqbins, length(f0))

    # Array{Cdouble,2} -> Array{Ptr{Cdouble}}
    caperiodicity = Array(Ptr{Cdouble}, size(aperiodicity, 2))
    ptrarray2d!(caperiodicity, aperiodicity)

    ccall((:D4C, libworld), Void,
          (Ptr{Cdouble}, Cint, Cint, Ptr{Cdouble}, Ptr{Cdouble}, Cint, Cint,
           Ptr{D4COption}, Ptr{Ptr{Cdouble}}),
          x, length(x), fs, timeaxis, f0, length(f0), fftsize, &opt,
          caperiodicity)

    # Array{Cdouble,2} <- Array{Ptr{Cdouble}}
    for i=1:length(f0), j=1:freqbins
        aperiodicity[j,i] = unsafe_load(caperiodicity[i], j)
    end

    aperiodicity
end

"""
$(SIGNATURES)

Synthesis synthesize the voice based on f0, spectrogram and
aperiodicity (not excitation signal.

**Parameters**

- `f0` : f0 contour
- `spectrogram` : Spectrogram estimated by CheapTrick
- `aperiodicity` : Aperiodicity spectrogram based on D4C
- `period` :  Temporal period used for the analysis
- `fs` : Sampling frequency
- `len` : Length of the output signal

**Returns**

- `y` : Calculated speech
"""
function synthesis(f0::StridedVector{Cdouble},
                   spectrogram::StridedMatrix{Cdouble},
                   aperiodicity::StridedMatrix{Cdouble},
                   period::Real, fs::Integer, len::Integer)
    fftsize = get_fftsize_for_cheaptrick(fs)

    # Array{Cdouble,2} -> Array{Ptr{Cdouble}}
    cspectrogram = Array(Ptr{Cdouble}, size(spectrogram, 2))
    ptrarray2d!(cspectrogram, spectrogram)

    caperiodicity = Array(Ptr{Cdouble}, size(aperiodicity, 2))
    ptrarray2d!(caperiodicity, aperiodicity)

    synthesized = Array(Cdouble, len)
    ccall((:Synthesis, libworld), Void,
          (Ptr{Cdouble}, Cint, Ptr{Ptr{Cdouble}}, Ptr{Ptr{Cdouble}},
           Cint, Cdouble, Cint, Cint, Ptr{Cdouble}),
          f0, length(f0), cspectrogram, caperiodicity, fftsize, period, fs, len,
          synthesized)

    synthesized
end

# matlabfunctions

"""
$(SIGNATURES)

inplace version of interp1

**Parameters**

- `x` : Input vector (Time axis)
- `y` : Values at x[n]
- `xi`: Required vector
- `yi` : Interpolated vector
"""
function interp1!(x::StridedVector{Cdouble},
                  y::StridedVector{Cdouble},
                  xi::StridedVector{Cdouble},
                  yi::StridedVector{Cdouble})
    @assert length(x) == length(y)
    @assert length(xi) == length(yi)
    ccall((:interp1, libworld), Void,
          (Ptr{Cdouble}, Ptr{Cdouble}, Cint, Ptr{Cdouble}, Cint, Ptr{Cdouble}),
          x, y, length(x), xi, length(xi), yi)
    yi
end

"""
$(SIGNATURES)

interp1 interpolates to find yi, the values of the underlying function Y
at the points in the vector or array xi. x must be a vector.
http://www.mathworks.co.jp/help/techdoc/ref/interp1.html

**Parameters**

- `x` : Input vector (Time axis)
- `y` : Values at x[n]
- `xi`: Required vector

**Returns**

- `yi` : Interpolated vector
"""
function interp1(x::StridedVector{Cdouble},
                 y::StridedVector{Cdouble},
                 xi::StridedVector{Cdouble})
    yi = similar(xi)
    interp1!(x, y, xi, yi)
end


text:
module LoggerFormatTest

using MiniLoggers
using ReTest

conts(s, sub) = match(Regex(sub), s) !== nothing

@testset "basic format" begin
    io = IOBuffer()
    logger = MiniLogger(io = io, minlevel = MiniLoggers.Debug,
                        format = "{timestamp}:{level}:{module}:{basename}:{filepath}:{line}:{group}:{module}:{id}:{message}")
    with_logger(logger) do
        @info "foo"
        s = String(take!(io))
        @test !conts(s, "timestamp")
        @test !conts(s, "level")
        @test !conts(s, "module")
        @test !conts(s, "basename")
        @test !conts(s, "filepath")
        @test !conts(s, "line")
        @test !conts(s, "group")
        @test !conts(s, "module")
        @test !conts(s, "id")
        @test conts(s, "foo")
        @test conts(s, "Info")
    end
end

@testset "squashing" begin
    io = IOBuffer()
    logger = MiniLogger(io = io, format = "{message}")
    with_logger(logger) do
        @info "foo\nbar"

        s = String(take!(io))
        @test s == "foo\tbar\n"
    end

    logger = MiniLogger(io = io, format = "{message}", message_mode = :notransformations)
    with_logger(logger) do
        @info "foo\nbar"

        s = String(take!(io))
        @test s == "foo\nbar\n"
    end
end

@testset "squashing values" begin
    io = IOBuffer()
    logger = MiniLogger(io = io, format = "{message}", squash_delimiter = "XXX")
    val = "foo\nbar"
    with_logger(logger) do
        @info "" x = val
        s = String(take!(io))

        @test s == "x = fooXXXbar\n"
    end

    val = "bar\n\rfoo"
    with_logger(logger) do
        @info "" x = val
        s = String(take!(io))

        @test s == "x = barXXXfoo\n"
    end
end

@testset "markdown support" begin
    io = IOBuffer()
    ioc = IOContext(io, :color => true)
    logger = MiniLogger(io = ioc, format = "{message}", message_mode = :markdown)
    with_logger(logger) do
        @info "**foo**"

        s = String(take!(io))
        @test s == "  \e[1mfoo\e[22m\n"
    end
end

@testset "badcaret" begin
    io = IOBuffer()
    logger = MiniLogger(io = io, format = "{message}")
    with_logger(logger) do
        @info "foo\r\nbar"

        s = String(take!(io))
        @test s == "foo\tbar\n"
    end
end

@testset "kwargs" begin
    io = IOBuffer()
    logger = MiniLogger(io = io, format = "{message}")
    x = 1
    y = "asd"

    with_logger(logger) do
        @info "values: " x y

        s = String(take!(io))
        @test s == "values: x = 1, y = asd\n"
    end
end

@testset "error catching" begin
    io = IOBuffer()
    logger = MiniLogger(ioerr = io, format = "{message}")

    with_logger(logger) do
        try
            error("ERROR")
        catch err
            @error (err, catch_backtrace())
        end

        s = String(take!(io))
        @test conts(s, "^ERROR\n *Stacktrace")

        try
            error("ERROR")
        catch err
            @error "foo " exception = (err, catch_backtrace())
        end

        s = String(take!(io))
        @test conts(s, "^foo exception = ERROR\n *Stacktrace")

        try
            error("ERROR")
        catch err
            @error err
        end
        s = String(take!(io))
        @test s == "ERROR\n"

        try
            error("ERROR")
        catch err
            @error catch_backtrace()
        end
        s = String(take!(io))
        @test conts(s, "^\n *Stacktrace")
    end
end

@testset "logstring colors" begin
    io = IOBuffer()
    ioc = IOContext(io, :color => true)
    logger = MiniLogger(io = ioc, ioerr = ioc, format = "{foo:cyan} {bar:bold} {level:func}:{message}", minlevel = MiniLoggers.Debug)

    with_logger(logger) do
        @debug "foo"
        s = String(take!(io))

        @test s == "\e[36mfoo\e[39m \e[0m\e[1mbar\e[22m \e[34mDebug\e[39m:foo\n"

        @info "foo"
        s = String(take!(io))

        @test s == "\e[36mfoo\e[39m \e[0m\e[1mbar\e[22m \e[36mInfo\e[39m:foo\n"

        @warn "foo"
        s = String(take!(io))

        @test s == "\e[36mfoo\e[39m \e[0m\e[1mbar\e[22m \e[33mWarn\e[39m:foo\n"

        @error "foo"
        s = String(take!(io))

        @test s == "\e[36mfoo\e[39m \e[0m\e[1mbar\e[22m \e[91mError\e[39m:foo\n"
    end
end

end # module


text:
function times2(x)
  return x*2
end

times2(3)

times2(-6)

max_stars_repo_path: Problem 1/problem1.jl

max_stars_repo_name: josiex/Coding-Problems

max_issues_repo_path: Problem 1/problem1.jl

max_issues_repo_name: josiex/Coding-Problems

max_forks_repo_path: Problem 1/problem1.jl

max_forks_repo_name: josiex/Coding-Problems


text:
export LSRK144, LSRK54, solve!, dostep!

function solve!(q, timeend, solver;
                after_step::Function = (x...) -> nothing,
                after_stage::Function = (x...) -> nothing)
  finalstep = false
  step = 0
  while true
    step += 1
    time = solver.time
    if time + solver.dt >= timeend
      solver.dt = timeend - time
      finalstep = true
    end
    dostep!(q, solver, after_stage)
    after_step(step, time, q)
    finalstep && break
  end
end

mutable struct LSRK{FT, AT, NS, RHS}
  dt::FT
  time::FT
  rhs!::RHS
  dq::AT
  rka::NTuple{NS, FT}
  rkb::NTuple{NS, FT}
  rkc::NTuple{NS, FT}

  function LSRK(rhs!, rka, rkb, rkc, q, dt, t0)
      FT = eltype(eltype(q))
      dq = similar(q)
      fill!(dq, zero(eltype(q)))
      AT = typeof(q)
      RHS = typeof(rhs!)
      new{FT, AT, length(rka), RHS}(FT(dt), FT(t0), rhs!, dq, rka, rkb, rkc)
  end
end

function dostep!(q, lsrk::LSRK, after_stage)
  @unpack rhs!, dq, rka, rkb, rkc, dt, time = lsrk
  for stage = 1:length(rka)
    stagetime = time + rkc[stage] * dt
    dq .*= rka[stage]
    rhs!(dq, q, stagetime)
    @. q += rkb[stage] * dt * dq
    after_stage(stagetime, q)
  end
  lsrk.time += dt
end

function LSRK144(rhs!, q, dt; t0 = 0)
  rka = (
     0,
    -0.7188012108672410,
    -0.7785331173421570,
    -0.0053282796654044,
    -0.8552979934029281,
    -3.9564138245774565,
    -1.5780575380587385,
    -2.0837094552574054,
    -0.7483334182761610,
    -0.7032861106563359,
     0.0013917096117681,
    -0.0932075369637460,
    -0.9514200470875948,
    -7.1151571693922548,
  )

  rkb = (
    0.0367762454319673,
    0.3136296607553959,
    0.1531848691869027,
    0.0030097086818182,
    0.3326293790646110,
    0.2440251405350864,
    0.3718879239592277,
    0.6204126221582444,
    0.1524043173028741,
    0.0760894927419266,
    0.0077604214040978,
    0.0024647284755382,
    0.0780348340049386,
    5.5059777270269628,
  )

  rkc = (
    0,
    0.0367762454319673,
    0.1249685262725025,
    0.2446177702277698,
    0.2476149531070420,
    0.2969311120382472,
    0.3978149645802642,
    0.5270854589440328,
    0.6981269994175695,
    0.8190890835352128,
    0.8527059887098624,
    0.8604711817462826,
    0.8627060376969976,
    0.8734213127600976,
  )
  LSRK(rhs!, rka, rkb, rkc, q, dt, t0)
end

function LSRK54(rhs!, q, dt; t0 = 0)
  rka = (
    (0),
    (-567301805773 // 1357537059087),
    (-2404267990393 // 2016746695238),
    (-3550918686646 // 2091501179385),
    (-1275806237668 // 842570457699),
  )

  rkb = (
    (1432997174477 // 9575080441755),
    (5161836677717 // 13612068292357),
    (1720146321549 // 2090206949498),
    (3134564353537 // 4481467310338),
    (2277821191437 // 14882151754819),
  )

  rkc = (
    (0),
    (1432997174477 // 9575080441755),
    (2526269341429 // 6820363962896),
    (2006345519317 // 3224310063776),
    (2802321613138 // 2924317926251),
  )
  LSRK(rhs!, rka, rkb, rkc, q, dt, t0)
end


text:
const LEAD_IN_LENGTH = 0x1c
const DECIMATE_MASK = 0b00100000
const convertible_fields = (:id, :name, :loc, :fs, :gain, :resp, :units, :src, :notes, :misc)
const nodalfields = (:id, :loc, :fs, :gain, :misc, :name, :notes, :resp, :src, :units, :t)
const tdms_dtos = round(Int64, d2u(DateTime("1904-01-01T00:00:00")))
const tdms_codes = Dict{UInt32, Type}(
  0x00000000 => UInt8,
  0x00000001 => Int8,
  0x00000002 => Int16,
  0x00000003 => Int32,
  0x00000004 => Int64,
  0x00000005 => UInt8,
  0x00000006 => UInt16,
  0x00000007 => UInt32,
  0x00000008 => UInt64,
  0x00000009 => Float32,
  0x0000000a => Float64,
  0x00000020 => Char,
  0x00000021 => Bool,
  0x00000044 => UInt64,   # *** convert to date string
)
const unindexed_fields = (:n, :ox, :oy, :oz, :data, :info, :x)

const TDMS = TDMSbuf(
    zero(UInt32),
    zero(UInt64),
    zero(UInt64),
    zero(UInt32),
    zero(Int64),
    zero(Float64),
    zero(Float64),
    zero(Float64),
    zero(Float64),
    "",
    Dict{String, Any}()
    )


text:
using Plots; pyplot()
using DeSolver
using LaTeXStrings

function test(u,t)
    du=sin(t)^2*u
    return du
end

function testexact(x)
    # exact solution of test with the initial condition u=1
    return exp(1/2*(x - sin(x)*cos(x)))
end

function wave(u,t)
    du= u[2]
    ddu = -u[1]
    return [du, ddu]
end

soltest4=odesolve(test, 1., (0.,10.), 1, method="RK4")
soltest2=odesolve(test, 1., (0.,10.), 1, method="RK2")
soltest5=odesolve(test, 1., (0.,10.), 0.25, method="RK4")

solwave4=odesolve(wave, [1., 0.], (0.,10.), 1, method="RK4")
solwave2=odesolve(wave, [1., 0.], (0.,10.), 1, method="RK2")
solwave5=odesolve(wave, [1., 0.], (0.,10.), 0.25, method="RK4")

display(plot(map(x->x[1],soltest4.t),map(x->x[1],soltest4.u), label=L"u'=sin(t)^2u"*"  with RK4 stepsize 1"))
display(plot!(map(x->x[1],soltest2.t),map(x->x[1],soltest2.u), label=L"u'=sin(t)^2u"*"  with RK2 stepsize 1"))
display(plot!(map(x->x[1],soltest5.t),map(x->x[1],soltest5.u), label=L"u'=sin(t)^2u"*"  with RK4 stepsize .25"))
display(plot!(testexact, label="Exact Solution of  "*L"u'=sin(t)^2u"))

display(plot(map(x->x[1],solwave4.t),map(x->x[1],solwave4.u), label=L"u''=-u"*"  with RK4 stepsize 1"))
display(plot!(map(x->x[1],solwave2.t),map(x->x[1],solwave2.u), label=L"u''=-u"*"  with RK2 stepsize 1"))
display(plot!(map(x->x[1],solwave5.t),map(x->x[1],solwave5.u), label=L"u''=-u"*"  with RK4 stepsize .25"))
display(plot!(cos, label = "Exact solution of  "*L"u''=-u"))

# gui()
# plot(p1)
# plot!(p4)


text:
# 2D Field tracing on a unstructured grid.

export trace

const Δ = 100000. # distance to the far field point
const ϵ = 1e-5 # small perturbation

"""
	 trace(mesh::SimpleMesh, vx, vy, xstart, ystart; maxIter=1000, maxLen=1000.)

2D stream tracing on unstructured quadrilateral and triangular mesh.
"""
function trace(mesh::SimpleMesh, vx, vy, xstart, ystart; maxIter=1000, maxLen=1000.)

   xStream = [fill(xs, maxIter) for xs in xstart]
   yStream = [fill(ys, maxIter) for ys in ystart]
   nIter = fill(maxIter, length(xStream))

   @inbounds for iS in 1:length(xStream)
      Pnow = Point(xStream[iS][1], yStream[iS][1])
      P⁺ = Point(0.0f0, 0.0f0)
      cellID = getCellID(mesh, Pnow)
      cellIDNew = 0

      for it = 1:maxIter
         Pfar = Pnow + Vec2f(vx[cellID]*Δ, vy[cellID]*Δ)
         element = getelement(mesh, cellID)

         if mesh[cellID] isa Quadrangle
            nEdge = 4
         elseif mesh[cellID] isa Triangle
            nEdge = 3
         end

         ray = Segment(Pnow, Pfar)
         # Loop over edges
         for i = 1:nEdge
            j = i % nEdge + 1
            P = Segment(element.vertices[i], element.vertices[j]) ∩ ray
            if P isa Point
               P⁺ = P + Vec2f(vx[cellID]*ϵ, vy[cellID]*ϵ)
               cellIDNew = getCellID(mesh, P⁺)
               break
            elseif P isa Segment
               P⁺ = element.vertices[j] + Vec2f(vx[cellID]*ϵ, vy[cellID]*ϵ)
               cellIDNew = getCellID(mesh, P⁺)
               break
            end
         end

         Pnow = P⁺

         xStream[iS][it+1] = Pnow.coords[1]
         yStream[iS][it+1] = Pnow.coords[2]
         nIter[iS] = it+1

         if cellIDNew == 0 # hit the boundary
            break
         else
            cellID = cellIDNew
         end
      end
   end

   @inbounds for i = 1:length(xStream)
      xStream[i] = xStream[i][1:nIter[i]]
      yStream[i] = yStream[i][1:nIter[i]]
   end

   xStream, yStream
end

"Return cell ID on the unstructured mesh."
function getCellID(mesh::SimpleMesh, point::Point2)
   for (i, element) in enumerate(elements(mesh))
      if point ∈ element
         return i
      end
   end
   return 0 # out of mesh boundary
end

"Return the `cellID`th element of the mesh."
function getelement(mesh, cellID)
   for (i, element) in enumerate(elements(mesh))
      if i == cellID
         return element
      end
   end
end


text:
module DICE

using JuMP;
#Ipopt is distributed under the EPL.
#We don't package it here though, and assume you have this set up on your system already.
using Ipopt;

include("Abstractions.jl")

#DICE Versions
abstract type Version end

#Configuration options & parameters, model settings and output
include("BaseTypes.jl")

#Scenarios
abstract type Scenario end
include("Scenarios.jl")

struct DICENarrative
    constants::Options
    parameters::Parameters
    model::JuMP.Model
    scenario::Scenario
    version::Version
    variables::Variables
    equations::Equations
    results::Results
end

function Base.show(io::IO, ::MIME"text/plain", dice::DICENarrative)
    println(io, "$(dice.scenario) scenario using $(dice.version).")
    show(io, dice.model)
end

function solve end
function options end

export solve, options

# A list of all linear solvers that Ipopt will accept.
# Note that DICE.jl is only tested against ma97 and mumps.
# We use this externally so users cannot accidentally send in nonsense strings,
# although the result is string based and used throughout as such.
@exported_enum ipoptLinearSolver ma27 ma57 ma77 ma86 ma97 pardiso wsmp mumps custom

# Test if system has HSL MA97 installed, fall back to MUMPS if not.
function selectLinearSolver(solver_name::ipoptLinearSolver = ma97)
    prob = Ipopt.CreateIpoptProblem(1,[1.],[1.],1,[1.],[1.],1,1,sum,sum,sum,sum,sum);
    Ipopt.AddIpoptStrOption(prob, "sb", "yes");
    Ipopt.AddIpoptIntOption(prob, "print_level", 0);
    solver_name_string = string(solver_name);
    # Initially, we must check that coinhsl is installed at all if we want a HSL solver.
    # If not, Ipopt will default to try and find any dynamically linked libhsl. If it
    # cannot find one it will hard panic and we can't capture that failure.
    if occursin("ma", solver_name_string)
        Ipopt.AddIpoptStrOption(prob, "linear_solver", "ma27");
        try
            # No HSL was found on the system, we return with fallback.
            return runLinearSolverCheck(prob, solver_name_string)
        catch
            #Continue, we have access to at least some version of HSL
        end
    end
    try
        # Outer try will fail if solver string is not in the list of
        # possible Ipopt solvers.
        # For now that's ma27, ma57, ma77, ma86, ma97, pardiso, wsmp, mumps, custom
        Ipopt.AddIpoptStrOption(prob, "linear_solver", solver_name_string);
        try
            # Inner try attempts to run the dummy program and will crash because
            # the dummy is malformed.
            # No error code will be returned if the solver is extant, and an Invalid_Option
            # if the library is not found.
            runLinearSolverCheck(prob, solver_name_string)
        catch
            # We can use the requested solver.
            solver_name_string
        end
    catch
        "mumps"
    end
end

function runLinearSolverCheck(prob::IpoptProblem, solver_name_string::String)
    result_code = Ipopt.IpoptSolve(prob);
    if Ipopt._STATUS_CODES[result_code] == :Invalid_Option
        @info "Unable to set linear_solver = $(solver_name_string), defaulting to MUMPS."
        return "mumps"
    else
        error("Attempts to identify linear solvers on system returned unexpected results.");
    end
end

# Include all version implementations
include("2013R.jl")
include("2016.jl")

end


text:
using MakieLayout
using AbstractPlotting
using Test

@testset "Layoutables constructors" begin
    scene, layout = layoutscene()
    ax = layout[1, 1] = LAxis(scene)
    cb = layout[1, 2] = LColorbar(scene)
    gl2 = layout[2, :] = MakieLayout.GridLayout()
    bu = gl2[1, 1] = LButton(scene)
    sl = gl2[1, 2] = LSlider(scene)

    scat = scatter!(ax, rand(10))
    le = gl2[1, 3] = LLegend(scene, [scat], ["scatter"])

    to = gl2[1, 4] = LToggle(scene)
    te = layout[0, :] = LText(scene, "A super title")
    me = layout[end+1, :] = LMenu(scene, options = ["one", "two", "three"])
end

max_stars_repo_name: jkrumbiegel/MakieLayout


text:
"""
    Path
A path is a sequence of node traversals that exist in a graph.
A path in a graph usually have flows.
"""
mutable struct Path
    nodelist::Vector
    keylist::Vector
end

Path() = Path([], [])

nodelist(p::Path) = p.nodelist
keylist(p::Path) = p.keylist

isempty(p::Path) = isempty(p.nodelist) && isempty(p.keylist)
==(p1::Path,p2::Path) = p1.nodelist == p2.nodelist && p1.keylist == p2.keylist

# deconstruct a path
"""
    ArcIter(p)
breaks down the arcs in a path.
"""
function ArcIter(p::Path)
    len = length(p.keylist)
    ulist = []; sizehint!(ulist, len)
    vlist = []; sizehint!(vlist, len)
    keylist = []; sizehint!(keylist, len)
    for i in eachindex(p.keylist)
        u = p.nodelist[i]; push!(ulist, u)
        v = p.nodelist[i + 1]; push!(vlist, v)
        key = p.keylist[i]; push!(keylist, key)
    end
    return ArcIter(len, ulist, vlist, keylist)
end

arcs(p::Path) = ArcIter(p)

"""
    Graph(p, flow, attrib)
creates a graph from the path with a fixed flow from beginning to end.
"""
function Graph(p::Path, flow=1, attrib=:flow)
    G = Graph()
    add_arcs_from!(G, collect(arcs(p)))
    for a in arcs(G)
        G[a, attrib] = flow
    end
    return Graph(G.nn, G.na, G.fadj, G.badj)
end

"""
    depth_first_search(g, u, v)
finds the shortest path (complete with keys) by first traversing and then assigning key.
"""
function depth_first_search(g::Graph, u, v, p=Path(), shortest=Path(), method=rand)
    P = Path(deepcopy(p.nodelist), deepcopy(p.keylist))
    push!(P.nodelist, u)
    if length(P.nodelist) > 1
        key = keys(g.fadj[last(p.nodelist)][u])
        push!(P.keylist, method(key))
    end
    # println("current path: $(P)")
    if u == v
        return P
    end
    for n in out_neighbor(g, u)
        if !(n in P.nodelist)
            if shortest.nodelist == [] || length(P.nodelist) < length(shortest.nodelist)
                new_p = depth_first_search(g, n, v, P, shortest)
                if new_p.nodelist != []
                    shortest = new_p
                end
            end
        end
    end
    return shortest
end

"""
    dijkstra_shortest_path(g, u, v)
finds the shortest path from u to v in g using the dijkstra algorithm.
"""
function dijkstra_shortest_path(g::Graph, u, v=nothing, w=:flow)
    pred = dictionary(k => [] for k in keys(g.fadj))
    dist = dictionary(k => Inf for k in keys(g.fadj)); H = PriorityQueue()
    dist[u] = 0; H[u] = 0
    while !isempty(H)
        u = dequeue!(H)
        d = dist[u]
        for v in out_neighbor(g, u)
            val, ind = findmin(map(p -> getindex(p, w), g.fadj[u][v]))
            alt = d + val
            if alt < dist[v]
                push!(pred[v], (node = u, key = ind))
                dist[v] = alt
                H[v] = alt
            end
        end
    end

    if isnothing(v)
        return dist, pred
    else
        S = Path()
        complete = false
        while !isempty(pred[v])
            push!(S.nodelist, v)
            push!(S.keylist, last(pred[v]).key)
            v = last(pred[v]).node
        end
        if !isempty(S)
            push!(S.nodelist, v)
            reverse!(S.nodelist)
            reverse!(S.keylist)
        end
        return S
    end
end


text:
module LocalHillResponses

using ForwardDiff

export local_sensitivity, respose_coefficient, local_hills

"""
    local_sensitivity(f, s)

Compute the local sensitivity (`df/ds`) of the function `f` for the singular input `s`.
"""
function local_sensitivity(f, s)
    return ForwardDiff.derivative(f, s)
end

"""
    respose_coefficient(f, s)

Compute the response coefficient (a.k.a. relative local sensitivity):

```math
R^{x}_{s} = \frac{dx}{ds} \frac{s}{x}
```

where ``R^{x}_{s}`` is the response coefficient of the input ``s`` and the output ``x = f(s)``.

# Warning

The result may be fluctuated due to the singularity point.
"""
function respose_coefficient(f, s, y = f(s))
    return local_sensitivity(f, s) * s / y
end

const relative_local_sensitivity = respose_coefficient

"""
Check if function `f(s)` is monotonous between `smin` and `smax`. A positive `direction` checks increasing function; A negative `direction` checks decreasing function. If no `direction` is given, checks either increasing or decreasing.
"""
function is_monotonous(f, smin, smax, direction; Δs = (smax - smin) * 1e-3)
    for s in smin:Δs:smax
        if ForwardDiff.derivative(f, s) * direction < 0
            return false
        end
    end
    return true
end

function is_monotonous(f, smin, smax; Δs = (smax - smin) * 1e-3)
    is_monotonous(f, smin, smax, +1; Δs = Δs) || is_monotonous(f, smin, smax, -1; Δs = Δs)
end

"""
Activated fraction for reponse `y` bounded by minimal response `ymin` and maximal response `ymax`.
"""
activated_fraction(y, ymin, ymax) = (y - ymin) / (ymax - ymin)

"""
Local Hill coefficient.
"""
function local_hills(f, s, ymin, ymax)
    responseCoef = respose_coefficient(f, s)
    actFrac = activated_fraction(y, ymin, ymax)
    nh = responseCoef / (one(actFrac) - actFrac)
    return nh
end

end # module


text:
# ---
# title: 242. Valid Anagram
# id: problem242
# author: zhwang
# date: 2022-02-15
# difficulty: Easy
# categories: Hash Table, Sort
# link: <https://leetcode.com/problems/valid-anagram/description/>
# hidden: true
# ---
#
# Given two strings _s_ and _t  _, write a function to determine if _t_ is an
# anagram of _s_.
#
# **Example 1:**
#
#
#
#     Input: _s_ = "anagram", _t_ = "nagaram"
#     Output: true
#
#
# **Example 2:**
#
#
#
#     Input: _s_ = "rat", _t_ = "car"
#     Output: false
#
#
# **Note:**
# You may assume the string contains only lowercase alphabets.
#
# **Follow up:**
# What if the inputs contain unicode characters? How would you adapt your
# solution to such case?
#
#
## @lc code=start
using LeetCode
using DataStructures
function is_anagram(s::String, t::String)::Bool
    length(s) != length(t) && return false
    count_s, count_t = DefaultDict{Char,Int}(0), DefaultDict{Char,Int}(0)
    for i in eachindex(s)
        count_s[s[i]] += 1
        count_t[t[i]] += 1
    end
    return all(count_s[c] == count_t[c] for c in 'a':'z')
end

## @lc code=end


text:
using Test
using SortingLab

x = [1,2,missing, 3]

@test isequal(sort(x), fsort(x))


text:
using Fermi
using Test
using Suppressor

molecules = ["water",   "ammonia",     "benzene", "ethanol", "formaldehyde", "glycine", "methane", "phosphaethene", "acetic_acid", "nitrogen"]
basis =     ["cc-pvtz", "aug-cc-pvdz", "6-31g",   "cc-pvdz", "6-31g*",       "sto-3g",  "cc-pvtz", "3-21g",         "6-311+g",     "cc-pcvdz"]

tol = 1E-8
@testset "Fermi" begin
    include("test_options.jl")
    include("test_output.jl")
    include("test_arrays.jl")
    include("test_contract.jl")
    include("test_geometry.jl")
    include("test_basisset.jl")
    include("test_integrals.jl")
    include("test_diis.jl")
    include("test_orbitals.jl")
    include("test_phycons.jl")
    include("test_RHF.jl")
    include("test_MP.jl")
    include("test_pT.jl")
    @reset
end


text:
using Plots
using Statistics
using StatsBase
import Base.length
# histograms

"""
    digitize(x, bins)

    Return the indices of the bins to which each value in input array belongs.
"""
function digitize(x::Vector{<:Real}, bins::LinRange{<:Real})
    return searchsortedlast.(Ref(bins), x)
end
# Notice that in this case the bins are expected to be sorted
function digitize(x::Vector{<:Real}, bins::Vector{<:Real})
    return searchsortedlast.(Ref(bins), x)
end

"""
    hist1d(x::Vector{T}, nbins::Integer, xmin::T=typemin(T), xmax::T=typemax(T))
    hist1d(x::Vector{T}, xs::String, nbins::Integer,
                    xmin::T=typemin(T), xmax::T=typemax(T); datap = true, fwl=true)
    hist1d(h::Histogram, xs::String; datap = true, markersize=3, fwl=false)

return a 1d histogram and its corresponding graphics (plots)
"""
function hist1d(x::Vector{<:Real}, bins::Vector{<:Real}, norm=false)
    h = fit(Histogram, x, bins)
    if norm
        h = StatsBase.normalize(h, mode=:density)
    end
    return h
end

function hist1d(x::Vector{T}, nbins::Integer,
                xmin::T=typemin(T), xmax::T=typemax(T), norm=false) where T
    xx   = in_range(x, xmin, xmax)
    bins = collect(LinRange(xmin, xmax, nbins + 1))
    h    = hist1d(xx, bins, norm)
    return h
end

"""
    edges(h::Histogram)

edges of the histogram
"""
function edges(h::Histogram, index::Int64=1)
    collect(h.edges[index])
end

"""
    centers(h::Histogram)

return the centres of the histogram bins.
"""
centers(h::Histogram) = centers(edges(h))

"""
    centers(edges::Vector{<:Real})

Calculate the bin centers from a vector of bin edges
"""
function centers(edges::Vector{<:Real})
    edges[1:end-1] + .-(@view(edges[2:end]), @view(edges[1:end-1])) / 2
end

"""
    hist_weights

return histogram weights of a set of data.
"""
function hist_weights(edges::Vector{<:Real})
  function get_weights(y::SubArray{<:Real})
    histo = fit(Histogram, y, edges)
    return histo.weights
  end
  return get_weights
end

"""
    hist1d(x::Vector{T}, xs::String, nbins::Integer,
           xmin::T=typemin(T), xmax::T=typemin(T);
           norm::Bool=false, datap::Bool=true,
           markersize::Int64=3, fwl::Bool=false,
           label::String="", legend::Bool=false) where T <: Real
return 1D histogram and plot axes for data x and number of bins nbins.
"""
function hist1d(x::Vector{T}, xs::String, nbins::Integer,
                xmin::T=typemin(T), xmax::T=typemax(T);
                norm::Bool=false, datap::Bool=true,
                markersize::Int64=3, fwl::Bool=false,
                label::String="", legend::Bool=false) where T <: Real

    h = hist1d(x, nbins, xmin, xmax, norm)
    return h, hist1d(h, xs, datap=datap, markersize=markersize, fwl=fwl,
                     label=label, legend=legend)
end

"""
    hist1d(h::Histogram, xs::String;
           datap::Bool=true, markersize::Int64=3, fwl::Bool=false,
           label::String="", legend::Bool=false)
Plot 1D histogram from Histogram object and options.
"""
function hist1d(h::Histogram, xs::String;
                datap::Bool=true, markersize::Int64=3, fwl::Bool=false,
                label::String="", legend::Bool=false)

    if datap
        xg = centers(h)
        yg = eltype(xg).(h.weights)
        p  = scatter(xg, yg, yerr=sqrt.(yg),
                     markersize=markersize, label=label, legend=legend)
        if fwl
            p = plot!(p, xg, yg, yerr=sqrt.(yg), fmt=:png,
                      linewidth=1, label=label, legend=legend)
        end
    else
        p = plot(h, xlabel=xs, fmt=:png, yl="frequency")
    end
    xlabel!(xs)
    ylabel!("frequency")

    return p
end

"""
    plot2hist1d(h1::Histogram, h2::Histogram, xs::String;
                markersize::Int64=2, norm::Bool=false)
Plot two 1D histograms on the same axes.
"""
function plot2hist1d(h1::Histogram, h2::Histogram, xs::String;
                     markersize::Int64=2, norm::Bool=false)

    if norm
        h1 = StatsBase.normalize(h1, mode=:density)
        h2 = StatsBase.normalize(h2, mode=:density)
    end

    xg1 = centers(h1)
    yg1 = eltype(xg1).(h1.weights)
    xg2 = centers(h2)
    yg2 = eltype(xg2).(h2.weights)

    p1 = scatter(     xg1, yg1, yerr=sqrt.(yg1), fmt=:png, markersize=markersize)
    p1 = scatter!(p1, xg2, yg2, yerr=sqrt.(yg2), fmt=:png, markersize=markersize)
    xlabel!(xs)
    ylabel!("frequency")

    return p1
end

"""
    hist2d(x::Vector{T}, y::Vector{T}, nbins::Integer,
           xl::String, yl::String,
           xmin::T=typemin(T), xmax::T=typemax(T),
           ymin::T=typemin(T), ymax::T=typemax(T)) where T <: Real

return a 2d histogram and its corresponding graphics (plots)
given data x, y, labels xl and yl and limits for each axis.
"""
function hist2d(x::Vector{T},y::Vector{T}, nbins::Integer,
                xl::String, yl::String,
                xmin::T=typemin(T), xmax::T=typemax(T),
                ymin::T=typemin(T), ymax::T=typemax(T); title="") where T <: Real
    function select_data(xval::T, yval::T)
        range_bound(xmin, xmax, OpenBound)(xval) .&
        range_bound(ymin, ymax, OpenBound)(yval)
    end
    mask = select_data.(x, y)
    data = (y[mask], x[mask])
    h    = fit(Histogram, data, nbins=nbins)
    ye   = centers(edges(h, 1))
    xe   = centers(edges(h, 2))
    hm   = heatmap(xe, ye, h.weights)
    xlabel!(xl)
    ylabel!(yl)
    if title != ""
        title!(title)
    end
    return h, hm
end

"""
    fmean_std

returns the weighted mean and std of an array given a minimum
bin weight for consideration.
"""
function fmean_std(x::Vector{<:Real}, min_prop::Float64=0.1)
    mask_func = weights -> weights .> min_prop * maximum(weights)
    function filt_mean(w::SubArray{<:Real})
        mask = mask_func(w)
        return [mean_and_std(x[mask], FrequencyWeights(w[mask]), corrected=true)]
    end
    return filt_mean
end

"""
    p1df(x, y, nbins; ybin_width, ymin, ymax, min_proportion)

return a profile DataFrame. This is a DF in which the variable y is
histogrammed as a function of the average of variable x in each bin.
The keyword arguments allow for a filtering in the y variable by
range and min_proportion of the maximum in the bin.
TODO: Protect against fine binning that results in zeros.
"""
function p1df(x::Vector{<:Real}, y::Vector{<:Real}, nbins::Integer;
              ybin_width::Real=0.1, ymin::Real=minimum(y), ymax::Real=maximum(y),
              min_proportion::Real=0.0)
    df           = DataFrame(:y => y)
    x_upper      = nextfloat(maximum(x))
    bin_edges    = LinRange(minimum(x), x_upper, nbins + 1)
    df[!, "bin"] = digitize(x, bin_edges)
    bin_centers  = centers(collect(bin_edges))
    bin_width    = .-(@view(bin_edges[2:end]), @view(bin_edges[1:end-1]))
    ## Bin in y so filtering can be done on bin content
    y_upper      = nextfloat(ymax)
    y_bins       = collect(ymin:ybin_width:y_upper)
    y_centers    = centers(y_bins)
    y_weights    = combine(groupby(df, :bin),
                           :y => hist_weights(y_bins) => :yweights,
                           ungroup = false)
    bin_stats    = fmean_std(y_centers, min_proportion)
    ndf = combine(y_weights,
                  :yweights => bin_stats => [:y_mean, :y_std])
    ndf = ndf[:, [:y_mean, :y_std]]
    ndf[!, :x_mean] = bin_centers
    ndf[!, :x_std]  = bin_width / 2

    p1 = plot(ndf.x_mean,ndf.y_mean, yerror=ndf.y_std, fmt = :png,
              shape = :circle, color = :black, legend=false)
    return ndf, p1
end

max_stars_repo_name: proscrite/LaserLab

max_issues_repo_name: proscrite/LaserLab

max_forks_repo_name: proscrite/LaserLab


text:
import Base: hash, parse

abstract type TxComponent end

mutable struct TxIn <: TxComponent
    prev_tx::Vector{UInt8}
    prev_index::UInt32
    script_sig::Script
    witness::Script
    sequence::UInt32
    TxIn(prev_tx, prev_index) = new(prev_tx, prev_index, Script(nothing), Script(nothing), 0xffffffff)
    TxIn(prev_tx, prev_index, script_sig::Script, sequence::Integer=0xffffffff) = new(prev_tx, prev_index, script_sig, Script(nothing), sequence)
    TxIn(prev_tx, prev_index, script_sig::Script, witness::Script, sequence::Integer=0xffffffff) = new(prev_tx, prev_index, script_sig, witness, sequence)
end

function show(io::IO, z::TxIn)
    print(io, "\n", bytes2hex(z.prev_tx), ":", z.prev_index, "\n", z.script_sig)
end

"""
Takes a byte stream and parses the tx_input at the start
return a TxIn object
"""
function TxIn(s::IOBuffer)
    prev_tx = read(s, 32)
    reverse!(prev_tx)
    bytes = read(s, 4)
    prev_index = to_int(bytes, little_endian=true)
    script_sig = scriptparse(s)
    readbytes!(s, bytes, 4)
    sequence = to_int(bytes, little_endian=true)
    return TxIn(prev_tx, prev_index, script_sig, sequence)
end

"""
    serialize(tx::TxIn) -> Vector{UInt8}

Returns the byte serialization of the transaction input
"""
function serialize(tx::TxIn)
    result = copy(tx.prev_tx)
    reverse!(result)
    append!(result, bytes(tx.prev_index, len=4, little_endian=true))
    append!(result, serialize(tx.script_sig))
    append!(result, bytes(tx.sequence, len=4, little_endian=true))
    return result
end

function get_tx(tx::TxIn; testnet::Bool=false)
    return get_tx(bytes2hex(tx.prev_tx), testnet=testnet)
end

"""
    value(txin::TxIn, testnet::Bool=false) ->

Get the outpoint value by looking up the tx hash
Returns the amount in satoshi
"""
function value(txin::TxIn, testnet::Bool=false)
    tx = get_tx(txin, testnet=testnet)
    return tx.tx_outs[txin.prev_index + 1].amount
end

"""
    script_pubkey(txin::TxIn, testnet::Bool=false)
    -> Script

Get the scriptPubKey by looking up the tx hash
Returns a Script object
"""
function script_pubkey(txin::TxIn, testnet::Bool=false)
    tx = get_tx(txin, testnet=testnet)
    return tx.tx_outs[txin.prev_index + 1].script_pubkey
end

struct TxOut <: TxComponent
    amount::UInt64
    script_pubkey::Script
    TxOut(amount, script_pubkey) = new(amount, script_pubkey)
end

function show(io::IO, z::TxOut)
    print(io, "\n", z.script_pubkey, "\namout (BTC) : ", z.amount / 100000000)
end

"""
     txoutparse(s::IOBuffer) -> TxOut

Takes a byte stream and parses the tx_output at the start
return a TxOut object
"""
function txoutparse(s::Base.GenericIOBuffer)
    bytes = UInt8[]
    readbytes!(s, bytes, 8)
    amount = to_int(bytes, little_endian=true)
    script_pubkey = scriptparse(s)
    return TxOut(amount, script_pubkey)
end

"""
    serialize(tx::TxOut) -> Vector{UInt8}

Returns the byte serialization of the transaction output
"""
function serialize(tx::TxOut)
    result = bytes(tx.amount, len=8, little_endian=true)
    append!(result, serialize(tx.script_pubkey))
    return result
end

mutable struct Tx <: TxComponent
    version::UInt32
    tx_ins::Vector{TxIn}
    tx_outs::Vector{TxOut}
    locktime::UInt32
    testnet::Bool
    segwit::Bool
    flag::UInt8
    _hash_prevouts
    _hash_sequence
    _hash_outputs
    Tx(version::Integer, tx_ins, tx_outs, locktime::Integer, testnet=false, segwit=false, flag=0x00) = new(UInt32(version), tx_ins, tx_outs, UInt32(locktime), testnet, segwit, flag, nothing, nothing, nothing)
end

function show(io::IO, z::Tx)
    print(io, "Transaction\n--------\nTestnet : ", z.testnet,
            "\nVersion : ", z.version,
            "\nLocktime : ", z.locktime,
            "\n--------\n",
            "\n", z.tx_ins,
            "\n--------\n",
            "\n", z.tx_outs)
end

function parse_legacy(s::Base.GenericIOBuffer, testnet::Bool=false)
    version = ltoh(reinterpret(UInt32, read(s, 4))[1])
    num_inputs = read_varint(s)
    inputs = []
    for i in 1:num_inputs
        input = TxIn(s)
        push!(inputs, input)
    end
    num_outputs = read_varint(s)
    outputs = []
    for i in 1:num_outputs
        output = txoutparse(s)
        push!(outputs, output)
    end
    locktime = ltoh(reinterpret(UInt32, read(s, 4))[1])
    return Tx(version, inputs, outputs, locktime, testnet)
end

function parse_segwit(s::Base.GenericIOBuffer, testnet::Bool=false)
    version = ltoh(reinterpret(UInt32, read(s, 4))[1])
    flag = read(s, 2)[2]
    num_inputs = read_varint(s)
    inputs = []
    for _ ∈ 1:num_inputs
        push!(inputs, TxIn(s))
    end
    num_outputs = read_varint(s)
    outputs = []
    for _ ∈ 1:num_outputs
        push!(outputs, txoutparse(s))
    end
    for tx_in ∈ inputs
        items = []
        num_items = read_varint(s)
        for _ in 1:num_items
            item_len = read_varint(s)
            if item_len == 0
                push!(items, 0)
            else
                push!(items, read(s, item_len))
            end
        end
        tx_in.witness.instructions = items
    end
    locktime = ltoh(reinterpret(UInt32, read(s, 4))[1])
    Tx(version, inputs, outputs, locktime,
    testnet, true, flag)
end

function parse(s::IOBuffer, testnet::Bool=false)::Tx
    if read(s, 5)[5] == 0x00
        f = parse_segwit
    else
        f = parse_legacy
    end
    seekstart(s)
    f(s, testnet)
end

function payload2tx(payload::Vector{UInt8})
    txparse(IOBuffer(payload))
end

"""
    serialize(tx::Tx) -> Vector{UInt8}

Returns the byte serialization of the transaction
"""
serialize(tx::Tx) = tx.segwit ? serialize_segwit(tx) : serialize_legacy(tx)

function serialize_legacy(tx::Tx)
    result = bytes(tx.version, len=4, little_endian=true)
    append!(result, encode_varint(length(tx.tx_ins)))
    for tx_in in tx.tx_ins
        append!(result, serialize(tx_in))
    end
    append!(result, encode_varint(length(tx.tx_outs)))
    for tx_out in tx.tx_outs
        append!(result, serialize(tx_out))
    end
    append!(result, bytes(tx.locktime, len=4, little_endian=true))
    return result
end

function serialize_segwit(tx::Tx)
    result = bytes(tx.version, len=4, little_endian=true)
    append!(result, [0x00, tx.flag])
    append!(result, encode_varint(length(tx.tx_ins)))
    for tx_in in tx.tx_ins
        append!(result, serialize(tx_in))
    end
    append!(result, encode_varint(length(tx.tx_outs)))
    for tx_out in tx.tx_outs
        append!(result, serialize(tx_out))
    end
    for tx_in in tx.tx_ins
        append!(result, UInt8(length(tx_in.witness.instructions)))
        for item in tx_in.witness.instructions
            if typeof(item) <: Vector
                append!(result, encode_varint(length(item)))
            end
            append!(result, item)
        end
    end
    append!(result, bytes(tx.locktime, len=4, little_endian=true))
    return result
end

"""
Binary hash of the legacy serialization
"""
function hash(tx::Tx)
    return reverse(hash256(serialize_legacy(tx)))
end

"""
    id(tx::Tx) -> String

Returns an hexadecimal string of the transaction hash
"""
function id(tx::Tx)
    return bytes2hex(hash(tx))
end

"""
    fee(tx::Tx) -> Integer

Returns the fee of this transaction in satoshi
"""
function fee(tx::Tx)
    input_sum, output_sum = 0, 0
    for tx_in in tx.tx_ins
        input_sum += value(tx_in, tx.testnet)
    end
    for tx_out in tx.tx_outs
        output_sum += tx_out.amount
    end
    return input_sum - output_sum
end

"""
    sig_hash(tx::Tx, input_index::Integer)::Vector{UInt8}

Returns the hash that needs to get signed for index input_index
"""
function sig_hash(tx::Tx, input_index::Integer, redeem_script::Union{Script,Nothing}=nothing)
    s = Vector(reinterpret(UInt8, [htol(tx.version)]))
    append!(s, encode_varint(length(tx.tx_ins)))

    i, script_sig = 0, Script(nothing)
    for tx_in in tx.tx_ins
        if i == input_index
            if redeem_script != nothing
                script_sig = redeem_script
            else
                script_sig = script_pubkey(tx_in, tx.testnet)
            end
        else
            script_sig = nothing
        end
        alt_tx_in = TxIn(tx_in.prev_tx,
                         tx_in.prev_index,
                         script_sig,
                         tx_in.sequence)

        append!(s, serialize(alt_tx_in))
    end
    append!(s, encode_varint(length(tx.tx_outs)))
    for tx_out in tx.tx_outs
        append!(s, serialize(tx_out))
    end
    append!(s, Vector(reinterpret(UInt8, [htol(tx.locktime)])))
    append!(s, Vector(reinterpret(UInt8, [htol(SIGHASH_ALL)])))
    return hash256(s)
end

function hash_prevouts(tx::Tx)
    if tx._hash_prevouts == nothing
        all_prevouts = UInt8[]
        all_sequence = UInt8[]
        for tx_in in tx.tx_ins
            append!(all_prevouts, reverse!(copy(tx_in.prev_tx)))
            append!(all_prevouts, reinterpret(UInt8, [htol(tx_in.prev_index)]))
            append!(all_sequence, reinterpret(UInt8, [htol(tx_in.sequence)]))
            tx._hash_prevouts = hash256(all_prevouts)
            tx._hash_sequence = hash256(all_sequence)
        end
    end
    return tx._hash_prevouts
end

function hash_sequence(tx::Tx)
    if tx._hash_sequence == nothing
        hash_prevouts(tx)
    end
    return tx._hash_sequence
end

function hash_outputs(tx::Tx)
    if tx._hash_outputs == nothing
        all_outputs = UInt8[]
        for tx_out in tx.tx_outs
            append!(all_outputs, serialize(tx_out))
        end
        tx._hash_outputs = hash256(all_outputs)
    end
    return tx._hash_outputs
end

"""
Returns the integer representation of the hash that needs to get
signed for index input_index
"""
function sig_hash_bip143(tx::Tx, input_index::Integer; redeem_script::Union{Script,Nothing}=nothing, witness_script::Union{Script,Nothing}=nothing)
    tx_in = tx.tx_ins[input_index+1]
    # per BIP143 spec
    s = Vector(reinterpret(UInt8, [htol(tx.version)]))
    append!(s, hash_prevouts(tx))
    append!(s, hash_sequence(tx))
    append!(s, reverse!(copy(tx_in.prev_tx)))
    append!(s, reinterpret(UInt8, [htol(tx_in.prev_index)]))

    if witness_script != nothing
        script_code = serialize(witness_script)
    elseif redeem_script != nothing
        script_code = serialize(p2pkh_script(redeem_script.instructions[2]))
    else
        script_code = serialize(p2pkh_script(script_pubkey(tx_in, tx.testnet).instructions[2]))
    end
    append!(s, script_code)
    append!(s, reinterpret(UInt8, [htol(value(tx_in, tx.testnet))]))
    append!(s, reinterpret(UInt8, [htol(tx_in.sequence)]))
    append!(s, hash_outputs(tx))
    append!(s, reinterpret(UInt8, [htol(tx.locktime)]))
    append!(s, reinterpret(UInt8, [htol(SIGHASH_ALL)]))

    return hash256(s)
end

"""
    verify(tx::Tx, input_index) -> Bool

Returns whether the input has a valid signature
"""
function verify(tx::Tx, input_index)
    tx_in = tx.tx_ins[input_index+1]
    script_pubkey_ = script_pubkey(tx_in, tx.testnet)
    if is_p2sh(script_pubkey_)
        raw_redeem = copy(tx_in.script_sig.instructions[end])
        length_ = UInt8(length(raw_redeem))
        pushfirst!(raw_redeem, length_)
        redeem_script = scriptparse(IOBuffer(raw_redeem))
        if is_p2wpkh(redeem_script)
            z = sig_hash_bip143(tx, input_index, redeem_script=redeem_script)
            witness = tx_in.witness
        elseif is_p2wsh(redeem_script)
            raw_witness = copy(tx_in.witness.instructions[end])
            length_ = encode_varint(length(raw_witness))
            prepend!(raw_witness, length_)
            witness_script = scriptparse(IOBuffer(raw_witness))
            z = sig_hash_bip143(tx, input_index, witness_script=witness_script)
            witness = tx_in.witness
        else
            z = sig_hash(tx, input_index, redeem_script)
            witness = nothing
        end
    else
        if is_p2wpkh(script_pubkey_)
            z = sig_hash_bip143(tx, input_index)
            witness = tx_in.witness
        elseif is_p2wsh(script_pubkey_)
            raw_witness = copy(tx_in.witness.instructions[end])
            length_ = encode_varint(length(raw_witness))
            prepend!(raw_witness, length_)
            witness_script = scriptparse(IOBuffer(raw_witness))
            z = sig_hash_bip143(tx, input_index, witness_script=witness_script)
            witness = tx_in.witness
        else
            z = sig_hash(tx, input_index)
            witness = nothing
        end
    end
    combined_script = Script(copy(tx_in.script_sig.instructions))
    append!(combined_script.instructions, script_pubkey(tx_in, tx.testnet).instructions)
    return evaluate(combined_script, to_int(z), witness)
end

"""
    verify(tx::Tx) -> Bool

Verify transaction `tx`
"""
function verify(tx::Tx)
    if fee(tx) < 0
        return false
    end
    for i in 1:length(tx.tx_ins)
        if !verify(tx, i - 1)
            return false
        end
    end
    return true
end

"""
Signs the input using the private key
"""
function txsigninput(tx::Tx, input_index::Integer, keypair::KeyPair)
    z = to_int(sig_hash(tx, input_index))
    sig = ECDSA.sign(keypair, z)
    txpushsignature(tx, input_index, z, sig, keypair.𝑄)
end

"""
Append Signature to the Script Pubkey of TxIn at index
"""
function txpushsignature(tx::Tx, input_index::Integer, z::Integer, sig::ECDSA.Signature, pubkey::Secp256k1.Point)
    der = serialize(sig)
    append!(der, bytes(SIGHASH_ALL))
    sec = serialize(pubkey)
    script_sig = Script([der, sec])
    tx.tx_ins[input_index + 1].script_sig = script_sig
    return verify(tx, input_index)
end

"""
Returns whether this transaction is a coinbase transaction or not
"""
function iscoinbase(tx::Tx)
    if length(tx.tx_ins) != 1
        return false
    end
    input = tx.tx_ins[1]
    if input.prev_tx != fill(0x00, 32) || input.prev_index != 0xffffffff
        return false
    end
    return true
end

"""
Returns the height of the block this coinbase transaction is in
Returns `nothing` if this transaction is not a coinbase transaction
"""
function coinbase_height(tx::Tx)
    if !iscoinbase(tx)
        return nothing
    end
    height_bytes = tx.tx_ins[1].script_sig.instructions[1]
    return to_int(height_bytes, little_endian=true)
end

@deprecate txparse(s::IOBuffer, testnet::Bool) parse(s::IOBuffer, testnet::Bool)::Tx
@deprecate txserialize(tx::Tx) serialize_legacy(tx::Tx)
@deprecate txoutserialize(tx::TxOut) serialize(tx::TxOut)
@deprecate txhash(tx::Tx) hash(tx::Tx)
@deprecate txinputverify(tx::Tx, input_index) verify(tx::Tx, input_index)
@deprecate txinvalue(txin::TxIn, testnet::Bool) value(txin::TxIn, testnet::Bool)
@deprecate txverify(tx::Tx) verify(tx::Tx)
@deprecate txid(tx::Tx) id(tx::Tx)
@deprecate txfee(tx::Tx) fee(tx::Tx)
@deprecate txsighash256(tx::Tx, input_index::Integer) sig_hash(tx::Tx, input_index::Integer)
@deprecate txin_scriptpubkey(txin::TxIn, testnet::Bool) script_pubkey(txin::TxIn, testnet::Bool)
@deprecate txinserialize(tx::TxIn) serialize(tx::TxIn)
@deprecate txin_fetchtx(tx::TxIn, testnet::Bool) fetch(tx::TxIn, testnet::Bool)
@deprecate txinparse(s::IOBuffer) TxIn(s::IOBuffer)


text:
####################### Graphs in the paper #######################
# Useful for testing whether the series generated by different solvers/algorithms match matlab's or not.

#Real GDP
GDP_plot=plot(1:p.N,[0 (sim_fun.GDP_Laspeyres[2:p.N]./sim_0.GDP .-1)']',title= "GDP Laspeyres baseline weighted 1 tau",legend=false)
savefig(GDP_plot,"GDP_Julia.png")

#Capital
Capital_plot=plot(1:p.N,[0 (sim_fun.K[2:p.N]./sim_0.K.-1)']', title="K baseline weighted 1 tau",legend=false)
savefig(Capital_plot,"Capital_Julia.png")

#Consumption
plot(1:p.N,[0 (sim_fun.C[2:p.N]./sim_0.C.-1)']',title="C baseline weighted 1 tau")

#Real Exports
plot(1:p.N,[0 (sim_fun.X_Laspeyres[2:p.N]./sim_0.PxYx.-1)']',title="X Laspeyres baseline weighted 1 tau")

#Real Domestic sales
plot(1:p.N,[0 (sim_fun.D_Laspeyres[2:p.N]./sim_0.PdYd.-1)']',title="D Laspeyres baseline weighted 1 tau")

#Real Exchange Rate
plot(1:p.N,[0 (sim_fun.ξ[2:p.N]./sim_0.ξ.-1)']',title="RER baseline weighted 1 tau")

#Real Wage
plot(1:p.N,[0 (sim_fun.w[2:p.N]./sim_0.w.-1)']',title=
"Wage baseline weighted 1 tau")

#Price of Capital
plot(1:p.N,[0 (sim_fun.Pk[2:p.N]./sim_0.Pk.-1)']',title="Pk baseline weighted 1 tau")

#Phi_h
plot(1:p.N,[0 (sim_fun.ϕ_h[2:p.N]./sim_0.ϕ_h.-1)']',title="Phih baseline weighted 1 tau")

#### MATLAB'S #####

#Real GDP
plot_gdp=plot(1:p.N,[0 (sim_fun_matlab.GDP_Laspeyres[2:p.N]./sim_0.GDP .-1)']',title= "GDP Laspeyres baseline weighted 1 tau")
savefig(plot_gdp,"GDP_Julia_MatlabGuesses2.png")
#Capital
plot_capital=plot(1:p.N,[0 (sim_fun_matlab.K[2:p.N]./sim_0.K.-1)']', title="K baseline weighted 1 tau")
savefig(plot_capital,"Capital_Julia_MatlabGuesses2.png")


text:
export Pade21

"""
    Pade21()

This uses a Padé-approximation of type (2,1) for prediction.
"""
struct Pade21 <: AbstractPredictor end
struct Pade21Cache{T, AV<:AbstractVector{T}} <: AbstractPredictorCache
    x²::Vector{T}
    x³::Vector{T}
    x_h::AV

    u::Vector{T}
    u₁::Vector{T}
    u₂::Vector{T}

    h₂::Float64
    h₃::Float64
end

function cache(::Pade21, H, x, ẋ, t)
    x², x³, x_h = copy(ẋ), copy(ẋ), copy(x)
    u = evaluate(H, x, t)
    u₁, u₂ = copy(u), copy(u)
    h₂ = nthroot(eps(), 2+2) # x² is the second order derivative and has a 2nd order approximation
    h₃ = nthroot(eps(), 3+2) # x³ is the third order derivative and has a 2nd order approximation
    Pade21Cache(x², x³, x_h, u, u₁, u₂, h₂, h₃)
end
@inline function g!(u, H, x, ẋ, t, h, x_h)
    @inbounds for i in eachindex(x)
        x_h[i] = x[i] + h * ẋ[i]
    end
    evaluate!(u, H, x_h, t + h)
end

@inline function g!(u, H, x, ẋ, ẍ, t, h, x_h)
    h² = 0.5 * h * h
    @inbounds for i in eachindex(x)
        x_h[i] = x[i] + h * ẋ[i] + h² * ẍ[i]
    end
    evaluate!(u, H, x_h, t + h)
end

function update!(cache::Pade21Cache, H, x, ẋ, t, Jac::Jacobian)
    # unpack stuff to make the rest easier to read
    u, u₁, u₂ = cache.u, cache.u₁, cache.u₂
    x_h, h₂, h₃ = cache.x_h, cache.h₂, cache.h₃
    x², x³ = cache.x², cache.x³

    h₂ *= max(abs(t), 1.0)
    # compute the second derivative at (x,t)
    g!(u₁, H, x, ẋ, t, h₂, x_h)
    g!(u₂, H, x, ẋ, t, -h₂, x_h)

    h² = h₂^2
    @inbounds for i in eachindex(u₁)
        u[i] = -(u₁[i] + u₂[i]) / h²
    end
    solve!(x², Jac, u)

    h₃ *= max(abs(t), 1.0)
    g!(u₁, H, x, ẋ, x², t, h₃, x_h)
    g!(u₂, H, x, ẋ, x², t, -h₃, x_h)

    h³ = h₃^3
    @inbounds for i in eachindex(u₁)
        u[i] = -3(u₁[i] - u₂[i]) / h³
    end
    solve!(x³, Jac, u)

    cache
end

function predict!(xnext, ::Pade21, cache::Pade21Cache, H::HomotopyWithCache, x, t, Δt, ẋ, Jac::Jacobian)
    x², x³ = cache.x², cache.x³
    @inbounds for i in eachindex(x)
        δ = 1 - Δt * x³[i] / (3 * x²[i])
        if isnan(δ) || iszero(δ)
            xnext[i] = x[i] + Δt * ẋ[i] + 0.5 * Δt^2 * x²[i]
        else
            xnext[i] = x[i] + Δt * ẋ[i] + (0.5 * Δt^2 * x²[i]) / δ
        end
    end
    nothing
end

order(::Pade21) = 4


text:
using TodoGrabber

function main()
    filenames = TodoGrabber.get_filenames("./tests")
    println("The filenames are:", filenames)

    todos = TodoGrabber.find_todos(filenames)
    TodoGrabber.write_todos(todos, "./tests/todostack.org")
    TodoGrabber.save_dict(todos, "serial.bin")
end

function maindb()
    if isfile("todo_cache.db")
        rm("todo_cache.db")
    end
    filenames = TodoGrabber.get_filenames("./tests")
    todos = TodoGrabber.find_todos(filenames)
    storetodos(todos)
end


text:
# This file was generated, do not modify it. # hide
perf_k = misclassification_rate(predict_mode(knn, Xtest), ytest)
perf_m = misclassification_rate(predict_mode(multi, Xtest), ytest)
println(rpad("KNN mcr:", 10), round(perf_k, sigdigits=3))
println(rpad("MNC mcr:", 10), round(perf_m, sigdigits=3))