text:
# Julia wrapper for header: /usr/include/scip/pub_heur.h
# Automatically generated using Clang.jl wrap_c

function SCIPheurComp(elem1, elem2)
    ccall((:SCIPheurComp, libscip), Cint, (Ptr{Cvoid}, Ptr{Cvoid}), elem1, elem2)
end

function SCIPheurCompName(elem1, elem2)
    ccall((:SCIPheurCompName, libscip), Cint, (Ptr{Cvoid}, Ptr{Cvoid}), elem1, elem2)
end

function SCIPheurGetData(heur)
    ccall((:SCIPheurGetData, libscip), Ptr{SCIP_HEURDATA}, (Ptr{SCIP_HEUR},), heur)
end

function SCIPheurSetData(heur, heurdata)
    ccall((:SCIPheurSetData, libscip), Cvoid, (Ptr{SCIP_HEUR}, Ptr{SCIP_HEURDATA}), heur, heurdata)
end

function SCIPheurGetName(heur)
    ccall((:SCIPheurGetName, libscip), Cstring, (Ptr{SCIP_HEUR},), heur)
end

function SCIPheurGetDesc(heur)
    ccall((:SCIPheurGetDesc, libscip), Cstring, (Ptr{SCIP_HEUR},), heur)
end

function SCIPheurGetDispchar(heur)
    ccall((:SCIPheurGetDispchar, libscip), UInt8, (Ptr{SCIP_HEUR},), heur)
end

function SCIPheurGetTimingmask(heur)
    ccall((:SCIPheurGetTimingmask, libscip), SCIP_HEURTIMING, (Ptr{SCIP_HEUR},), heur)
end

function SCIPheurSetTimingmask(heur, timingmask)
    ccall((:SCIPheurSetTimingmask, libscip), Cvoid, (Ptr{SCIP_HEUR}, SCIP_HEURTIMING), heur, timingmask)
end

function SCIPheurUsesSubscip(heur)
    ccall((:SCIPheurUsesSubscip, libscip), UInt32, (Ptr{SCIP_HEUR},), heur)
end

function SCIPheurGetPriority(heur)
    ccall((:SCIPheurGetPriority, libscip), Cint, (Ptr{SCIP_HEUR},), heur)
end

function SCIPheurGetFreq(heur)
    ccall((:SCIPheurGetFreq, libscip), Cint, (Ptr{SCIP_HEUR},), heur)
end

function SCIPheurSetFreq(heur, freq)
    ccall((:SCIPheurSetFreq, libscip), Cvoid, (Ptr{SCIP_HEUR}, Cint), heur, freq)
end

function SCIPheurGetFreqofs(heur)
    ccall((:SCIPheurGetFreqofs, libscip), Cint, (Ptr{SCIP_HEUR},), heur)
end

function SCIPheurGetMaxdepth(heur)
    ccall((:SCIPheurGetMaxdepth, libscip), Cint, (Ptr{SCIP_HEUR},), heur)
end

function SCIPheurGetNCalls(heur)
    ccall((:SCIPheurGetNCalls, libscip), Clonglong, (Ptr{SCIP_HEUR},), heur)
end

function SCIPheurGetNSolsFound(heur)
    ccall((:SCIPheurGetNSolsFound, libscip), Clonglong, (Ptr{SCIP_HEUR},), heur)
end

function SCIPheurGetNBestSolsFound(heur)
    ccall((:SCIPheurGetNBestSolsFound, libscip), Clonglong, (Ptr{SCIP_HEUR},), heur)
end

function SCIPheurIsInitialized(heur)
    ccall((:SCIPheurIsInitialized, libscip), UInt32, (Ptr{SCIP_HEUR},), heur)
end

function SCIPheurGetSetupTime(heur)
    ccall((:SCIPheurGetSetupTime, libscip), Cdouble, (Ptr{SCIP_HEUR},), heur)
end

function SCIPheurGetTime(heur)
    ccall((:SCIPheurGetTime, libscip), Cdouble, (Ptr{SCIP_HEUR},), heur)
end

function SCIPheurGetDivesets(heur)
    ccall((:SCIPheurGetDivesets, libscip), Ptr{Ptr{SCIP_DIVESET}}, (Ptr{SCIP_HEUR},), heur)
end

function SCIPheurGetNDivesets(heur)
    ccall((:SCIPheurGetNDivesets, libscip), Cint, (Ptr{SCIP_HEUR},), heur)
end

function SCIPdivesetGetHeur(diveset)
    ccall((:SCIPdivesetGetHeur, libscip), Ptr{SCIP_HEUR}, (Ptr{SCIP_DIVESET},), diveset)
end

function SCIPdivesetGetWorkSolution(diveset)
    ccall((:SCIPdivesetGetWorkSolution, libscip), Ptr{SCIP_SOL}, (Ptr{SCIP_DIVESET},), diveset)
end

function SCIPdivesetSetWorkSolution(diveset, sol)
    ccall((:SCIPdivesetSetWorkSolution, libscip), Cvoid, (Ptr{SCIP_DIVESET}, Ptr{SCIP_SOL}), diveset, sol)
end

function SCIPdivesetGetName(diveset)
    ccall((:SCIPdivesetGetName, libscip), Cstring, (Ptr{SCIP_DIVESET},), diveset)
end

function SCIPdivesetGetMinRelDepth(diveset)
    ccall((:SCIPdivesetGetMinRelDepth, libscip), Cdouble, (Ptr{SCIP_DIVESET},), diveset)
end

function SCIPdivesetGetMaxRelDepth(diveset)
    ccall((:SCIPdivesetGetMaxRelDepth, libscip), Cdouble, (Ptr{SCIP_DIVESET},), diveset)
end

function SCIPdivesetGetSolSuccess(diveset)
    ccall((:SCIPdivesetGetSolSuccess, libscip), Clonglong, (Ptr{SCIP_DIVESET},), diveset)
end

function SCIPdivesetGetNCalls(diveset)
    ccall((:SCIPdivesetGetNCalls, libscip), Cint, (Ptr{SCIP_DIVESET},), diveset)
end

function SCIPdivesetGetNSolutionCalls(diveset)
    ccall((:SCIPdivesetGetNSolutionCalls, libscip), Cint, (Ptr{SCIP_DIVESET},), diveset)
end

function SCIPdivesetGetMinDepth(diveset)
    ccall((:SCIPdivesetGetMinDepth, libscip), Cint, (Ptr{SCIP_DIVESET},), diveset)
end

function SCIPdivesetGetMaxDepth(diveset)
    ccall((:SCIPdivesetGetMaxDepth, libscip), Cint, (Ptr{SCIP_DIVESET},), diveset)
end

function SCIPdivesetGetAvgDepth(diveset)
    ccall((:SCIPdivesetGetAvgDepth, libscip), Cdouble, (Ptr{SCIP_DIVESET},), diveset)
end

function SCIPdivesetGetMinSolutionDepth(diveset)
    ccall((:SCIPdivesetGetMinSolutionDepth, libscip), Cint, (Ptr{SCIP_DIVESET},), diveset)
end

function SCIPdivesetGetMaxSolutionDepth(diveset)
    ccall((:SCIPdivesetGetMaxSolutionDepth, libscip), Cint, (Ptr{SCIP_DIVESET},), diveset)
end

function SCIPdivesetGetAvgSolutionDepth(diveset)
    ccall((:SCIPdivesetGetAvgSolutionDepth, libscip), Cdouble, (Ptr{SCIP_DIVESET},), diveset)
end

function SCIPdivesetGetNLPIterations(diveset)
    ccall((:SCIPdivesetGetNLPIterations, libscip), Clonglong, (Ptr{SCIP_DIVESET},), diveset)
end

function SCIPdivesetGetNProbingNodes(diveset)
    ccall((:SCIPdivesetGetNProbingNodes, libscip), Clonglong, (Ptr{SCIP_DIVESET},), diveset)
end

function SCIPdivesetGetNBacktracks(diveset)
    ccall((:SCIPdivesetGetNBacktracks, libscip), Clonglong, (Ptr{SCIP_DIVESET},), diveset)
end

function SCIPdivesetGetNConflicts(diveset)
    ccall((:SCIPdivesetGetNConflicts, libscip), Clonglong, (Ptr{SCIP_DIVESET},), diveset)
end

function SCIPdivesetGetNSols(diveset)
    ccall((:SCIPdivesetGetNSols, libscip), Clonglong, (Ptr{SCIP_DIVESET},), diveset)
end

function SCIPdivesetGetMaxLPIterQuot(diveset)
    ccall((:SCIPdivesetGetMaxLPIterQuot, libscip), Cdouble, (Ptr{SCIP_DIVESET},), diveset)
end

function SCIPdivesetGetMaxLPIterOffset(diveset)
    ccall((:SCIPdivesetGetMaxLPIterOffset, libscip), Cint, (Ptr{SCIP_DIVESET},), diveset)
end

function SCIPdivesetGetUbQuotNoSol(diveset)
    ccall((:SCIPdivesetGetUbQuotNoSol, libscip), Cdouble, (Ptr{SCIP_DIVESET},), diveset)
end

function SCIPdivesetGetAvgQuotNoSol(diveset)
    ccall((:SCIPdivesetGetAvgQuotNoSol, libscip), Cdouble, (Ptr{SCIP_DIVESET},), diveset)
end

function SCIPdivesetGetUbQuot(diveset)
    ccall((:SCIPdivesetGetUbQuot, libscip), Cdouble, (Ptr{SCIP_DIVESET},), diveset)
end

function SCIPdivesetGetAvgQuot(diveset)
    ccall((:SCIPdivesetGetAvgQuot, libscip), Cdouble, (Ptr{SCIP_DIVESET},), diveset)
end

function SCIPdivesetUseBacktrack(diveset)
    ccall((:SCIPdivesetUseBacktrack, libscip), UInt32, (Ptr{SCIP_DIVESET},), diveset)
end

function SCIPdivesetGetLPSolveFreq(diveset)
    ccall((:SCIPdivesetGetLPSolveFreq, libscip), Cint, (Ptr{SCIP_DIVESET},), diveset)
end

function SCIPdivesetGetLPResolveDomChgQuot(diveset)
    ccall((:SCIPdivesetGetLPResolveDomChgQuot, libscip), Cdouble, (Ptr{SCIP_DIVESET},), diveset)
end

function SCIPdivesetUseOnlyLPBranchcands(diveset)
    ccall((:SCIPdivesetUseOnlyLPBranchcands, libscip), UInt32, (Ptr{SCIP_DIVESET},), diveset)
end

function SCIPdivesetSupportsType(diveset, divetype)
    ccall((:SCIPdivesetSupportsType, libscip), UInt32, (Ptr{SCIP_DIVESET}, SCIP_DIVETYPE), diveset, divetype)
end

function SCIPdivesetGetRandnumgen(diveset)
    ccall((:SCIPdivesetGetRandnumgen, libscip), Ptr{SCIP_RANDNUMGEN}, (Ptr{SCIP_DIVESET},), diveset)
end

function SCIPvariablegraphBreadthFirst(scip, vargraph, startvars, nstartvars, distances, maxdistance, maxvars, maxbinintvars)
    ccall((:SCIPvariablegraphBreadthFirst, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_VGRAPH}, Ptr{Ptr{SCIP_VAR}}, Cint, Ptr{Cint}, Cint, Cint, Cint), scip, vargraph, startvars, nstartvars, distances, maxdistance, maxvars, maxbinintvars)
end

function SCIPvariableGraphCreate(scip, vargraph, relaxdenseconss, relaxdensity, nrelaxedconstraints)
    ccall((:SCIPvariableGraphCreate, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_VGRAPH}}, UInt32, Cdouble, Ptr{Cint}), scip, vargraph, relaxdenseconss, relaxdensity, nrelaxedconstraints)
end

function SCIPvariableGraphFree(scip, vargraph)
    ccall((:SCIPvariableGraphFree, libscip), Cvoid, (Ptr{SCIP_}, Ptr{Ptr{SCIP_VGRAPH}}), scip, vargraph)
end


text:
## general data and market data =================================
seed = 2            ## seed for random generator
n_scen = 1_000_000  ## number of monte carlo scenarios
α = 0.99            ## confidence level for expected shortfall
coc_rate = 0.06     ## cost of capital rate for Mindestbetrag
spot = [0.00155, 0.00013, -0.00011]  ## spot rates as of 2012
stock_increase = 0.05   ## expected rel. increase of stock index

## risk factors =================================================

## CHF spot rates
x0_spot = [0.0, 0.0, 0.0]  ## as of 2012
h_spot = [0.01, 0.01, 0.01]
add_spot = [true, true, true]
index_spot = collect(1:3) # Arg. must always be of the form 1:T
## stock index
x0_stock = 1.0
h_stock = 0.1
add_stock = false
index_stock = [4]
## mortality
x0_mort = 1.0
h_mort = 0.1
add_mort = false
index_mort = [5]

## data provied by FINMA (for 2012) -----------------------------
## components link to index_spot, index_stock, index_mort

## standard deviations of risk factors
σ = [0.0060292643195745,         # spot_chf[1]
     0.00605990680656078,        # spot_chf[2]
     0.00633152336734695,        # spot_chf[3]
     0.150517788474455,          # stock
     0.05 ]                      # q

## correlations for risk factors
##       spot_chf[1] spot_chf[2] spot_chf[3] stock       mort
corr =  [1.000000000 0.721560179 0.545556472 0.402239567 0;
         0.721560179 1.000000000 0.953194547 0.433388849 0;
         0.545556472 0.953194547 1.000000000 0.413041139 0;
         0.402239567 0.433388849 0.413041139 1.000000000 0;
         0.000000000 0.000000000 0.000000000 0.000000000 1]

## stress scenarios

## data frame that describes the stresst scenarios to be added
## n:       number of scenarios
## name:    stress scenario name / description
## target:  true: include, false: do not include
## prob:    probability for this stress scenario occuring
## Δx:      effect of this scenario on risk factor
## Δrtk:    efect of this scenario on RTK (to be calculated)
stress =
  Stress(8,
         ["SZ01: Equity drop -60%",
          "SZ03: Stock mkt crash (87)",
          "SZ04: Nikkei crash (89/90)",
          "SZ05: European currency crisis (92)",
          "SZ06: US interest crisis (94)",
          "SZ07: Russian crisis/LTCM (98)",
          "SZ08: Stock mkt crash (00/01)",
          "SZ11: Financial crisis (08)"],
         [true, true, true, true, true, true, true, true],
         [0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001],
         [ 0.00000  0.00000 0.00000 -0.6000 0;
          -0.00155 -0.00013 0.00000 -0.2323 0;
           0.01563  0.01098 0.01177 -0.2643 0;
          -0.00155 -0.00013 0.00000 -0.0580 0;
           0.01109  0.01406 0.01509 -0.1852 0;
          -0.00155 -0.00013 0.00000 -0.2841 0;
          -0.00155 -0.00013 0.00000 -0.3567 0;
          -0.00155 -0.00013 0.00000 -0.3881 0],
         zeros(Float64, 8))

## portfolio data  ==============================================
invest =
  DataFrame(maturity = [NaN, 1, 2, 3],
            index = [4, 1, 2, 3],
            kind = ["stock", "zerobond", "zerobond", "zerobond"],
            nominal = [300.0, 0.0, 300.0, 500.0])

n_ins = 10                         ## number insurance contracts
ins_sum = 100.0                    ## insured sum per contract
qx_flat = 0.02                     ## mortality probability
qx = [qx_flat, qx_flat, qx_flat]   ## mortality probabilities
B_PX  =[0.0, 0.0, n_ins * ins_sum] ## benefit profile (survival)
T = length(B_PX)                   ## Time until run off

ext: jl


text:
using DrWatson
@quickactivate
using ArgParse
using GenerativeAD
import StatsBase: fit!, predict
using StatsBase
using BSON, FileIO
using Flux
using GenerativeModels
using DistributionsAD
using ValueHistories

s = ArgParseSettings()
@add_arg_table! s begin
	"modelname"
		default = "vae"
		arg_type = String
		help = "model name"
	"datatype"
		default = "tabular"
		arg_type = String
		help = "tabular or image"
	"dataset"
		default = "iris"
		arg_type = String
		help = "dataset"
	"--seed"
		default = nothing
		help = "if specified, only results for a given seed will be recomputed"
	"--anomaly_class"
		default = nothing
		help = "if specified, only results for a given anomaly class will be recomputed"
end
parsed_args = parse_args(ARGS, s)
@unpack dataset, datatype, modelname, seed, anomaly_class = parsed_args

masterpath = datadir("experiments/$(datatype)/$(modelname)/$(dataset)")
files = GenerativeAD.Evaluation.collect_files(masterpath)
mfiles = filter(f->occursin("model", f), files)
if seed != nothing
	filter!(x->occursin("/seed=$seed/", x), mfiles)
end
if anomaly_class != nothing
	filter!(x->occursin("/ac=$(anomaly_class)/", x), mfiles)
end

kld_batched(m,x,batchsize) =
	vcat(map(y->vec(Base.invokelatest(GenerativeModels.kl_divergence,
			condition(m.encoder,y), m.prior)), Flux.Data.DataLoader(x, batchsize=batchsize))...)
kld_batched_gpu(m,x,batchsize) =
	vcat(map(y-> cpu(vec(Base.invokelatest(GenerativeModels.kl_divergence,
			Base.invokelatest(condition, m.encoder, gpu(Array(y))), m.prior))), Flux.Data.DataLoader(x, batchsize=batchsize))...)

function save_elbo_score(f::String, data, seed::Int, ac=nothing)
	# get model
	savepath = dirname(f)
	mdata = load(f)
	model = mdata["model"]

	# setup entries to be saved
	save_entries = (
		modelname = modelname,
		fit_t = mdata["fit_t"],
		history = mdata["history"],
		dataset = dataset,
		npars = sum(map(p->length(p), Flux.params(model))),
		model = nothing,
		seed = seed
		)
	save_entries = (ac == nothing) ? save_entries : merge(save_entries, (anomaly_class=ac,))

	# test if the file alread exists
	parameters = merge(mdata["parameters"], (L=100, score = "elbo",))
	savef = joinpath(savepath, savename(parameters, "bson", digits=5))
	if !isfile(savef)
		@info "computing sample score for $f"

		# now load the sampled reconstruction probability score
		namepars = DrWatson.parse_savename(f)[2]
		sf = filter(x->
			occursin("init_seed=$(namepars["init_seed"])", x) &&
			occursin("score=reconstruction-sampled", x),
			readdir(savepath))
		if length(sf)==0
			@info "sampled reconstruction score for $f not found!"
			return
		end
		sf = joinpath(savepath, sf[1])
		score_data = load(sf)
		rec_scores = (score_data[:tr_scores], score_data[:val_scores], score_data[:tst_scores])

		# compute the full scores
		score_fun(x) = (ac == nothing) ?
			kld_batched(model, x, 512) : kld_batched_gpu(gpu(model), x, 512)
		tr_data, val_data, tst_data = data

		# extract scores
		tr_scores, tr_eval_t, _, _, _ = @timed score_fun(tr_data[1])
		val_scores, val_eval_t, _, _, _ = @timed score_fun(val_data[1])
		tst_scores, tst_eval_t, _, _, _ = @timed score_fun(tst_data[1])

		result = (
			parameters = parameters,
			tr_scores = tr_scores .+ rec_scores[1],
			tr_labels = tr_data[2],
			tr_eval_t = tr_eval_t + score_data[:tr_eval_t],
			val_scores = val_scores .+ rec_scores[2],
			val_labels = val_data[2],
			val_eval_t = val_eval_t + score_data[:val_eval_t],
			tst_scores = tst_scores .+ rec_scores[3],
			tst_labels = tst_data[2],
			tst_eval_t = tst_eval_t + score_data[:tst_eval_t]
			)

		result = Dict{Symbol, Any}([sym=>val for (sym,val) in pairs(merge(result, save_entries))]) # this has to be a Dict
		tagsave(savef, result, safe = true)
		(@info "Results saved to $savef")
	end
end

for f in mfiles
	# get data
	savepath = dirname(f)
	local seed = parse(Int, replace(basename(savepath), "seed=" => ""))
	ac = occursin("ac=", savepath) ? parse(Int, replace(basename(dirname(savepath)), "ac=" => "")) : nothing
	data = (ac == nothing) ?
		GenerativeAD.load_data(dataset, seed=seed) :
		GenerativeAD.load_data(dataset, seed=seed, anomaly_class_ind=ac)

	# compute and save the score
	save_elbo_score(f, data, seed, ac)
end
@info "DONE"


text:
using ..Ocean: kinematic_stress

"""
    ocean_velocity_boundary_state!(::NumericalFluxFirstOrder, ::Impenetrable{NoSlip}, ::HBModel)

apply no slip boundary condition for velocity
sets reflective ghost point
"""
function ocean_velocity_boundary_state!(
    nf::NumericalFluxFirstOrder,
    bc_velocity::Impenetrable{NoSlip},
    ocean,
    Q⁺,
    A⁺,
    n⁻,
    Q⁻,
    A⁻,
    t,
)
    Q⁺.u = -Q⁻.u
    A⁺.w = -A⁻.w

    return nothing
end

"""
    ocean_velocity_boundary_state!(::NumericalFluxGradient, ::Impenetrable{NoSlip}, ::HBModel)

apply no slip boundary condition for velocity
set numerical flux to zero for u
"""
function ocean_velocity_boundary_state!(
    nf::NumericalFluxGradient,
    bc_velocity::Impenetrable{NoSlip},
    ocean,
    Q⁺,
    A⁺,
    n⁻,
    Q⁻,
    A⁻,
    t,
)
    FT = eltype(Q⁺)
    Q⁺.u = SVector(-zero(FT), -zero(FT))
    A⁺.w = -zero(FT)

    return nothing
end

"""
    ocean_velocity_boundary_state!(::NumericalFluxSecondOrder, ::Impenetrable{NoSlip}, ::HBModel)

apply no slip boundary condition for velocity
sets ghost point to have no numerical flux on the boundary for u
"""
@inline function ocean_velocity_boundary_state!(
    nf::NumericalFluxSecondOrder,
    bc_velocity::Impenetrable{NoSlip},
    ocean,
    Q⁺,
    D⁺,
    A⁺,
    n⁻,
    Q⁻,
    D⁻,
    A⁻,
    t,
)
    Q⁺.u = -Q⁻.u
    A⁺.w = -A⁻.w
    D⁺.ν∇u = D⁻.ν∇u

    return nothing
end

"""
    ocean_velocity_boundary_state!(::NumericalFluxFirstOrder, ::Impenetrable{FreeSlip}, ::HBModel)

apply free slip boundary condition for velocity
sets reflective ghost point
"""
function ocean_velocity_boundary_state!(
    nf::NumericalFluxFirstOrder,
    bc_velocity::Impenetrable{FreeSlip},
    ocean,
    Q⁺,
    A⁺,
    n⁻,
    Q⁻,
    A⁻,
    t,
)
    v⁻ = @SVector [Q⁻.u[1], Q⁻.u[2], A⁻.w]
    v⁺ = v⁻ - 2 * n⁻ ⋅ v⁻ .* SVector(n⁻)
    Q⁺.u = @SVector [v⁺[1], v⁺[2]]
    A⁺.w = v⁺[3]

    return nothing
end

"""
    ocean_velocity_boundary_state!(::NumericalFluxGradient, ::Impenetrable{FreeSlip}, ::HBModel)

apply free slip boundary condition for velocity
sets non-reflective ghost point
"""
function ocean_velocity_boundary_state!(
    nf::NumericalFluxGradient,
    bc_velocity::Impenetrable{FreeSlip},
    ocean,
    Q⁺,
    A⁺,
    n⁻,
    Q⁻,
    A⁻,
    t,
)
    v⁻ = @SVector [Q⁻.u[1], Q⁻.u[2], A⁻.w]
    v⁺ = v⁻ - n⁻ ⋅ v⁻ .* SVector(n⁻)
    Q⁺.u = @SVector [v⁺[1], v⁺[2]]
    A⁺.w = v⁺[3]

    return nothing
end

"""
    ocean_velocity_normal_boundary_flux_second_order!(::NumericalFluxSecondOrder, ::Impenetrable{FreeSlip}, ::HBModel)

apply free slip boundary condition for velocity
apply zero numerical flux in the normal direction
"""
function ocean_velocity_boundary_state!(
    nf::NumericalFluxSecondOrder,
    bc_velocity::Impenetrable{FreeSlip},
    ocean,
    Q⁺,
    D⁺,
    A⁺,
    n⁻,
    Q⁻,
    D⁻,
    A⁻,
    t,
)
    Q⁺.u = Q⁻.u
    A⁺.w = A⁻.w
    D⁺.ν∇u = n⁻ * (@SVector [-0, -0])'

    return nothing
end

"""
    ocean_velocity_boundary_state!(::Union{NumericalFluxFirstOrder, NumericalFluxGradient}, ::Penetrable{FreeSlip}, ::HBModel)

apply free slip boundary condition for velocity
sets non-reflective ghost point
"""
function ocean_velocity_boundary_state!(
    nf::Union{NumericalFluxFirstOrder, NumericalFluxGradient},
    bc_velocity::Penetrable{FreeSlip},
    ocean,
    args...,
)
    return nothing
end

"""
    ocean_velocity_boundary_state!(::NumericalFluxSecondOrder, ::Penetrable{FreeSlip}, ::HBModel)

apply free slip boundary condition for velocity
sets non-reflective ghost point
"""
function ocean_velocity_boundary_state!(
    nf::NumericalFluxSecondOrder,
    bc_velocity::Penetrable{FreeSlip},
    ocean,
    Q⁺,
    D⁺,
    A⁺,
    n⁻,
    Q⁻,
    D⁻,
    A⁻,
    t,
)
    Q⁺.u = Q⁻.u
    A⁺.w = A⁻.w
    D⁺.ν∇u = n⁻ * (@SVector [-0, -0])'

    return nothing
end

"""
    ocean_velocity_boundary_state!(::Union{NumericalFluxFirstOrder, NumericalFluxGradient}, ::Impenetrable{KinematicStress}, ::HBModel)

apply kinematic stress boundary condition for velocity
applies free slip conditions for first-order and gradient fluxes
"""
function ocean_velocity_boundary_state!(
    nf::Union{NumericalFluxFirstOrder, NumericalFluxGradient},
    bc_velocity::Impenetrable{<:KinematicStress},
    ocean,
    args...,
)
    return ocean_velocity_boundary_state!(
        nf,
        Impenetrable(FreeSlip()),
        ocean,
        args...,
    )
end

"""
    ocean_velocity_boundary_state!(::NumericalFluxSecondOrder, ::Impenetrable{KinematicStress}, ::HBModel)

apply kinematic stress boundary condition for velocity
sets ghost point to have specified flux on the boundary for ν∇u
"""
@inline function ocean_velocity_boundary_state!(
    nf::NumericalFluxSecondOrder,
    bc_velocity::Impenetrable{<:KinematicStress},
    ocean,
    Q⁺,
    D⁺,
    A⁺,
    n⁻,
    Q⁻,
    D⁻,
    A⁻,
    t,
)
    Q⁺.u = Q⁻.u
    D⁺.ν∇u =
        n⁻ * kinematic_stress(ocean.problem, A⁻.y, ocean.ρₒ, bc_velocity.drag)'

    return nothing
end

"""
    ocean_velocity_boundary_state!(::Union{NumericalFluxFirstOrder, NumericalFluxGradient}, ::Penetrable{KinematicStress}, ::HBModel)

apply kinematic stress boundary condition for velocity
applies free slip conditions for first-order and gradient fluxes
"""
function ocean_velocity_boundary_state!(
    nf::Union{NumericalFluxFirstOrder, NumericalFluxGradient},
    bc_velocity::Penetrable{<:KinematicStress},
    ocean,
    args...,
)
    return ocean_velocity_boundary_state!(
        nf,
        Penetrable(FreeSlip()),
        ocean,
        args...,
    )
end

"""
    ocean_velocity_boundary_state!(::NumericalFluxSecondOrder, ::Penetrable{KinematicStress}, ::HBModel)

apply kinematic stress boundary condition for velocity
sets ghost point to have specified flux on the boundary for ν∇u
"""
@inline function ocean_velocity_boundary_state!(
    nf::NumericalFluxSecondOrder,
    bc_velocity::Penetrable{<:KinematicStress},
    ocean,
    Q⁺,
    D⁺,
    A⁺,
    n⁻,
    Q⁻,
    D⁻,
    A⁻,
    t,
)
    Q⁺.u = Q⁻.u
    D⁺.ν∇u =
        n⁻ *
        kinematic_stress(
            ocean.problem,
            A⁻.y,
            ocean.ρₒ,
            bc_velocity.drag.stress,
        )'

    return nothing
end


text:
using Documenter, MshReader

makedocs( sitename="MshReader" )


text:
module Gates

using QSpice.BitOps
using QSpice.State

export identity, hadamard, not, cnot, ccnot,
       swap, cswap, sqrtswap, phaseshift, cphase,
       paulix, pauliy, pauliz,
       unitary, choose1, superposition, join, split, qft,
       measure, partialmeasure, probe

function identity(state::QuantumState)
    return copy(state)
end

function superposition(state1::QuantumState, state2::QuantumState)
    return fromstates(state1, state2)
end

function join(state1::QuantumState, state2::QuantumState)
    return QuantumState(vcat(state1.vector, state2.vector), state1.bits + state2.bits)
end

function split(state::QuantumState, lastbit::Int)
    lpos = 2^lastbit
    return QuantumState(state[1:lpos], lastbit), QuantumState(state[lpos:end], state.bits - lastbit)
end

# Performs the following mapping:
# |0> to (|0> + |1>) / sqrt(2)
# |1> to (|0> - |1>) / sqrt(2)
function hadamard(state::QuantumState, bit::Int)
    newstate = Vector{Complex{Float64}}(length(state))
    bit = state.bits - bit
    invsqrt2 = 1.0 / sqrt(2)

    for i = 1:length(state)
        basis  = i - 1
        mapped0 = clearbit(basis, bit)
        mapped1 = setbit(basis, bit)

        if iszero(basis, bit)
            newstate[i] = invsqrt2 * (state[mapped0 + 1] + state[mapped1 + 1])
        else
            newstate[i] = invsqrt2 * (state[mapped0 + 1] - state[mapped1 + 1])
        end
    end

    return QuantumState(newstate, state.bits)
end

# Also called the Pauli-X gate. Performs the following mapping:
# |0> to |1>
# |1> to |0>
function not(state::QuantumState, bit::Int)
    newstate = Vector{Complex{Float64}}(length(state))
    bit = state.bits - bit

    for i = 1:length(state)
        basis = i - 1
        mapped = flipbit(basis, bit)
        newstate[i] = state[mapped + 1]
    end

    return QuantumState(newstate, state.bits)
end

# Performs the following mapping: |x, y> to |x, x + y mod2>
function cnot(state::QuantumState, ctrl::Int, flip::Int)
    newstate = Vector{Complex{Float64}}(length(state))

    ctrl = state.bits - ctrl
    flip = state.bits - flip

    for i = 1:length(state)
        basis = i - 1
        if iszero(basis, ctrl)
            newstate[i] = state[i]
        else
            mapped = flipbit(basis, flip)
            newstate[i] = state[mapped + 1]
        end
    end

    return QuantumState(newstate, state.bits)
end

# Performs the following mapping: |x, y, z> to |x, y, xy + z mod2>
function ccnot(state::QuantumState, ctrl1::Int, ctrl2::Int, flip::Int)
    newstate = Vector{Complex{Float64}}(length(state))

    ctrl1 = state.bits - ctrl1
    ctrl2 = state.bits - ctrl2
    flip  = state.bits - flip

    for i = 1:length(state)
        basis = i - 1
        if !iszero(basis, ctrl1) && !iszero(basis, ctrl2)
            mapped = flipbit(basis, flip)
            newstate[i] = state[mapped + 1]
        else
            newstate[i] = state[i]
        end
    end

    return QuantumState(newstate, state.bits)
end

# Performs the following mapping: |x, y> to |y, x>
function swap(state::QuantumState, x::Int, y::Int)
    newstate = Vector{Complex{Float64}}(length(state))

    x = state.bits - x
    y = state.bits - y

    for i = 1:length(state)
        basis = i - 1
        mapped = swapbit(basis, x, y)
        newstate[i] = state[mapped + 1]
    end

    return QuantumState(newstate, state.bits)
end

# Given the state |x, y, z> performs the following mapping:
# If x equals to 0
#     |x, y, z> to |x, y, z>
# If x equals to 1
#     |x, y, z> to |x, z, y>
function cswap(state::QuantumState, ctrl::Int, x::Int, y::Int)
    newstate = Vector{Complex{Float64}}(length(state))

    ctrl = state.bits - ctrl
    x = state.bits - x
    y = state.bits - y

    for i = 1:length(state)
        basis = i - 1
        if iszero(basis, ctrl)
            newstate[i] = state[i]
        else
            mapped = swapbit(basis, x, y)
            newstate[i] = state[mapped + 1]
        end
    end

    return QuantumState(newstate, state.bits)
end

# Performs the following mapping:
# |00> to |00>
# |01> to (0.5 + 0.5i) * |01> + (0.5 - 0.5i) * |10>
# |10> to (0.5 - 0.5i) * |01> + (0.5 + 0.5i) * |10>
# |11> to |11>
# This more generally translates to the following, less efficient
# implementation:
# |x, y> to (0.5 + 0.5i) * |x, y> + (0.5 - 0.5i) * |y, x>
function sqrtswap(state::QuantumState, x::Int, y::Int)
    newstate = Vector{Complex{Float64}}(length(state))

    x = state.bits - x
    y = state.bits - y

    mulorig = 0.5 + 0.5im
    mulswap = 0.5 - 0.5im

    for i = 1:length(state)
        basis = i - 1
        mapped = swapbit(basis, x, y)
        if basis == mapped
            newstate[i] = state[i]
        else
            newstate[i] = mulorig * state[i] + mulswap * state[mapped + 1]
        end
    end

    return QuantumState(newstate, state.bits)
end

# Given theta, it performs the following mapping:
# |0> to |0>
# |1> to e^(i * theta) * |1>
function phaseshift(state::QuantumState, bit::Int, theta::Float64)
    newstate = Vector{Complex{Float64}}(length(state))
    bit = state.bits - bit
    mul = exp(im * theta)

    for i = 1:length(state)
        basis = i - 1
        if iszero(basis, bit)
            newstate[i] = state[i]
        else
            newstate[i] = mul * state[i]
        end
    end

    return QuantumState(newstate, state.bits)
end

function cphase(state::QuantumState, control::Int, bit::Int, theta::Float64)
    newstate = Vector{Complex{Float64}}(length(state))
    control = state.bits - bit
    bit = state.bits - bit
    mul = exp(im * theta)

    for i = 1:length(state)
        basis = i - 1
        if iszero(basis, bit) || iszero(basis, control)
            newstate[i] = state[i]
        else
            newstate[i] = mul * state[i]
        end
    end

    return QuantumState(newstate, state.bits)
end

paulix(state::QuantumState, bit::Int) = not(state, bit)

# Performs the following mapping:
# |0> to i|1>
# |1> to -i|0>
function pauliy(state::QuantumState, bit::Int)
    newstate = Vector{Complex{Float64}}(length(state))
    bit = state.bits - bit

    for i = 1:length(state)
        basis = i - 1
        if iszero(basis, bit)
            mapped = setbit(basis, bit)
            newstate[i] =  im * state[mapped + 1]
        else
            mapped = clearbit(basis, bit)
            newstate[i] = -im * state[mapped + 1]
        end
    end

    return QuantumState(newstate, state.bits)
end

# A special case of the phaseshift gate. Performs the following mapping:
# |0> to |0>
# |1> to -|1>
function pauliz(state::QuantumState, bit::Int)
    newstate = Vector{Complex{Float64}}(length(state))
    bit = state.bits - bit

    for i = 1:length(state)
        basis = i - 1
        if iszero(basis, bit)
            newstate[i] =  state[i]
        else
            newstate[i] = -state[i]
        end
    end

    return QuantumState(newstate, state.bits)
end

function unitary(state::QuantumState, matrix::Array)
    if length(state.vector) != size(matrix, 1) || length(state.vector) != size(matrix, 2)
        error("The unitary matrix dimensions don't match the quantum state vector's");
    end

    return QuantumState(matrix * state.vector, state.bits)
end

function unitary(state::QuantumState, matrix::Array, x0::Int, xs::Int...)
    bits = [x0, xs...]

    if 2^length(bits) != size(matrix, 1) || 2^length(bits) != size(matrix, 2)
        error("The unitary matrix does not map enough bits")
    end

    newstate = state

    # First swap the bits to the first N positions
    for i = 1:length(bits)
        newstate = swap(newstate, i, bits[i])

        # Track swaps that would effect later qubits. An example:
        # Given a 3 bit state with an operation defined on bits 3 and 1 in that order,
        # let the bits be the x, y and x, their indices 1 2 and 3. We want to achieve
        # the following operation: x, y, z -> z, x, y
        # Step 0: Initial state: x, y, z [remaining bits: 3, 1]
        # Step 1: Swap the bit at position 1 with the first bit of the operation:
        #         x, y, z -> z, y, x [remaining bits: 1]
        # Step 2: Update the remaing bits. Replace every reference to position 1
        #         with the position it's being swapped to: [remaining bits: 3]
        # Step 3: Proceed to Step 1 with the updated remaining list
        for j = i + 1:length(bits)
            if bits[j] == i
                bits[j] = bits[i]
            end
        end
    end

    # Perform the unitary operation
    if length(state) == size(matrix, 1)
        newstate.vector = matrix * newstate.vector
    else
        gatematrix = kron(matrix, speye(2^(state.bits - length(bits))))
        newstate.vector = gatematrix * newstate.vector
    end

    # The bits variable stores the actual operations that were performed to achieve
    # the current bit order. Performing these in reverse order restores the correct
    # qubit order
    sb = length(bits) + 1
    reverse!(bits)
    for i = 1:length(bits)
        newstate = swap(newstate, sb - i, bits[i])
    end

    return newstate
end

function qft(state::QuantumState)
    omega = exp(2 * pi * im / state.bits)
    matrix = [omega^((i - 1) * (j - 1)) for i = 1:length(state), j = 1:length(state)] .* (1.0 / sqrt(length(state)))
    newstate = matrix * state.vector
    return QuantumState(newstate, state.bits)
end

function qft(state::QuantumState, bit::Int, bits::Int...)
    partialsize = 2^length([bit, bits...])
    omega = exp(2 * pi * im / state.bits)
    matrix = [omega^((i - 1) * (j - 1)) for i = 1:partialsize, j = 1:partialsize] .* (1.0 / sqrt(partialsize))
    return unitary(state, matrix, bit, bits...)
end

# Chooses one of the given gates, based on the decimal representation of the classical bit state given in
# the bits argument. It performs the following mapping:
# |x> to gates[decimal(bits) + 1](|x>)
# Where the +1 adjusts the result to Julia's 1-based indexing system
function choose1(state::QuantumState, bits::Vector{Int}, gates::Vector{Vector{Tuple{Function, Vector{Any}}}})
    index = todecimal(bits) + 1
    chain = gates[index]

    newstate::QuantumState = chain[1][1](state, chain[1][2]...)
    for i = 2:length(chain)
        newstate = chain[i][1](newstate, chain[i][2]...)
    end
    return newstate
end

# Selects a random basis vector based on the probabilities given in the argument.
# It uses the same algorithm as the one used in Fitness Proportionate Selection.
# TODO(gustorn): See if the stochastic acceptance variant performs better for large
#                quantum states
function randbasis{T<:Real}(probabilities::Vector{T})
    psum = sum(probabilities)

    target = rand() * psum
    for (i, v) in enumerate(probabilities)
        target -= v
        if (target <= 0)
            return i - 1
        end
    end

    # In case of floating point rounding error, return the element
    # with the otherwise highest probability
    return indmax(probabilities) - 1
end

# Measures the quantum state, collapsing it and returning both the resulting
# posterior state (with the appropriate basis vector having a probability of 1.0)
# and the result of the measurement.
function measure(state::QuantumState)
    probabilities = map(abs2, state.vector)::Vector{Float64}
    basis = randbasis(probabilities)

    # Measure the quantum state
    measurement = reverse(digits(basis, 2, state.bits))

    # Collapse the state
    newstate = fill(0.0 + 0im, length(state))
    newstate[basis + 1] = 1.0 + 0im
    posterior = QuantumState(newstate, state.bits)

    return posterior, measurement
end

# Performs partial measurement on the given bit of the quantum state
# in the computational basis
function partialmeasure(state::QuantumState, bit::Int)
    if state.bits == 1
        return measure(state)
    end

    posterior::Vector{Complex{Float64}} = fill(0.0 + 0im, 2^(state.bits - 1))
    bit = state.bits - bit

    # First collect the coefficients for the probabilities of the qubit being |0>
    prob0 = 0.0
    for i = 1:length(state)
        basis = i - 1
        if iszero(basis, bit)
            prob0 += abs2(state[i])
        end
    end

    # And measure the bit.
    measurement = rand() <= prob0 ? 0 : 1

    # Next collect the coefficients for the basis vectors of the
    # posterior states
    for i = 1:length(state)
        basis = i - 1
        if getbit(basis, bit) == measurement
            postbasis = removebit(basis, bit)
            posterior[postbasis + 1] += state[i]
        end
    end

    # Now normalize the new state
    # TODO(gustorn): see if the error from FP precision has any visible
    # effect on the posterior states
    n = sqrt(sumabs2(posterior))
    posterior = posterior .* (1.0 / n)

    return QuantumState(posterior, state.bits - 1), [measurement]
end

function partialmeasure(state::QuantumState, x0::Int, x1::Int, xs::Int...)
    # First construct the array of bits we want to measure.
    bits = [x0, x1, xs...]
    bits = unique(bits)
    nbits  = length(bits)

    if nbits == state.bits
        return measure(state)
    end

    map!(b -> state.bits - b, bits)
    sort!(bits)

    nposterior = state.bits - nbits
    probabilities::Vector{Float64} = fill(0.0, 2^nbits)
    posterior::Vector{Complex{Float64}} = fill(0.0 + 0im, 2^nposterior)

    # First build up the probabilities for the bits we're going to measure
    for i = 1:length(state)
        basis = i - 1
        measuredbasis = 0
        shiftby = 0

        for bit in bits
            measuredbasis = measuredbasis | (getbit(basis, bit) << shiftby)
            shiftby += 1
        end

        probabilities[measuredbasis + 1] += abs2(state[i])
    end

    # Measure it.
    measurement = digits(randbasis(probabilities), 2, nbits)

    # Now build up the posterior state
    for basis = 0:length(state) - 1
        # Check if the current basis vector matches the measurement at the specified bits
        basismatches = true
        for (bit, result) in zip(bits, measurement)
            if getbit(basis, bit) != result
                basismatches = false
                break
            end
        end

        # We only need to bother with the coefficients if the measurement matches the
        # partial basis vector
        if basismatches
            nremoved = 0
            postindex = basis
            for bit in bits
                postindex = removebit(postindex, bit - nremoved)
                nremoved += 1
            end

            # Add the coefficients from the old state to the new one
            posterior[postindex + 1] += state[basis + 1]
        end
    end

    # Normalize the new state
    n = sqrt(sumabs2(posterior))
    posterior = posterior .* (1 / n)

    # Match the bit order of the regular `measure` function
    reverse!(measurement)
    return QuantumState(posterior, state.bits - length(bits)), measurement
end

# Prints the current state of the qubit in terms of basis vectors, coefficients
# and probabilities. Useful for debugging and visualization. Guaranteed not to disturb
# the quantum state
function probe(state::QuantumState, name::AbstractString)
    println(name, "\n", state)
    return state
end

# The same as probe, but without the need for an explicit name
function probe(state::Vector{Int}, name::AbstractString)
    println(name, "\n", state, "\n")
    return state
end

end

max_stars_repo_name: leios/quantum

max_issues_repo_name: leios/quantum

max_forks_repo_name: leios/quantum


text:
#! /usr/bin/julia

# Rosetta Code, Loops/Nested

M = [rand(1:20) for i in 1:5, j in 1:10]
R, C = size(M)

println("The full matrix is:")
println(M, "\n")

println("Find the first 20:")
for i in 1:R, j in 1:C
    n = M[i,j]
    @printf "%4d" n
    if n == 20
        println()
        break
    elseif j == C
        println()
    end
end

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
using NewtonsMethod
using Test

@testset "NewtonsMethod.jl" begin
    # Write your own tests here.
    @test newtonroot(x -> x^2 , 1.0; tolerance = 1E-7, maxiter = 1) == nothing

    # several @test for the root of a known function
    @test newtonroot(x -> x^2 , x -> 2*x; tolerance = 1E-7, maxiter =1000) ≈≈ 5.960464477539063e-8
    @test newtonroot(x->x^3 , x->3*x^2, tolerance = 1E-7, maxiter =100) ≈≈ 1.3565659025724944e-7

    #tests of those roots using the automatic differentiation version of the function
    @test newtonroot(x -> x^2 , 1.0; tolerance = 1E-7, maxiter =1000) ≈≈ 5.960464477539063e-8
    @test newtonroot(x -> x^3 , 2.0, tolerance = 1E-7, maxiter =100) ≈≈ 1.808754536763326e-7

    #test of non-convergence for a function without a root (e.g. f(x)=2+x2 )
    @test newtonroot(x-> 2+x^2 , 1.0; tolerance = 1E-7, maxiter =1000) ≈≈ 5.960464477539063e-8

    #test to ensure that the maxiter is working (e.g. what happens if you call maxiter = 5
    @test newtonroot(x -> x^2 , 1.0; tolerance = 1E-7, maxiter = 5) == nothing

    #test to ensure that tol is working
    @test newtonroot(x -> x^2 , 1.0; tolerance = 1E-2, maxiter = 5) == nothing
end


text:
using LinearAlgebra
using Flux
using CUDA
import Flux.gpu, Flux.cpu

includet("./solvers.jl")
includet("./fixed_point_layer_general.jl")
includet("./output_layer.jl")

# Implicit, direct parametrizations
abstract type implicit_param end

# Explicit types that we simulate
mutable struct explicit_ff_cell
    ϕ
    A
    B1
    B2
    C1
    D11
    D12
    bx
    bv
end

Flux.trainable(L::explicit_ff_cell) = [L.A, L.B1, L.B2, L.C1, L.D11, L.D12, L.bx, L.bv]

## Stable REN cell
mutable struct implicit_ff_cell{T} <: implicit_param
    ϕ
    V::Union{Matrix{T},CuMatrix{T}}
    S_1::Union{Matrix{T},CuMatrix{T}}
    B2::Union{Matrix{T},CuMatrix{T}}
    D12::Union{Matrix{T},CuMatrix{T}}
    bx::Union{Vector{T},CuVector{T}}
    bv::Union{Vector{T},CuVector{T}}
    ϵ::T
end
Flux.trainable(L::implicit_ff_cell) = [L.V, L.S_1, L.B2, L.D12, L.bx, L.bv]

function Flux.gpu(M::implicit_ff_cell{T}) where T
    if T != Float32
        println("Moving type: ", T, " to gpu may not be supported. Try Float32!")
    end
    return implicit_ff_cell{T}(M.ϕ, gpu(M.V), gpu(M.S_1), gpu(M.B2),
                             gpu(M.D12), gpu(M.bx), gpu(M.bv), M.ϵ)
end

function Flux.cpu(M::implicit_ff_cell{T}) where T
    return implicit_ff_cell{T}(M.ϕ, cpu(M.V), cpu(M.S_1), cpu(M.B2),
                             cpu(M.D12), cpu(M.bx), cpu(M.bv), M.ϵ)
end

function implicit_ff_cell{T}(nu, nx, nv; nl=relu, ϵ=0.01, bx_scale=0.0, bv_scale=1.0) where T
    glorot_normal(n, m) = convert.(T, randn(n, m) / sqrt(n + m))
    E = Matrix{T}(I, nx, nx)
    F = Matrix{T}(I, nx, nx)
    P = Matrix{T}(I, nx, nx)
    B1 = zeros(T, nx, nv)
    B2 = glorot_normal(nx, nu)

    C1 = zeros(T, nv, nx)
    D11 = glorot_normal(nv, nv)
    D12 = zeros(T, nv, nu)

    Λ = I  #Choose sum abs values to make diagonally dominant
    H22 = 2Λ - D11 - D11'
    Htild = [(E + E' - P) -C1' F';
             -C1 H22 B1'
             F B1  P] + ϵ * I

    S_1 = glorot_normal(nx, nx)

    V = Matrix{T}(cholesky(Htild).U) # H = V'*V

    bv = convert.(T, bv_scale * randn(nv) / sqrt(nv))
    bx = convert.(T, bx_scale * rand(nx) / sqrt(nx))

    return implicit_ff_cell{T}(nl, V, S_1, B2, D12, bx, bv, ϵ)
end
# constructor if type is not specified
implicit_ff_cell(nu, nx, nv; nl=relu, ϵ=0.01, bx_scale=0.0, bv_scale=1.0) = implicit_ff_cell{Float64}(nu, nx, nv; nl=nl, ϵ=ϵ, bx_scale=bx_scale, bv_scale=bv_scale)

function sample_ff_ren(nu, nx, nv; nl=relu, ϵ=0.01, bx_scale=0.0, bv_scale=1.0)
    # glorot_normal(n, m) = randn(n, m) / sqrt(n + m)
    B2 = randn(nx, nu) / sqrt(nx+nu)

    D12 = randn(nv, nu) / sqrt(nv+nu)
    V = 2 * randn(2nx + nv, 2nx + nv) / sqrt(2nx + nv)

    bv = bv_scale * randn(nv) / sqrt(nv)
    bx = bx_scale * rand(nx) / sqrt(nx)

    S_1 = randn(nx, nx) / sqrt(nx+nx)

    return implicit_ff_cell(nl, V, S_1, B2, D12, bx, bv, ϵ)
end

function explicit(model::implicit_ff_cell)
    nx = size(model.B2, 1)
    nu = size(model.B2, 2)
    nv = size(model.D12, 1)

    H = model.V' * model.V + model.ϵ * I

    # For some reason taking view doesn't work with CUDA
    H11 = H[1:nx, 1:nx]
    H22 = H[ nx + 1:nx + nv, nx + 1:nx + nv]
    H33 = H[ nx + nv + 1:2nx + nv, nx + nv + 1:2nx + nv]
    H21 = H[nx + 1:nx + nv, 1:nx]
    H31 = H[ nx + nv + 1:2nx + nv, 1:nx]
    H32 = H[ nx + nv + 1:2nx + nv, nx + 1:nx + nv]

    # Implicit model parameters
    S_1 = (model.S_1 - model.S_1') / 2

    B2 = model.B2
    D12 = model.D12

    P = H33
    E = (H11 + P + S_1) / 2
    F = H31

    # equilibrium network stuff
    C1 = - (H21)
    B1 = H32
    Λᵢ = diag(H22) / 2
    D11 = -tril(H22, -1)  # strictly lower triangular part

    # Construct explicit rnn model.
    # Use \bb to differentiate from implicit model
    𝔸 = E \ F
    𝔹_1 = E \ B1
    𝔹_2 = E \ B2

    ℂ_1 = (1 ./ Λᵢ) .* C1
    𝔻_11 = (1 ./ Λᵢ) .* D11
    𝔻_12 = (1 ./ Λᵢ) .* D12

    bx = model.bx
    bv = model.bv

    return explicit_ff_cell(model.ϕ, 𝔸, 𝔹_1, 𝔹_2, ℂ_1, 𝔻_11, 𝔻_12, bx, bv)
end

function init_state(model::implicit_ff_cell{T}, batches) where T
    nx = size(model.B2, 1)
    return typeof(model.V)(zeros(T, nx, batches))
end

function (exp_cell::explicit_ff_cell)(xt, ut)
    b = exp_cell.C1 * xt + exp_cell.D12 * ut .+ exp_cell.bv
    wt = tril_eq_layer(exp_cell.ϕ, exp_cell.D11, b)
    xn = exp_cell.A * xt + exp_cell.B1 * wt + exp_cell.B2 * ut .+ exp_cell.bx
    return xn, (xt, wt)
end

function simulate(exp_cell::explicit_ff_cell, x0, ut)
    recurrent = Flux.Recur(exp_cell, x0)
    return unzip(recurrent.(ut))
end

function (implicit_cell::implicit_ff_cell)(xt, ut)
    exp_cell = explicit(implicit_cell)
    return exp_cell(xt, ut)
end

function simulate(cell::implicit_ff_cell, x0, ut)
    exp_cell = explicit(cell)
    return simulate(exp_cell, x0, ut)
end

# Dissipative feedforward RNN code

mutable struct dissipative_ff_rnn{T}
    nu
    nx
    nv
    ny
    implicit_cell::implicit_ff_cell{T}
    output::output{T}
    Q::Matrix{T}
    S::Matrix{T}
    R::Matrix{T}
end

Flux.trainable(model::dissipative_ff_rnn) = (Flux.trainable(model.implicit_cell)..., Flux.trainable(model.output)...)

# Constructors
function dissipative_ff_rnn{T}(nu, nx, nv, ny, Q, S, R; nl=relu, ϵ=0.001, bx_scale=0.0, bv_scale=1.0) where T
    cell_params = implicit_ff_cell{T}(nu, nx, nv; nl=nl, ϵ, bx_scale, bv_scale)
    output_params = output{T}(nu, nx, nv, ny)
    return dissipative_ff_rnn{T}(nu, nx, nv, ny, cell_params, output_params, Q, S, R)
end

function init_state(model::dissipative_ff_rnn{T}, batches) where T
    return typeof(model.implicit_cell.V)(zeros(T, model.nx, batches))
end

# Common forms of dissipativity
function bounded_ff_rnn(T::DataType, nu, nx, nv, ny, γ;
                        nl=relu, ϵ=0.001, bx_scale=0.0, bv_scale=1.0)
    R = Matrix{T}(γ * I, nu, nu)
    S = zeros(T, nu, ny)
    Q = Matrix{T}(-I / γ, ny, ny)
    return dissipative_ff_rnn{T}(nu, nx, nv, ny, Q, S, R; nl=relu, ϵ=0.001, bx_scale=0.0, bv_scale=1.0)
end

function stable_ff_rnn(T::DataType, nu, nx, nv, ny;
                        nl=relu, ϵ=0.001, bx_scale=0.0, bv_scale=1.0)
    return bounded_ff_rnn(T, nu, nx, nv, ny, Inf; nl=nl, ϵ=ϵ, bx_scale=bx_scale, bv_scale=bv_scale)
end

function passive_ff_rnn(T::DataType, nu, nx, nv, ny;
                     η=0.0, nl=relu, ϵ=0.001, bx_scale=0.0, bv_scale=1.0)

    R = zeros(T, nu, nu)
    S = Matrix(I / 2, nu, ny)
    Q = Matrix(-η * I, nu, nu)
    return dissipative_ff_rnn(nu, nx, nv, ny, Q, S, R; nl=relu, ϵ=0.001, bx_scale=0.0, bv_scale=1.0)
end

# default type is Float64
bounded_ff_rnn(nu, nx, nv, ny, γ; nl=relu, ϵ=0.001, bx_scale=0.0, bv_scale=1.0) =
                bounded_ff_rnn(Float64, nu, nx, nv, ny, γ;
                            nl=nl, ϵ=ϵ, bx_scale=bx_scale, bv_scale=bv_scale)

stable_ff_rnn(nu, nx, nv, ny; nl=relu, ϵ=0.001, bx_scale=0.0, bv_scale=1.0) =
                stable_ff_rnn(Float64, nu, nx, nv, ny;
                                nl=nl, ϵ=ϵ, bx_scale=bx_scale, bv_scale=bv_scale)

passive_ff_rnn(nu, nx, nv, ny; η=0.0, nl=relu, ϵ=0.001, bx_scale=0.0, bv_scale=1.0) =
                passive_ff_rnn(Float64, nu, nx, nv, ny;
                                        η=η, nl=nl, ϵ=ϵ, bx_scale=bx_scale, bv_scale=bv_scale)

# Construct and explicit
function explicit(model::dissipative_ff_rnn)
    nx = model.nx
    nu = model.nu
    ny = model.ny
    nv = model.nv

    # dissipation parameter
    Q = model.Q
    S = model.S
    R = model.R

    # Implicit model parameters
    S_1 = (model.implicit_cell.S_1 - model.implicit_cell.S_1') / 2

    C2 = model.output.C2
    D21 = model.output.D21
    D22 = model.output.D22

    B2 = model.implicit_cell.B2
    D12 = model.implicit_cell.D12

    # RHS of dissipation inequality
    Γ1 = [C2'; D21'; zeros(nx, ny)] * Q * [C2 D21 zeros(ny, nx)]  # possibly transpose zeros
    Γ2 = [(C2' * S'); (D21' * S' - D12); B2] * inv(R) * [(S * C2) (S * D21 - D12') B2']

    H = model.implicit_cell.V' * model.implicit_cell.V - Γ1 + Γ2 + model.implicit_cell.ϵ * I
    H11 = H[1:nx, 1:nx]
    H22 = H[ nx + 1:nx + nv, nx + 1:nx + nv]
    H33 = H[ nx + nv + 1:2nx + nv, nx + nv + 1:2nx + nv]
    H21 = H[nx + 1:nx + nv, 1:nx]
    H31 = H[ nx + nv + 1:2nx + nv, 1:nx]
    H32 = H[ nx + nv + 1:2nx + nv, nx + 1:nx + nv]

    # extract parameters from implicit rnn
    Λᵢ = diag(H22) / 2
    D11 = -tril(H22, -1)
    P = H33
    E = (H11 + P + S_1) / 2
    F = H31
    C1 = -H21
    B1 = H32

    # Construct explicit rnn model. Use \bb font to differentiate
    A = E \ F
    𝔹_1 = E \ B1
    𝔹_2 = E \ B2

    ℂ_1 = (1 ./ Λᵢ) .* C1
    𝔻_11 = (1 ./ Λᵢ) .* D11
    𝔻_12 = (1 ./ Λᵢ) .* D12

    bx = model.implicit_cell.bx
    bv = model.implicit_cell.bv

    return explicit_ff_cell(model.implicit_cell.ϕ, A, 𝔹_1, 𝔹_2, ℂ_1, 𝔻_11, 𝔻_12, bx, bv)
end

function (model::dissipative_ff_rnn)(x0, ut)
    exp_cell = explicit(model)
    xt, wt = simulate(exp_cell, x0, ut)
    return model.output.(xt, wt, ut), xt[end]
end

function check_lmi(model::dissipative_ff_rnn)
    nx = model.nx
    nu = model.nu
    ny = model.ny
    nv = model.nv

    # dissipation parameter
    Q = model.Q
    S = model.S
    R = model.R

    # Implicit model parameters
    S_1 = (model.implicit_cell.S_1 - model.implicit_cell.S_1') / 2

    C2 = model.output.C2
    D21 = model.output.D21
    D22 = model.output.D22

    B2 = model.implicit_cell.B2
    D12 = model.implicit_cell.D12

    # RHS of dissipation inequality
    Γ1 = [C2'; D21'; zeros(nx, ny)] * Q * [C2 D21 zeros(ny, nx)]  # possibly transpose zeros
    Γ2 = [(C2' * S'); (D21' * S' - D12); B2] * inv(R) * [(S * C2) (S * D21 - D12') B2']

    H = model.implicit_cell.V' * model.implicit_cell.V - Γ1 + Γ2 + model.implicit_cell.ϵ * I
    H11 = H[1:nx, 1:nx]
    H22 = H[ nx + 1:nx + nv, nx + 1:nx + nv]
    H33 = H[ nx + nv + 1:2nx + nv, nx + nv + 1:2nx + nv]
    H21 = H[nx + 1:nx + nv, 1:nx]
    H31 = H[ nx + nv + 1:2nx + nv, 1:nx]
    H32 = H[ nx + nv + 1:2nx + nv, nx + 1:nx + nv]

    # extract parameters from implicit rnn
    Λᵢ = diag(H22) / 2
    Λ = Diagonal(Λᵢ)
    D11 = -tril(H22, -1)
    P = H33
    E = (H11 + P + S_1) / 2
    F = H31
    C1 = -H21
    B1 = H32

    # Construct the LMI
    lmi = [(E + E' - P + C2' * Q * C2) (C2' * Q * D21 - C1') (C2' * S') F';
           (D21' * Q * C2 - C1) (2Λ - D11 - D11' + D21' * Q * D21) (D21' * S' - D12) B1';
           (S * C2) (S * D21 - D12') R B2';
           F B1 B2 P]
end

function test_ff_REN()

    device = cpu
    # Check equilibrium solve accuracy
    nu = 5
    model = implicit_ff_cell{Float64}(5, 5, 10; nl=tanh) |> device
    nPoints = 100

    x0 = init_state(model, nPoints) |> device
    utrain = randn(nu, nPoints) |> device
    ytrain = sin.(utrain)

    exp_cell = explicit(model)

    xn, (xt, wt) = exp_cell(x0, utrain)

    b = exp_cell.C1 * x0 + exp_cell.D12 * utrain .+ exp_cell.bv
    err = norm(wt - exp_cell.ϕ.(exp_cell.D11 * wt + b))

    if err > 1E-8
        println("Forward solve doesnt seem to be working well")
    end

    # check gradient calculation via FD
    sample_loss(y1, y2) = norm(y1 - y2)^2
    L() = mean(sample_loss.(ytrain, model(x0, utrain)[1]))
    fd_test_grads(L, Flux.params(model))  # param gradients
    fd_test_grads(L, Flux.Params([x0, utrain]))  # input gradients

    # Try simple noiseless regression task
    model = implicit_ff_cell{Float32}(1, 1, 40; nl=tanh, ϵ=1.0) |> device

    nPoints = 50
    utrain = 2 * randn(1, nPoints) |> device
    ytrain = (utrain.^3)
    x0 = init_state(model, nPoints)

    loss() = norm(ytrain - model(x0, utrain)[1])
    ps = Flux.Params(Flux.trainable(model))

    opt = Flux.Optimise.ADAM(1E-3)
    for k in 1:2000
        train_loss, back = Zygote.pullback(loss, ps)

        # calculate gradients and update loss
        ∇J = back(one(train_loss))
        update!(opt, ps, ∇J)

        printfmt("Iteration: {1:2d}\tTraining loss: {2:1.2E}\n", k, train_loss)
    end

    exp_cell = explicit(model)
    ϕ = exp_cell.ϕ
    W = exp_cell.D11
    b = exp_cell.C1 * x0 + exp_cell.D12 * utrain .+ exp_cell.bv

    xn, (xt, wt) = model(x0, utrain)
    err = norm(wt - ϕ.(W * wt + b))
    zeq = wt

    # Run to test if gradient is wrong
    fd_test_grads(loss, Flux.params(model))  # param gradients
    fd_test_grads(loss, Flux.Params([x0, utrain]))  # input gradients

    yest, (xt, wt) = model(x0, utrain)

    plot(utrain', ytrain'; seriestype=:scatter)
    plot!(utrain', yest'; seriestype=:scatter)

end

max_stars_repo_name: imanchester/REN

max_issues_repo_name: imanchester/REN

max_forks_repo_name: imanchester/REN


text:
# Julia 0.6 for now

# Based heavily on CUFFT.jl/gen/wrap_cufft.jl
# https://github.com/JuliaArchive/CUFFT.jl/blob/master/gen/wrap_cufft.jl

using Clang, Clang.cindex
# We'll replace some of the wrap functions, and we need to call the following:
import Clang.wrap_c: repr_jl, rep_type, rep_args, name_safe

include("../src/find_ros2.jl")

clang_includes=[
    ROS_INCLUDE_PATH,
    "/usr/include",
    "/usr/lib/llvm-6.0/lib/clang/6.0.0/include",    # stdbool.h
]

function header_wrapped(header_file, cursor_name)
    @show header_file, cursor_name
    return startswith(header_file, ROS_INCLUDE_PATH) &&
        (startswith(cursor_name, "rcl_") || startswith(cursor_name, "rmw_") || startswith(cursor_name, "rcutils_"))
end

# TODO: figure out the Void -> Cvoid conversion to eliminate deprecation warnings

context = wrap_c.init(
    output_file="gen_librcl.jl",
    common_file="gen_librcl_h.jl",
    clang_includes=clang_includes,
    header_wrapped = (header_file, cursor_name) -> startswith(cursor_name, ROS_INCLUDE_PATH),
    header_library = header_file -> "lib$(dirname(relpath(header_file, ROS_INCLUDE_PATH)))",
)

context.headers = [
    joinpath(ROS_INCLUDE_PATH, "rcutils", "allocator.h"),
    joinpath(ROS_INCLUDE_PATH, "rcutils", "error_handling.h"),
    joinpath(ROS_INCLUDE_PATH, "rcl", "node.h"),
    joinpath(ROS_INCLUDE_PATH, "rcl", "rcl.h"),
    ]
run(context)


text:
module ScriptUtilities

using Test, Pkg, Suppressor, Dates

export print_test_results, capture_test

"""
    print_test_results( test_func;
        params=[] :: Vector{Any},
        stdout=true :: Bool,
        dir="./" :: String
    ) :: Bool

Prints the results the test results and meta data to `.txt.` file or `STDOUT`.
If results are printed to `.txt` file, then the file name is `"test_func_datetime.txt"`
where `"datetime"` is a stringified datetime.
This method returns a Boolean value which is `true` if the test passes and `false`
otherwise.

Arguments:
* `test_func` - a generic method that runs a `@testset`.
* `params` - optional keyword argument, which contains the parameters to run `test_func`.
* `stdout` - optional keyword argument, if `true` print results to `STDOUT` instead of `.txt` file.
* `dir` - optional keyword argument, directory to which `.txt` file is saved.

Results are printed in the form:

```
Version :    # SignalingDimension.jl version used

Status `~/.julia/dev/SignalingDimension.jl/script/Project.toml`
  [96aab1c2] SignalingDimension v0.1.1 [`..`]

Test Completed : 2021-02-26T10:17:26        # Datetime of test completion
Elapsed Time : 16.799114946                     # Time elapsed during test
Test Method : verify_maximum_likelihood_facet    # Test method used
Test Args : (30)                                # Arguments used for test method
Test Pass : true                                # All tests pass

Test Summary:                                                     | Pass  Total
Testing all maximum likelihood facets of size `(30, 30)` or less. |  406    406
```

!!! warning "Do not run within @testset"
    Running within a @testset prevents the test results from being properly captured.
    If this method must be run in a testset, then `test_func` not contain a testset
    and instead print results to STDOUT and throw an exception on failure.
"""
function print_test_results(test_func; params=[] :: Vector{Any}, stdout=true :: Bool, dir="./" :: String) :: Bool
    name = string(test_func)
    datetime = replace(string(now()), r"\.\d+$" => "")
    datetime_file = replace(datetime, ":" => "-")

    filename = join([dir, "/", name, "_", datetime_file, ".txt"])

    elapsed_time = @elapsed (pass, results) = capture_test(test_func, params=params)

    version = @capture_out Pkg.status("SignalingDimension")

    results_string = join([
        "Version :\n\n", version,
        "\nTest Completed : ", datetime,
        "\nElapsed Time : ", elapsed_time,
        "\nTest Method : ", name,
        "\nTest Args : ", "(" * join(params, ", ") * ")",
        "\nTest Pass : ", string(pass),
        "\n\n", results
    ])

    if !stdout
        open(filename, "w") do io
            println(io, results_string)
        end
    else
        println(results_string)
    end

    return pass
end

"""
    capture_test(test_func; params=[] :: Vector{Any}) :: Tuple{Bool, String}

Captures the `STDOUT` and pass staus of the provided `test_func`.
A tuple is returned as `(pass, results)` where `pass` is a boolean indicating if
the test passedand `results` is a string containing the `STDOUT` from `test_func`.

Arguments:
* `test_func` - a generic method that runs a `@testset`.
* `params` - optional keyword argument, which contains the parameters to run `test_func`.

!!! warning "Do not run within @testset"
    Running within a @testset prevents the test results from being properly captured.
    If this method must be run in a testset, then `test_func` not contain a testset
    and instead print results to STDOUT and throw an exception on failure.
"""
function capture_test(test_func; params=[] :: Vector{Any}) :: Tuple{Bool, String}
    pass = nothing
    results = @capture_out try
        test_func(params...); pass = true
    catch err
        println(err); pass = false
    end

    (pass, results)
end

end  # module ScriptUtilities


text:
type Music
    ptr::Ptr{Void}

    function Music(ptr::Ptr{Void})
        m = new(ptr)
        finalizer(m, destroy)
        m
    end
end

function Music(filename::AbstractString)
    Music(ccall((:sfMusic_createFromFile, libcsfml_audio), Ptr{Void}, (Ptr{Cchar},), filename))
end

function destroy(music::Music)
    ccall((:sfMusic_destroy, libcsfml_audio), Void, (Ptr{Void},), music.ptr)
end

function set_loop(music::Music, loop::Bool)
    ccall((:sfMusic_setLoop, libcsfml_audio), Void, (Ptr{Void}, Int32,), music.ptr, loop)
end

function get_loop(music::Music)
    return Bool(ccall((:sfMusic_getLoop, libcsfml_audio), Int32, (Ptr{Void},), music.ptr))
end

function play(music::Music)
    ccall((:sfMusic_play, libcsfml_audio), Void, (Ptr{Void},), music.ptr)
end

function pause(music::Music)
    ccall((:sfMusic_pause, libcsfml_audio), Void, (Ptr{Void},), music.ptr)
end

function stop(music::Music)
    ccall((:sfMusic_pause, libcsfml_audio), Void, (Ptr{Void},), music.ptr)
end

function get_duration(music::Music)
    return ccall((:sfMusic_getDuration, libcsfml_audio), Time, (Ptr{Void},), music.ptr)
end

function get_channelcount(music::Music)
    return Real(ccall((:sfMusic_getChannelCount, libcsfml_audio), UInt32, (Ptr{Void},), music.ptr))
end

function get_samplerate(music::Music)
    return Real(ccall((:sfMusic_getSampleRate, libcsfml_audio), UInt32, (Ptr{Void},), music.ptr))
end

function set_pitch(music::Music, pitch::Real)
    ccall((:sfMusic_setPitch, libcsfml_audio), Void, (Ptr{Void}, Cfloat,), music.ptr, pitch)
end

function set_volume(music::Music, volume::Real)
    ccall((:sfMusic_setVolume, libcsfml_audio), Void, (Ptr{Void}, Cfloat,), music.ptr, volume)
end

function get_pitch(music::Music)
    return Real(ccall((:sfMusic_getPitch, libcsfml_audio), Cfloat, (Ptr{Void},), music.ptr))
end

function get_volume(music::Music)
    return Real(ccall((:sfMusic_getVolume, libcsfml_audio), Cfloat, (Ptr{Void},), music.ptr))
end

function get_status(music::Music)
    return ccall((:sfMusic_getStatus, libcsfml_audio), Int32, (Ptr{Void},), music.ptr)
end


text:
@testset "Numerical integral" begin
    gl1 = BasicBSpline.GaussLegendre(1)
    gl2 = BasicBSpline.GaussLegendre(2)
    gl3 = BasicBSpline.GaussLegendre(3)
    gl4 = BasicBSpline.GaussLegendre(4)

    I = 2..8
    coef = rand(10)

    f(n,x) = sum(coef[i]*x^(i-1) for i in 1:n)
    F(n,x) = sum(coef[i]*x^(i)/i for i in 1:n)

    @test F(1,maximum(I)) - F(1,minimum(I)) ≈ BasicBSpline.integrate(x->f(1,x), I, gl1)
    @test F(2,maximum(I)) - F(2,minimum(I)) ≈ BasicBSpline.integrate(x->f(2,x), I, gl1)
    @test F(3,maximum(I)) - F(3,minimum(I)) ≈ BasicBSpline.integrate(x->f(3,x), I, gl2)
    @test F(4,maximum(I)) - F(4,minimum(I)) ≈ BasicBSpline.integrate(x->f(4,x), I, gl2)
    @test F(5,maximum(I)) - F(5,minimum(I)) ≈ BasicBSpline.integrate(x->f(5,x), I, gl3)
    @test F(6,maximum(I)) - F(6,minimum(I)) ≈ BasicBSpline.integrate(x->f(6,x), I, gl3)
    @test F(7,maximum(I)) - F(7,minimum(I)) ≈ BasicBSpline.integrate(x->f(7,x), I, gl4)
    @test F(8,maximum(I)) - F(8,minimum(I)) ≈ BasicBSpline.integrate(x->f(8,x), I, gl4)
end


text:
using Documenter
using HyrodiumHome

# English docs
makedocs(;
    modules=[HyrodiumHome],
    authors="hyrodium <hyrodium@gmail.com>",
    repo="https://github.com/hyrodium/hyrodium.github.io/blob/{commit}{path}#L{line}",
    sitename="Hyrodium Home",
    doctest=false,
    format=Documenter.HTML(;
        prettyurls=true,
        canonical="https://hyrodium.github.io",
        assets=String[],
        lang="en",
    ),
    source="src_en",
    build="build_en",
    pages=[
        "Home" => "index.md",
        "PDF documents" => "pdf.md",
        "History" => "history.md",
        "SNS and etc." => "sns.md",
        "GitHub pages" => "githubpages.md",
    ],
)

# Japanese docs
makedocs(;
    modules=[HyrodiumHome],
    authors="hyrodium <hyrodium@gmail.com>",
    repo="https://github.com/hyrodium/hyrodium.github.io/blob/{commit}{path}#L{line}",
    sitename="Hyrodium Home",
    doctest=false,
    format=Documenter.HTML(;
        prettyurls=true,
        canonical="https://hyrodium.github.io",
        assets=String[],
        lang="ja",
    ),
    source="src_ja",
    build="build_ja",
    pages=[
        "Home" => "index.md",
        "PDF資料" => "pdf.md",
        "出展・発表歴" => "history.md",
        "SNS・連絡先" => "sns.md",
        "GitHub pages" => "githubpages.md",
    ],
)

for (root, dirs, files) in (walkdir("docs/build_en") ∪ walkdir("docs/build_ja"))
    for file in files
        if file == "index.html"
            path_html = joinpath(root, file)
            script = read(path_html, String)
            script = replace(script, ">Version<" => ">Language / 言語<")
            write(path_html, script)
        end
    end
end

# 404 docs
makedocs(;
    modules=[HyrodiumHome],
    authors="hyrodium <hyrodium@gmail.com>",
    repo="https://github.com/hyrodium/hyrodium.github.io/blob/{commit}{path}#L{line}",
    sitename="Hyrodium Home",
    doctest=false,
    format=Documenter.HTML(;
        prettyurls=true,
        canonical="https://hyrodium.github.io",
        assets=String[],
        lang="en",
    ),
    source="src_404",
    build="build_404",
    pages=[
        "404" => "index.md",
    ],
)
mv("docs/build_404/index.html", "docs/build_404/404.html")

script = read("docs/build_404/404.html", String)
script = replace(script, "assets/" => "https://hyrodium.github.io/assets/")
script = replace(script, "../versions.js" => "https://hyrodium.github.io/versions.js")
script = replace(script, ">Version<" => ">Language / 言語<")
write("docs/build_404/404.html", script)
write("docs/build_404/index.html", """<!DOCTYPE html><meta http-equiv="refresh" content="0; url=https://hyrodium.github.io/en" />""")

deploydocs(;
    target="build_404",
    repo="github.com/hyrodium/hyrodium.github.io",
    versions=nothing
)

deploydocs(;
    target="build_en",
    devurl="en",
    repo="github.com/hyrodium/hyrodium.github.io"
)

deploydocs(;
    target="build_ja",
    devurl="ja",
    repo="github.com/hyrodium/hyrodium.github.io"
)


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

name = "Libtask"
version = v"1.0.0"

# Collection of sources required to build Libtask
sources = [
    "https://github.com/TuringLang/Libtask.git" =>
    "fbe2feb0dab2f9d3d59f7710ff95c0cb2e61a835",

]

# Bash recipe for building across all platforms
script = raw"""
if [ $target = "arm-linux-gnueabihf" ]; then
cd $WORKSPACE/srcdir
wget "https://julialang-s3.julialang.org/bin/linux/armv7l/1.0/julia-1.0.0-linux-armv7l.tar.gz"
tar xzvf julia-1.0.0-linux-armv7l.tar.gz
rm *.tar.gz
mv julia* julia
LIBS="`pwd`/julia/lib"
LIBSJL="`pwd`/julia/lib/julia"
INCLUDES="`pwd`/julia/include/julia"
cd Libtask
gcc -O2 -shared -std=gnu99 -I$INCLUDES -DJULIA_ENABLE_THREADING=1 -fPIC -L$LIBSJL -L$LIBS -Wl,--export-dynamic -ljulia task.c -o libtask.so
mkdir $prefix/lib
mv libtask.so $prefix/lib
exit
fi

if [ $target = "i686-linux-gnu" ]; then
cd $WORKSPACE/srcdir
wget "https://julialang-s3.julialang.org/bin/linux/x86/1.0/julia-1.0.0-linux-i686.tar.gz"
tar xzvf julia-1.0.0-linux-i686.tar.gz
rm *.tar.gz
mv julia* julia
LIBS="`pwd`/julia/lib"
LIBSJL="`pwd`/julia/lib/julia"
INCLUDES="`pwd`/julia/include/julia"
cd Libtask
gcc -march=pentium4 -O2 -shared -std=gnu99 -I$INCLUDES -DJULIA_ENABLE_THREADING=1 -fPIC -L$LIBSJL -L$LIBS -Wl,--export-dynamic -ljulia task.c -o libtask.so
mkdir $prefix/lib
mv libtask.so $prefix/lib
exit
fi

if [ $target = "x86_64-linux-gnu" ]; then
cd $WORKSPACE/srcdir
wget "https://julialang-s3.julialang.org/bin/linux/x64/1.0/julia-1.0.0-linux-x86_64.tar.gz"
tar xzvf julia-1.0.0-linux-x86_64.tar.gz
rm *.tar.gz
mv julia* julia
LIBS="`pwd`/julia/lib"
LIBSJL="`pwd`/julia/lib/julia"
INCLUDES="`pwd`/julia/include/julia"
cd Libtask
gcc -O2 -shared -std=gnu99 -I$INCLUDES -DJULIA_ENABLE_THREADING=1 -fPIC -L$LIBSJL -L$LIBS -Wl,--export-dynamic -ljulia task.c -o libtask.so
mkdir $prefix/lib
mv libtask.so $prefix/lib
exit
fi

if [ $target = "x86_64-w64-mingw32" ]; then
cd $WORKSPACE/srcdir
wget "http://mlg.eng.cam.ac.uk/hong/julia-1.0.0-win64.tar.gz"
tar xzvf julia-1.0.0-win64.tar.gz
rm *.tar.gz
mv julia* julia
LIBS="`pwd`/julia/bin"
INCLUDES="`pwd`/julia/include/julia"
cd Libtask
gcc -O2 -shared -std=gnu99 -I$INCLUDES -DJULIA_ENABLE_THREADING=1 -fPIC -L$LIBS -Wl,--export-all-symbols -ljulia -lopenlibm task.c -o libtask.dll
mkdir $prefix/bin
mv libtask.dll $prefix/bin
exit

fi

if [ $target = "i686-w64-mingw32" ]; then
cd $WORKSPACE/srcdir
wget "http://mlg.eng.cam.ac.uk/hong/julia-1.0.0-win32.tar.gz"
tar xzvf julia-1.0.0-win32.tar.gz
rm *.tar.gz
mv julia* julia
LIBS="`pwd`/julia/bin"
INCLUDES="`pwd`/julia/include/julia"
cd Libtask
gcc -march=pentium4 -O2 -shared -std=gnu99 -I$INCLUDES -DJULIA_ENABLE_THREADING=1 -fPIC -L$LIBS -Wl,--export-all-symbols -ljulia -lopenlibm task.c -o libtask.dll
mkdir $prefix/bin
mv libtask.dll $prefix/bin
exit

fi

if [ $target = "x86_64-apple-darwin14" ]; then
cd $WORKSPACE/srcdir
wget "http://mlg.eng.cam.ac.uk/hong/julia-1.0.0-mac64.tar.gz"
tar xzvf julia-1.0.0-mac64.tar.gz
rm *.tar.gz
mv julia* julia
LIBS="`pwd`/julia/Contents/Resources/julia/lib"
LIBSJL="`pwd`/julia/Contents/Resources/julia/lib/julia"
INCLUDES="`pwd`/julia/Contents/Resources/julia/include/julia"
cd Libtask
gcc -O2 -shared -std=gnu99 -I$INCLUDES -DJULIA_ENABLE_THREADING=1 -fPIC -L$LIBSJL -L$LIBS -ljulia task.c -o libtask.dylib
mkdir $prefix/lib
mv libtask.dylib $prefix/lib
exit

fi

"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:i686, :glibc),
    Linux(:x86_64, :glibc),
    Linux(:armv7l, :glibc, :eabihf),
    MacOS(:x86_64),
    Windows(:i686),
    Windows(:x86_64)
]

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libtask", :libtask)
]

# Dependencies that must be installed before this package can be built
dependencies = [

]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)


text:
# ------------------------------------------------------------------------------------------
# # Julia is fast
# (Originally from https://juliabox.com under tutorials/intro-to-julia/short-
# version/05.Julia_is_fast.ipynb)
#
# Very often, benchmarks are used to compare languages.  These benchmarks can lead to long
# discussions, first as to exactly what is being benchmarked and secondly what explains the
# differences.  These simple questions can sometimes get more complicated than you at first
# might imagine.
#
# The purpose of this notebook is for you to see a simple benchmark for yourself.  One can
# read the notebook and see what happened on the author's Macbook Pro with a 4-core Intel
# Core I7, or run the notebook yourself.
#
# (This material began life as a wonderful lecture by Steven Johnson at MIT:
# https://github.com/stevengj/18S096/blob/master/lectures/lecture1/Boxes-and-
# registers.ipynb.)
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# # Outline of this notebook
#
# - Define the sum function
# - Implementations & benchmarking of sum in...
#     - C (hand-written)
#     - python (built-in)
#     - python (numpy)
#     - python (hand-written)
#     - Julia (built-in)
#     - Julia (hand-written)
#     - Julia (hand-written with SIMD)
# - Summary of benchmarks
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# # `sum`: An easy enough function to understand
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# Consider the  **sum** function `sum(a)`, which computes
# $$
# \mathrm{sum}(a) = \sum_{i=1}^n a_i,
# $$
# where $n$ is the length of `a`.
# ------------------------------------------------------------------------------------------

a = rand(10^7) # 1D vector of random numbers, uniform on [0,1)
a'

sum(a)

# ------------------------------------------------------------------------------------------
# The expected result is 5 * 10^6, since the mean of each entry is 0.5
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# # Benchmarking a few ways in a few languages
# ------------------------------------------------------------------------------------------

using Pkg
for p in ("BenchmarkTools","Plots","PyCall","Conda")
    haskey(Pkg.installed(),p) || Pkg.add(p)
end

@time sum(a)

@time sum(a)

@time sum(a)

# ------------------------------------------------------------------------------------------
# The `@time` macro can yield noisy results, so it's not our best choice for benchmarking!
#
# Luckily, Julia has a `BenchmarkTools.jl` package to make benchmarking easy and accurate:
# ------------------------------------------------------------------------------------------

using BenchmarkTools

# ------------------------------------------------------------------------------------------
# #  1. The C language
#
# C is often considered the gold standard: difficult on the human, nice for the machine.
# Getting within a factor of 2 of C is often satisfying. Nonetheless, even within C, there
# are many kinds of optimizations possible that a naive C writer may or may not get the
# advantage of.
#
# The current author does not speak C, so he does not read the cell below, but is happy to
# know that you can put C code in a Julia session, compile it, and run it. Note that the
# `"""` wrap a multi-line string.
# ------------------------------------------------------------------------------------------

using Libdl

C_code = """
#include <stddef.h>
double c_sum(size_t n, double *X) {
    double s = 0.0;
    for (size_t i = 0; i < n; ++i) {
        s += X[i];
    }
    return s;
}
"""

const Clib = tempname()   # make a temporary file
# compile to a shared library by piping C_code to gcc
# (works only if you have gcc installed):
open(`gcc -std=c99 -fPIC -O3 -msse3 -xc -shared -o $(Clib * "." * Libdl.dlext) -`, "w") do f
    print(f, C_code)
end
# define a Julia function that calls the C function:
c_sum(X::Array{Float64}) = ccall(("c_sum", Clib), Float64, (Csize_t, Ptr{Float64}), length(X), X)

c_sum(a)

c_sum(a) ≈ sum(a) # type \approx and then <TAB> to get the ≈ symbolb

c_sum(a) - sum(a)

≈  # alias for the `isapprox` function

?isapprox

# ------------------------------------------------------------------------------------------
# We can now benchmark the C code directly from Julia:
# ------------------------------------------------------------------------------------------

c_bench = @benchmark c_sum($a)

println("C: Fastest time was $(minimum(c_bench.times) / 1e6) msec")

d = Dict()  # a "dictionary", i.e. an associative array
d["C"] = minimum(c_bench.times) / 1e6  # in milliseconds
d

using Plots, Statistics
default(fmt = :png)

t = c_bench.times / 1e6 # times in milliseconds
m, σ = minimum(t), std(t)

histogram(t, bins=500,
    xlim=(m - 0.01, m + σ),
    xlabel="milliseconds", ylabel="count", label="")

# ------------------------------------------------------------------------------------------
# # 2. Python's built in `sum`
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# The `PyCall` package provides a Julia interface to Python:
# ------------------------------------------------------------------------------------------

using PyCall

# Call a low-level PyCall function to get a Python list, because
# by default PyCall will convert to a NumPy array instead (we benchmark NumPy below):

apy_list = PyCall.array2py(a)

# get the Python built-in "sum" function:
pysum = pybuiltin("sum")

pysum(a)

pysum(a) ≈ sum(a)

py_list_bench = @benchmark $pysum($apy_list)

d["Python built-in"] = minimum(py_list_bench.times) / 1e6
d

# ------------------------------------------------------------------------------------------
# # 3. Python: `numpy`
#
# ## Takes advantage of hardware "SIMD", but only works when it works.
#
# `numpy` is an optimized C library, callable from Python.
# It may be installed within Julia as follows:
# ------------------------------------------------------------------------------------------

using Conda

numpy_sum = pyimport("numpy")."sum"
apy_numpy = PyObject(a) # converts to a numpy array by default

numpy_sum(apy_list) # python thing

py_numpy_bench = @benchmark $numpy_sum($apy_numpy)

numpy_sum(apy_list) ≈ sum(a)

d["Python numpy"] = minimum(py_numpy_bench.times) / 1e6
d

# ------------------------------------------------------------------------------------------
# # 4. Python, hand-written
# ------------------------------------------------------------------------------------------

py"""
def py_sum(A):
    s = 0.0
    for a in A:
        s += a
    return s
"""

sum_py = py"py_sum"

sum_py(apy_list)

py_hand = @benchmark $sum_py($apy_list)

sum_py(apy_list) ≈ sum(a)

d["Python hand-written"] = minimum(py_hand.times) / 1e6
d

# ------------------------------------------------------------------------------------------
# # 5. Julia (built-in)
#
# ## Written directly in Julia, not in C!
# ------------------------------------------------------------------------------------------

@which sum(a)

sum(a)

j_bench = @benchmark sum($a)

d["Julia built-in"] = minimum(j_bench.times) / 1e6
d

# ------------------------------------------------------------------------------------------
# # 6. Julia (hand-written)
# ------------------------------------------------------------------------------------------

function mysum(A)
    s = 0.0 # s = zero(eltype(a))
    for a in A
        s += a
    end
    s
end

mysum(a)

j_bench_hand = @benchmark mysum($a)

d["Julia hand-written"] = minimum(j_bench_hand.times) / 1e6
d

# ------------------------------------------------------------------------------------------
# # 7. Julia (hand-written w. simd)
# ------------------------------------------------------------------------------------------

function mysum_simd(A)
    s = 0.0 # s = zero(eltype(A))
    @simd for a in A
        s += a
    end
    s
end

mysum_simd(a)

j_bench_hand_simd = @benchmark mysum_simd($a)

d["Julia hand-written simd"] = minimum(j_bench_hand_simd.times) / 1e6
d

# ------------------------------------------------------------------------------------------
# # Summary
# ------------------------------------------------------------------------------------------

for (key, value) in sort(collect(d), by=last)
    println(rpad(key, 25, "."), lpad(round(value, digits=1), 6, "."))
end

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
using Documenter, CatmullClark

makedocs(sitename = "CatmullClark Module Documentation", format = Documenter.HTML(prettyurls = false))


text:
"""
Structure for MH options
"""
struct MHOptions <: Options
    samplesteps::Int64
    annealsteps::Int64
    warmupsteps::Int64
    maxtime::Float64
    temp::Float64
    tempanneal::Float64
end
"""
Structure for MH results
"""
struct Fit <: Results
    param::Array
    ll::Array
    parml::Array
    llml::Float64
    ppd::Array
    pwaic::Array
    prior::Float64
    accept::Int
    total::Int
end

"""
Structure for parameter statistics results
"""
struct Stats
    meanparam::Array
    stdparam::Array
    medparam::Array
    madparam::Array
    qparam::Array
    corparam::Array
    covparam::Array
    covlogparam::Array
end

"""
run_mh(data,model,options)

Run Metropolis-Hastings MCMC algorithm and compute statistics of results
data is a HistogramData structure
and must have a datahistogram function
model is a StochasticGR model with a logprior function
model and data must have a likelihoodfn function
"""
function run_mh(data::HistogramData,model::StochasticGRmodel,options::MHOptions)
    fit,waic = metropolis_hastings(data,model,options)
    if options.samplesteps > 0
        stats = compute_stats(fit.param)
    else
        stats = 0
    end
    return fit, stats, waic
end
"""
run_mh(data,model,options,nchains)

"""
function run_mh(data::HistogramData,model::StochasticGRmodel,options::MHOptions,nchains)
    if nchains == 1
        return run_mh(data,model,options)
    else
        sd = run_chains(data,model,options,nchains)
        chain = extract_chain(sd)
        waic = pooled_waic(chain)
        fit = merge_fit(chain)
        stats = compute_stats(fit.param)
        return fit, stats, waic
    end
end
"""
run_chains(data,model,options,nchains)

"""
function run_chains(data,model,options,nchains)
    sd = Array{Future,1}(undef,nchains)
    for chain in 1:nchains
        sd[chain] = @spawn metropolis_hastings(data,model,options)
    end
    return sd
end

"""
metropolis_hastings(param,data,model,options)

Metropolis-Hastings MCMC algorithm
param = array of parameters to be fit
model, data, and options are structures

stochastic kinetic transcription models
Data can include ON and OFF MS2/PP7 reporter time distributions from live cell recordings,
scRNA or FISH mRNA data, and burst correlations between alleles
"""
function metropolis_hastings(data,model,options)
    param,d = initial_proposal(model)
    ll,predictions = loglikelihood(param,data,model)
    totalsteps = options.warmupsteps + options.samplesteps + options.annealsteps
    if options.annealsteps > 0
        param,parml,ll,llml,predictions,temp = anneal(predictions,param,param,ll,ll,d,model.proposal,data,model,options.annealsteps,options.temp,options.tempanneal,time(),options.maxtime*options.annealsteps/totalsteps)
    else
        temp = options.temp
    end
    if options.warmupsteps > 0
        param,parml,ll,llml,d,proposalcv,predictions = warmup(predictions,param,param,ll,ll,d,model.proposal,data,model,options.warmupsteps,temp,time(),options.maxtime*options.warmupsteps/totalsteps)
    else
        parml = param
        llml = ll
        proposalcv = model.proposal
    end
    fit=sample(predictions,param,parml,ll,llml,d,proposalcv,data,model,options.samplesteps,temp,time(),options.maxtime*options.samplesteps/totalsteps)
    waic = compute_waic(fit.ppd,fit.pwaic,data)
    return fit, waic
end
"""
function anneal(predictions,param,parml,ll,llml,d,proposalcv,data,model,samplesteps,temp,t1,maxtime)

"""
function anneal(predictions,param,parml,ll,llml,d,proposalcv,data,model,samplesteps,temp,tempanneal,t1,maxtime)
    parout = Array{Float64,2}(undef,length(param),samplesteps)
    prior = logprior(param,model)
    step = 0
    annealrate = 3/samplesteps
    anneal = 1 - annealrate  # annealing rate with time constant of samplesteps/3
    while  step < samplesteps && time() - t1 < maxtime
        step += 1
        _,predictions,param,ll,prior = mhstep(predictions,param,ll,prior,d,proposalcv,model,data,tempanneal)
        if ll < llml
            llml,parml = ll,param
        end
        parout[:,step] = param
        tempanneal = anneal * tempanneal + annealrate * temp
    end
    return param,parml,ll,llml,predictions,temp
end

"""
warmup(predictions,param,rml,ll,llml,d,sigma,data,model,samplesteps,temp,t1,maxtime)

"""
function warmup(predictions,param,parml,ll,llml,d,proposalcv,data,model,samplesteps,temp,t1,maxtime)
    parout = Array{Float64,2}(undef,length(param),samplesteps)
    prior = logprior(param,model)
    step = 0
    while  step < samplesteps && time() - t1 < maxtime
        step += 1
        _,predictions,param,ll,prior = mhstep(predictions,param,ll,prior,d,proposalcv,model,data,temp)
        if ll < llml
            llml,parml = ll,param
        end
        parout[:,step] = param
    end
    covlogparam = cov(log.(parout[:,1:step]'))
    if isposdef((2.4)^2 / length(param)*covlogparam)
        d=proposal(param,covlogparam)
        proposalcv = covlogparam
    end
    return param,parml,ll,llml,d,proposalcv,predictions
end
"""
sample(predictions,param,parml,ll,llml,d,proposalcv,data,model,samplesteps,temp,t1,maxtime)

"""
function sample(predictions,param,parml,ll,llml,d,proposalcv,data,model,samplesteps,temp,t1,maxtime)
    llout = Array{Float64,1}(undef,samplesteps)
    pwaic = (0,log.(max.(predictions,eps(Float64))),zeros(length(predictions)))
    ppd = zeros(length(predictions))
    accepttotal = 0
    parout = Array{Float64,2}(undef,length(param),samplesteps)
    prior = logprior(param,model)
    step = 0
    while  step < samplesteps && time() - t1 < maxtime
        step += 1
        accept,predictions,param,ll,prior = mhstep(predictions,param,ll,prior,d,proposalcv,model,data,temp)
        if ll < llml
            llml,parml = ll,param
        end
        parout[:,step] = param
        llout[step] = ll
        accepttotal += accept
        pwaic = update_waic(ppd,pwaic,predictions)
    end
    pwaic = step > 1 ? pwaic[3]/(step -1) :  pwaic[3]
    ppd /= step
    Fit(parout[:,1:step],llout[1:step],parml,llml,ppd,pwaic,prior,accepttotal,step)
end

"""
mhstep(predictions,param,ll,prior,d,sigma,model,data,temp)

ll is negative log likelihood
"""
function mhstep(predictions,param,ll,prior,d,proposalcv,model,data,temp)
    paramt,dt = proposal(d,proposalcv)
    priort = logprior(paramt,model)
    llt,predictionst = loglikelihood(paramt,data,model)
    if rand() < exp((ll + prior - llt - priort + mhfactor(param,d,paramt,dt))/temp)
        return 1,predictionst,paramt,llt,priort,dt
    else
        return 0,predictions,param,ll,prior,d
    end
end
"""
computewaic(ppd::Array{T},pwaic::Array{T},data) where {T}

Compute WAIC and SE of WAIC using ppd and pwaic computed in metropolis_hastings()
"""
function compute_waic(ppd::Array{T},pwaic::Array{T},data) where {T}
    hist = datahistogram(data)
    se = sqrt(sum(hist)*(var(ppd,weights(hist),corrected=false)+var(pwaic,weights(hist),corrected=false)))
    lppd = hist'*log.(max.(ppd,eps(T)))
    pwaic = hist'*pwaic
    return -2*(lppd-pwaic), 2*se
end
"""
aic(fit)

AIC
"""
aic(fit::Fit) = 2*length(fit.parml) + 2*fit.llml
aic(nparams::Int,llml) = 2*nparams + 2*llml

"""
initial_proposal(model)
return parameters to be fitted and an initial proposal distribution

"""
function initial_proposal(model)
    param = get_param(model)
    d=proposal(param,model.proposal)
    return param,d
end
"""
proposal(d,cv)
return rand(d) and proposal distribution for cv (vector or covariance)

"""
function proposal(d::MultivariateDistribution,cv::Matrix)
    param = rand(d)
    return param, proposal(param,cv)
end
function proposal(d::Product,cv)
    param = rand(d)
    return param, proposal(param,cv)
end
"""
proposal(param,cv)
return proposal distribution specified by location and scale

"""
# function proposal(param,cv::Vector)
#     d = Array{Truncated{Normal{Float64},Continuous}}(undef,0)
#     for i in eachindex(param)
#         push!(d,truncated(Normal(param[i],cv[i]*param[i]),0.,1000.))
#     end
#     product_distribution(d)
# end
function proposal(param::Vector,cv::Float64)
    d = Vector{LogNormal{Float64}}(undef,0)
    for i in eachindex(param)
        push!(d,LogNormal(log(max(param[i],1e-100))-.5*log(1+cv^2),sqrt(log(1+cv^2))))
    end
    product_distribution(d)
end

function proposal(param::Vector,cv::Vector)
    d = Vector{LogNormal{Float64}}(undef,0)
    for i in eachindex(param)
        push!(d,LogNormal(log(max(param[i],1e-100))-.5*log(1+cv[i]^2),sqrt(log(1+cv[i]^2))))
    end
    product_distribution(d)
end

function proposal(param::Vector,cov::Matrix)
    c = (2.4)^2 / length(param)
    if isposdef(c*cov)
        return MvLogNormal(log.(max.(param,eps(Float64))) - .5*c*diag(cov),c*cov)
    else
        return MvLogNormal(log.(max.(param,eps(Float64))) - .5*c*diag(cov),sqrt.(abs.(diag(c*cov))))
    end
end

"""
mhfactor(param,d,paramt,dt)
Metropolis-Hastings correction factor for asymmetric proposal distribution
"""
mhfactor(param,d,paramt,dt) = logpdf(dt,param)-logpdf(d,paramt) #pdf(dt,param)/pdf(d,paramt)

"""
update_waic(ppd,pwaic,predictions)

"""
function update_waic(ppd,pwaic,predictions)
    ppd .+= predictions
    var_update(pwaic,log.(max.(predictions,eps(Float64))))
end
# Functions for parallel computing on multiple chains
"""
extract_chain(sdspawn)

"""
function extract_chain(sdspawn)
    chain = Array{Tuple,1}(undef,length(sdspawn))
    for i in eachindex(sdspawn)
        chain[i] = fetch(sdspawn[i])
    end
    return chain
end
"""
collate_fit(chain)
collate fits from multiple metropolis_hastings chains
into and array of Fit structures

"""
function collate_fit(chain)
    fits = Array{Fit,1}(undef,length(chain))
    for i in eachindex(chain)
        fits[i] = chain[i][1]
    end
    return fits
end
"""
collate_waic(chain)
collate waic results from multiple chains
into a 2D array
"""
function collate_waic(chain)
    waic = Array{Float64,2}(undef,length(chain),3)
    for i in eachindex(chain)
        waic[i,1] = chain[i][2][1]
        waic[i,2] = chain[i][2][2]
        waic[i,3] = chain[i][1].total
    end
    return waic
    # mean(waic,dims=1),median(waic,dims=1), mad(waic[:,1],normalize=false),waic
end
function pooled_waic(chain)
    waic = collate_waic(chain)
    return pooled_mean(waic[:,1],waic[:,3]), pooled_std(waic[:,2],waic[:,3])
end

function merge_param(fits::Vector)
    param = fits[1].param
    for i in 2:length(fits)
        param = [param fits[i].param]
    end
    return param
end
function merge_ll(fits::Vector)
    ll = fits[1].ll
    for i in 2:length(fits)
        ll = [ll; fits[i].ll]
    end
    return ll
end

merge_fit(chain::Array{Tuple,1}) = merge_fit(collate_fit(chain))

function merge_fit(fits::Array{Fit,1})
    param = merge_param(fits)
    ll = merge_ll(fits)
    parml,llml = find_ml(fits)
    ppd = fits[1].ppd
    pwaic = fits[1].pwaic
    prior = fits[1].prior
    accept = fits[1].accept
    total = fits[1].total
    for i in 2:length(fits)
        accept += fits[i].accept
        total += fits[i].total
        prior += fits[i].prior
        ppd += fits[i].ppd
        pwaic += fits[i].pwaic
    end
    Fit(param,ll,parml,llml,ppd/length(fits),pwaic/length(fits),prior/length(fits),accept,total)
end

"""
compute_stats(fit::Fit)

Compute mean, std, median, mad, quantiles and correlations, covariances of parameters
"""
function compute_stats(param::Array{Float64,2})
    meanparam = mean(param,dims=2)
    stdparam = std(param,dims=2)
    corparam = cor(param')
    covparam = cov(param')
    covlogparam = cov(log.(param'))
    medparam = median(param,dims=2)
    np = size(param,1)
    madparam = Array{Float64,1}(undef,np)
    qparam = Array{Array,1}(undef,np)
    for i in 1:np
        madparam[i] = mad(param[i,:],normalize=false)
        qparam[i] = quantile(param[i,:],[.025;.5;.975])
    end
    Stats(meanparam,stdparam,medparam,madparam,qparam,corparam,covparam,covlogparam)
end

"""
find_ml(fits)
find the maximum likelihood out of multiple chains
"""
function find_ml(fits::Array)
    llml = Array{Float64,1}(undef,length(fits))
    for i in eachindex(fits)
        llml[i] = fits[i].llml
    end
    return fits[argmin(llml)].parml, llml[argmin(llml)]
end

"""
loglikelihood(param,data,model)

Calls likelihoodfn and datahistogram
provided for each data and model type
"""
function loglikelihood(param,data,model)
    predictions = likelihoodfn(param,data,model)
    hist = datahistogram(data)
    return crossentropy(predictions,hist),predictions
end
"""
crossentropy(predictions::Array{T},data::Array{T}) where {T}

compute crosscentropy between data histogram and likelilhood
"""
function crossentropy(predictions::Array{T1},hist::Array{T2}) where {T1,T2}
        lltot = hist' * log.(max.(predictions,eps(T1)))
        isfinite(lltot) ? -lltot : Inf
end
# function crossentropy(predictions::Array{T1},hist::Array{T2}) where {T1,T2}
#     if minimum(predictions) > -0.00000001
#         lltot = hist' * log.(normalize_histogram(max.(predictions,eps(T1))))
#         return isfinite(lltot) ? -lltot : Inf
#     else
#         return Inf
#     end
# end
"""
hist_entropy(hist)

Compute entropy of a normalized histogram
"""
function hist_entropy(hist::Array{Float64,1})
    -hist' * (log.(normalize_histogram(max.(hist,eps(Float64)))))
end
function hist_entropy(hist::Array{Array,1})
    l = 0
    for h in hist
        l +=hist_entropy(h)
    end
    return l
end
"""
deviance(predictions,data)
deviance(fit,data,model)

Deviance
"""
# deviance(predictions::Array,data::HistogramData) = -2*datahistogram(data)'*(log.(max.(predictions,eps(Float64))) - log.(max.(datapdf(data),eps(Float64))))
deviance(predictions::Array,data::HistogramData) = deviance(predictions,datapdf(data))

# (crossentropy(predictions,datahistogram(data)) - hist_entropy(datahistogram(data)))
function deviance(fit::Fit,data::HistogramData,model)
    h=likelihoodfn(fit.parml,data,model)
    deviance(h,data)
end

function deviance(data::HistogramData,model::StochasticGRmodel)
    h = likelihoodfn(model.rates[model.fittedparam],data,model)
    deviance(h,data)
end

deviance(predictions::Array,hist::Array) = 2*hist'*(log.(max.(hist,eps(Float64)))-log.(max.(predictions,eps(Float64))))

max_stars_repo_name: nih-niddk-mbs/StochasticGene

max_issues_repo_name: nih-niddk-mbs/StochasticGene

max_forks_repo_name: nih-niddk-mbs/StochasticGene


text:
using ExecutableSpecifications.Gherkin: parsefeature, issuccessful

@testset "Feature              " begin
    @testset "Feature description" begin
        @testset "Read feature description; Description matches input" begin
            text = """
            Feature: This is a feature
            """

            result = parsefeature(text)

            @test issuccessful(result)
            feature = result.value
            @test feature.header.description == "This is a feature"
        end

        @testset "Read another feature description; Description matches input" begin
            text = """
            Feature: This is another feature
            """

            result = parsefeature(text)

            @test issuccessful(result)
            feature = result.value
            @test feature.header.description == "This is another feature"
        end

        @testset "Read long feature description" begin
            text = """
            Feature: This is another feature
              This is the long description.
              It contains several lines.
            """

            result = parsefeature(text)

            @test issuccessful(result)
            feature = result.value
            @test "This is the long description." in feature.header.long_description
            @test "It contains several lines." in feature.header.long_description
        end

        @testset "Scenarios are not part of the feature description" begin
            text = """
            Feature: This is another feature
                This is the long description.
                It contains several lines.

                Scenario: Some scenario
                    Given a precondition
            """

            result = parsefeature(text)

            @test issuccessful(result)
            feature = result.value
            @test ("Given a precondition" in feature.header.long_description) == false
        end
    end

    @testset "Read scenarios" begin
        @testset "Feature has one scenario; one scenarios is parsed" begin
            text = """
            Feature: This feature has one scenario

                Scenario: This is one scenario
                    Given a precondition
            """

            result = parsefeature(text)

            @test issuccessful(result)
            feature = result.value
            @test length(feature.scenarios) == 1
        end

        @testset "Feature has two scenarios; two scenarios are parsed" begin
            text = """
            Feature: This feature has one scenario

                Scenario: This is one scenario
                    Given a precondition

                Scenario: This is a second scenario
                    Given a precondition
            """

            result = parsefeature(text)

            @test issuccessful(result)
            feature = result.value
            @test length(feature.scenarios) == 2
        end

        @testset "Feature has one scenario; The description is read from the scenario" begin
            text = """
            Feature: This feature has one scenario

                Scenario: This is one scenario
                    Given a precondition
            """

            result = parsefeature(text)

            @test issuccessful(result)
            feature = result.value
            @test feature.scenarios[1].description == "This is one scenario"
        end

        @testset "Feature has two scenarios; two scenarios are parsed" begin
            text = """
            Feature: This feature has one scenario

                Scenario: This is one scenario
                    Given a precondition

                Scenario: This is a second scenario
                    Given a precondition
            """

            result = parsefeature(text)

            @test issuccessful(result)
            feature = result.value
            @test feature.scenarios[1].description == "This is one scenario"
            @test feature.scenarios[2].description == "This is a second scenario"
        end

        @testset "Scenario with three steps; The parsed scenario has three steps" begin
            text = """
            Feature: This feature has one scenario

                Scenario: This is one scenario
                    Given a precondition
                    When an action is performed
                    Then some postcondition holds
            """

            result = parsefeature(text)

            @test issuccessful(result)
            feature = result.value
            @test length(feature.scenarios[1].steps) == 3
        end

        @testset "Scenario with one step; The parsed scenario has one step" begin
            text = """
            Feature: This feature has one scenario

                Scenario: This is one scenario
                    Given a precondition
            """

            result = parsefeature(text)

            @test issuccessful(result)
            feature = result.value
            @test length(feature.scenarios[1].steps) == 1
        end

        @testset "Feature has a scenario outline; The feature scenarios list has one element" begin
            text = """
            Feature: This feature has one scenario

                Scenario Outline: This is one scenario outline
                    Given a precondition with field <Foo>

                Examples:
                    | Foo |
                    | 1   |
                    | 2   |
            """

            result = parsefeature(text)

            @test issuccessful(result)
            feature = result.value
            @test length(feature.scenarios) == 1
        end

        @testset "Feature has a scenario outline and a normal scenario; Two scenarios are parsed" begin
            text = """
            Feature: This feature has one scenario

                Scenario Outline: This is one scenario outline
                    Given a precondition with field <Foo>

                Examples:
                    | Foo |
                    | 1   |
                    | 2   |

                Scenario: A normal scenario
                    Given some precondition
            """

            result = parsefeature(text)

            @test issuccessful(result)
            feature = result.value
            @test length(feature.scenarios) == 2
        end
    end

    @testset "Robustness" begin
        @testset "Many empty lines before scenario; Empty lines are ignored" begin
            text = """
            Feature: This feature has many empty lines between scenarios

                Scenario: This is one scenario
                    Given a precondition

                Scenario: This is another scenario
                    Given a precondition
            """

            result = parsefeature(text)

            @test issuccessful(result)
            feature = result.value
            @test length(feature.scenarios) == 2
        end
    end

    @testset "Malformed features" begin
        @testset "Scenario found before feature; Parse fails with feature expected" begin
            text = """
                Scenario: This is one scenario
                    Given a precondition

            Feature: This feature has one scenario

                Scenario: This is a second scenario
                    Given a precondition
            """

            result = parsefeature(text)

            @test issuccessful(result) == false
            @test result.reason == :unexpected_construct
            @test result.expected == :feature
            @test result.actual == :scenario
        end
    end
end


text:
module App

export start_app

using Web.Backend.Overview: constructor_overview
importall Bukdu

type OverviewController <: ApplicationController
    conn::Conn
end

function get_file(filename)
    loaded_file = join(readlines(open("../src/web/frontend/$filename")), "")
    function get_loaded_file(void)
        return loaded_file
    end
end

function start_app()
    overview_builder = constructor_overview()
    function overview(c::OverviewController)
        q = c[:params]
        # x, y = map(v -> parse(Int, v), (q[:x], q[:y]))
        # render(JSON, x + 2*y)
        return overview_builder()
    end

    Router() do
        get("/", OverviewController, overview)
    end

    Endpoint() do
        plug(Plug.Logger)
        plug(Plug.Static, at="/static", from="../src/web/frontend/static")
        plug(Router)
    end

    Bukdu.start(8080, getaddrinfo("0.0.0.0"))
end

end # module

max_stars_repo_name: Gollor/JuliaAutomation

max_issues_repo_name: Gollor/JuliaAutomation

max_forks_repo_name: Gollor/JuliaAutomation