text:
module clcircularplatecl_examples
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsDeforLinear.AlgoDeforLinearModule
using FinEtools.MeshExportModule
using Statistics: mean

# Clamped square plate with concentrated force
#  Data listed in the Simo 1990 paper 'A class of... '
#  Analytical solution for the vertical deflection under the load
#    Zienkiewicz, Taylor, The finite element method, fifth edition, volume 2,
#     analyt_sol = 3*(1-nu^2)*Magnitude*R^2/(4*pi*E*thickness^3);
E = 1e7;
nu = 0.3;
Magnitude = 10;
R = 100.0;
thickness = 2.5;
tolerance = 0.0001*thickness;
@show analyt_sol = 3*(1-nu^2)*Magnitude*R^2/(4*pi*E*thickness^3);

function clcircularplatecl_h8_full()
	nt = 1
	for nperradius in [2, 4, 8]
		nt=nt+1;
    	fens, fes  = Q4circlen(R, nperradius)
    	fens, fes = H8extrudeQ4(fens, fes, nt, (x, k) -> [x[1],x[2],k*thickness/nt])

    	MR = DeforModelRed3D
    	material = MatDeforElastIso(MR, E, nu)
    	femm = FEMMDeforLinear(MR, IntegDomain(fes, GaussRule(3, 2)), material)

    	boundaryfes = meshboundary(fes);
    	topl = selectelem(fens,boundaryfes, box =  [-Inf Inf -Inf Inf thickness thickness], inflate =   tolerance);
    	botl = selectelem(fens,boundaryfes, box =  [-Inf Inf -Inf Inf 0.0 0.0], inflate =   tolerance);
    	x0l = selectelem(fens,boundaryfes, box =  [0.0 0.0 -Inf Inf 0.0 thickness], inflate =   tolerance);
    	y0l = selectelem(fens,boundaryfes, box =  [-Inf Inf 0.0 0.0 0.0 thickness], inflate =   tolerance);
    	cyll = setdiff(1:count(boundaryfes),topl,botl,x0l,y0l);

    	geom = NodalField(fens.xyz)
    	u = NodalField(zeros(size(fens.xyz,1), 3)) # displacement field

    	x0nl = connectednodes(subset(boundaryfes, x0l))
    	setebc!(u,x0nl,true,1,0.0)

    	y0nl = connectednodes(subset(boundaryfes, y0l))
    	setebc!(u,y0nl,true,2,0.0)

    	cylnl = connectednodes(subset(boundaryfes, cyll))
    	setebc!(u,cylnl,true,1,0.0)
    	setebc!(u,cylnl,true,2,0.0)
    	setebc!(u,cylnl,true,3,0.0)

    	applyebc!(u)
    	numberdofs!(u)

    	cnl = selectnode(fens; box = [0 0 0 0 0 thickness], inflate  =  tolerance)

    	nfemm = FEMMBase(IntegDomain(FESetP1(reshape(cnl, length(cnl), 1)), PointRule()))
  		F = distribloads(nfemm, geom, u, ForceIntensity([0; 0; Magnitude/4/length(cnl)]), 3)

  		associategeometry!(femm, geom)

  		K = stiffness(femm, geom, u)

  		scattersysvec!(u, K\F)

  		u0z = mean(u.values[cnl, 3]);
  		println("Deflection under the load: $(round((u0z / analyt_sol)* 100000)/100000*100) %")

    	# File =  "clcircularplatecl_1_$(nperradius).vtk"
    	# vtkexportmesh(File, fens, fes; vectors = [("u", u.values)])
    	# @async run(`"paraview.exe" $File`)

    end

    true
end # clcircularplatecl_h8_full

function clcircularplatecl_h8_uri()
	nt = 1
	for nperradius in [2, 4, 8]
		nt=nt+1;
		fens, fes  = Q4circlen(R, nperradius)
		fens, fes = H8extrudeQ4(fens, fes, nt, (x, k) -> [x[1],x[2],k*thickness/nt])

    	MR = DeforModelRed3D
    	material = MatDeforElastIso(MR, E, nu)
    	femm = FEMMDeforLinear(MR, IntegDomain(fes, GaussRule(3, 1)), material)

    	boundaryfes = meshboundary(fes);
    	topl = selectelem(fens,boundaryfes, box =  [-Inf Inf -Inf Inf thickness thickness], inflate =   tolerance);
    	botl = selectelem(fens,boundaryfes, box =  [-Inf Inf -Inf Inf 0.0 0.0], inflate =   tolerance);
    	x0l = selectelem(fens,boundaryfes, box =  [0.0 0.0 -Inf Inf 0.0 thickness], inflate =   tolerance);
    	y0l = selectelem(fens,boundaryfes, box =  [-Inf Inf 0.0 0.0 0.0 thickness], inflate =   tolerance);
    	cyll = setdiff(1:count(boundaryfes),topl,botl,x0l,y0l);

    	geom = NodalField(fens.xyz)
    	u = NodalField(zeros(size(fens.xyz,1), 3)) # displacement field

    	x0nl = connectednodes(subset(boundaryfes, x0l))
    	setebc!(u,x0nl,true,1,0.0)

    	y0nl = connectednodes(subset(boundaryfes, y0l))
    	setebc!(u,y0nl,true,2,0.0)

    	cylnl = connectednodes(subset(boundaryfes, cyll))
    	setebc!(u,cylnl,true,1,0.0)
    	setebc!(u,cylnl,true,2,0.0)
    	setebc!(u,cylnl,true,3,0.0)

    	applyebc!(u)
    	numberdofs!(u)

    	cnl = selectnode(fens; box = [0 0 0 0 0 thickness], inflate  =  tolerance)

    	nfemm = FEMMBase(IntegDomain(FESetP1(reshape(cnl, length(cnl), 1)), PointRule()))
  		F = distribloads(nfemm, geom, u, ForceIntensity([0; 0; Magnitude/4/length(cnl)]), 3)

  		associategeometry!(femm, geom)

  		K = stiffness(femm, geom, u)

  		scattersysvec!(u, K\F)

  		u0z = mean(u.values[cnl, 3]);
  		println("Deflection under the load: $(round((u0z / analyt_sol)* 100000)/100000*100) %")

    	# File =  "clcircularplatecl_1_$(nperradius).vtk"
    	# vtkexportmesh(File, fens, fes; vectors = [("u", u.values)])
    	# @async run(`"paraview.exe" $File`)

    end

    true
end # clcircularplatecl_h8_uri

function clcircularplatecl_h8_ms()
    nt = 1
    for nperradius in [2, 4, 8]
    	nt=nt+1;
    	fens, fes  = Q4circlen(R, nperradius)
    	fens, fes = H8extrudeQ4(fens, fes, nt, (x, k) -> [x[1],x[2],k*thickness/nt])

    	MR = DeforModelRed3D
    	material = MatDeforElastIso(MR, E, nu)
    	femm = FEMMDeforLinearMSH8(MR, IntegDomain(fes, GaussRule(3, 2)), material)

    	boundaryfes = meshboundary(fes);
    	topl = selectelem(fens,boundaryfes, box =  [-Inf Inf -Inf Inf thickness thickness], inflate =   tolerance);
    	botl = selectelem(fens,boundaryfes, box =  [-Inf Inf -Inf Inf 0.0 0.0], inflate =   tolerance);
    	x0l = selectelem(fens,boundaryfes, box =  [0.0 0.0 -Inf Inf 0.0 thickness], inflate =   tolerance);
    	y0l = selectelem(fens,boundaryfes, box =  [-Inf Inf 0.0 0.0 0.0 thickness], inflate =   tolerance);
    	cyll = setdiff(1:count(boundaryfes),topl,botl,x0l,y0l);

    	geom = NodalField(fens.xyz)
    	u = NodalField(zeros(size(fens.xyz,1), 3)) # displacement field

    	x0nl = connectednodes(subset(boundaryfes, x0l))
    	setebc!(u,x0nl,true,1,0.0)

    	y0nl = connectednodes(subset(boundaryfes, y0l))
    	setebc!(u,y0nl,true,2,0.0)

    	cylnl = connectednodes(subset(boundaryfes, cyll))
    	setebc!(u,cylnl,true,1,0.0)
    	setebc!(u,cylnl,true,2,0.0)
    	setebc!(u,cylnl,true,3,0.0)

    	applyebc!(u)
    	numberdofs!(u)

    	cnl = selectnode(fens; box = [0 0 0 0 0 thickness], inflate  =  tolerance)

    	nfemm = FEMMBase(IntegDomain(FESetP1(reshape(cnl, length(cnl), 1)), PointRule()))
  		F = distribloads(nfemm, geom, u, ForceIntensity([0; 0; Magnitude/4/length(cnl)]), 3)

  		associategeometry!(femm, geom)
  		@show minimum(femm.phis), maximum(femm.phis)

  		K = stiffness(femm, geom, u)

  		scattersysvec!(u, K\F)

  		u0z = mean(u.values[cnl, 3]);
  		println("Deflection under the load: $(round((u0z / analyt_sol)* 100000)/100000*100) %")

    	# File =  "clcircularplatecl_1_$(nperradius).vtk"
    	# vtkexportmesh(File, fens, fes; vectors = [("u", u.values)])
    	# @async run(`"paraview.exe" $File`)

    end

    true
end # clcircularplatecl_h8_ms

function clcircularplatecl_h8u_ms()
    nt = 1
    for nperradius in [2, 4, 8]
    	nt=nt+1;
    	fens,fes = T4quartercyln(R, thickness, nperradius, nt)
    	fens,fes = T4toH8(fens,fes)

    	MR = DeforModelRed3D
    	material = MatDeforElastIso(MR, E, nu)
    	femm = FEMMDeforLinearMSH8(MR, IntegDomain(fes, GaussRule(3, 2)), material)

    	boundaryfes = meshboundary(fes);
    	topl = selectelem(fens,boundaryfes, box =  [-Inf Inf -Inf Inf thickness thickness], inflate =   tolerance);
    	botl = selectelem(fens,boundaryfes, box =  [-Inf Inf -Inf Inf 0.0 0.0], inflate =   tolerance);
    	x0l = selectelem(fens,boundaryfes, box =  [0.0 0.0 -Inf Inf 0.0 thickness], inflate =   tolerance);
    	y0l = selectelem(fens,boundaryfes, box =  [-Inf Inf 0.0 0.0 0.0 thickness], inflate =   tolerance);
    	cyll = setdiff(1:count(boundaryfes),topl,botl,x0l,y0l);

    	geom = NodalField(fens.xyz)
    	u = NodalField(zeros(size(fens.xyz,1), 3)) # displacement field

    	x0nl = connectednodes(subset(boundaryfes, x0l))
    	setebc!(u,x0nl,true,1,0.0)

    	y0nl = connectednodes(subset(boundaryfes, y0l))
    	setebc!(u,y0nl,true,2,0.0)

    	cylnl = connectednodes(subset(boundaryfes, cyll))
    	setebc!(u,cylnl,true,1,0.0)
    	setebc!(u,cylnl,true,2,0.0)
    	setebc!(u,cylnl,true,3,0.0)

    	applyebc!(u)
    	numberdofs!(u)

    	cnl = selectnode(fens; box = [0 0 0 0 0 thickness], inflate  =  tolerance)

    	nfemm = FEMMBase(IntegDomain(FESetP1(reshape(cnl, length(cnl), 1)), PointRule()))
  		F = distribloads(nfemm, geom, u, ForceIntensity([0; 0; Magnitude/4/length(cnl)]), 3)

  		associategeometry!(femm, geom)
  		@show minimum(femm.phis), maximum(femm.phis)

  		K = stiffness(femm, geom, u)

  		scattersysvec!(u, K\F)

  		u0z = mean(u.values[cnl, 3]);
  		println("Deflection under the load: $(round((u0z / analyt_sol)* 100000)/100000*100) %")

    	File =  "clcircularplatecl_h8u_ms_$(nperradius).vtk"
    	vtkexportmesh(File, fens, fes; vectors = [("u", u.values)])
    	@async run(`"paraview.exe" $File`)

    end

    true
end # clcircularplatecl_h8u_ms

function clcircularplatecl_h8_export()
    nt = 1
    for nperradius in [2, 4, 8]
    	nt=nt+1;
    	fens, fes  = Q4circlen(R, nperradius)
    	fens, fes = H8extrudeQ4(fens, fes, nt, (x, k) -> [x[1],x[2],k*thickness/nt])

    	MR = DeforModelRed3D
    	material = MatDeforElastIso(MR, E, nu)
    	femm = FEMMDeforLinearMSH8(MR, IntegDomain(fes, GaussRule(3, 2)), material)

    	boundaryfes = meshboundary(fes);
    	topl = selectelem(fens,boundaryfes, box =  [-Inf Inf -Inf Inf thickness thickness], inflate =   tolerance);
    	botl = selectelem(fens,boundaryfes, box =  [-Inf Inf -Inf Inf 0.0 0.0], inflate =   tolerance);
    	x0l = selectelem(fens,boundaryfes, box =  [0.0 0.0 -Inf Inf 0.0 thickness], inflate =   tolerance);
    	y0l = selectelem(fens,boundaryfes, box =  [-Inf Inf 0.0 0.0 0.0 thickness], inflate =   tolerance);
    	cyll = setdiff(1:count(boundaryfes),topl,botl,x0l,y0l);

    	geom = NodalField(fens.xyz)
    	u = NodalField(zeros(size(fens.xyz,1), 3)) # displacement field

    	x0nl = connectednodes(subset(boundaryfes, x0l))
    	setebc!(u,x0nl,true,1,0.0)

    	y0nl = connectednodes(subset(boundaryfes, y0l))
    	setebc!(u,y0nl,true,2,0.0)

    	cylnl = connectednodes(subset(boundaryfes, cyll))
    	setebc!(u,cylnl,true,1,0.0)
    	setebc!(u,cylnl,true,2,0.0)
    	setebc!(u,cylnl,true,3,0.0)

    	applyebc!(u)
    	numberdofs!(u)

    	cnl = selectnode(fens; box = [0 0 0 0 0 thickness], inflate  =  tolerance)

    	AE = AbaqusExporter("clcircularplatecl_h8_export_$(nperradius)");
    	HEADING(AE, "Clamped square plate with concentrated force");
    	PART(AE, "part1");
    	END_PART(AE);
    	ASSEMBLY(AE, "ASSEM1");
    	INSTANCE(AE, "INSTNC1", "PART1");
    	NODE(AE, fens.xyz);
    	ELEMENT(AE, "c3d8rh", "AllElements", 1, connasarray(femm.integdomain.fes))
    	NSET_NSET(AE, "cnl", cnl)
    	NSET_NSET(AE, "x0nl", x0nl)
    	NSET_NSET(AE, "y0nl", y0nl)
    	NSET_NSET(AE, "cylnl", cylnl)
    	ORIENTATION(AE, "GlobalOrientation", vec([1. 0 0]), vec([0 1. 0]));
    	SOLID_SECTION(AE, "elasticity", "GlobalOrientation", "AllElements", "Hourglassctl");
    	END_INSTANCE(AE);
    	END_ASSEMBLY(AE);
    	MATERIAL(AE, "elasticity")
    	ELASTIC(AE, E, nu)
    	SECTION_CONTROLS(AE, "Hourglassctl", "HOURGLASS=ENHANCED")
    	STEP_PERTURBATION_STATIC(AE)
    	BOUNDARY(AE, "ASSEM1.INSTNC1.x0nl", 1)
    	BOUNDARY(AE, "ASSEM1.INSTNC1.y0nl", 2)
    	BOUNDARY(AE, "ASSEM1.INSTNC1.cylnl", 1)
    	BOUNDARY(AE, "ASSEM1.INSTNC1.cylnl", 2)
    	BOUNDARY(AE, "ASSEM1.INSTNC1.cylnl", 3)
    	CLOAD(AE, "ASSEM1.INSTNC1.cnl", 3, Magnitude/4/length(cnl))
    	END_STEP(AE)
    	close(AE)

    	# File =  "clcircularplatecl_h8_export_$(nperradius).vtk"
    	# vtkexportmesh(File, fens, fes; vectors = [("u", u.values)])
    	# @async run(`"paraview.exe" $File`)

    end

    true
end # clcircularplatecl_h8_export

#     true
# end # twisted_beam_export

function allrun()
    println("#####################################################")
    println("# clcircularplatecl_h8_full ")
    clcircularplatecl_h8_full()
    println("#####################################################")
    println("# clcircularplatecl_h8_uri ")
    clcircularplatecl_h8_uri()
    println("#####################################################")
    println("# clcircularplatecl_h8_ms ")
    clcircularplatecl_h8_ms()
    println("#####################################################")
    println("# clcircularplatecl_h8_export ")
    clcircularplatecl_h8_export()
    return true
end # function allrun

end # module clcircularplatecl_examples


text:
module PowerFlowData

using DocStringExtensions
using InlineStrings: InlineString1, InlineString3, InlineString15
using Parsers: Parsers, xparse, checkdelim!
using Parsers: codes, eof, invalid, invaliddelimiter, newline, valueok, peekbyte
using PrettyTables: pretty_table
using Tables

export parse_network
export Network
export CaseID, Buses, Buses30, Buses33, Branches, Branches30, Branches33
export FixedShunts, Loads, Generators, Transformers, AreaInterchanges
export SwitchedShunts, SwitchedShunts30, SwitchedShunts33
export TwoTerminalDCLines, TwoTerminalDCLines30, TwoTerminalDCLines33
export VSCDCLines, ImpedanceCorrections
export MultiTerminalDCLines, MultiTerminalDCLine, DCLineID, DCLineID30, DCLineID33
export ACConverters, DCBuses, DCLinks
export MultiSectionLineGroups, MultiSectionLineGroups30, MultiSectionLineGroups33
export Zones, InterAreaTransfers, Owners
export FACTSDevices, FACTSDevices30, FACTSDevices33

include("debug.jl")
include("types.jl")
include("parsing.jl")

end  # module

ext: jl


text:
module Juniper

using JSON
using LinearAlgebra
using Random
const JUNIPER_RNG = MersenneTwister(1)

using Distributed
using Statistics

import MutableArithmetics
const MA = MutableArithmetics

using MathOptInterface
const MOI = MathOptInterface
const MOIU = MOI.Utilities

# functions
const SVF = MOI.SingleVariable
const SAF = MOI.ScalarAffineFunction{Float64}
const SQF = MOI.ScalarQuadraticFunction{Float64}
const VECTOR = MOI.VectorOfVariables

include("types.jl")

function Base.show(io::IO, opts::SolverOptions)
    longest_field_name = maximum([length(string(fname)) for fname in fieldnames(SolverOptions)])+2
    for name in fieldnames(SolverOptions)
        sname = string(name)
        pname = sname*repeat(" ", longest_field_name-length(sname))
        if getfield(opts,name) === nothing
            println(io, pname, ": NA")
        else
            println(io, pname, ": ", getfield(opts,name))
        end
    end
end

include("util.jl")
include("filter.jl")
include("printing.jl")
include("solver.jl")
include("model.jl")
include("BnBTree.jl")
include("MOI_wrapper/MOI_wrapper.jl")

end


text:
using MacroTools: @q

# TODO: `pure` functions should probably have type signatures, to avoid being
# over-general. For now this makes things convenient.

function purem(P, f)
  quote
    Mjolnir.abstract(::$(esc(P)), ::AType{typeof($(esc(f)))}, args::Const...) =
      Const($(esc(f))(map(arg -> arg.value, args)...))
    Mjolnir.abstract(::$(esc(P)), ::AType{typeof($(esc(f)))}, args...) =
      Core.Compiler.return_type($(esc(f)), Tuple{widen.(args)...})
  end
end

macro pure(P, fs)
  if isexpr(fs, :tuple)
    Expr(:block, [purem(P, f) for f in fs.args]...)
  else
    purem(P, fs)
  end
end

named(arg) = isexpr(arg, :(::)) && length(arg.args) == 1 ? :($(gensym())::$(arg.args[1])) : arg

typeless(x) = MacroTools.postwalk(x -> isexpr(x, :(::), :kw) ? x.args[1] : x, x)
isvararg(x) = isexpr(x, :(::)) && namify(x.args[2]) == :Vararg

wraptype(x) = namify(x) in (:Partial, :Const, :AType, :Shape) ? x : :(Mjolnir.AType{<:$x})
wraptypes(x) = MacroTools.postwalk(x -> isexpr(x, :(::)) ? Expr(:(::), x.args[1], wraptype(x.args[2])) : x, x)

const_kw(x) = MacroTools.postwalk(x -> isexpr(x, :kw) ? Expr(:kw, x.args[1], :($Const($(x.args[2])))) : x, x)

unwrap_kw(kw::Const) = kw.value

function abstractm(ex, P, abstract)
  @capture(shortdef(ex), (name_(args__) = body_) |
                         (name_(args__) where {Ts__} = body_)) || error("Need a function definition")
  kw = length(args) > 1 && isexpr(args[1], :parameters) ? esc(popfirst!(args)) : nothing
  kw = const_kw(kw)
  f, T = isexpr(name, :(::)) ?
    (length(name.args) == 1 ? (esc(gensym()), esc(name.args[1])) : esc.(name.args)) :
    (esc(gensym()), :(Core.Typeof($(esc(name)))))
  Ts == nothing && (Ts = [])
  args = named.(args)
  argnames = Any[typeless(arg) for arg in args]
  !isempty(args) && isvararg(args[end]) && (argnames[end] = :($(argnames[end])...,))
  args = wraptypes.(args)
  args = esc.(args)
  argnames = esc.(argnames)
  Ts = esc.(Ts)
  fargs = kw == nothing ? [:($f::AType{<:$T}), args...] : [kw, :($f::AType{<:$T}), args...]
  func = @q f($(fargs...)) where $(Ts...) = $(esc(body))
  quote
    $func
    function Mjolnir.$abstract(::$(esc(P)), $f::AType{$T}, $(args...)) where $(Ts...)
      f($f, $(argnames...))
    end
    function Mjolnir.$abstract(::$(esc(P)), ::AType{Core.kwftype($T)}, kw, $f::AType{<:$T}, $(args...)) where $(Ts...)
      f($f, $(argnames...); unwrap_kw(kw)...)
    end
    nothing
  end
end

macro abstract(P, ex)
  abstractm(ex, P, :abstract)
end

macro partial(P, ex)
  abstractm(ex, P, :partial)
end


text:
# NOTE: there's a `kmeans!` function that updates centers, maybe a candidate
# for the `update` machinery. Same for `kmedoids!`
# NOTE: if the prediction is done on the original array, just the assignment
# should be returned, unclear what's the best way of doing this.

module Clustering_

import MLJBase
import MLJBase: @mlj_model
using ScientificTypes

import ..Clustering # strange sytax for lazy-loading

using Distances
using LinearAlgebra: norm

const C = Clustering

const KMeansDescription =
    """
    K-Means algorithm: find K centroids corresponding to K clusters in the data.
    """

const KMedoidsDescription =
    """
    K-Medoids algorithm: find K centroids corresponding to K clusters in the data.
    Unlike K-Means, the centroids are found among data points themselves."
    """

const KMFields =
    """
    ## Keywords

    * `k=3`     : number of centroids
    * `metric`  : distance metric to use
    """

"""
KMeans(; kwargs...)

$KMeansDescription

$KMFields

See also the [package documentation](http://juliastats.github.io/Clustering.jl/latest/kmeans.html).
"""
@mlj_model mutable struct KMeans <: MLJBase.Unsupervised
    k::Int = 3::(_ ≥ 2)
    metric::SemiMetric = SqEuclidean()
end

"""
KMedoids(; kwargs...)

$KMedoidsDescription

$KMFields

See also the [package documentation](http://juliastats.github.io/Clustering.jl/latest/kmedoids.html).
"""
@mlj_model mutable struct KMedoids <: MLJBase.Unsupervised
    k::Int = 3::(_ ≥ 2)
    metric::SemiMetric = SqEuclidean()
end

function MLJBase.fit(model::KMeans
                   , verbosity::Int
                   , X)
    # NOTE: using transpose here to get a LinearAlgebra.Transpose object which Kmeans can handle
    Xarray = transpose(MLJBase.matrix(X))

    result    = C.kmeans(Xarray, model.k; distance=model.metric)
    fitresult = result.centers # centers (p x k)
    cache     = nothing
    report    = (assignments=result.assignments,) # size n

    return fitresult, cache, report
end

MLJBase.fitted_params(::KMeans, fitresult) = (centers=fitresult,)

function MLJBase.transform(model::KMeans
                         , fitresult
                         , X)
    # pairwise distance from samples to centers
    X̃ = pairwise(model.metric, transpose(MLJBase.matrix(X)), fitresult, dims=2)
    return MLJBase.table(X̃, prototype=X)
end

function MLJBase.fit(model::KMedoids
                   , verbosity::Int
                   , X)

    # NOTE: using transpose=true will materialize the transpose (~ permutedims), KMedoids
    # does not yet accept LinearAlgebra.Transpose
    Xarray = MLJBase.matrix(X, transpose=true)
    # cost matrix: all the pairwise distances
    Carray    = pairwise(model.metric, Xarray, dims=2) # n x n
    result    = C.kmedoids(Carray, model.k)
    fitresult = view(Xarray, :, result.medoids) # medoids
    cache     = nothing
    report    = (assignments=result.assignments,) # size n

    return fitresult, cache, report
end

MLJBase.fitted_params(::KMedoids, fitresult) = (medoids=fitresult,)

function MLJBase.transform(model::KMedoids
                         , fitresult
                         , X)
    # pairwise distance from samples to medoids
    X̃ = pairwise(model.metric, MLJBase.matrix(X, transpose=true), fitresult, dims=2)
    return MLJBase.table(X̃, prototype=X)
end

####
#### Predict methods
####

function MLJBase.predict(model::Union{KMeans,KMedoids}, fitresult, Xnew)

    Xarray = MLJBase.matrix(Xnew)
    (n, p), k = size(Xarray), model.k

    pred = zeros(Int, n)
    @inbounds for i ∈ 1:n
        minv = Inf
        for j ∈ 1:k
            curv    = evaluate(model.metric,
                                  view(Xarray, i, :), view(fitresult, :, j))
            P       = curv < minv
            pred[i] =    j * P + pred[i] * !P # if P is true --> j
            minv    = curv * P +    minv * !P # if P is true --> curvalue
        end
    end
    return pred
end

####
#### METADATA
####

import ..metadata_pkg, ..metadata_model

metadata_pkg.((KMeans, KMedoids),
    name="Clustering",
    uuid="aaaa29a8-35af-508c-8bc3-b662a17a0fe5",
    url="https://github.com/JuliaStats/Clustering.jl",
    julia=true,
    license="MIT",
    is_wrapper=false
    )

metadata_model(KMeans,
    input=MLJBase.Table(MLJBase.Continuous),
    output=MLJBase.Table(MLJBase.Continuous),
    weights=false,
    descr=KMeansDescription
    )

metadata_model(KMedoids,
    input=MLJBase.Table(MLJBase.Continuous),
    output=MLJBase.Table(MLJBase.Continuous),
    weights=false,
    descr=KMedoidsDescription
    )

end # module


text:
Base.promote_rule(::Type{Decimal}, ::Type{<:Real}) = Decimal

# override definitions in Base
Base.promote_rule(::Type{BigFloat}, ::Type{Decimal}) = Decimal
Base.promote_rule(::Type{BigInt}, ::Type{Decimal}) = Decimal

# Addition
# To add, convert both decimals to the same exponent.
# (If the exponents are different, use the smaller exponent
# to make sure we're adding integers.)
function +(x::Decimal, y::Decimal)
    cx = (-1)^x.s * x.c * 10^max(x.q - y.q, 0)
    cy = (-1)^y.s * y.c * 10^max(y.q - x.q, 0)
    s = (abs(cx) > abs(cy)) ? x.s : y.s
    c = BigInt(cx) + BigInt(cy)
    normalize(Decimal(s, abs(c), min(x.q, y.q)))
end
@deprecate (+)(x::Number, y::Array{Decimal}) x .+ y
@deprecate (+)(x::Array{Decimal}, y::Number) x .+ y
@deprecate (+)(x::Array{<:Number}, y::Array{Decimal}) x .+ y
@deprecate (+)(x::Array{Decimal}, y::Array{<:Number}) x .+ y

# Negation
-(x::Decimal) = Decimal((x.s == 1) ? 0 : 1, x.c, x.q)

# Subtraction
-(x::Decimal, y::Decimal) = +(x, -y)
@deprecate (-)(x::Number, y::Array{Decimal}) x .- y
@deprecate (-)(x::Array{Decimal}, y::Number) x .- y
@deprecate (-)(x::Array{<:Number}, y::Array{Decimal}) x .- y
@deprecate (-)(x::Array{Decimal}, y::Array{<:Number}) x .- y

# Multiplication
function *(x::Decimal, y::Decimal)
    s = (x.s == y.s) ? 0 : 1
    normalize(Decimal(s, BigInt(x.c) * BigInt(y.c), x.q + y.q))
end

# Inversion
function Base.inv(x::Decimal)
    str = string(x)
    if str[1] == '-'
        str = str[2:end]
    end
    b = ('.' in str) ? length(split(str, '.')[1]) : 0
    c = round(BigInt(10)^(-x.q + DIGITS) / x.c)
    q = (x.q < 0) ? 1 - b - DIGITS : -b - DIGITS
    normalize(Decimal(x.s, c, q))
end
@deprecate inv(x::Array{Decimal}) map(inv, x)

# Division
/(x::Decimal, y::Decimal) = x * inv(y)

# TODO exponentiation


text:
module RPRMakie

using Makie
using RadeonProRender
using GeometryBasics
using Colors
using FileIO

const RPR = RadeonProRender
using Makie: colorbuffer

const NUM_ITERATIONS = Ref(200)
const RENDER_RESOURCE = Ref(RPR.RPR_CREATION_FLAGS_ENABLE_GPU0)
const RENDER_PLUGIN = Ref(RPR.Tahoe)

include("scene.jl")
include("lines.jl")
include("meshes.jl")
include("volume.jl")

"""
    RPRMakie.activate!(;
        iterations=200,
        resource=RPR.RPR_CREATION_FLAGS_ENABLE_GPU0,
        plugin=RPR.Tahoe)

- iterations: Iterations of light simulations. The more iterations, the less noisy the picture becomes, but higher numbers take much longer. For e.g. `iterations=10` one should expect rendering to take a couple of seconds but it will produce pretty noisy output. 200 is a good middle ground taking around ~30s on an old nvidia 1060. For highest quality output, numbers above 500 are required.
- resource: GPU or CPU to use. Multiple GPUs and CPUs can be used together by using `&` (e.g. `RPR.RPR_CREATION_FLAGS_ENABLE_GPU0 & RPR.RPR_CREATION_FLAGS_ENABLE_CPU`)
- plugin:
    * `RPR.Tahoe`, the legacy RadeonProRender backend. It's the most stable, but doesn't have all new features (e.g. the `RPR.MatX` material), and may be slower than others
    * `RPR.Northstar`, the new rewritten backend, faster and optimized for many iterations. Single iterations are much slower, so less usable for interactive display. Sometimes, Northstar just produces black, jiggly objects. It's not clear yet, if that's just a bug, or the result of using an unsupported/deprecated feature. Switch to `Tahoe` if that happens.
    * `RPR.Hybrid`: Vulkan backend, fit for real time rendering, using AMDs and NVIDIAs new hardware accelerated ray tracing. Doesn't work reliably yet and only works with `RPR.Uber` material.
    * `RPR.HybridPro`: The same as Hybrid, but works only for Radeon GPUs, using AMDs own hardware acceleration API.

"""
function activate!(; iterations=200, resource=RENDER_RESOURCE[], plugin=RENDER_PLUGIN[])
    NUM_ITERATIONS[] = iterations
    RENDER_RESOURCE[] = resource
    RENDER_PLUGIN[] = plugin
    b = RPRBackend()
    Makie.register_backend!(b)
    Makie.current_backend[] = b
    Makie.inline!(true)
    return
end

function __init__()
    activate!()
    return
end

for name in names(Makie; all=true)
    if Base.isexported(Makie, name)
        @eval using Makie: $(name)
        @eval export $(name)
    end
end

export RPRScreen, RPR, colorbuffer

end


text:
using BGEN
using Random
using MLJBase
using MLJLinearModels
using TMLE
using JLD2

genotypesfile = joinpath("data", "genotypes.csv")
confoundersfile = joinpath("data", "confounders.csv")
queryfile = joinpath("config", "query.toml")
continuous_phenotypefile = joinpath("data", "continuous_phenotypes.csv")
binary_phenotypefile = joinpath("data", "binary_phenotypes.csv")
phenotypelist_file = joinpath("data", "phen_list.csv")

function build_query_file(;path=queryfile, threshold=0.9, het="AG")
    queries = Dict(
        "threshold" => threshold,
        "SNPS" => Dict(
            "RSID_10"  => BGEN.datadir("example.8bits.bgen"),
            "RSID_100" => BGEN.datadir("example.8bits.bgen")),
        "QUERY_1" => Dict(
            "RSID_10"  => "GG -> "*het,
            "RSID_100" => "GG -> "*het),
        "QUERY_2" => Dict(
            "RSID_10"  => "GG -> "*het,
            "RSID_100" => "GG -> AA"),
        "QUERY_3" => Dict(
            "RSID_10"  => "TT -> "*het,
            "RSID_100" => "GG -> AA")
    )
    open(path, "w") do io
        TOML.print(io, queries)
    end
end

function build_ate_query_file(;path=queryfile, threshold=0.9)
    queries = Dict(
        "threshold" => threshold,
        "SNPS" => Dict("RSID_10"  => BGEN.datadir("example.8bits.bgen")),
        "QUERY_1" => Dict("RSID_10"  => "GG -> AG"),
        "QUERY_2" => Dict("RSID_10"  => "GG -> AA")
    )
    open(path, "w") do io
        TOML.print(io, queries)
    end
end

function clean_hdf5estimates_files(prefix)
    rm(string(prefix, "_batch_1_Real.hdf5"))
    rm(string(prefix, "_batch_1_Bool.hdf5"))
    rm(string(prefix, "_batch_2_Bool.hdf5"))
end

function build_results_files(grm_ids, prefix)
    rng = Xoshiro(0)
    n = size(grm_ids, 1)
    T = (t₁=categorical(rand(rng, [0, 1], n)),)
    W = (w₁=rand(rng, n), w₂=rand(rng, n))
    height = rand(rng, n)
    bmi = 2convert(Vector{Float64}, T.t₁) + 0.2rand(rng, n)
    cancer = categorical(rand(rng, [0,1], n))
    query_1 = Query(case=(t₁=1,), control=(t₁=0,), name="QUERY_1")
    query_2 = Query(case=(t₁=0,), control=(t₁=1,), name="QUERY_2")

    # Continuous single batch
    cont_path = string(prefix, "_batch_1_Real.hdf5")
    tmle_reg = TMLEstimator(
            LinearRegressor(),
            LogisticClassifier(),
            query_1,
            query_2
    )
    TMLE.fit(tmle_reg, T, W, (height=height, bmi=bmi);
        verbosity=0,
        cache=false,
        callbacks=[TMLEEpistasis.JLD2Saver(cont_path, true)])
    jldopen(cont_path, "a") do io
            io["SAMPLE_IDS"] = Dict(
                "height" => string.(grm_ids.SAMPLE_ID),
                "bmi" => string.(grm_ids.SAMPLE_ID)
            )
    end
    # Binary single batch
    bin_path = string(prefix, "_batch_1_Bool.hdf5")
    tmle_bin = TMLEstimator(
        LogisticClassifier(),
        LogisticClassifier(),
        query_1,
        query_2)
    TMLE.fit(tmle_bin, T, W, (cancer=cancer,);
        verbosity=0,
        cache=false,
        callbacks=[TMLEEpistasis.JLD2Saver(bin_path, true)])
    jldopen(bin_path, "a") do io
        io["SAMPLE_IDS"] = Dict(
            "cancer" => string.(grm_ids.SAMPLE_ID)
        )
    end

    # Empty file that may be output by the tmle procedure
    jldopen(string(prefix, "_batch_2_Bool.hdf5"), "a") do io
        JLD2.Group(io, "TMLEREPORTS")
    end
end

function test_queries(queries, expected_queries)
    for (i, query) in enumerate(queries)
        @test query.case == expected_queries[i].case
        @test query.control == expected_queries[i].control
        @test query.name == expected_queries[i].name
    end
end

function _test_wiped_data(mach)
    @test mach.data == ()
    @test :data ∉ mach.cache
    @test mach.args == ()
end

function test_data_has_been_removed(tmle_mach)
    _test_wiped_data(tmle_mach)

    for submach in tmle_mach.report.machines
        _test_wiped_data(submach)
    end

    for submach in report(tmle_mach).G.machines
        _test_wiped_data(submach)
    end

    for submach in report(tmle_mach).Q̅.machines
        _test_wiped_data(submach)
    end
end


text:
function n_bar_equation(cur_reactor::AbstractReactor, cur_value::Number)
  cur_G_T = cur_value

  cur_G_T /= K_n(cur_reactor)

  cur_equation = Equation(-2, 0, 1, cur_G_T)

  cur_equation
end


text:
# Standalone https://github.com/JuliaLang/julia/pull/18777 for use
# with earlier versions of Julia

import Base: first, last, start, next, done, step, convert, promote_rule,
             show, size, isempty, length, minimum, maximum,
             ctranspose, transpose, copy, getindex, intersect, findin, vcat,
             reverse, issorted, sort, sort!, sortperm, sum, mean, median,
             in, map, float, big

@static if VERSION < v"0.6.0-dev.1614"
    import Base: ==, -, +, .+, .-, .*, ./
else
    import Base: ==, -, +, *, /, broadcast
end

@static if VERSION < v"0.6.0-dev.2618"
    const typclause = ""
    const whereclause = ""
else
    const typclause = "{T}"
    const whereclause = "where T"
end

include_string(
"""
immutable LinSpace{T} <: Range{T}
    start::T
    stop::T
    len::Int
    lendiv::Int

    function LinSpace$(typclause)(start,stop,len) $(whereclause)
        len >= 0 || error("linspace(\$start, \$stop, \$len): negative length")
        if len == 1
            start == stop || error("linspace(\$start, \$stop, \$len): endpoints differ")
            return new(start, stop, 1, 1)
        end
        new(start,stop,len,max(len-1,1))
    end
end
""")

function LinSpace(start, stop, len::Integer)
    T = typeof((stop-start)/len)
    LinSpace{T}(start, stop, len)
end

linspace(start, stop, len::Real=50) = LinSpace(start, stop, Int(len))

# range(start, step, len) = linspace(start, fma(len-1, step, start), len)
range(start, step, len) = linspace(start, start+(len-1)*step, len)

function show(io::IO, r::LinSpace)
    print(io, "Ranges.linspace(")
    show(io, first(r))
    print(io, ',')
    show(io, last(r))
    print(io, ',')
    show(io, length(r))
    print(io, ')')
end

logspace(start::Real, stop::Real, n::Integer=50) = 10.^linspace(start, stop, n)

## interface implementations

isempty(r::LinSpace) = length(r) == 0

step(r::LinSpace) = (last(r)-first(r))/r.lendiv

length(r::LinSpace) = r.len

first(r::LinSpace) = r.start

last(r::LinSpace) = r.stop

start(r::LinSpace) = 1
done(r::LinSpace, i::Int) = length(r) < i
@inline function next(r::LinSpace, i::Int)
    unsafe_getindex(r, i), i+1
end

@inline function getindex(r::LinSpace, i::Integer)
    @boundscheck checkbounds(r, i)
    unsafe_getindex(r, i)
end

# This is separate to make it useful even when running with --check-bounds=yes
function unsafe_getindex(r::LinSpace, i::Integer)
    d = r.lendiv
    j, a, b = ifelse(2i >= length(r), (i-1, r.start, r.stop), (length(r)-i, r.stop, r.start))
    lerpi(j, d, a, b)
end

@inline function lerpi{T}(j::Integer, d::Integer, a::T, b::T)
    t = j/d
    # computes (1-t)*a + t*b
    # T(fma(t, b, fma(-t, a, a)))
    # Until issues of precision, etc, are worked out in
    # https://github.com/JuliaLang/julia/pull/18777, for safety use
    # the simplest approach
    convert(T, (1-t)*a + t*b)
end

@inline function getindex{T}(r::LinSpace{T}, s::OrdinalRange)
    @boundscheck checkbounds(r, s)
    sl = length(s)
    ifirst = first(s)
    ilast = last(s)
    vfirst = unsafe_getindex(r, ifirst)
    vlast  = unsafe_getindex(r, ilast)
    return linspace(vfirst, vlast, sl)
end

=={T<:LinSpace}(r::T, s::T) = (first(r) == first(s)) & (length(r) == length(s)) & (last(r) == last(s))

-(r::LinSpace) = LinSpace(-r.start, -r.stop, r.len)

@static if VERSION < v"0.6.0-dev.1614"
    function .+(x::Real, r::LinSpace)
        LinSpace(x + r.start, x + r.stop, r.len)
    end
    function .+(x::Number, r::LinSpace)
        LinSpace(x + r.start, x + r.stop, r.len)
    end
    function .-(x::Real, r::LinSpace)
        LinSpace(x - r.start, x - r.stop, r.len)
    end
    function .-(x::Number, r::LinSpace)
        LinSpace(x - r.start, x - r.stop, r.len)
    end
    function .-(r::LinSpace, x::Real)
        LinSpace(r.start - x, r.stop - x, r.len)
    end
    function .-(r::LinSpace, x::Number)
        LinSpace(r.start - x, r.stop - x, r.len)
    end

    .*(x::Real, r::LinSpace)     = LinSpace(x * r.start, x * r.stop, r.len)
    .*(r::LinSpace, x::Real)     = x .* r

    ./(r::LinSpace, x::Real)     = LinSpace(r.start / x, r.stop / x, r.len)
else
    broadcast(::typeof(+), x::Real, r::LinSpace) =   LinSpace(x + r.start, x + r.stop, r.len)
    broadcast(::typeof(+), x::Number, r::LinSpace) = LinSpace(x + r.start, x + r.stop, r.len)
    broadcast(::typeof(-), x::Real, r::LinSpace) =   LinSpace(x - r.start, x - r.stop, r.len)
    broadcast(::typeof(-), x::Number, r::LinSpace) = LinSpace(x - r.start, x - r.stop, r.len)
    broadcast(::typeof(-), r::LinSpace, x::Real) =   LinSpace(r.start - x, r.stop - x, r.len)
    broadcast(::typeof(-), r::LinSpace, x::Number) = LinSpace(r.start - x, r.stop - x, r.len)
    broadcast(::typeof(*), x::Real, r::LinSpace) =   LinSpace(x * r.start, x * r.stop, r.len)
    broadcast(::typeof(*), r::LinSpace, x::Real) =   x .* r
    broadcast(::typeof(/), r::LinSpace, x::Real) =   LinSpace(r.start / x, r.stop / x, r.len)
end

promote_rule{T1,T2}(::Type{LinSpace{T1}},::Type{LinSpace{T2}}) =
    LinSpace{promote_type(T1,T2)}
convert{T<:AbstractFloat}(::Type{LinSpace{T}}, r::LinSpace{T}) = r
convert{T<:AbstractFloat}(::Type{LinSpace{T}}, r::LinSpace) =
    LinSpace{T}(r.start, r.stop, r.len)

promote_rule{F,OR<:OrdinalRange}(::Type{LinSpace{F}}, ::Type{OR}) =
    LinSpace{promote_type(F,eltype(OR))}
convert{T<:AbstractFloat}(::Type{LinSpace{T}}, r::OrdinalRange) =
    linspace(convert(T, first(r)), convert(T, last(r)), length(r))
convert{T}(::Type{LinSpace}, r::OrdinalRange{T}) =
    convert(LinSpace{typeof(float(first(r)))}, r)

# Promote FloatRange to LinSpace
@static if VERSION < v"0.6.0-dev.2376"
    promote_rule{F,OR<:FloatRange}(::Type{LinSpace{F}}, ::Type{OR}) =
        LinSpace{promote_type(F,eltype(OR))}
    convert{T<:AbstractFloat}(::Type{LinSpace{T}}, r::FloatRange) =
        linspace(convert(T, first(r)), convert(T, last(r)), length(r))
    convert{T<:AbstractFloat}(::Type{LinSpace}, r::FloatRange{T}) =
        convert(LinSpace{T}, r)
else
    include_string(
"""
    promote_rule{F,OR<:StepRangeLen{<:AbstractFloat}}(::Type{LinSpace{F}}, ::Type{OR}) =
        LinSpace{promote_type(F,eltype(OR))}
    convert(::Type{LinSpace{T}}, r::StepRangeLen{<:AbstractFloat}) where {T<:AbstractFloat} =
        linspace(convert(T, first(r)), convert(T, last(r)), length(r))
    convert(::Type{LinSpace}, r::StepRangeLen{T}) where {T<:AbstractFloat} =
        convert(LinSpace{T}, r)
""")
end

reverse(r::LinSpace)     = LinSpace(r.stop, r.start, length(r))

# Extra methods
function map{T<:AbstractFloat}(::Type{T}, r::LinSpace)
    LinSpace(T(r.start), T(r.stop), length(r))
end

float(r::LinSpace) = LinSpace(float(r.start), float(r.stop), length(r))
big(r::LinSpace) = LinSpace(big(r.start), big(r.stop), length(r))

function lerpi(j::Integer, d::Integer, A::AbstractArray, B::AbstractArray)
    broadcast((a,b) -> lerpi(j, d, a, b), A, B)
end
function lerpi(j::Integer, d::Integer, a::BigFloat, b::BigFloat)
    t = BigFloat(j)/d
    fma(t, b, fma(-t, a, a))
end
function lerpi(j::Integer, d::Integer, a::Rational, b::Rational)
    ((d-j)*a)/d + (j*b)/d
end

function show(io::IO, ::MIME"text/plain", r::LinSpace)
    # show for linspace, e.g.
    # linspace(1,3,7)
    # 7-element LinSpace{Float64}:
    #   1.0,1.33333,1.66667,2.0,2.33333,2.66667,3.0
    print(io, summary(r))
    if !isempty(r)
        println(io, ":")
        Base.print_range(io, r)
    end
end


text:
function export_wav(data)

end


text:
using FLOWMath
using Test
import ForwardDiff
import FiniteDiff
using LinearAlgebra: diag

@testset "FLOWMath.jl" begin

# ------ trapz --------
# tests from matlab trapz docs
x = 1:5
y = [1 4 9 16 25]
z = trapz(x, y)
@test z == 42.0

x = range(0.0, stop=pi+1e-15, step=pi/100)
y = sin.(x)
z = trapz(x, y)
@test isapprox(z, 1.999835503887444, atol=1e-15)

# -------------------------

# ------ Brent's method ------

f(x) = x^2 - 1.0

xstar, _ = brent(f, -2.0, 0)
@test xstar == -1.0
xstar, _ = brent(f, 0.0, 2)
@test xstar == 1.0

f(x) = x^3 - 1
xstar, _ = brent(f, 0, 3)
@test xstar == 1.0

f(x) = sin(x)
atol = 2e-12
xstar, _ = brent(f, 1, 4, atol=atol)
@test isapprox(xstar, pi, atol=atol)

atol = 1e-15
xstar, _ = brent(f, 1, 4, atol=atol)
@test isapprox(xstar, pi, atol=atol)

# -------------------------

# ------ abs_smooth ---------
x = -3.0
delta_x = 0.0
y = abs_smooth(x, delta_x)
@test y == 3.0

x = 3.0
delta_x = 0.1
y = abs_smooth(x, delta_x)
@test y == 3.0

x = -3.0
delta_x = 0.1
y = abs_smooth(x, delta_x)
@test y == 3.0

x = 0
delta_x = 0.1
y = abs_smooth(x, delta_x)
@test y == 0.05

x = 0.1
delta_x = 0.1
y = abs_smooth(x, delta_x)
@test y == 0.1

x = -0.1
delta_x = 0.1
y = abs_smooth(x, delta_x)
@test y == 0.1

x = 0.05
delta_x = 0.1
y = abs_smooth(x, delta_x)
@test y == 0.0625

x = -0.05
delta_x = 0.1
y = abs_smooth(x, delta_x)
@test y == 0.0625
# -------------------------

# ------ ksmax ---------
x = [0.0, 0.0]
x_max_smooth = ksmax(x)
@test isapprox(x_max_smooth, 0.013862943611198907)

# overflow
x = [1e6, 1e6]
x_max_smooth = ksmax(x)
@test isapprox(x_max_smooth, 1.0000000138629436e6)

# underflow
x = [-1e6, -1e6]
x_max_smooth = ksmax(x)
@test isapprox(x_max_smooth, -999999.9861370564)

# hardness
x = [0.0, 0.0]
hardness = 100.0
x_max_smooth = ksmax(x, hardness)
@test isapprox(x_max_smooth, 0.006931471805599453)

# ------ ksmin ---------
x = [0.0, 0.0]
x_max_smooth = ksmin(x)
@test isapprox(x_max_smooth, -0.013862943611198907)

# overflow
x = [1e6, 1e6]
x_max_smooth = ksmin(x)
@test isapprox(x_max_smooth, 999999.9861370564)

# underflow
x = [-1e6, -1e6]
x_max_smooth = ksmin(x)
@test isapprox(x_max_smooth, -1.0000000138629436e6)

# hardness
x = [0.0, 0.0]
hardness = 100.0
x_max_smooth = ksmin(x, hardness)
@test isapprox(x_max_smooth, -0.006931471805599453)

# -------------------------

# ------ ksmax_adaptive ---------
x = [0.0, 0.0]
x_max_smooth = ksmax_adaptive(x)
@test isapprox(x_max_smooth, 0.0008325546111576975)

# overflow
x = [1e6, 1e6]
x_max_smooth = ksmax_adaptive(x)
@test isapprox(x_max_smooth, 1.0000000008325547e6)

# underflow
x = [-1e6, -1e6]
x_max_smooth = ksmax_adaptive(x)
@test isapprox(x_max_smooth, -999999.9991674453)

# hardness
x = [0.0, 0.0]
hardness = 100.0
x_max_smooth = ksmax_adaptive(x, hardness)
@test isapprox(x_max_smooth, 0.0008325546111576975)

x = [-0.1, 0.0]
hardness = 100.0
x_max_smooth = ksmax_adaptive(x, hardness)
@test isapprox(x_max_smooth, 4.5398899216870535e-7)

# tolerance
x = [0.0, 0.0]
tol = 1e-3
x_max_smooth = ksmax_adaptive(x, tol=tol)
@test isapprox(x_max_smooth, 0.013862943611198907)

# smoothing_fraction
x = [-0.165, 0.0]
smoothing_fraction = 0.1
x_max_smooth = ksmax_adaptive(x, smoothing_fraction=smoothing_fraction)
@test isapprox(x_max_smooth, 5.261253727654178e-6)

x = [-0.165, 0.0]
smoothing_fraction = 0.2
x_max_smooth = ksmax_adaptive(x, smoothing_fraction=smoothing_fraction)
@test isapprox(x_max_smooth, 5.2856933329025475e-6)

# ------ ksmin_adaptive ---------
x = [0.0, 0.0]
x_max_smooth = ksmin_adaptive(x)
@test isapprox(x_max_smooth, -0.0008325546111576975)

# overflow
x = [1e6, 1e6]
x_max_smooth = ksmin_adaptive(x)
@test isapprox(x_max_smooth, 999999.9991674453)

# underflow
x = [-1e6, -1e6]
x_max_smooth = ksmin_adaptive(x)
@test isapprox(x_max_smooth, -1.0000000008325547e6)

# hardness
x = [0.0, 0.0]
hardness = 100.0
x_max_smooth = ksmin_adaptive(x, hardness)
@test isapprox(x_max_smooth, -0.0008325546111576975)

x = [0.1, 0.0]
hardness = 100.0
x_max_smooth = ksmin_adaptive(x, hardness)
@test isapprox(x_max_smooth, -4.5398899216870535e-7)

# tolerance
x = [0.0, 0.0]
tol = 1e-3
x_max_smooth = ksmin_adaptive(x, tol=tol)
@test isapprox(x_max_smooth, -0.013862943611198907)

# smoothing_fraction
x = [0.165, 0.0]
smoothing_fraction = 0.1
x_max_smooth = ksmin_adaptive(x, smoothing_fraction=smoothing_fraction)
@test isapprox(x_max_smooth, -5.261253727654178e-6)

x = [0.165, 0.0]
smoothing_fraction = 0.2
x_max_smooth = ksmin_adaptive(x, smoothing_fraction=smoothing_fraction)
@test isapprox(x_max_smooth, -5.2856933329025475e-6)

# -------------------------

# ------ sigmoid ---------
x = -0.5
y = sigmoid(x)
@test isapprox(y, 0.37754066879814546)

x = 0.0
y = sigmoid(x)
@test isapprox(y, 0.5)

x = 0.5
y = sigmoid(x)
@test isapprox(y, 0.6224593312018546)

# overflow
x = 1e6
y = sigmoid(x)
@test isapprox(y, 1.0)

# underflow
x = -1e6
y = sigmoid(x)
@test isapprox(y, 0.0)

# -------------------------

# ------ sigmoid_blend ---------
x = -0.1
f1x = x
f2x = x^2
xt = 0.0
y = sigmoid_blend(f1x, f2x, x, xt)
@test isapprox(y, -0.09926378639832867)

x = 0.0
f1x = x
f2x = x^2
xt = 0.0
y = sigmoid_blend(f1x, f2x, x, xt)
@test isapprox(y, 0.0)

x = 0.1
f1x = x
f2x = x^2
xt = 0.0
y = sigmoid_blend(f1x, f2x, x, xt)
@test isapprox(y, 0.010602356583185632)

# overflow
x = 1e6
f1x = x
f2x = x^2
xt = 0.0
y = sigmoid_blend(f1x, f2x, x, xt)
@test isapprox(y, f2x)

# underflow
x = -1e6
f1x = x
f2x = x^2
xt = 0.0
y = sigmoid_blend(f1x, f2x, x, xt)
@test isapprox(y, f1x)

# hardness
x = -0.1
f1x = x
f2x = x^2
xt = 0.0
hardness = 100
y = sigmoid_blend(f1x, f2x, x, xt, hardness)
@test isapprox(y, -0.09999500623444274)

x = 0.0
f1x = x
f2x = x^2
xt = 0.0
hardness = 100
y = sigmoid_blend(f1x, f2x, x, xt, hardness)
@test isapprox(y, 0.0)

x = 0.1
f1x = x
f2x = x^2
xt = 0.0
hardness = 100
y = sigmoid_blend(f1x, f2x, x, xt, hardness)
@test isapprox(y, 0.010004085808183225)

# vectorized
x = -0.25:0.05:0.25
f1x = x
f2x = x.^2
xt = 0.0
hardness = 100
y = sigmoid_blend.(f1x, f2x, x, xt, hardness)
ytest = [-0.24999999999566003,
         -0.19999999950532316,
         -0.14999994723186585,
         -0.09999500623444274,
         -0.04964862532647505,
          0.0,
          0.0028179104189035298,
          0.010004085808183225,
          0.022500039002533945,
          0.040000000329784596,
          0.062500000002604]
@test isapprox(y, ytest)

# -------------------------

# ------ cubic_blend ---------
x = -3.0
f1x = x
f2x = x^2
xt = 0.0
delta_x = 0.1
y = cubic_blend(f1x, f2x, x, xt, delta_x)
@test y == -3.0

x = 3.0
f1x = x
f2x = x^2
xt = 0.0
delta_x = 0.1
y = cubic_blend(f1x, f2x, x, xt, delta_x)
@test y == 9.0

x = 0.0
f1x = x
f2x = x^2
xt = 0.0
delta_x = 0.1
y = cubic_blend(f1x, f2x, x, xt, delta_x)
@test y == 0.0

x = -0.1
f1x = x
f2x = x^2
xt = 0.0
delta_x = 0.1
y = cubic_blend(f1x, f2x, x, xt, delta_x)
@test y == f1x

x = 0.1
f1x = x
f2x = x^2
xt = 0.0
delta_x = 0.1
y = cubic_blend(f1x, f2x, x, xt, delta_x)
@test y == f2x

x = -0.05
f1x = x
f2x = x^2
xt = 0.0
delta_x = 0.1
y = cubic_blend(f1x, f2x, x, xt, delta_x)
@test isapprox(y, -0.041796875000000004)

x = 0.05
f1x = x
f2x = x^2
xt = 0.0
delta_x = 0.1
y = cubic_blend(f1x, f2x, x, xt, delta_x)
@test isapprox(y, 0.009921875000000004)

# vectorized
x = -0.25:0.05:0.25
f1x = x
f2x = x.^2
xt = 0.0
delta_x = 0.1
y = cubic_blend.(f1x, f2x, x, xt, delta_x)
ytest = [-0.25,
         -0.2,
         -0.15,
         -0.1,
         -0.041796875000000004,
          0.0,
          0.009921875000000004,
          0.010000000000000002,
          0.0225,
          0.04000000000000001,
          0.0625]
@test isapprox(y, ytest)

# -------------------------
# ------ quintic_blend ---------
x = -3.0
f1x = x
f2x = x^2
xt = 0.0
delta_x = 0.1
y = quintic_blend(f1x, f2x, x, xt, delta_x)
@test y == -3.0

x = 3.0
f1x = x
f2x = x^2
xt = 0.0
delta_x = 0.1
y = quintic_blend(f1x, f2x, x, xt, delta_x)
@test y == 9.0

x = 0.0
f1x = x
f2x = x^2
xt = 0.0
delta_x = 0.1
y = quintic_blend(f1x, f2x, x, xt, delta_x)
@test y == 0.0

x = -0.1
f1x = x
f2x = x^2
xt = 0.0
delta_x = 0.1
y = quintic_blend(f1x, f2x, x, xt, delta_x)
@test y == f1x

x = 0.1
f1x = x
f2x = x^2
xt = 0.0
delta_x = 0.1
y = quintic_blend(f1x, f2x, x, xt, delta_x)
@test y == f2x

x = -0.05
f1x = x
f2x = x^2
xt = 0.0
delta_x = 0.1
y = quintic_blend(f1x, f2x, x, xt, delta_x)
@test isapprox(y, -0.044565429687500005)

x = 0.05
f1x = x
f2x = x^2
xt = 0.0
delta_x = 0.1
y = quintic_blend(f1x, f2x, x, xt, delta_x)
@test isapprox(y, 0.007416992187499999)

# vectorized
x = -0.25:0.05:0.25
f1x = x
f2x = x.^2
xt = 0.0
delta_x = 0.1
y = quintic_blend.(f1x, f2x, x, xt, delta_x)
ytest = [ -0.25,
          -0.2,
          -0.15,
          -0.1,
          -0.044565429687500005,
           0.0,
           0.007416992187499999,
           0.010000000000000002,
           0.0225,
           0.04000000000000001,
           0.0625]
@test isapprox(y, ytest)
# -------------------------

# ---------- akima interpolate -------------

x = 0:pi/4:2*pi
y = sin.(x)
xpt = 0:pi/16:2*pi

ypt = akima(x, y, xpt)

# FIXME: this is not an independent test.  right now it is just testing against refactoring.
ytest = [0.0, 0.21394356492851746, 0.40088834764831843, 0.5641656518405971, 0.7071067811865475, 0.8281808021479407, 0.920495128834866, 0.9793385864009951, 1.0, 0.9873810649285174, 0.9419417382415921, 0.8523082377305078, 0.7071067811865476, 0.5303300858899105, 0.35355339059327384, 0.17677669529663714, 1.2246467991473532e-16, -0.17677669529663698, -0.35355339059327373, -0.5303300858899104, -0.7071067811865475, -0.8523082377305078, -0.941941738241592, -0.9873810649285175, -1.0, -0.9793385864009951, -0.9204951288348657, -0.8281808021479408, -0.7071067811865477, -0.5641656518405971, -0.400888347648319, -0.21394356492851788, 0.0]

@test isapprox(ypt, ytest, atol=1e-8)

spl = Akima(x, y)

dydx = derivative(spl, pi/16)

@test isapprox(dydx, 1.0180260961104746, atol=1e-12)

dydx = gradient(spl, xpt)

dydx_test = [1.1640128599466306, 1.0180260961104746, 0.889005522603991, 0.7769511394271804, 0.6818629465800423, 0.5473879346340718, 0.3889191062220746, 0.20645646134405082, 1.8419755006770065e-16, -0.13818978335121798, -0.33430576382780836, -0.5883479414297708, -0.900316316157106, -0.9003163161571062, -0.9003163161571062, -0.9003163161571061, -0.9003163161571061, -0.9003163161571064, -0.9003163161571064, -0.9003163161571065, -0.9003163161571061, -0.5883479414297709, -0.33430576382780913, -0.13818978335121862, -1.8419755006770062e-16, 0.20645646134405057, 0.3889191062220747, 0.5473879346340714, 0.6818629465800421, 0.7769511394271802, 0.8890055226039906, 1.0180260961104746, 1.1640128599466313]

@test isapprox(dydx, dydx_test, atol=1e-12)

# differentiate with dual numbers

dydx = ForwardDiff.derivative(spl, pi/16)
@test isapprox(dydx, 1.0180260961104746, atol=1e-12)

J = ForwardDiff.jacobian(spl, xpt)
dydx = diag(J)
@test isapprox(dydx, dydx_test, atol=1e-12)

# differentiate coordinates of spline

wrapper(y) = Akima(x, y, 0.1)(xpt)
J = ForwardDiff.jacobian(wrapper, y)
J2 = FiniteDiff.finite_difference_jacobian(wrapper, y)
@test maximum(abs.(J - J2)) < 1e-6

wrapper2(y) = Akima(x, y, 0.1)(xpt)
J = ForwardDiff.jacobian(wrapper2, x)
J2 = FiniteDiff.finite_difference_jacobian(wrapper2, x)
@test maximum(abs.(J - J2)) < 1e-6

# ---------------------------

# ----- linear interpolation ------
xvec = [1.0, 2.0, 4.0, 5.0]
yvec = [2.0, 3.0, 5.0, 8.0]

y = linear(xvec, yvec, 1.0)
@test y == 2.0
y = linear(xvec, yvec, 1.5)
@test y == 2.5
y = linear(xvec, yvec, 3.0)
@test y == 4.0
y = linear(xvec, yvec, 4.5)
@test y == 6.5
y = linear(xvec, yvec, 5.0)
@test y == 8.0

y = linear(xvec, yvec, [1.0, 1.5, 3.0, 4.5, 5.0])
@test y == [2.0, 2.5, 4.0, 6.5, 8.0]

dydx = derivative(xvec, yvec, 1.0)
@test dydx == 1.0
dydx = derivative(xvec, yvec, 1.5)
@test dydx == 1.0
dydx = derivative(xvec, yvec, 3.0)
@test dydx == 1.0
dydx = derivative(xvec, yvec, 4.5)
@test dydx == 3.0
dydx = derivative(xvec, yvec, 5.0)
@test dydx == 3.0

dydx = gradient(xvec, yvec, [1.0, 1.5, 3.0, 4.5, 5.0])
@test dydx == [1.0, 1.0, 1.0, 3.0, 3.0]
# ---------------------------

end


text:
const DEBUG = true

"""
    @dbgassert(expr, [message]) -> true

If `DEBUG = true`, throws an error with an optionally specified `message` if the
given code expression evaluates to `false`, otherwise returns `true`.
"""
macro dbgassert(expr, msgs...)
    if DEBUG
        msg_str = isempty(msgs) ? string(expr) : msgs[1]
        esc(:($expr || throw(AssertionError($msg_str))))
    else
        true
    end
end

"""
    @dbgasserteq(left, right, [message]) -> true

If `DEBUG = true`, throws an error with an optionally specified `message` if the
given `left` code expression does not evaluates to the same value as the given
`right` expression, otherwise returns `true`.
"""
macro dbgasserteq(left, right, msgs...)
    if DEBUG
        msg_str = isempty(msgs) ? "$(left) == $(right)" : msgs[1]
        esc(:($left == $right || throw(AssertionError($msg_str))))
    else
        true
    end
end

max_stars_repo_name: lars-frogner/Whirl

max_issues_repo_name: lars-frogner/Whirl

max_forks_repo_name: lars-frogner/Whirl


text:
module Datax
using  Requires, Formatting

export datax,@datax

# TODO: use printf

"""
```julia
datax(...)
```
Print the arguments to a file readable by pgfkeys (Best for use with the `datax` LaTeX
package). A string will be printed as is, a number will be wrapped in `siunitx`'s
`\\num`, and a tuple `(value::Number, unit::String)` or a `Unitful` quantity will be wrapped in `\\SI`.
If the argument is a `Tuple` and the first argument is a string, it is used as a format
string for the rest of the tuple.
There's also a macro form, `@datax`, which reuses the variable names from the script.

# Examples
```julia
datax(a=2,b=1.24,c="hi",d=(24,"\\meter"),e=15u"kg/s^2")
```
Save the given variables in siunitx form
"""
function datax(;kwargs...)
    datax(Dict(kwargs))
end

function datax(d::Dict{Symbol,<:Any};filename::String="data.tex",format::String="%.4g")
    open(filename,"w") do f
        println(f,"% Autogenerated by Datax.jl, will be overwritten")
        for (k,v) in d
            print(f,"\\pgfkeyssetvalue{/datax/$k}{")
            printdata(f,v,format)
            print(f,"}\n")
        end
    end
end

"""
```julia
@datax ...
```
Print the arguments to a file readable by pgfkeys (Best for use with the `datax` LaTeX package)
Like the `datax()` function, but use the variables' names.

# Examples
```julia
a=2;
b=3.2u"m"
c="hi"
@datax a b c
```
"""
macro datax(args...)
    esc(datax_helper(args...))
end
function datax_helper(args...)
    names = Expr(:vect, QuoteNode.(args)...)
    values = Expr(:vect, args...)
    quote
        datax(Dict($names .=> $values))
    end
end

printdata(f::IO,v::String,fmt::String) = print(f,v)
printdata(f::IO,v::Number,fmt::String) = print(f,"\\num{"*sprintf1(fmt,v)*"}")
printdata(f::IO,v::AbstractArray{<:Number},fmt::String) = print(f,"\\numlist{"*prod(sprintf1.(fmt,v).*";")*"}")
printdata(f::IO,v::Tuple{<:Number,String},fmt::String) = print(f,"\\SI{"*sprintf1(fmt,v[1])*"}{"*v[2]*"}")
printdata(f::IO,v::Tuple{AbstractArray{<:Number},String},fmt::String) = print(f,"\\SIlist{"*prod(sprintf1.(fmt,v[1]).*";")*"}{"*v[2]*"}")
printdata(f::IO,v::Tuple{String,Vararg},fmt::String) = printdata(f,length(v)>2 ? v[2:end] : v[2], v[1])

function __init__()
    @require Unitful="1986cc42-f94f-5a68-af5c-568840ba703d" begin

        function printdata(f::IO,v::Unitful.Quantity,fmt::String)
            print(f,"\\SI{"*sprintf1(fmt,v.val)*"}{")
            siunitxprint(f,Unitful.unit(v))
            print(f,"}")
        end
        function printdata(f::IO,v::AbstractArray{T},fmt::String) where T<:Unitful.Quantity
            print(f,"\\SIlist{"*prod(sprintf1.(fmt,Unitful.ustrip(v)).*";")*"}{")
            siunitxprint(f,Unitful.unit(T))
            print(f,"}")
        end

        function siunitxprint(f::IO,u::Unitful.FreeUnits)
            prefixes = Dict(
                            -24 => "\\yocto",
                            -21 => "\\zepto",
                            -18 => "\\atto",
                            -15 => "\\femto",
                            -12 => "\\pico",
                            -9  => "\\nano",
                            -6  => "\\micro",
                            -3  => "\\milli",
                            -2  => "\\centi",
                            -1  => "\\deci",
                            0   => "",
                            1   => "\\deka",
                            2   => "\\hecto",
                            3   => "\\kilo",
                            6   => "\\mega",
                            9   => "\\giga",
                            12  => "\\tera",
                            15  => "\\peta",
                            18  => "\\exa",
                            21  => "\\zetta",
                            24  => "\\yotta"
                           )
            for p in typeof(u).parameters[1]
                p.power<0 && print(f,"\\per")
                iszero(p.tens) || print(f,prefixes[p.tens])
                print(f,"\\$(lowercase(String(Unitful.name(p))))")
                abs(p.power)==1 || print(f,"\\tothe$(Int64(abs(p.power)))")
            end
        end

    end
end

                                                           end # module


text:
using Ipopt
using Statistics
using StatsBase

# Smoothing function for tampering
function smoothstep(x, N)
    if x >= 1
        return(1)
    end

    if x <= 0
        return(0)
    end

    f = 0

    for i in 0:N
        f = f .+ binomial(N + i, i) * binomial(2 * N + 1, N - i) * (- x) ^ i
    end

    f = x ^ (N+1) * f

    if 0 < x < 1
        return(f)
    end
end

# Tampering of the properties
function tamp!(Q,n)
    nz = size(Q,1)
    nx = size(Q,2)

    for i in 1:nz
        for j in 1:nx
            lat = max( n - i, i - (nx - n) )
            vert = max( n - j, j - (nz - n) )

            dist = max(lat, vert)

            Q[i,j] = Q[i,j] * smoothstep(1 - (dist+1) / n, 3)
        end
    end
end

# Generic objective function evaluation
function function_eval(x, paE, x0, c_y, Cinv; layered=false, log_normal_param=false)
    nz = size(paE.mgrid[1],1)
    nx = size(paE.mgrid[2],1)

#     If layered model, then build 2D model
    if layered
        X = zeros(nz,nx)

        for i in 1:nz
            X[i,:] .= x[i]
        end
    else
        X = copy(x)
    end

#     If log normal parametrization, get physical model
    if log_normal_param
        X = exp.(X)
    end

#     If non log normal parametrization, nothing more to do
    if !log_normal_param

    end

    data_loss = 1 / c_y * func(X,paE)

    if layered
        prior_loss = (x .- x0)' * Cinv * (x .- x0)
    else
        prior_loss = (reshape(x,nz*nx) .- reshape(x0,nz*nx))' * Cinv * (reshape(x,nz*nx) .- reshape(x0,nz*nx))
    end

    return (data_loss + prior_loss)
end

# Generic gradient eval
function gradient_eval(x, grad_f, paE, x0, c_y, Cinv; layered=true, log_normal_param=false)
    nz = size(paE.mgrid[1],1)
    nx = size(paE.mgrid[2],1)

#     If layered model, then build 2D model
    if layered
        X = zeros(nz,nx)

        for i in 1:nz
            X[i,:] .= x[i]
        end
    else
        X = copy(x)
    end

#     If log normal parametrization, get physical model
    if log_normal_param
        X = exp.(X)
    end

    mod!(paE,reshape(X,nz*nx),FGσ())

    grad_tot = 1 / c_y * reshape(paE.g, (nz * nx,1))[:,1]

    if log_normal_param
        grad_tot = grad_tot .* reshape(X, nz * nx)
    end

    if layered
        proj = zeros(nz, nz * nx)

        for i in 1:nz
            proj_l = zeros(nz, nx)
            proj_l[i,:] .= 1.

            proj[i,:] .= reshape(proj_l, nz * nx)
        end

        grad_tot = proj * grad_tot .+ 2 * Cinv * (x .- x0)
    else
        grad_tot = grad_tot .+ 2 * Cinv * (reshape(x,nz*nx) .- reshape(x0,nz*nx))
    end

    grad_f .= grad_tot
end

# Optimization routine
function solve_inverse_problem(paE, x0, c_y, Cinv; layered=true, log_normal_param=false)
    nz = size(paE.mgrid[1],1)
    nx = size(paE.mgrid[2],1)

    if layered
        x = copy(paE.σ[:,1])
    else
        x = copy(paE.σ)
    end

    f(x) = function_eval(x, paE, x0, c_y, Cinv; layered=layered, log_normal_param=log_normal_param)
    gradf(x,grad_f) = gradient_eval(x, grad_f, paE, x0, c_y, Cinv; layered=layered, log_normal_param=log_normal_param)

    function void_g(x, g)

    end

    function void_g_jac(x, mode, rows, cols, values)

    end

    if log_normal_param
        prob = createProblem(nz, -10. .* ones(nz), 10. .* ones(nz), 0,
                    Array{Float64}(undef,0), Array{Float64}(undef,0), 0, 0,
                    f, void_g, gradf, void_g_jac, nothing)
    else
        prob = createProblem(nz, 0. .* ones(nz), 100. .* ones(nz), 0,
                    Array{Float64}(undef,0), Array{Float64}(undef,0), 0, 0,
                    f, void_g, gradf, void_g_jac, nothing)
    end

    prob.x = copy(x0)

    addOption(prob, "hessian_approximation", "limited-memory")
    addOption(prob, "max_iter", 30)
    addOption(prob, "limited_memory_max_history", 30)

    solveProblem(prob)

    return prob.x
end


text:
"""
    OutNeighborVector <: AbstractVector

A structure for iterating over the out neighbors in a graph for a certain vertex.
"""
struct OutNeighborVector{V, G <: LG.AbstractGraph{V}} <: AbstractVector{V}
    graph::G
    vertex::V
end

"""
    SimpleEdgeVector <: AbstractVector

A structure for iterating over the edges of a graph
"""
struct SimpleEdgeVector{V, G <: LG.AbstractGraph{V}} <: AbstractVector{LG.Edge{V}}
    graph::G
end


text:
# Copyright 2017 Eric S. Tellez
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

using Random

export montecarlo, kfolds

function montecarlo(evalfun::Function, X::Vector{ItemType}, y::Vector{LabelType}; runs::Int=3, trainratio=0.7, testratio=0.3) where {ItemType, LabelType}
    Z = zip(X, y) |> collect
    R = Float64[]
    for i in 1:runs
        shuffle!(Z)
        sep = floor(Int, length(Z) * trainratio)
        last = min(length(Z), sep + 1 + floor(Int, length(Z) * testratio))
        train = @view Z[1:sep]
        test = @view Z[1+sep:last]
        s = evalfun([z[1] for z in train], [z[2] for z in train], [z[1] for z in test], [z[2] for z in test])
        push!(R, s)
    end

    return R
end

function kfolds(evalfun::Function, X::Vector{ItemType}, y::Vector{LabelType}; folds::Int=3, shuffle=true) where {ItemType, LabelType}
    Z = zip(X, y) |> collect
    if shuffle
        shuffle!(Z)
    end
    R = Float64[]
    bsize = floor(Int, length(Z) / folds)
    begin
        test = @view Z[1:bsize]
        train = @view Z[bsize+1:end]
        s = evalfun([z[1] for z in train], [z[2] for z in train], [z[1] for z in test], [z[2] for z in test])
        push!(R, s)
    end

    for i in 1:folds-1
        train = vcat(Z[1:bsize*i], Z[bsize*(i+1)+1:end])
        test = @view Z[bsize*i+1:bsize*(i+1)]
        s = evalfun([z[1] for z in train], [z[2] for z in train], [z[1] for z in test], [z[2] for z in test])
        push!(R, s)
    end

    return R
end


text:
"""
This function returns an empty vector of existing projects,
if data for exsting projects does not exist.
"""
function create_project_existing(projectdata::Nothing,
                                simulation_data::AgentSimulationData,
                                sys_UC::Union{Nothing, PSY.System},
                                investor_name::String,
                                investor_dir::String,
                                scenario_names::Vector{String})
    return Project{Existing}[]
end

"""
This function returns a vector of existing projects for the investor.
"""
function create_project_existing(projectdata::DataFrames.DataFrame,
                                simulation_data::AgentSimulationData,
                                sys_UC::Nothing,
                                investor_name::String,
                                investor_dir::String,
                                scenario_names::Vector{String})

    project_existing = Vector{Project{Existing}}(undef, size(projectdata,1))

    for i in 1:size(projectdata,1)
        project_existing[i] = create_project(projectdata[i,:],
                                           simulation_data,
                                           investor_name,
                                           investor_dir,
                                           scenario_names,
                                           false)
        println("Created project $(get_name(project_existing[i])) for $(investor_name)")
    end

    return project_existing
end

"""
This function returns a vector of existing projects for the investor if PSY data exists.
"""
function create_project_existing(projectdata::DataFrames.DataFrame,
                                simulation_data::AgentSimulationData,
                                sys_UC::PSY.System,
                                investor_name::String,
                                investor_dir::String,
                                scenario_names::Vector{String})

    existing_gens = PSY.get_components(PSY.Generator, sys_UC)
    existing_storage = PSY.get_components(PSY.Storage, sys_UC)

    existing_techs = union(existing_gens, existing_storage)

    project_existing = Vector{Project{Existing}}(undef, DataFrames.nrow(projectdata))

    for i in 1:DataFrames.nrow(projectdata)

        device = PSY.get_components_by_name(PSY.Device, sys_UC, projectdata[i, "GEN_UID"])

        project_existing[i] = create_project(projectdata[i, :],
                                           device[1],
                                           simulation_data,
                                           investor_name,
                                           investor_dir,
                                           scenario_names,
                                           false)

        println("Created project $(get_name(project_existing[i])) for $(investor_name)")
    end

    return project_existing
end

"""
This function returns an empty vector of option projects,
if data for option projects does not exist.
"""
function create_project_options(projectdata::Nothing,
                                   simulation_data::AgentSimulationData,
                                   investor_name::String,
                                   investor_dir::String,
                                   scenario_names::Vector{String})
    return Project{Option}[]
end

"""
This function returns a vector of project options for the investor.
"""
function create_project_options(projectdata::DataFrames.DataFrame,
                                   simulation_data::AgentSimulationData,
                                   investor_name::String,
                                   investor_dir::String,
                                   scenario_names::Vector{String})
    project_option = Vector{Project{Option}}(undef, size(projectdata,1))
    for i in 1:size(projectdata,1)
        project_option[i] = create_project(projectdata[i,:],
                                         simulation_data,
                                         investor_name,
                                         investor_dir,
                                         scenario_names,
                                         true)

        println("Created project $(get_name(project_option[i])) for $(investor_name)")
    end
        return project_option
end


text:
#Julia Parallel Processing

addprocs(1)
srand(111)
r = remotecall(rand, 2, 3, 4)
s = @spawnat 2 1 .+ fetch(r)
fetch(s)

max_stars_repo_name: PacktPublishing/Learning-Jupyter

max_issues_repo_name: PacktPublishing/Learning-Jupyter

max_forks_repo_name: PacktPublishing/Learning-Jupyter


text:
"""
Provides a dynamically typed singly-linked list implementation, `List` similar to that
available in most lisp dialects. This module provides basic Julia functionality to
manipulate such lists, exported under both the Julia and Scheme names where they differ
(such as `Base.first` vs. `car`.
"""
module Lists

import FunctionalCollections: append
using Base.Iterators

const car = Base.first

export Cons, List, isnil, ispair, car, cdr, caar, cadr, cddr, nil, append, ++, cons, list,
       islist, Nil

"""
An object representing nothing, or an empty list.

The singleton instance of this type is called `nil`. This type is isomorphic to, and very
similar, to `Nothing`. However, it is often useful to distinguish the `nil` used within many
lisp dialects from a true `Nothing`, because `nil` is iterable and represents an empty list.
"""
struct Nil end
const nil = Nil.instance

"""
An object which is essentially a pair of two objects.

The two objects are referred to, for historical reasons, as `car` and `cdr`.These
abbreviations are not semantically relevant today, so can generally be thought of as the
head object and the tail object. For `Cons` objects which are (proper) lists, the `car`
(head) object will be the first element of the list, and the `cdr` (tail) object will be a
list representing the remaining elements.
"""
struct Cons
    car
    cdr
end
const cons = Cons

car(α::Cons) = α.car
"""
The tail of a `List`.

Note that the tail need not be a list (or even a `List`) if the list is improper.
"""
cdr(α::Cons) = α.cdr
caar(α::Cons) = car(car(α))
cadr(α::Cons) = car(cdr(α))
cddr(α::Cons) = cdr(cdr(α))

Base.map(f, ::Nil, ::Nil...) = nil
Base.map(f, α::Cons, βs::Cons...) =
    Cons(f(car(α), map(car, βs)...), map(f, cdr(α), map(cdr, βs)...))
Base.filter(p, ::Nil) = nil
Base.filter(p, α::Cons) = let β = filter(p, cdr(α))
    p(car(α)) ? Cons(car(α), β) : β
end

Base.:(==)(α::Cons, β::Cons) = car(α) == car(β) && cdr(α) == cdr(β)

"""
An immutable singly linked list, which may be improper.

A list is either `Nil` (empty) or `Cons`, where the first element is an arbitrary object and
the second element is a list. The `List` type, however, includes all instances of `Cons`,
including those for which the second element is not a list. Such `List` instances are called
improper lists, which are allowed because they are used in many lisp dialects.
"""
const List = Union{Cons, Nil}

List() = nil
List(x, xs...) = Cons(x, List(xs...))
const list = List

isnil(::Nil) = true
isnil(::Cons) = false

ispair(α::List) = !isnil(α)

"""
Return `true` if the provided `List` is in fact a list, i.e., it is not improper.

This is the case if the list is a `Nil` (empty list) or a `Cons` with a proper list as its
tail.
"""
islist(::Nil) = true
islist(α::Cons) = islist(cdr(α))

Base.iterate(α::List, β=α) = isnil(β) ? nothing : (car(β), cdr(β))

Base.length(::Nil) = 0
Base.length(α::Cons) = 1 + length(cdr(α))

Base.getindex(α::Nil, b) = throw(BoundsError(α, b))
Base.getindex(α::Cons, b) = b == 1 ? car(α) : cdr(α)[b - 1]

Base.convert(::Type{List}, xs::List) = xs
Base.convert(::Type{List}, xs) = List(xs...)
append(::Nil, β::List) = β
append(α::Cons, β::List) = Cons(car(α), append(cdr(α), β))
append(α::List, β::List, γ::List, γs::List...) = append(append(α, β), γ, γs...)
const (++) = append

include("broadcast.jl")
include("show.jl")

end  # Lists module


text:
facts("solving MIP with mixintprog") do
    solver = SCIPSolver("display/verblevel", 0)

    # integer knapsack problem
    sol = mixintprog(-[5,3,2,7,4],Float64[2 8 4 2 5],'<',10,:Int,0,1,solver)
    @fact sol.status --> :Optimal
    @fact sol.objval --> roughly(-16.0)
    @fact sol.sol --> roughly([1.0, 0.0, 0.0, 1.0, 1.0])
end


text:
include("../src/datastructs_module.jl")

module ttest
using ..datastructs

# still doesnt work correctly
#
# probably the rotations are wrong

rbtree  = datastructs.RedBlackTree(1)
n1      = datastructs.RedBlackTreeNode(2)
n2      = datastructs.RedBlackTreeNode(5)
n3      = datastructs.RedBlackTreeNode(7)
n4      = datastructs.RedBlackTreeNode(8)
n5      = datastructs.RedBlackTreeNode(11)
n6      = datastructs.RedBlackTreeNode(14)
n7      = datastructs.RedBlackTreeNode(15)
n8      = datastructs.RedBlackTreeNode(4)

insert!(rbtree, n6)
insert!(rbtree, n7)

insert!(rbtree, n1)
insert!(rbtree, n2)
insert!(rbtree, n3)
insert!(rbtree, n4)
insert!(rbtree, n5)

insert!(rbtree, n8)
datastructs.tree_view(rbtree)

# datastructs.tree_view(rbtree)

end

#=
module binary_tree_test
using ..datastructs

btree = datastructs.BinaryTree{Int}(10)
level1_left  = datastructs.BinaryTreeNode(5)
level1_right = datastructs.BinaryTreeNode(15)

level12_left  = datastructs.BinaryTreeNode(3)
level12_right = datastructs.BinaryTreeNode(7)

level22_left  = datastructs.BinaryTreeNode(12)
level22_right = datastructs.BinaryTreeNode(20)

level312_left  = datastructs.BinaryTreeNode(4)
level312_right = datastructs.BinaryTreeNode(2)

insert!(btree, level1_left)
insert!(btree, level1_right)

insert!(btree, level12_left)
insert!(btree, level12_right)

insert!(btree, level22_left)
insert!(btree, level22_right)

insert!(btree, level312_left)
insert!(btree, level312_right)

datastructs.tree_view(btree)

datastructs.successor(ttest.btree.node.left.right).data

end
=#

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
overflows(::Type{T}) where {T} = true
overflows(::Type{BigInt}) = false
overflowval(::Type{T}) where {T <: Integer} = div(typemax(T) - T(9), T(10))
# if we eventually support non-base 10
# overflowval(::Type{T}, base) where {T <: Integer} = div(typemax(T) - base + 1, base)

@inline function typeparser(::Type{T}, source, pos, len, b, code, options) where {T <: Integer}
    x = zero(T)
    neg = false
    # start actual int parsing
    neg = b == UInt8('-')
    if neg || b == UInt8('+')
        pos += 1
        incr!(source)
    end
    if eof(source, pos, len)
        code |= INVALID | EOF
        @goto done
    end
    b = peekbyte(source, pos) - UInt8('0')
    if b > 0x09
        # character isn't a digit, INVALID value
        code |= INVALID
        @goto done
    end
    while true
        x = T(10) * x + unsafe_trunc(T, b)
        pos += 1
        incr!(source)
        if eof(source, pos, len)
            x = ifelse(neg, -x, x)
            code |= OK | EOF
            @goto done
        end
        b = peekbyte(source, pos) - UInt8('0')
        if b > 0x09
            # detected a non-digit, time to bail on value parsing
            x = ifelse(neg, -x, x)
            code |= OK
            @goto done
        end
        overflows(T) && x > overflowval(T) && break
    end
    # extra loop because we got too close to overflowing while parsing digits
    x = ifelse(neg, -x, x)
    while true
        x, ov_mul = Base.mul_with_overflow(x, T(10))
        x, ov_add = Base.add_with_overflow(x, ifelse(neg, -T(b), T(b)))
        if ov_mul | ov_add
            # we overflowed, mark as OVERFLOW
            code |= OVERFLOW
            # in this case, we know b is a digit that caused us
            # to overflow, so we don't really want to consider it for
            # a close quote character or delimiter, so if b is the last character
            # let's just bail as eof
            pos += 1
            incr!(source)
            if eof(source, pos, len)
                code |= EOF
            end
            @goto done
        end
        pos += 1
        incr!(source)
        if eof(source, pos, len)
            code |= OK | EOF
            @goto done
        end
        b = peekbyte(source, pos) - UInt8('0')
        if b > 0x09
            code |= OK
            @goto done
        end
    end

@label done
    return x, code, pos
end


text:
using Test
using AnalogDataAssimilation

include("test_model_forecasting.jl")
include("test_lorenz96.jl")
include("test_analog_forecasting.jl")
include("test_time_series.jl")
include("test_generate_data.jl")


text:
"""
    StochInt

StochInt is a collection of functions about stochastic integral. It consists of:
1. simu_BM(): generate samples of Brownian motion (discretized) which is a drift-diffusion process
2. simu_OU(): generate samples of Ornstein-Uhlenbeck process (discretized) which is a mean-reverting process
3. solve_SS_MarkovChain(): solve steady state (stationary distribution) for a given Markov chain (of course, discrete)
4. solve_SS_MarkovChain_FUCK(): a fucking savage algorithm to compute steady state, through computing P^n for enough n degrees until a tolerance is met
5. simu_MarkovChain(): generate samples of a Markov Chain, requires initial state (integer from 1,...,k; where k is the number of all states )

And there are some non-export underlying functions:
1. randDiscrete(): generate a random integer from a given discrete distribution, e.g. [0.1, 0.2, 0.4, 0.3] --> 2 from {1,2,3,4}

There are some notes:
1. I do not provide functions for Geometric Brownian Motion (GBM), because it can always be converted from a Brownian motion by exponentialization;

References:
1. "Stokey, N. L. . (2008). The economics of inaction: stochastic control models with fixed costs."
2. "Press, William H., Flannery, B. P., Teukolsky, S. A., & Vetterling, W. T. (1988). Numerical recipes in C: the art of scientific computing. Numerical Recipes in C The Art of Scientific Computing."

If any question or issue, please email to:
{tianhao.zhao@emory.edu}, Tianhao Zhao (GitHub: Clpr)
"""
module StochInt
    export randDiscrete, simu_BM, simu_OU, simu_MarkovChain, solve_SS_MarkovChain_FUCK

# ----------
"""
    randDiscrete( pVec::Vector{Float64} )

generates a random integer based on a given discrete distribution.
e.g. input a pVec = [0.3, 0.4, 0.3], select one element from {1,2,3},
where {1,2,3} correspond to probabilities [0.3,0.4,0.3] in order.

Receives:
1. pVec::Vector{Float64}, discrete distribution function
Returns:
1. Res::Int, a random integer

Timing:
about 0.0015 seconds (10.00 k allocations: 1.068 MiB) for F = [0.1,0.2,0.3,0.4] per 10000 times;
the testing is in main namespace, which is much slower than called by other functions in practice.

Mathematics:
for a discrete distribution like F = {1,...,k} -> [p1,...,pk] where sum(F) = 1.0,
firstly generate a uniform random number R from [0,1] uniform distribution,
then try to locate R between two probabilities pi & pj;
"""
randDiscrete( pVec::Vector{Float64} ) = begin
    # validation
    @assert( isapprox(1.0, sum(pVec)) , "requires a discrete distribution whose added probabilities is 1, but received $(sum(pVec))" )
    # generate
    local kState = length(pVec)
    local pVecSum = insert!( cumsum(pVec), 1, -1E-8 )  # convert to CDF & add a first minor negative number (to avoid the branch for exact 0)
    # to avoid the branch for exact 1
    pVecSum[end] += 1E-8
    local RndUnif = rand()  # generate a random number from [0,1] uniform distribution
    local Res = 1::Int # result

    while Res < kState  # totally compare T times (including in [0.0, p1])
        if pVecSum[Res] <= RndUnif < pVecSum[Res + 1]
            break
        end
        Res += 1
    end
    return Res::Int
end # end function randDiscrete
# ----------
"""
    simu_BM( T::Int ; N::Int = 1, mu::Float64 = 0.0, sigma::Float64 = 1.0 )

approximates a Brownian motion `` dX_t = \\mu dt + \\sigma dW_t `` with parameter ``( \\mu, \\sigma^2 )``,
where ``X_t`` is a Brownian motion with **constant** drift ``\\mu`` and **constant** diffusion ``\\sigma``, ``W_t`` is a standard Gaussian distribution;

Receives:
1. T::Int, how many periods to simulate (e.g. years)
2. N::Int, how many points during one period, by default 1
3. mu::Float64, drift term, a constant, by default 0.0
4. sigma::Float64, diffusion term, a constant, by default 1.0
Returns:
1. Vector{Float64} = Array{Float64,1}, a simulated path, with length T * N, starting from X0 = 0

According to (Stokey, 2008) Ch 2.4, I use ``X_{t+\\Delta t} = X_{t} + \\Delta X`` approximate the defined Brownian motion,
where ``\\Delta X = \\begin{cases}  h, P(\\Delta X = h) = p  \\\\  -h, P(\\Delta X = -h) = 1 - p     \\end{cases}`` is a random variable.
We have ``h = \\begin{cases} \\frac{1}{2}( 1 + \\frac{\\mu\\sqrt{\\Delta t}}{\\sqrt{ \\sigma^2 + \\mu^2 \\Delta t }} )   \\\\  \\approx  \\frac{1}{2}( 1 + \\frac{\\mu}{\\sigma}\\sqrt{\\Delta t} )  ,  \\Delta t \\ll \\frac{\\sigma^2}{\\mu^2}  \\end{cases}``,
and `` h = \\sigma \\sqrt{\\Delta t} ``;
"""
function simu_BM( T::Int ; N::Int = 1, mu::Float64 = 0.0, sigma::Float64 = 1.0 )
    # decalre
    local DeltaT = 1 / N  # step length
    local Res = [0.0]  # result, use push!() then to acquire higher performance
    local RndShock = rand(Float64, T * N - 1 )  # generating random numbers first; NOTE: minus 1 because it is increments
    # NOTE: I do not use conventional approximation ``pThreshold = \\frac{1}{2}()`` to allow user-friendly simulations of daily stock returns which has DeltaT = 1 (usually larger than ``\\sigma^2/\\mu^2``)
    local pThreshold = 0.5 * ( 1.0 + mu * sqrt(DeltaT) / sqrt( sigma^2 + mu^2 * DeltaT) )  # probability threshold to RISE hMove::Float64, for each discretized delta-t period
    local hRise = sigma * sqrt(DeltaT)  # movement size, if pThreshold happens, rise hRise; otherwise, fall hRise
    # generate Delta X
    for t in 1:(T * N - 1)
        if RndShock[t] <= pThreshold
            push!( Res, Res[t] + hRise )  # rise case
        else
            push!( Res, Res[t] - hRise )  # fall case
        end
    end
    # return
    return Res::Vector{Float64}
end  # end function simu_BM
# -----------
"""
    simu_OU( T::Int ; N::Int = 1, mu::Float64 = 0.0, sigma::Float64 = 1.0 )

generates simulated series following a Ornstein-Uhlenbeck process with parameter ``(\\mu, \\sigma^2)``,
where ``\\mu`` is drift term, ``\\sigma`` is diffusion term;
According to (Stokey, 2008), we use algorithms similar to that in simu_BM()::Function to approximate
`` dX_t = -\\mu dt + \\sigma dW_t ``, where ``W_t`` is a Wiener process.

Receives:
1. T::Int, how many periods to simulate (e.g. years)
2. N::Int, how many points during one period, by default 1
3. mu::Float64, drift term, a constant, by default 0.0
4. sigma::Float64, diffusion term, a constant, by default 1.0
Returns:
1. Vector{Float64} = Array{Float64,1}, a simulated path, with length T * N, starting from X0 = 0

Mathematics:
1. ``h = \\sigma \\sqrt{\\Delta t}``
2. `` p(X_t) = 0.5 - \\frac{\\alpha X_t}{ 2\\sigma } \\sqrt{\\Delta t} ``
3. `` \\Delta X_t = \\begin{cases} h, \\text{Prob}=p(X_t)  \\\\ -h, \\text{Prob} = 1 - p(X_t)  \\end{cases} ``
3. then, `` X_{t+\\Delta t} = X_t + \\Delta X_t ``

"""
function simu_OU( T::Int ; N::Int = 1, mu::Float64 = 0.0, sigma::Float64 = 1.0 )
    # decalre
    local DeltaT = 1 / N  # step length
    local Res = [0.0]  # result, use push!() then to acquire higher performance
    local RndShock = rand(Float64, T * N - 1 )  # generating random numbers first; NOTE: minus 1 because it is increments
    # NOTE: different from BM & GBM, the pThreshold::Float64 for a OU process may depend on X_t, which means we have to dynamically compute it
    local hRise = sigma * sqrt(DeltaT)  # movement size, if pThreshold happens, rise hRise; otherwise, fall hRise
    # generate Delta X
    for t in 1:( T * N - 1 )
        if RndShock[t] <= (  0.5 - mu * Res[t] * sqrt(DeltaT) / 2 / sigma  )  # this is pThreshold
            push!( Res, Res[t] + hRise )
        else
            push!( Res, Res[t] - hRise )
        end
    end
    return Res::Vector{Float64}
end  # end function simu_OU
# -----------------
# NOTE: need debug (math is wrong)
# """
#     solve_SS_MarkovChain( Pmat::Array{Float64,2} )
#
# solves the steady state (stationary distribution) for a given Markov chain,
# where the Markov chain is denoted with a one-step matrix ``P``.
# In matrix P, element x{row,coloumn} indicates the one-step transition probability
# from state row to state column;
#
# Receives:
# 1. Pmat::Array{T,2} where T, one-step transition probability matrix
# Returns:
# 1. Res::Vector{Float64}, vector of steady state (probabilities), states are sorted along with the order of matrix P's rows
#
# Notes:
# 1. we just solve the equations `` \\vec{\\pi} \\cdot \\bm{P} = \\vec{\\pi} ``
# 2. following conventions, we add a new all-one column to matrix P to constrain the steady state: `` \\sum^k \\pi_k = 1 ``
# 2. only least validations (e.g. inversibility of matrix P) to acquire better performance
# 3. without loss of generality, we only provide a Float64 version function
#
# Mathematics (how to solve `` \\vec{\\pi} \\cdot \\bm{P} = \\vec{\\pi} ``):
# 1. write element-level equation for a specific state j in total k states: `` \\pi_1 p_{1,j} + \\dots + \\pi_jp_{j,j} + \\dots + \\pi_{k} p_{k,j} = \\pi_j ``
# 2. move ``\\pi_j`` to the left side: `` \\pi_1 p_{1,j} + \\dots + \\pi_j(p_{j,j}-1) + \\dots + \\pi_{k} p_{k,j} = 0  ``
# 3. substitute the constraint of ``\\pi_j``: ``\\sum^k_j \\pi_j = 1``
# 4. get: `` \\pi_1 p_{1,j} + \\dots + (1-\\sum^{i\\neq j} \\pi_i)(p_{j,j}-1) + \\dots + \\pi_{k} p_{k,j} = 0   ``
# 5. simplify and get: `` \\pi_1 (p_{1,j} - p_{j,j} + 1) + \\dots + \\pi_{j} 0 + \\dots + \\pi_k(p_{k,j} - p_{j,j} + 1) = 1 - p_{j,j}  ``
# 6. write the equations in a linear algebra form: `` \\pi_{1\\times k} A_{k \\times k} = b_{1 \\times k} ``
# 7. to compute matrix A, firtst compute a matrix A': `` A'_{k\\times k} = P_{k\\times k} + 1_{k\\times k} - 1_{k\\times 1}\\cdot \\text{diag}(P_{k\\times k})_{1\\times k} ``
# 8. then matrix A is: `` A_{k\\times k} = A'_{k\\times k} - \\text{diag}(A')_{k\\times k} ``
# 9. matrix b is: `` b_{1\\times k} = 1_{1\\times k} - \\text{diag}(P_{k\\times k})_{1\\times k} ``
# 10. please note, index j traversals on rows, and index k traversals on columns
# """
# function solve_SS_MarkovChain( Pmat::Array{Float64,2} )
#     # declare
#     local kState = size(Pmat)  # number of states
#     # assertion
#     @assert( kState[1] == kState[2], " requires an n*n Array{Float64,2} matrix but received a $(kState) $(typeof(Pmat)) " )
#     # construct matrix A & matrix b
#     local Amat = copy(Pmat); local bmat = zeros(Float64, 1, kState[1])
#     for idxj in 1:kState[1]  # loop on rows/states
#         # construct matrix b
#         bmat[1,idxj] = 1.0 - Amat[idxj,idxj]
#         # construct matrix A'
#         for idxk in 1:kState[2]  # loop on columns (though A is a n*n matrix, we use [1],[2] to distinguish rows & columns)
#             Amat[idxj,idxk] += 1.0 - Amat[idxj,idxj]
#         end
#         # then set diag(A) as 0 (equals to A = A' - diag(A'))
#         Amat[idxj,idxj] = 0.0
#     end
#     # solve the linear equation system
#     local Res::Vector{Float64} = ( bmat * inv(Amat) )[1,:]  # squeeze to 1-D vector
#
#     return Res::Vector{Float64}
# end # end function solve_SS_MarkovChain
# --------------------
"""
    solve_SS_MarkovChain_FUCK( Pmat::Array{Float64,2} ; atol::Float64 = 1E-6, maxiter::Int = 100 )

solves the steady state of a given Markov chain with one-step transition probability matrix Pmat::Array{Float64,2},
where element p_{i,j} denotes the probability from state i to state j.
uses a fucking savage algorithm: computing P*P*P.... until converges.

Receives:
1. Pmat::Array{Float64,2}, one-step transition probability matrix
2. atol::Float64, tolerance to converge
3. maxiter::Int, maximum iteration times to stop
Retruns:
1. Res::Vector{Float64}, result vector (steady state distribution)

(Actually, fuckingly naive but very simple & high performance!)
"""
function solve_SS_MarkovChain_FUCK( Pmat::Array{Float64,2} ; atol::Float64 = 1E-6, maxiter::Int = 100 )
    # assertion (one-step transition matrix definition)
    local kState = size(Pmat)  # number of states
    local Pmat2 = copy(Pmat)  # temp variable
    @assert( kState[1] == kState[2], " requires an n*n Array{Float64,2} matrix but received a $(kState) $(typeof(Pmat)) " )
    @assert( all( abs.(sum(Pmat, dims=2) .- 1.0) .< 1E-6 ), "not a one-step transition matrix, requires the sum of every row is one!")
    # decalre
    local Res = zeros(Float64,kState[1]);  # result, also temp variable
    local flag_exit = true;  # convergence flag, false to stop loops
    local counter = 1;  # counter
    # fuckingly compute
    while flag_exit
        Res[:] .= Pmat2[1,:]  # current value
        Pmat2[:,:] *= Pmat2 * Pmat  # update
        if counter >= maxiter | isapprox(Pmat2[1,:], Res, atol = atol)  # if converged
            flag_exit = false
        end
        counter += 1
    end
    return Res::Vector{Float64}
end  # end function solve_SS_MarkovChain_FUCK
# ----------------
"""
    simu_MarkovChain( Pmat::Matrix{Float64}, InitState::Int, T::Int )

generates a random Markov Chain when given one-step transition matrix, an initial state and periods of time to simulate.

Receives:
1. Pmat::Matrix{Float64}, one-step transition probability matrix, where P{i,j} means the probability from state i to state j
2. InitState::Int, initial state to start, from 1 to size(Pmat)[1] (row/column number of Pmat)
3. T::Int, time periods to simulate, a finite positive integer
Returns:
1. Res::Vector{Int}, a simulated Markov Chain with length T
Dependency [no external]:
1. StochInt.randDiscrete(), generate a random integer from a given discrete distribution

"""
function simu_MarkovChain( Pmat::Matrix{Float64}, InitState::Int, T::Int )
    # declare
    local kState = size(Pmat)  # size of one-step transition matrix
    local Res = [InitState]  # result, a Vector{Int}
    # validation (one-step transition matrix definition)
    @assert( kState[1] == kState[2], " requires an n*n Array{Float64,2} matrix but received a $(kState) $(typeof(Pmat)) " )
    @assert( all( abs.(sum(Pmat, dims=2) .- 1.0) .< 1E-6 ), "not a one-step transition matrix, requires the sum of every row is one!")
    # generate
    for t in 1:T-1 # do not use 2:T to avoid bounding error
        # select a discrete distribution from the transition matrix according to current state
        # then randomly generate a new state
        push!(Res, randDiscrete( Pmat[Res[end],:] ) )  # [:] auto squeezes the row-vector to a one-dimension vector
    end
    return Res::Vector{Int}
end # end function simu_MarkovChain
# -------------

end  # end module: StochInt

max_stars_repo_name: Clpr/NumToolLibs

max_issues_repo_name: Clpr/NumToolLibs

max_forks_repo_name: Clpr/NumToolLibs


text:
# misc dofhandler unit tests
@testset "dofs" begin

# set up a test DofHandler
grid = generate_grid(Triangle, (10, 10))
dh = DofHandler(grid)
push!(dh, :u, 2, Lagrange{2,RefTetrahedron,2}())
push!(dh, :p, 1, Lagrange{2,RefTetrahedron,1}())
close!(dh)

# renumber!
perm = randperm(ndofs(dh))
iperm = invperm(perm)
dofs = copy(dh.cell_dofs)
Ferrite.renumber!(Ferrite.renumber!(dh, perm), iperm)
@test dofs == dh.cell_dofs

# dof_range
@test (@inferred dof_range(dh, :u)) == 1:12
@test (@inferred dof_range(dh, :p)) == 13:15
# dof_range for FieldHandler (use with MixedDofHandler)
ip = Lagrange{2, RefTetrahedron, 1}()
field_u = Field(:u, ip, 2)
field_c = Field(:c, ip, 1)
fh = FieldHandler([field_u, field_c], Set(1:getncells(grid)))
@test dof_range(fh, :u) == 1:6
@test dof_range(fh, :c) == 7:9

end # testset

@testset "Dofs for Line2" begin

nodes = [Node{2,Float64}(Vec(0.0,0.0)), Node{2,Float64}(Vec(1.0,1.0)), Node{2,Float64}(Vec(2.0,0.0))]
cells = [Line2D((1,2)), Line2D((2,3))]
grid = Grid(cells,nodes)

#2d line with 1st order 1d interpolation
dh = DofHandler(grid)
push!(dh, :x, 2)
close!(dh)

@test celldofs(dh,1) == [1,2,3,4]
@test celldofs(dh,2) == [3,4,5,6]

#2d line with 2nd order 1d interpolation
dh = DofHandler(grid)
push!(dh, :x, 2, Lagrange{1,RefCube,2}())
close!(dh)

@test celldofs(dh,1) == [1,2,3,4,5,6]
@test celldofs(dh,2) == [3,4,7,8,9,10]

#3d line with 2nd order 1d interpolation
dh = DofHandler(grid)
push!(dh, :u, 3, Lagrange{1,RefCube,2}())
push!(dh, :θ, 3, Lagrange{1,RefCube,2}())
close!(dh)

@test celldofs(dh,1) == collect(1:18)
@test celldofs(dh,2) == [4,5,6,19,20,21,22,23,24,    # u
                         13,14,15,25,26,27,28,29,30] # θ
end

@testset "Dofs for quad in 3d (shell)" begin

nodes = [Node{3,Float64}(Vec(0.0,0.0,0.0)), Node{3,Float64}(Vec(1.0,0.0,0.0)),
            Node{3,Float64}(Vec(1.0,1.0,0.0)), Node{3,Float64}(Vec(0.0,1.0,0.0)),
            Node{3,Float64}(Vec(2.0,0.0,0.0)), Node{3,Float64}(Vec(2.0,2.0,0.0))]

cells = [Quadrilateral3D((1,2,3,4)), Quadrilateral3D((2,5,6,3))]
grid = Grid(cells,nodes)

#3d quad with 1st order 2d interpolation
dh = DofHandler(grid)
push!(dh, :u, 3, Lagrange{2,RefCube,1}())
push!(dh, :θ, 3, Lagrange{2,RefCube,1}())
close!(dh)

@test celldofs(dh,1) == collect(1:24)
@test celldofs(dh,2) == [4,5,6,25,26,27,28,29,30,7,8,9, # u
                         16,17,18,31,32,33,34,35,36,19,20,21]# θ

#3d quads with two quadratic interpolations fields
#Only 1 dim per field for simplicity...
dh = DofHandler(grid)
push!(dh, :u, 1, Lagrange{2,RefCube,2}())
push!(dh, :θ, 1, Lagrange{2,RefCube,2}())
close!(dh)

@test celldofs(dh,1) == collect(1:18)
@test celldofs(dh,2) == [2, 19, 20, 3, 21, 22, 23, 6, 24, 11, 25, 26, 12, 27, 28, 29, 15, 30]

# test reshape_to_nodes
## DofHandler
mesh = generate_grid(Quadrilateral, (1,1))
dh = DofHandler(mesh)
push!(dh, :v, 2)
push!(dh, :s, 1)
close!(dh)

u = [1.1, 1.2, 2.1, 2.2, 4.1, 4.2, 3.1, 3.2, 1.3, 2.3, 4.3, 3.3]

s_nodes = reshape_to_nodes(dh, u, :s)
@test s_nodes ≈ [i+0.3 for i=1:4]'
v_nodes = reshape_to_nodes(dh, u, :v)
@test v_nodes ≈ [i==3 ? 0.0 : j+i/10 for i=1:3, j=1:4]
end


text:
function getBirkhoffStepSize(x_star,x,y)
    return minimum((x_star-x).*y - (y.-1));
end