text:
# This file holds the function to upload the initial poulation snapshot.

export uploadSnapshot

function uploadSnapshot( mpSim::ManpowerSimulation, snapName::String,
    idCol::Int, colsToImport::Dict{Int,String} = Dict{Int, String}(),
    catalogueName::String = ""; recCol::Int = -1, isRecDate::Bool = true,
    ageCol::Int = -2, defaultRecAge::Float64 = 18.0, isBirthDate::Bool = true,
    stateCol::Int = -3 )::Void

    tmpSnapName = snapName
    tmpSnapName *= endswith( snapName, ".xlsx" ) ? "" : ".xlsx"

    if !ispath( tmpSnapName )
        warn( "The file `$tmpSnapName` does not exist. Cannot upload snapshot." )
        return
    end  # if !ispath( tmpSnapName )

    # Check if required columns (id, age, etc...) are all different
    colsToIgnore = [ idCol, recCol, ageCol, stateCol ]

    if length( unique( colsToIgnore ) ) < length( colsToIgnore )
        warn( "Some of the required columns (id, age, tenure, state) are the same. Cannot upload snapshot." )
        return
    end  # if length( unique( colsToIgnore ) ) < length( colsToIgnore )

    verifyAttributes( colsToImport, colsToIgnore, mpSim )

    # Make sure the name of the ID column is the one stored in the
    #   simulation configuration.
    if idCol > 0
        colsToImport[ idCol ] = mpSim.idKey
    end  # if idCol > 0

    # Add recruitment time information to the database.
    if recCol > 0
        colsToImport[ recCol ] = "timeEntered"
    end  # if recCol > 0

    # Add recruitment age information to the database.
    if ageCol > 0
        colsToImport[ ageCol ] = "ageAtRecruitment"
    end  # if ageCol > 0

    dataColsToGrab = collect( keys( colsToImport ) )
    systemAttrs = [ mpSim.idKey, "timeEntered", "ageAtRecruitment",
        "expectedRetirementTime" ]

    XLSX.openxlsx( tmpSnapName ) do xf
        # Get the personnel data.
        dataSheet = xf[ 1 ]
        nEntries, nCols = XLSX.size( XLSX.get_dimension( dataSheet ) )
        nEntries -= 1
        excessCols = filter( ii -> ii > nCols, dataColsToGrab )

        if nEntries == 0
            warn( "Snapeshot file '$snapName' empty." )
            return
        end  # if nEffectiveEntries == 0

        # Purge non-existent columns.
        if !isempty( excessCols )
            filter!( ii -> ii <= nCols, dataColsToGrab )
            warn( "Column(s) $(join( excessCols, ", " )) do(es) not exist in snapshot. Not importing those columns." )
        end  # if !isempty( excessCols )

        contents = XLSX.getdata( dataSheet )[ 2:end, dataColsToGrab ]
        colNames = map( ii -> colsToImport[ ii ], dataColsToGrab )
        idColIndex = findfirst( dataColsToGrab .== idCol )

        # Add generated id if no ID column is present in the retrieved data.
        if ( idCol <= 0 ) || ( idCol > nCols )
            contents = hcat( contents,
                "Sim" .* string.( collect( 1:nEntries ) ) )
            push!( colNames, mpSim.idKey )
        else
            uniqueIndices = indexin( unique( contents[ :, idColIndex ] ),
                contents[ :, idColIndex ] )
            contents = contents[ uniqueIndices, : ]
        end  # if ( idCol <= 0 ) || ...

        nEffectiveEntries = size( contents )[ 1 ]
        println( nEntries - nEffectiveEntries,
            " entry/ies with duplicate IDs. Retaining only one of each." )

        if catalogueName != ""
            tmpCatalogueName = endswith( catalogueName, ".xlsx" ) ?
                catalogueName : catalogueName * ".xlsx"
            isEntryOkay = validateContents( contents,
                colNames[ vcat( 1:(idColIndex-1), (idColIndex+1):end ) ],
                mpSim.idKey, tmpCatalogueName )
            contents = contents[ isEntryOkay, : ]
            nEffectiveEntries = sum( isEntryOkay )

            if nEffectiveEntries == 0
                warn( "No valid personnel entries left in snapshot file '$snapName'." )
                return
            end  # if nEffectiveEntries == 0
        end  # if catalogueName != ""

        recColIndex = length( colNames ) + 1

        # Adjust recruitment time information.
        if ( recCol <= 0 ) || ( recCol > nCols )
            contents = hcat( contents, zeros( Float64, nEffectiveEntries ) )
            push!( colNames, "timeEntered" )
        else
            recColIndex = findfirst( colName -> colName == "timeEntered",
                colNames )

            if isRecDate
                contents[ :, recColIndex ] = map( recDate ->
                    -( mpSim.simStartDate - recDate ).value,
                    contents[ :, recColIndex ] ) ./ ( 365.0 / 12.0 )
            else
                contents[ :, recColIndex ] *= -12.0
            end  # if isRecDate
        end  # if ( recCol <= 0 ) || ...

        ageColIndex = length( colNames ) + 1

        # Adjust recruitment age information.
        if ( ageCol <= 0 ) || ( ageCol > nCols )
            contents = hcat( contents, ones( Float64, nEffectiveEntries ) .*
                ( 12.0 * defaultRecAge ) )
            push!( colNames, "ageAtRecruitment" )
        else
            ageColIndex = findfirst( colName -> colName == "ageAtRecruitment",
                colNames )

            # Determine ages in months.
            if isBirthDate
                contents[ :, ageColIndex ] = map( birthDate ->
                    ( mpSim.simStartDate - birthDate ).value,
                    contents[ :, ageColIndex ] ) ./ ( 365.0 / 12.0 )
            else
                contents[ :, ageColIndex ] *= 12.0
            end  # if isBirthDate

            contents[ :, ageColIndex ] += contents[ :, recColIndex ]
        end  # if ( ageCol <= 0 ) || ...

        # Determine retirement ages/times.
        push!( colNames, "expectedRetirementTime" )
        contents = hcat( contents, zeros( Float64, nEffectiveEntries ) )

        # Retrieve time of last transition.
        lastTransTime = - 12.0 .* ones( nEntries )

        if 0 < stateCol <= nCols
            lastTransTime = map( transDate ->
            -( mpSim.simStartDate - transDate ).value,
                XLSX.getdata( dataSheet )[ 2:end, stateCol ] ) ./
                ( 365.0 / 12.0 )
        end  # if 0 < stateCol <= nCols

        # Discover states of personnel members.
        attrList = filter( attrName -> attrName ∉ systemAttrs, colNames )
        attrIndices = map( attrName -> findfirst( attrName .== colNames ),
            attrList )
        attrVals = Dict{String, Any}()
        transCmd = "('" .* contents[ :, idColIndex ] .* "', " .*
            string.( contents[ :, recColIndex ] ) .* ", 'snapshot', 'active')"
        stateList = Dict{String, Vector{State}}()
        attritionSchemes = Vector{String}( nEffectiveEntries )
        attritionTimes = Vector{Union{String, Float64}}( nEffectiveEntries )

        for ii in 1:nEffectiveEntries
            for jj in attrIndices
                attrVals[ colNames[ jj ] ] = contents[ ii, jj ]
            end  # for jj in attrIndices

            # Identify all initial states the person belongs to and add entry info to
            #   each of those states.
            persStates = collect( Iterators.filter(
                state -> isPersonnelOfState( attrVals,
                    mpSim.stateList[ state ] ), keys( mpSim.stateList ) ) )
                # XXX Iterators.filter is needed to avoid deprecation warnings.

            id = contents[ ii, idColIndex ]

            # Check if each person can be assigned to exactly one state.
            if mpSim.isWellDefined
                if isempty( persStates )
                    mpSim.isWellDefined = false
                    warn( "Person read that couldn't be assigned to any state. Please check if all required attributes to determine state have been read." )
                elseif length( persStates ) > 1
                    mpSim.isWellDefined = false
                    warn( "Person read that can be assigned to multiple states. Please check system configuration for consistency." )
                end  # if isempty( initPersStates )
            end  # if mpSim.isWellDefined

            # Determine person's expected retirement age.
            stateRetAge = isempty( persStates ) ? 0 :
                mpSim.stateList[ persStates[ 1 ] ].stateRetAge
            contents[ ii, end ] = computeExpectedRetirementTime(
                mpSim, mpSim.retirementScheme, contents[ ii, ageColIndex ],
                stateRetAge, contents[ ii, recColIndex ] )
            attrScheme = mpSim.defaultAttritionScheme
            transTime = lastTransTime[ ii ]
            timeOfAttr = transTime

            # Add person to state list.
            for stateName in persStates
                state = mpSim.stateList[ stateName ]
                state.inStateSince[ id ] = transTime
                state.isLockedForTransition[ id ] = false
                attrScheme = state.attrScheme
                push!( transCmd,
                    "('$id', $transTime, 'snapshot', '$stateName')" )
            end  # for state in persStates

            while isa( timeOfAttr, Real ) && ( timeOfAttr <= 0.0 )
                timeOfAttr = generateTimeOfAttrition( attrScheme,
                    transTime )
            end  # while isa( timeOfAttr, Real )

            attritionSchemes[ ii ] = attrScheme.name
            attritionTimes[ ii ] = timeOfAttr
        end  # for ii in 1:nEffectiveEntries

        map!( dataPoint -> isa( dataPoint, String ) ? "'" * dataPoint * "'" :
            dataPoint, contents, contents )

        contents = hcat( contents, attritionTimes,
            "'" .* attritionSchemes .* "'" )

        # Inject data into database.
        persCmd = map( ii -> "(" * join( contents[ ii, : ], ", " ) *
            ", 'active')", 1:nEffectiveEntries )
        persCmd = "INSERT INTO $(mpSim.personnelDBname)
            ('$(join( colNames, "', '" ))', expectedAttritionTime, attritionScheme, status)
            VALUES $(join( persCmd, ", " ))"
        SQLite.execute!( mpSim.simDB, persCmd )

        # Generate the history database.
        isAttrFixed = Dict{String, Bool}()
        fixedAttrList = filter( attr -> attr.isFixed,
            vcat( mpSim.initAttrList, mpSim.otherAttrList ) )
        fixedAttrList = vcat( map( attr -> attr.name, fixedAttrList ),
            systemAttrs )
        varColNames = filter( attrName -> attrName ∉ fixedAttrList, colNames )

        for attrName in varColNames
            attrIndex = findfirst( attrName .== colNames )
            histCmd = "(" .* contents[ :, idColIndex ] .* ", '" .* attrName .*
                "', " .* string.( lastTransTime ) .* ", " .*
                contents[ :, attrIndex ] .* ")"
                # XXX -12.0 is chosen since actual information on last attribute
                #   change isn't available, and for consistency when tenure
                #   information is missing.
            histCmd = "INSERT INTO $(mpSim.historyDBname)
                ($(mpSim.idKey), attribute, timeIndex, strValue)
                VALUES $(join( histCmd, ", " ))"
            SQLite.execute!( mpSim.simDB, histCmd )
        end  # for attrName in varColNames

        # Generate transition database.
        transCmd = "INSERT INTO $(mpSim.transitionDBname)
            ($(mpSim.idKey), timeIndex, transition, endState )
            VALUES $(join( transCmd, ", " ))"
        SQLite.execute!( mpSim.simDB, transCmd )
        mpSim.personnelSize = nEffectiveEntries
        mpSim.resultSize = nEffectiveEntries

        # Statistics of entered data.
        println( "Entered $nEffectiveEntries of $nEntries persons ($(100.0 * nEffectiveEntries / nEntries)%) in snapshot into database." )

        # for ii in eachindex( colNames )
        #     if colNames[ ii ] ∉ systemAttrs
        #         counts = StatsBase.countmap( contents[ :, ii ] )
        #         println( "Attribute '$(colNames[ ii ])':" )
        #         foreach( val -> println( "   $val: $(counts[ val ])  ($(100.0 * counts[ val ] / nEffectiveEntries )%)" ), keys( counts ) )
        #     end  # if colNames( ii ) ∉ [ mpSim.idKey, ...
        # end  # for ii in eachindex( colNames )
    end  # XLSX.openxlsx( tmpSnapName ) do xf

    return

end  # uploadSnapshot( mpSim, snapName, idCol, colsToImport, catalogueName,
    #    recCol, isRecDate, ageCol, defaultRecAge, isBirthDate, stateCol )

# ==============================================================================
# Non-exported methods.
# ==============================================================================

function verifyAttributes( colsToImport::Dict{Int, String},
    colsToIgnore::Vector{Int}, mpSim::ManpowerSimulation )::Void

    # Import only columns which have been defined as attributes.
    attrNameList = map( attr -> attr.name,
        vcat( mpSim.initAttrList, mpSim.otherAttrList ) )
    colsNotToImport = Vector{Int}()

    for ii in keys( colsToImport )
        if ( ii ∉ colsToIgnore ) && ( colsToImport[ ii ] ∉ attrNameList )
            push!( colsNotToImport, ii )
        end  # if ii != idcol
    end  # for ii in keys( colsToImport )

    for ii in colsNotToImport
        delete!( colsToImport, ii )
    end  # for ii in colsNotToImport

    return

end  # verifyAttributes( colsToImport, colsToIgnore, mpSim )

function validateContents( dataMatrix::Array, attrNames::Vector{String},
    idKey::String, catalogueName::String )::Vector{Bool}

    nEntries = size( dataMatrix )[ 1 ]
    isEntryOkay = trues( nEntries )

    XLSX.openxlsx( catalogueName ) do xf
        dataSheet = xf[ "General" ]
        nAttrs = dataSheet[ "B5" ]
        dataSheet = xf[ "Attributes" ]
        catalogueAttrNames = dataSheet[ "A2:A$(nAttrs + 1)" ]

        for ii in eachindex( attrNames )
            attrName = attrNames[ ii ]

            if attrName ∉ [ idKey, "timeEntered", "ageAtRecruitment" ]
                attrIndex = findfirst( catalogueAttrNames .== attrName )
                nVals = dataSheet[ "E$(attrIndex + 1)" ]
                vals = dataSheet[ XLSX.CellRange( attrIndex + 1, 6,
                    attrIndex + 1, 5 + nVals ) ]
                isEntryOkay .&= map( val -> val ∈ vals, dataMatrix[ :, ii ] )
            end  # if attrName != idKey
        end  # for attrName in attrNames
    end  # XLSX.openxlsx( catalogueName ) do xf

    println( nEntries - sum( isEntryOkay ), " entries contained attribute values not defined in the catalogue." )

    return isEntryOkay

end  # validateContents( dataMatrix, attrNames, idKey, catalogueName )

function readSnapshot( mpSim::ManpowerSimulation )::Void

    tStart = now()

    XLSX.openxlsx( mpSim.parFileName ) do xf
        if !XLSX.hassheet( xf, "Snapshot" )
            warn( "Configuration file has no sheet named 'Snapshot'. Not uploading initial population." )
            return
        end  # if !hassheet( xf, "Snapshot" )

        snapSheet = xf[ "Snapshot" ]

        if snapSheet[ "B3" ] == "NO"
            return
        end  # if snapSheet[ "B3" ] == "NO"

        snapName = snapSheet[ "B4" ]
        snapName *= endswith( snapName, ".xlsx" ) ? "" : ".xlsx"
        snapName = joinpath( dirname( mpSim.parFileName ), snapName )

        if !ispath( snapName )
            warn( "Snapshot file '$snapName' does not exist. Not uploading initial population." )
        end  # if !ispath( snapName )

        # Get special column parameters.
        idCol = snapSheet[ "B5" ]
        idCol = isa( idCol, Missings.Missing ) ? 0 : idCol
        recCol = snapSheet[ "B6" ]
        recCol = isa( recCol, Missings.Missing ) ? -1 : recCol
        ageCol = snapSheet[ "B8" ]
        ageCol = isa( ageCol, Missings.Missing ) ? -2 : ageCol
        stateCol = snapSheet[ "B10" ]
        stateCol = isa( stateCol, Missings.Missing ) ? -3 : stateCol
        isRecDate = snapSheet[ "B7" ] == "YES"
        isBirthDate = snapSheet[ "B9" ] == "YES"

        # Get columns to import.
        nCols = snapSheet[ "B13" ]
        colsToImport = Dict{Int, String}()
        foreach( ii -> colsToImport[ snapSheet[ "A$(ii+15)" ] ] =
            snapSheet[ "B$(ii+15)" ], 1:nCols )

        uploadSnapshot( mpSim, snapName, idCol, colsToImport, mpSim.catFileName,
            recCol = recCol, isRecDate = isRecDate, ageCol = ageCol,
            isBirthDate = isBirthDate, stateCol = stateCol )
    end  # XLSX.openXLSX( mpSim.parFileName) do xf

    println( "Processing initial population snapshot took ",
        ( now() - tStart ).value / 1000, " seconds." )

    return

end  # readSnapshot( mpSim )


text:
#=
composite_types:
- Julia version: 1.5.0
- Author: shisa
- Date: 2020-08-09
=#

module ComposTypes
    # struct
    struct Point
        x::Int
        y::Int
    end

    mutable struct PointMut
        x::Int
        y::Int
    end

    # methods
    function get_distance(p::Point)
        return sqrt(p.x^2 + p.y^2)
    end

    function main()
        p = Point(3, 5)
        println("defined struct Point p = ", p)
        println("p.x = $(p.x), p.y = $(p.y)")
        println("Distance = $(get_distance(p))")
        p_mut = PointMut(3, 5)
        println("defined mutable struct Point p_mut = ", p_mut)
        println("p_mut.x = $(p_mut.x), p_mut.y = $(p_mut.y)")
        println("set values to p_mut: x = 10, y = 20")
        p_mut.x = 10
        p_mut.y = 20
        println("p_mut.x = $(p_mut.x), p_mut.y = $(p_mut.y)")
    end
end

if abspath(PROGRAM_FILE) == @__FILE__
    using .ComposTypes
    ComposTypes.main()
end

ext: jl

max_stars_repo_name: ShisatoYano/JuliaPractice

max_issues_repo_name: ShisatoYano/JuliaPractice

max_forks_repo_name: ShisatoYano/JuliaPractice


text:
using CSV, Shapefile
using JSServe, WGLMakie, AbstractPlotting, Markdown
using JSServe.DOM
using JSServe: styled_slider
using AbstractPlotting.MakieLayout

read_data(path) = open(path) do io
    # skip first line since it contains an invalid comment
    readuntil(io, "\n")
    CSV.File(io, delim=';')
end

function normalize_names(x)
    no_umlaut = replace(x, r"[üöä]" => (x)-> Dict("ü" => "ue", "ö"=>"ae", "ä"=>"oe")[x])
    bb = replace(no_umlaut, "Brandenburg/Berlin" => "Berlin")
    if any(y-> occursin(y, bb), ("Hamburg", "Bremen", "Niedersachsen"))
        return "Niedersachsen"
    end
    return bb
end

data_dir = joinpath(@__DIR__, "data")
mkdir(data_dir)
data_temp_path = joinpath(data_dir, "temperature.txt")
data_prec_path = joinpath(data_dir, "precipitation.txt")
download("https://opendata.dwd.de/climate_environment/CDC/regional_averages_DE/annual/air_temperature_mean/regional_averages_tm_year.txt", data_temp_path)
download("https://opendata.dwd.de/climate_environment/CDC/regional_averages_DE/annual/precipitation/regional_averages_rr_year.txt", data_prec_path)

temp_data = read_data(data_temp_path)
prec_data = read_data(data_prec_path)

url_gadm = "https://biogeo.ucdavis.edu/data/gadm3.6/shp/gadm36_DEU_shp.zip"
shape_path = joinpath(data_dir, "gadm36_DEU_1.shp")
shape_zip = joinpath(data_dir, "gadm36_DEU_1.zip")
download(url_gadm, shape_zip)
run(`unzip $shape_zip -d $data_dir`)
tbl_states = Shapefile.Table(shape_path)

states = map(normalize_names, tbl_states.NAME_1)
states_geom_map = Dict(zip(map(normalize_names, tbl_states.NAME_1), tbl_states.NAME_1))
states_sym = propertynames(temp_data)[3:end]
statenames_to_sym = Dict(zip(normalize_names.(string.(states_sym)), states_sym))

states_to_plot = filter(states) do state
    haskey(statenames_to_sym, state)
end

function data_extrema(data)
    combinations = [(s, i) for s in states_to_plot, i in 1:length(data)]
    return extrema(map(combinations) do (state, idx)
        return data[idx][statenames_to_sym[state]]
    end)
end

function visualize_data(scene, data, cmap, colorrange, s)
    colormap = to_colormap(cmap)
    colors = map(s) do idx
        map(states_to_plot) do state
            datapoint = data[idx][statenames_to_sym[state]]
            return AbstractPlotting.interpolated_getindex(colormap, datapoint, colorrange)
        end
    end

    polys = map(states_to_plot) do state
        idx = findfirst(x-> x.NAME_1 == states_geom_map[state], tbl_states)
        return tbl_states[idx].Geometry
    end
    # Plot each state individually, so we can update one color per state
    for i in 1:length(states_to_plot)
        poly!(scene, polys[i], color=lift(x-> x[i], colors))
    end
    return scene
end
fig = Figure()
ax = fig[1,1] = Axis(fig)

fig = Figure(resolution = (900, 700))
to_plot = (
    (1, "Temperatur", temp_data, :heat),
    (2, "Regen", prec_data, :blues))
s = Observable(1)
foreach(to_plot) do (i, title, data, cmap)
    colorrange = data_extrema(data)
    ax = fig[1, i] = Axis(fig; title=title)
    hidedecorations!(ax, grid = false)
    visualize_data(ax, data, cmap, colorrange, s)
    layout[2, i] = Colorbar(fig, colormap=cmap, limits=colorrange, vertical=false, height=30, flipaxisposition=false, labelvisible=false,
    ticklabelalign=(:center, :top), ticksize=5)
    return
end
display(fig)

function handler(session, req)
    s = JSServe.Slider(1:length(temp_data))

    scene, layout = layoutscene(20, resolution = (900, 700))
    to_plot = (
        (1, "Temperatur", temp_data, :heat),
        (2, "Regen", prec_data, :blues))

    foreach(to_plot) do (i, title, data, cmap)
        colorrange = data_extrema(data)
        ax = layout[1, i] = LAxis(scene; title=title)
        hidedecorations!(ax, grid = false)
        visualize_data(ax, data, cmap, colorrange, s)
        layout[2, i] = LColorbar(scene, colormap=cmap, limits=colorrange, vertical=false, height=30, flipaxisposition=false, labelvisible=false,
        ticklabelalign=(:center, :top), ticksize=5)
        return
    end

    year = map(idx-> temp_data[idx].Jahr, s)

    markdown = md"""
    # Durchschnitts Temperatur per Bundesland

    ## Reise vom Jahr 1881 bis 2019
    $(styled_slider(s, year))

    $(scene)

    [Quelle: Deutscher Wetterdienst](https://opendata.dwd.de/climate_environment/CDC/regional_averages_DE/annual/)

    [Quellcode für Visualisierung](https://github.com/SimonDanisch/JSServe.jl/blob/master/examples/german_clima.jl)
    """

    dom = DOM.div(JSServe.MarkdownCSS, JSServe.TailwindCSS, JSServe.Styling, markdown)
    # return dom
    return JSServe.record_state_map(session, dom).dom
end

# Either export standalone
JSServe.export_standalone(handler, "dev/WGLDemos/german_heat", clear_folder=true)

# Or serve!
# app = JSServe.Server(handler, "0.0.0.0", 8082)


text:
using StaticArrays, Test
using OrdinaryDiffEq

u0 = fill(zero(MVector{2,Float64}), 2)
u0[1] = ones(MVector{2,Float64}) .+ 1
f = (du,u,p,t) -> du .= u
ode = ODEProblem(f, u0, (0.,1.))
sol = solve(ode, Euler(), dt=1.e-2)
sol = solve(ode, Tsit5())

u0 = fill(zero(SVector{2,Float64}), 2) .+ 1
u0[1] = ones(SVector{2,Float64}) .+ 1
ode = ODEProblem(f, u0, (0.,1.))
sol = solve(ode, Euler(), dt=1.e-2)
sol = solve(ode, Tsit5())

sol = solve(ode, SSPRK22(), dt=1.e-2)

u0 = ones(MVector{2,Float64})
ode = ODEProblem(f, u0, (0.,1.))
sol = solve(ode, Euler(), dt=1.e-2)
sol = solve(ode, Tsit5(), dt=1.e-2)

u0 = ones(SVector{2,Float64})
f = (u,p,t) -> u
ode = ODEProblem(f, u0, (0.,1.))
sol = solve(ode, Euler(), dt=1.e-2)
sol = solve(ode, ImplicitEuler(nlsolve=NLAnderson()), dt=1.e-2)
sol = solve(ode, Tsit5(), dt=1.e-2)


text:
using Test
using Printf: @printf
using LinearAlgebra: tr
using QuantumLattices.Essentials.Frameworks
using QuantumLattices.Essentials.Terms: Parameters
using QuantumLattices.Interfaces: id,register!
using QuantumLattices.Prerequisites: Float
import QuantumLattices.Interfaces: prepare!,run!,update!
import QuantumLattices.Essentials.Frameworks: dependences

struct FEngine <: Engine end
Base.repr(::FEngine)="FEngine"
Base.show(io::IO,::FEngine)=@printf io "FEngine"
struct FApp <: App end

@testset "Assignment" begin
    @test FApp()==FApp()
    @test isequal(FApp(),FApp())
    @test update!(FApp())==FApp()

    assign=Assignment(:FApp,FApp(),(t=1.0,U=8.0),dependences=(:FApp1,:FApp2))
    @test deepcopy(assign)==assign
    @test isequal(deepcopy(assign),assign)
    @test assign|>valtype==assign|>typeof|>valtype==Any
    @test assign|>id==assign|>typeof|>id==:FApp
    update!(assign,t=2.0)
    @test assign.parameters==(t=2.0,U=8.0)
end

@testset "Algorithm" begin
    @test FEngine()==FEngine()
    @test isequal(FEngine(),FEngine())
    @test update!(FEngine())==FEngine()

    alg=Algorithm("Alg",FEngine(),parameters=(t=1.0,U=8.0))
    @test string(alg)==repr(alg)=="Alg_FEngine_1.0_8.0"
    @test repr(alg,(:U,))=="Alg_FEngine_1.0"
    register!(alg,:FApp1,FApp(),parameters=(U=5.0,))
    register!(alg,:FApp2,FApp(),parameters=(U=6.0,),dependences=(:FApp1,))
    @test get(alg,Val(:FApp1))==get(alg,:FApp1)==Assignment(:FApp1,FApp(),(t=1.0,U=5.0),virgin=false)
    @test get(alg,Val(:FApp2))==get(alg,:FApp2)==Assignment(:FApp2,FApp(),(t=1.0,U=6.0),dependences=(:FApp1,),virgin=false)
    @test prepare!(alg,get(alg,:FApp1))===nothing
    @test run!(alg,get(alg,:FApp1))===nothing
    @test dependences(alg,get(alg,:FApp1))==()
    @test dependences(alg,get(alg,:FApp2))==(:FApp1,)
    @test dependences(alg,get(alg,:FApp2),(:FApp1,))==()
end

mutable struct VCA <: Engine
    t::Float
    U::Float
end
function update!(vca::VCA;kwargs...)
    vca.t=get(kwargs,:t,vca.t)
    vca.U=get(kwargs,:U,vca.U)
    return vca
end
Parameters(vca::VCA)=Parameters{(:t,:U)}(vca.t,vca.U)
gf(alg::Algorithm{VCA})=get(alg,Val(:_VCAGF_)).data

mutable struct GF <: App
    dim::Int
    count::Int
end

mutable struct DOS <: App
    mu::Float
end
update!(eb::DOS;kwargs...)=(eb.mu=get(kwargs,:mu,eb.mu);eb)

dependences(alg::Algorithm{VCA},assign::Assignment{GF},::Tuple{}=())=assign.dependences
dependences(alg::Algorithm{VCA},assign::Assignment,::Tuple{}=())=(:_VCAGF_,assign.dependences...)
prepare!(alg::Algorithm{VCA},assign::Assignment{GF})=assign.virgin && (assign.data=Matrix{valtype(assign)|>eltype}(undef,assign.app.dim,assign.app.dim))
run!(alg::Algorithm{VCA},assign::Assignment{GF})=(assign.app.count+=1; assign.data[:,:].=alg.engine.t+alg.engine.U)

function run!(alg::Algorithm{VCA},assign::Assignment{DOS})
    rundependences!(alg,assign)
    assign.data=tr(gf(alg))
end

@testset "Framework" begin
    engine=VCA(1.0,8.0)
    gf=Assignment(:_VCAGF_,GF(4,0),Parameters(engine),data=Matrix{Complex{Float}}(undef,0,0))
    vca=Algorithm("Test",engine,assignments=(gf,))
    register!(vca,:DOS,DOS(-3.5),parameters=(U=7.0,),map=params::Parameters->Parameters{(:t,:U,:mu)}(params.t,params.U,-params.U/2))
    dos=get(vca,:DOS)
    @test dos.data==32.0
    @test dos.app.mu==-3.5
    update!(dos,U=6.0)
    run!(vca,:DOS,false)
    @test dos.data==28.0
    @test dos.app.mu==-3.0
end


text:
struct NpyPickler{PROTO} <: AbstractPickle
  memo::Memo
  stack::PickleStack
  mt::HierarchicalTable
end

function NpyPickler(proto=DEFAULT_PROTO, memo=Dict())
  mt = HierarchicalTable()
  mt["numpy.core.multiarray._reconstruct"] = np_multiarray_reconstruct
  mt["numpy.dtype"] = np_dtype
  mt["numpy.core.multiarray.scalar"] = np_scalar
  mt["__build__.Pickle.NpyDtype"] = build_npydtype
  mt["__build__.Pickle.NpyArrayPlaceholder"] = build_nparray
  mt["scipy.sparse.csr.csr_matrix"] = sparse_matrix_reconstruct
  mt["__build__.Pickle.SpMatrixPlaceholder"] = build_spmatrix

  return Pickler{proto}(Memo(memo), PickleStack(), mt)
end

npyload(f) = load(NpyPickler(), f)

struct NpyArrayPlaceholder end

function np_multiarray_reconstruct(subtype, shape, dtype)
    @assert subtype.head == Symbol("numpy.ndarray")
    @assert isempty(subtype.args)
    @assert shape == (0,)
    @assert dtype == b"b" || dtype == "b"
    return NpyArrayPlaceholder()
end

struct SpMatrixPlaceholder end

sparse_matrix_reconstruct() = SpMatrixPlaceholder()

struct NpyDtype{T}
    little_endian::Bool
    dstring::String
    align::Bool
    copy::Bool
end

Base.eltype(::NpyDtype{T}) where T = T

function npy_typechar_to_jltype(t, n)
    if t in ("?", "b", "B")
        @assert isempty(n) || n == "1"
    end
    n = tryparse(Int, n)
    n = (isnothing(n) ? 4 : n) * 8
    @assert n in (8, 16, 32, 64, 128)

    if t == "?"
        return Bool
    elseif t == "b"
        return Int8
    elseif t == "B"
        return UInt8
    elseif t == "i"
        if n == 8
            return Int8
        elseif n == 16
            return Int16
        elseif n == 32
            return Int32
        elseif n == 64
            return Int64
        elseif n == 128
            return Int128
        else
            error("unsupport length $n for $t")
        end
    elseif t == "u"
        if n == 8
            return UInt8
        elseif n == 16
            return UInt16
        elseif n == 32
            return UInt32
        elseif n == 64
            return UInt64
        elseif n == 128
            return UInt128
        else
            error("unsupport length $n for $t")
        end
    elseif t == "f"
        if n == 16
            return Float16
        elseif n == 32
            return Float32
        elseif n == 64
            return Float64
        else
            error("unsupport length $n for $t")
        end
    else
        error("unsupport type $t")
    end
end

struct NString{N} end

slen(::NString{N}) where N = N

# TODO: handle dtype comprehensively
# https://github.com/numpy/numpy/blob/eeef9d4646103c3b1afd3085f1393f2b3f9575b2/numpy/core/src/multiarray/descriptor.c#L2442
function np_dtype(obj, align, copy)
    align, copy = Bool(align), Bool(copy)
    @assert !align "structure dtype disallow"
    @assert copy
    m = match(r"^([<=>])?([?bBiufU])(\d*)$", obj)
    if isnothing(m)
        @warn "unsupported dtype $obj: consider file an issue"
        return Defer(Symbol("numpy.dtype"), obj, align, copy)
    end

    ei, t, n = m.captures
    # '>': big, '<': little, '=': hardware-native
    islittle = ei == ">" ? false : ei == "<" ? true : islittle_endian()
    if t == "U"
        n = tryparse(Int, n)
        T = NString{isnothing(n) ? 1 : n}()
    else
        T = npy_typechar_to_jltype(t, n)
    end
    return NpyDtype{T}(islittle, obj, align, copy)
end

function build_npydtype(npydtype, state)
    metadata = length(state) == 9 ? state[9] : nothing
    ver, ei, sub_descrip, _names, _fields, elsize, alignment, flags = state

    T = eltype(npydtype)
    @assert isnothing(sub_descrip)
    @assert isnothing(_names)
    @assert isnothing(_fields)
    if T isa NString
        n = slen(T)
        @assert elsize == 4n
        @assert alignment == 4
        @assert flags == 8
    else
        @assert elsize == -1
        @assert alignment == -1
        @assert flags == 0
    end
    # '>': big, '<': little, '=': hardware-native
    islittle = ei == ">" ? false : ei == "<" ? true : islittle_endian()
    return NpyDtype{T}(islittle, npydtype.dstring, npydtype.align, npydtype.copy)
end

c2f(arr, shape) = PermutedDimsArray(reshape(arr, reverse(shape)), reverse(ntuple(identity, length(shape))))
c2f(arr, shape, n) = c2f(arr, (shape..., n))

function nstring(cs)
    i = findfirst(isequal('\0'), cs)
    return isnothing(i) ? String(cs) : String(@view(cs[1:i-1]))
end

# https://github.com/numpy/numpy/blob/6568c6b022e12ab6d71e7548314009ced6ccabe9/numpy/core/src/multiarray/methods.c#L1711
# TODO: support picklebuffer (Pickle v5)
function build_nparray(_, args)
    ver, shp, dtype, is_column_maj, data = args
    if dtype isa Defer
        return Defer(Symbol("build.nparray"), args)
    end
    T = eltype(dtype)

    data = data isa String ? codeunits(data) : data # old numpy use string instead of bytes

    if T isa NString
        n = slen(T)
        _data = reinterpret(UInt32, data)
        _arr = dtype.little_endian ? Char.(Base.ltoh.(_data)) : Char.(Base.ntoh.(_data))
        arr = is_column_maj ? reshape(_arr, n, shp...) : c2f(_arr, shp, n)
        return reshape(mapslices(nstring, arr; dims=ndims(arr)), shp)
    else
        _data = reinterpret(T, data)
        _arr = dtype.little_endian ? Base.ltoh.(_data) : Base.ntoh.(_data)
        arr = is_column_maj ? reshape(_arr, shp) : c2f(_arr, shp)
        return collect(arr)
    end
end

function np_scalar(dtype, data)
    T = eltype(dtype)
    if T isa NString
        n = slen(T)
        _data = reinterpret(UInt32, data)
        _arr = dtype.little_endian ? Char.(Base.ltoh.(_data)) : Char.(Base.ntoh.(_data))
        return String(_arr)
    else
        _data = reinterpret(T, data)
        _arr = dtype.little_endian ? Base.ltoh.(_data) : Base.ntoh.(_data)
        return _arr[]
    end
end

function build_spmatrix(_, args)
    shape = args["_shape"]
    nzval, colptr, rowval = csr_to_csc(shape..., args["data"], args["indptr"] .+ 1, args["indices"] .+ 1)
    return SparseArrays.SparseMatrixCSC(shape..., colptr, rowval, nzval)
end


text:
# This script times the results from using sparse arrays for caching; sparse Jacobians
# of μ, ξ, and 𝒱 ; and sparse Jacobians for calls to nlsolve.
using RiskAdjustedLinearizations, LinearAlgebra, SparseArrays
using BenchmarkTools, Test, SparseDiffTools

# Settings
define_functions = true
time_methods     = true
algorithm        = :relaxation # Note only relaxation works for sparse differentiation.
N_approx         = 10          # Number of periods ahead used for forward-difference equations

if define_functions
    include(joinpath(dirname(@__FILE__), "..", "nk_with_capital", "nk_with_capital.jl"))
end

# Set up

## Instantiate object
m_nk = NKCapital(; N_approx = N_approx) # create parameters
m = nk_capital(m_nk; sparse_arrays = true, sparse_jacobian = [:μ, :ξ])
zinit = deepcopy(m.z)
yinit = deepcopy(m.y)
Ψinit = deepcopy(m.Ψ)

## Solve for steady state once and update sparsity pattern
solve!(m; algorithm = algorithm, verbose = :none)

sparsity = Dict()
colorvec = Dict()
sparsity[:J𝒱] = sparse(m[:JV])
colorvec[:J𝒱] = isempty(sparsity[:J𝒱].nzval) ? ones(Int64, size(sparsity[:J𝒱], 2)) : matrix_colors(sparsity[:J𝒱])
update_sparsity_pattern!(m, [:𝒱]; sparsity = sparsity, colorvec = colorvec)

## Solve w/sparse array caching; sparse differentiation of Jacobians of
## μ, ξ, and 𝒱 ; and sparse differentiation of the objective functions in `nlsolve`
jac_cache = preallocate_jac_cache(m, algorithm)
update!(m, zinit, yinit, Ψinit)
solve!(m; algorithm = algorithm, sparse_jacobian = true, jac_cache = jac_cache)

if time_methods
    m_dense = nk_capital(m_nk)

    @info "Timing solve! with varying degrees of sparsiy"

    println("Dense Array Caches and Dense Jacobians")
    @btime begin
        update!(m_dense, zinit, yinit, Ψinit)
        solve!(m_dense; algorithm = algorithm, verbose = :none)
    end
    # ~ 2.48 s

    println("Sparse Array Caches and Sparse Jacobians for Equilibrium Functions")
    @btime begin
        update!(m, zinit, yinit, Ψinit)
        solve!(m; algorithm = algorithm, verbose = :none)
    end
    # ~ 2.37 s

    println("Sparse Jacobians for nlsolve")
    @btime begin
        update!(m_dense, zinit, yinit, Ψinit)
        solve!(m_dense; algorithm = algorithm, sparse_jacobian = true,
               jac_cache = jac_cache, verbose = :none)
    end
    # ~ 0.85 s

    println("Sparse Array Caches, Sparse Jacobians for Equilibrium Functions, and Sparse Jacobians for nlsolve")
    @btime begin
        update!(m, zinit, yinit, Ψinit)
        solve!(m; algorithm = algorithm, sparse_jacobian = true, jac_cache = jac_cache, verbose = :none)
    end
    # ~ 0.9s

    @test m_dense.z ≈ m.z
    @test m_dense.y ≈ m.y
    @test m_dense.Ψ ≈ m.Ψ
end


text:
using JLD2
export read_JLD22Dict, save_CorrData2JLD2, save_Dict2JLD2, save_Array2JLD2

"""
    read_JLD22Dict(inFile::Array{JLD2.JLDFile,1}, tstamplist::AbstractArray, times::Array{Int64,1})

    read data structures from JLD2 into a list of dictionaries
"""
function read_JLD22Dict(inFile::Array{JLD2.JLDFile{JLD2.MmapIO},1}, tstamplist::AbstractArray)
    # create a list of dictionaries containing input data for each time step
    # this assumes inFile is in the same order as tstamplist
    dsets = [Dict("$ts/$stn" => @views inFile[i]["$ts/$stn"] for stn in keys(inFile[i][ts])) for (i, ts) in enumerate(tstamplist)]
    return dsets
end

"""
    read_JLD22Dict(inFile, tstamplist::AbstractArray, times::Array{Int64,1})

    read data structures from JLD2 into a list of dictionaries
"""
function read_JLD22Dict(inFile::JLD2.JLDFile, tstamplist::AbstractArray)
    # create a list of dictionaries containing input data for each time step
    dsets = [Dict("$ts/$stn" => inFile["$ts/$stn"] for stn in keys(inFile[ts])) for ts in tstamplist]
    return dsets
end

"""
    save_CorrData2JLD2(foname::String, varname::String, CD::CorrData)

    save CorrData structure to JLD2
"""
function save_CorrData2JLD2(foname::String, varname::String, CD::CorrData)
    file = jldopen(foname, "a+")
    file[varname] = CD
    JLD2.close(file)
end

"""
    save_Dict2JLD2(foname::String, varname::String, D::Dict)

    save Dictionary to JLD2
"""
function save_Dict2JLD2(foname::String, varname::String, D::Dict)
    file = jldopen(foname, "r+")
    file[varname] = D
    JLD2.close(file)
end

"""
    save_Dict2JLD2(foname::String, varname::String, D::Dict)

    save Dictionary to JLD2
"""
function save_Array2JLD2(foname::String, varname::String, A::AbstractArray)
    file = jldopen(foname, "a+")
    file[varname] = A
    JLD2.close(file)
end

max_stars_repo_name: kura-okubo/SeisXcorrelation

max_issues_repo_name: kura-okubo/SeisXcorrelation

max_forks_repo_name: kura-okubo/SeisXcorrelation


text:
main_folder = "$(dirname(@__FILE__))/.."

include("$main_folder/src/functions/utils/flatten.jl")
include("$main_folder/src/functions/utils/delete_method.jl")
include("$main_folder/src/functions/utils/get_all_files.jl")
include("$main_folder/src/functions/utils/get_all_symbols.jl")
include("$main_folder/src/functions/utils/include_all_files.jl")
include("$main_folder/src/macros/export_all_files.jl")

loaded_folders = [ "config/initializers" , "src" , "lib/tasks" ]

for loaded_folder in loaded_folders
  loaded_folder = "$main_folder/$loaded_folder"

  cur_package_name = "Julz"
  Julz.include_all_files(loaded_folder, package_name=cur_package_name)
  @eval Julz.@export_all_files $loaded_folder $cur_package_name
end

function load_input(raw_input, is_file_input=false, can_be_missing=false, is_init_load=false)
  if is_file_input
    file_path = "$main_folder/$raw_input"

    if !isfile(file_path)
      file_path = "$main_folder/lib/input_decks/$raw_input"
    end

    if !isfile(file_path)
      if can_be_missing ; return ; end
    end

    open(file_path) do file
      file_lines = map(x -> split(x, "#")[1], readlines(file))
      file_lines = map(x -> strip(x), file_lines)

      raw_input = join(file_lines, ";")

      botched_chars = [ "(" , "," , "[" ]

      for cur_char in botched_chars
        raw_input = replace(raw_input, "$cur_char;", "$cur_char")
      end

      raw_input *= "; nothing ;"
    end
  end

  old_symbols = get_all_symbols(Julz)
  eval( parse(raw_input) )
  new_symbols = setdiff(get_all_symbols(Julz), old_symbols)

  is_init_load && ( new_symbols = [] )

  for cur_symbol in setdiff(new_symbols, old_symbols)
    cur_expression = quote
      $(cur_symbol) = $(getfield(Julz, cur_symbol))
    end

    eval(Main, cur_expression)
  end
end

print_workspace() = Julz.print_workspace(Julz)

defaults_file_name = "defaults.jl"
load_input(defaults_file_name, true)

scratch_file_name = "tmp/scratch.jl"
Julz.load_input(scratch_file_name, true, true, true)


text:
module node_partitioning_examples
using FinEtools
using PyCall

function node_partitioning_1()

    H = 100. # strip width
    a = 10. # radius of the hole
    L = 200. # length of the strip
    nL = 15
    nH = 10
    nR = 50
    fens,fes = Q4elliphole(a, a, L/2, H/2, nL, nH, nR)

    npartitions = 4
    partitioning = nodepartitioning(fens, npartitions)
    partitionnumbers = unique(partitioning)

    @pyimport matplotlib.pyplot as plt
    fig = plt.figure()
    ax = plt.axes()
    for ixxxx = 1:length(partitionnumbers)
        i1 = find(x -> x == partitionnumbers[ixxxx], partitioning)
        ax[:plot](vec(fens.xyz[i1, 1]), vec(fens.xyz[i1, 2]), linestyle="none", marker=:o)
    end
    plt.show()

end # node_partitioning_1

function node_partitioning_2()
    a = 10. # radius of the hole
    nC = 20
    nR = 4
    fens,fes = Q4annulus(a, 1.5*a, nR, nC, 2*pi)

    npartitions = 8
    partitioning = nodepartitioning(fens, npartitions)
    partitionnumbers = unique(partitioning)

    @pyimport matplotlib.pyplot as plt
    plt.style[:use]("seaborn-whitegrid")
    fig = plt.figure()
    ax = plt.axes()
    for ixxxx = 1:length(partitionnumbers)
        i1 = find(x -> x == partitionnumbers[ixxxx], partitioning)
        ax[:plot](vec(fens.xyz[i1, 1]), vec(fens.xyz[i1, 2]), linestyle="none", marker=:o)
    end
    ax[:set_xlabel]("x")
    ax[:set_ylabel]("y")
    plt.axis("equal")
    plt.show()

end # node_partitioning_2

function node_partitioning_3()
    H = 100. # strip width
    a = 10. # radius of the hole
    L = 200. # length of the strip
    nL = 15
    nH = 10
    nR = 50
    fens,fes = Q4elliphole(a, a, L/2, H/2, nL, nH, nR)

    npartitions = 4
    partitioning = nodepartitioning(fens, npartitions)
    partitionnumbers = unique(partitioning)

    # Render the  points as dots ("Points" in paraview)
    Points = FESetP1(reshape(collect(1:count(fens)), count(fens), 1))
    File = "Q4elliphole_mesh.vtk"
    vtkexportmesh(File, fens, Points;
    scalars=[("Partition", vec(partitioning))])
    @async run(`"paraview.exe" $File`)

end # node_partitioning_3

function node_partitioning_4()
    H = 30. # strip width
    R = 10. # radius of the hole
    L = 20. # length of the strip
    nL = 15
    nH = 10
    nR = 5
    fens,fes = H8block(L, H, R, nL, nH, nR)

    npartitions = 16
    partitioning = nodepartitioning(fens, npartitions)
    partitionnumbers = unique(partitioning)
    @test norm(sort(partitionnumbers) - sort(1:npartitions)) < 1.e-5

    # Visualize partitioning
    for gp = partitionnumbers
      groupnodes = findall(k -> k == gp, partitioning)
      File =  "partition-nodes-$(gp).vtk"
      vtkexportmesh(File, fens, FESetP1(reshape(groupnodes, length(groupnodes), 1)))
    end
    File =  "partition-mesh.vtk"
    vtkexportmesh(File, fens, fes)
    @async run(`"paraview.exe" $File`)
end

function allrun()
    println("#####################################################")
    println("# node_partitioning_1 ")
    node_partitioning_1()
    println("#####################################################")
    println("# node_partitioning_2 ")
    node_partitioning_2()
    println("#####################################################")
    println("# node_partitioning_3 ")
    node_partitioning_3()
    println("#####################################################")
    println("# node_partitioning_4 ")
    node_partitioning_4()
    return true
end # function allrun

end # module node_partitioning_examples


text:
using FEMSolids
using Test

@testset "FEMSolids.jl" begin
    include("test_element_routine.jl")
    include("test_adaptivity.jl")
end


text:
"""
	grdvolume(cmd0::String="", arg1=[], kwargs...)

Reads one 2-D grid and returns xyz-triplets.

Full option list at [`grdvolume`](http://gmt.soest.hawaii.edu/doc/latest/grdvolume.html)

Parameters
----------

- **C** : **contour** : -- Str or List --   Flags = cval or low/high/delta or rlow/high or rcval

    Find area, volume and mean height (volume/area) inside the cval contour.
    [`-C`](http://gmt.soest.hawaii.edu/doc/latest/grdvolume.html#c)
- **L** : **base_level** : -- Number --           Flags = base

    Also add in the volume from the level of the contour down to base [Default base is contour].
    [`-L`](http://gmt.soest.hawaii.edu/doc/latest/grdvolume.html#l)
- $(GMT.opt_R)
- **S** : **unit** : -- Str --              Flags = e|f|k|M|n|u

    For geographical grids, append a unit from e|f|k|M|n|u [Default is meter (e)].
    [`-S`](http://gmt.soest.hawaii.edu/doc/latest/grdvolume.html#s)
- **T** : -- Str --                         Flags = [c|h]

    Determine the single contour that maximized the average height (= volume/area).
    [`-T`](http://gmt.soest.hawaii.edu/doc/latest/grdvolume.html#t)
- $(GMT.opt_V)
- **Z** : **scale** : -- Str or List --     Flags = fact[/shift]

    Optionally subtract shift before scaling data by fact. [Default is no scaling].
    [`-Z`](http://gmt.soest.hawaii.edu/doc/latest/grdvolume.html#z)
- $(GMT.opt_f)
- $(GMT.opt_o)
"""
function grdvolume(cmd0::String="", arg1=[]; kwargs...)

	length(kwargs) == 0 && occursin(" -", cmd0) && return monolitic("grdvolume", cmd0, arg1)	# Speedy mode

	d = KW(kwargs)

	cmd, = parse_R("", d)
	cmd = parse_V_params(cmd, d)
	cmd, = parse_f(cmd, d)
	cmd, = parse_o(cmd, d)

	cmd = add_opt(cmd, 'C', d, [:C :contour])
	cmd = add_opt(cmd, 'L', d, [:L :base_level])
	cmd = add_opt(cmd, 'S', d, [:S :unit])
	cmd = add_opt(cmd, 'T', d, [:T])
	cmd = add_opt(cmd, 'Z', d, [:Z :scale])

	cmd, got_fname, arg1 = find_data(d, cmd0, cmd, 1, arg1)
	return common_grd(d, cmd, got_fname, 1, "grdvolume", arg1)		# Finish build cmd and run it
end

# ---------------------------------------------------------------------------------------------------
grdvolume(arg1=[], cmd0::String=""; kw...) = grdvolume(cmd0, arg1; kw...)


text:
using SnakeAI
using Test

# @testset "SnakeAI.jl" begin
#     # Write your tests here.
# end

@testset "Points.jl" begin

    p1 = Point(100, 200)
    p2 = Point(300, 400)
    p3 = Point(300, 400)

    @test p1.x == 100 && p1.y == 200
    @test p1 != p2
    @test p2 == p3

end


text:
###############################################################################
# GRAPHICS
###############################################################################
"""
Cartesian Coordinate.

$(TYPEDEF)
$(TYPEDFIELDS)
"""
struct Coordinate{T <: Number}
    x::T
    y::T
end
Coordinate() = Coordinate{Int}(0, 0)
Coordinate(x::T) where T<:Number = Coordinate{T}(x, zero(x))
Base.eltype(c::Coordinate{T}) where {T<:Number} = T

#-------------------
"""
Fill attributes as strings.

$(TYPEDEF)
$(TYPEDFIELDS)
"""
struct Fill
    color::Maybe{String}
    image::Maybe{String}
    gradient_color::Maybe{String}
    gradient_rotation::Maybe{String}
end
function Fill(; color=nothing, image=nothing, gradient_color=nothing, gradient_rotation=nothing)
    Fill(color, image, gradient_color, gradient_rotation )
end

#-------------------
"""
Font attributes as strings.

$(TYPEDEF)
$(TYPEDFIELDS)
"""
struct Font
    family    ::Maybe{String}
    style     ::Maybe{String}
    weight    ::Maybe{String}
    size      ::Maybe{String}
    align     ::Maybe{String}
    rotation  ::Maybe{String}
    decoration::Maybe{String}
end
function Font(; family=nothing, style=nothing, weight=nothing,
              size=nothing, align=nothing, rotation=nothing, decoration=nothing)
    Font(family, style, weight, size, align, rotation, decoration)
end

#-------------------
"""
Line attributes as strings.

$(TYPEDEF)
$(TYPEDFIELDS)
"""
struct Line
    color::Maybe{String}
    shape::Maybe{String}
    style::Maybe{String}
    width::Maybe{String}
end
function Line(; color=nothing, shape=nothing, style=nothing, width=nothing)
    Line(color, shape, style, width)
end

#-------------------
"""
PNML Graphics can be attached to many parts of PNML models.

$(TYPEDEF)
$(TYPEDFIELDS)
"""
struct Graphics #{COORD,FILL,FONT,LINE}
    dimension::Maybe{Coordinate}
    fill::Maybe{Fill}
    font::Maybe{Font}
    line::Maybe{Line}
    offset::Maybe{Coordinate}
    position::Maybe{Vector{Coordinate}}
end

function Graphics(;dim=nothing, fill=nothing, font=nothing,
                  line=nothing, offset=nothing, position=nothing)
    Graphics(dim, fill, font, line, offset, position)
end


text:
"""
    PurifiedState{ST<:DualValuedState} <: DualValuedComposedState

State composed by two identical states ST, used to represent states of a
density matrix.
"""
mutable struct PurifiedState{ST<:FiniteBasisState} <: FiniteBasisState
    n_sys::Int # number of states
    n_add::Int # number of states
    σ_sys::ST # left
    σ_add::ST # right

    space_dim::Int
end

## Constructors
PurifiedState(sys::ST, add::ST) where {ST<:FiniteBasisState} =
    PurifiedState(nsites(sys),   nsites(add),
                  deepcopy(sys), deepcopy(add),
                  spacedimension(sys)*spacedimension(add))

PurifiedState{ST}(n, add, i_σ=0) where ST =
    set!(PurifiedState(ST(n, 0), ST(add, 0)), i_σ)

# Property Accessors
@inline spacedimension(state::PurifiedState) = state.space_dim
@inline nsites(state::PurifiedState) = nsites(state.σ_sys) + nsites(state.σ_add)
@inline toint(state::PurifiedState) = _toint(add(state), sys(state)) #toint(state.σ_row, state.σ_col)
@inline index(state::PurifiedState) = index(sys(state))
@inline flipped(a::PurifiedState, b::PurifiedState) = (flipped(a.σ_sys, b.σ_sys), flipped(a.σ_add, b.σ_add))
@inline index_to_int(state::PurifiedState, id) = @error("Undefined: index is on sys, add is undefined.")
@inline index_to_int(state::PurifiedState, id, add::FiniteBasisState) =
        index_to_int(state, id, toint(add))
@inline index_to_int(state::PurifiedState, id, add) = add*spacedimension(sys(state)) + index_to_int(sys(state), id)
@inline eltype(state::PurifiedState) = eltype(sys(state))

# custom accessor
sys(v::PurifiedState) = v.σ_sys
add(v::PurifiedState) = v.σ_add
config(v::PurifiedState) = (config(v.σ_sys), config(v.σ_add))

# checs
same_basis(v::PurifiedState{ST}, v2::PurifiedState{ST2}) where {ST, ST2} =
    ST==ST2 && same_basis(sys(v), sys(v2)) && same_basis(add(v), add(v2))

function flipat!(rng::AbstractRNG, v::PurifiedState, i::Int)
    i > v.n_sys ? flipat!(v.σ_add, i-v.n_sys) : flipat!(v.σ_sys, i)
end

function setat!(v::PurifiedState, i::Int, val)
    i > v.n_sys ? setat!(v.σ_add, i-v.n_sys, val) : setat!(v.σ_sys, i, val)
end

set_index!(v::PurifiedState, i::Integer, a) = set!(v, index_to_int(v, i, a))
function set!(v::PurifiedState, i::Integer)
    i_add, i_sys = divrem(i, spacedimension(v.σ_sys))
    set!(v.σ_sys, i_sys)
    set!(v.σ_add, i_add)

    v
end
set!(v::PurifiedState, sys, add) = set!(sys(v), row) && set!(add(v), col) && v

function rand!(rng::AbstractRNG, state::PurifiedState)
    rand!(rng, state.σ_sys)
    rand!(rng, state.σ_add)
end

set_sys!(state::PurifiedState, i::Integer) = set!(state.σ_sys, i)
set_add!(state::PurifiedState, i::Integer) = set!(state.σ_add, i)

# -------------- Base.show extension for nice printing -------------- #
Base.show(io::IO, mm::MIME"text/plain", bs::PurifiedState) =
    print(io, "PurifiedState($(bs.n_sys), $(bs.n_add)) : \n\tRow: ", bs.σ_sys,"\n\tCol: ", bs.σ_add, "\n")

Base.show(io::IO, bs::PurifiedState) = print(io, "(",bs.σ_sys,", ", bs.σ_add, ")")


text: const ASSET_FINGERPRINT = "968fc8acb3cb5030a26f"

max_stars_repo_name: tetsugps/julia

max_issues_repo_name: tetsugps/julia

max_forks_repo_name: tetsugps/julia


text:
using MOOSE

function doIt()
    startRootLog!(main_perf_log)

    # Create the Mesh
    mesh = buildSquare(0, 1, 0, 1, 500, 500)

    # Create the System to hold the equations
    diffusion_system = System{Float64}(mesh)

    # Add a variable to solve for
    u = addVariable!(diffusion_system, "u")

    # Apply the Laplacian operator to the variable
    addKernel!(diffusion_system, Diffusion(u))

    # Add in a Convection operator with a velocity vector
    addKernel!(diffusion_system, Convection(u, ContMechTensors.Vec{2}((10.,0.))))

    # u = 0 on the Left
    addBC!(diffusion_system, DirichletBC(u, [4], 0.0))

    # u = 1 on the Right
    addBC!(diffusion_system, DirichletBC(u, [2], 1.0))

    # Initialize the system of equations
    initialize!(diffusion_system)

    # Create a solver and solve
    solver = PetscImplicitSolver(diffusion_system)

    # Initialize the solver so that we can call it multiple times
    initialize!(solver)

    solve!(solver)

    # Output
    out = VTKOutput()
    output(out, solver, "convection_diffusion_out")

    stopRootLog!(main_perf_log)
end

doIt()
clear!(main_perf_log)
doIt()


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

using JuliaFEM
using JuliaFEM.Testing

#=
In [36]: C = Matrix([[0], [30], [15]]) # node coordinates
In [37]: A = Matrix([P.subs({x: C[i,0]}).T for i in range(len(P))])
In [38]: N = P.T*A.inv()
In [39]: Me = integrate(N.T*N, (x, 0, 30))
In [40]: De = diag(*integrate(N, (x, 0, 30)))
In [41]: Me
Out[41]:
Matrix([
[ 4, -1,  2],
[-1,  4,  2],
[ 2,  2, 16]])
In [42]: De
Out[42]:
Matrix([
[5, 0,  0],
[0, 5,  0],
[0, 0, 20]])
=#

@testset "dirichlet problem in 1 dimension" begin
    element = Element(Seg2, [1, 2])
    element["geometry"] = Vector{Float64}[[0.0, 0.0], [6.0, 0.0]]
    element["temperature 1"] = 0.0
    p1 = Problem(Dirichlet, "test problem 1", 1, "temperature")
    p1.properties.dual_basis = false
    p2 = Problem(Dirichlet, "test problem 2", 1, "temperature")
    p2.properties.dual_basis = true
    assemble!(p1, element)
    assemble!(p2, element)
    C1 = full(p1.assembly.C1)
    C2 = full(p1.assembly.C2)
    @test isapprox(C1, C2)
    @test isapprox(C1, [2.0 1.0; 1.0 2.0])
    C1 = full(p2.assembly.C1)
    C2 = full(p2.assembly.C2)
    @test isapprox(C1, C2)
    @test isapprox(C1, [3.0 0.0; 0.0 3.0])

    element = Element(Seg3, [1, 2, 3])
    element["geometry"] = Vector{Float64}[[0.0, 0.0], [30.0, 0.0], [15.0, 0.0]]
    element["temperature 1"] = 0.0
    p1 = Problem(Dirichlet, "quadratic 1", 1, "temperature")
    p1.properties.dual_basis = false
    p2 = Problem(Dirichlet, "quadratic 1", 1, "temperature")
    p2.properties.dual_basis = true
    assemble!(p1, element)
    assemble!(p2, element)
    C1 = full(p1.assembly.C1)
    C2 = full(p1.assembly.C2)
    @test isapprox(C1, C2)
    @test isapprox(C1, [4.0 -1.0 2.0; -1.0 4.0 2.0; 2.0 2.0 16.0])
    C1 = full(p2.assembly.C1)
    C2 = full(p2.assembly.C2)
    @test isapprox(C1, C2)
    @test isapprox(C1, [5.0 0.0 0.0; 0.0 5.0 0.0; 0.0 0.0 20.0])
end

#=
@testset "dirichlet problem using tri3 surface element" begin
    element = Tri3([1, 2, 3])
    element["geometry"] = Node[[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0]]
    element["temperature"] = 0.0
    problem = Problem(Dirichlet, "test problem", 1, "temperature")
    push!(problem, element)
    assemble!(problem, 0.0)
    C1 = full(problem.assembly.C1)
    C2 = full(problem.assembly.C2)
    @test isapprox(C1, C2)
    @test isapprox(C1, 1/24*[2 1 1; 1 2 1; 1 1 2])
end

@testset "dirichlet problem in 2 dimensions" begin
    element = Seg2([1, 2])
    element["geometry"] = Node[[1.0, 1.0], [0.0, 1.0]]
    element["displacement 1"] = 0.0
    element["displacement 2"] = 0.0
    problem = Problem(Dirichlet, "test problem", 2, "displacement")
    push!(problem, element)
    assemble!(problem, 0.0)
    C1 = full(problem.assembly.C1)
    C2 = full(problem.assembly.C2)
    g = full(problem.assembly.g)
    @test isapprox(C1, C2)
    C1_expected = 1/6*[2 0 1 0; 0 2 0 1; 1 0 2 0; 0 1 0 2]
    @test isapprox(C1, C1_expected)
    @test isapprox(g, [0.0, 0.0, 0.0, 0.0])
end

@testset "dirichlet problem in 2 dimensions, with 1 dof fixed" begin
    element = Seg2([1, 2])
    element["geometry"] = Node[[1.0, 1.0], [0.0, 1.0]]
    element["displacement 2"] = 0.0
    problem = Problem(Dirichlet, "test problem", 2, "displacement")
    push!(problem, element)
    assemble!(problem, 0.0)
    C1 = full(problem.assembly.C1)
    C2 = full(problem.assembly.C2)
    g = full(problem.assembly.g)
    @test isapprox(C1, C2)
    C1_expected = 1/6*[
        0 0 0 0
        0 2 0 1
        0 0 0 0
        0 1 0 2]
    @test isapprox(C1, C1_expected)
    @test isapprox(g, [0.0, 0.0, 0.0, 0.0])
end
=#

@testset "test analytical boundary condition" begin
    X = Dict{Int64, Vector{Float64}}(
        1 => [0.0, 0.0],
        2 => [1.0, 0.0])
    element = Element(Seg2, [1, 2])
    update!(element, "geometry", X)
    update!(element, "displacement 1", 0.0)
    f(xi, time) = begin
        info("function called at xi = $xi, time = $time")
        X = element("geometry", xi, time)
        info("geometry at xi, X = $X")
        val = X[1]*time
        info("result for field at xi = $val")
        return val
    end
    update!(element, "displacement 2", f)
    p = Problem(Dirichlet, "test boundary", 2, "displacement")
    push!(p, element)
    assemble!(p, 0.0)
    g1 = full(p.assembly.g, 4, 1)
    @test isapprox(g1, [0.0, 0.0, 0.0, 0.0])
    empty!(p.assembly)
    assemble!(p, 1.0)
    g2 = full(p.assembly.g, 4, 1)
    C2 = full(p.assembly.C2, 4, 4)
    u = C2 \ g2
    info("u = $u")
    @test isapprox(u, [0.0, 0.0, 0.0, 1.0])
end


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, Pkg

name = "prrte"
version = v"2.0.0"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/openpmix/prrte/releases/download/v$(version)/prte-$(version).tar.bz2", "9f4abc0b1410e0fa74ed7b00cfea496aa06172e12433c6f2864d11b534becc25")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd prte-*
./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target} \
    --enable-shared \
    --with-libevent=${prefix} \
    --with-hwloc=${prefix} \
    --with-pmix=${prefix} \
    --without-tests-examples \
    --disable-man-pages
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Platform("i686", "linux"; libc = "glibc"),
    Platform("x86_64", "linux"; libc = "glibc"),
    Platform("aarch64", "linux"; libc = "glibc"),
    Platform("armv7l", "linux"; call_abi = "eabihf", libc = "glibc"),
    Platform("powerpc64le", "linux"; libc = "glibc"),
    Platform("i686", "linux"; libc = "musl"),
    Platform("x86_64", "linux"; libc = "musl"),
    Platform("aarch64", "linux"; libc = "musl"),
    Platform("armv7l", "linux"; call_abi = "eabihf", libc = "musl"),
    Platform("x86_64", "macos"; )
]

# The products that we will ensure are always built
products = [
    ExecutableProduct("prte", :prte)
    ExecutableProduct("prun", :prun)
    ExecutableProduct("prte_info", :prte_info)
    ExecutableProduct("prterun", :prterun)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="libevent_jll", uuid="1080aeaf-3a6a-583e-a51c-c537b09f60ec"))
    Dependency(PackageSpec(name="Hwloc_jll", uuid="e33a78d0-f292-5ffc-b300-72abe9b543c8"))
    Dependency(PackageSpec(name="PMIx_jll", uuid="32165bc3-0280-59bc-8c0b-c33b6203efab"), compat="4.1.0")
]

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

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
## gradlogpdf ##

# TODO incorrect, NO_FIELDS is wrong
# function CRC.rrule(::typeof(logpdf), d::UnivariateDistribution, x)
    # v = logpdf(d, x)
    # pullback(y) = (NO_FIELDS, d. @thunk(gradlogpdf(d, x)))
    # return (v, pullback)
# end

function CRC.rrule(::Type{D}, θ) where {D <: Exponential}
    function D_pullback(dy)
        return (NO_FIELDS, dy.θ)
    end
    return D(θ), D_pullback
end

function CRC.rrule(::Type{D}, a, b) where {D <: Union{Uniform, Arcsine}}
    function D_pullback(dy)
        return (NO_FIELDS, dy.a, dy.b)
    end
    return D(a, b), D_pullback
end

function CRC.rrule(::Type{D}, α, β) where {D <: Union{Beta, BetaPrime}}
    function D_pullback(dy)
        return (NO_FIELDS, dy.α, dy.β)
    end
    return D(α, β), D_pullback
end

function CRC.rrule(::Type{D}, µ, σ) where {D <: Union{Normal, LogNormal}}
    function D_pullback(dy)
        return (NO_FIELDS, dy.µ, dy.σ)
    end
    return D(µ, σ), D_pullback
end

function CRC.rrule(::Type{<:Weibull}, α, θ)
    function weibull_pullback(dy)
        return (NO_FIELDS, dy.α, dy.θ)
    end
    return Weibull(α, θ), weibull_pullback
end


text:
"""
	map function for parallelized processing of PGV and PGA

2021.1.4 Kurama Okubo
"""

using CSV, DataFrames

"""
	map_seisPGVPGA(startid, InputDict::Dict)

Download 3 components seismogram and compute PGV, PGA and noise strength d_rms.
"""
function map_seisPGVPGA(startid, InputDict::Dict)

	download_time_unit = InputDict["download_time_unit"]
	fodir              = InputDict["fodir"]
	starttimelist = InputDict["starttimelist"]

	#SeisIO getdata option
	!haskey(InputDict, "get_data_opt") && (InputDict["get_data_opt"] = [true, true, true, true, true])#: [unscale, demean, detrend, taper, ungap]
	!haskey(InputDict, "sampling_frequency") && (InputDict["sampling_frequency"] = false) # default margin: 5 minutes
	!haskey(InputDict, "download_margin") && (InputDict["download_margin"] = 5 * 60) # default margin: 5 minutes

	#show progress
	if starttimelist[startid][end-8:end] == "T00:00:00"
	    println("start downloading $(starttimelist[startid])")
	end

	request_src_chanks = InputDict["requeststations"]

	for src in keys(request_src_chanks)
	    # replace source name
	    if src=="IRISDMC"
			reqsrc="IRIS"
		else
			reqsrc=src
		end

	    #---download data---#
	    requeststrs = get_requeststr(request_src_chanks[src], 1) #  request single station by station

	    for requeststr in requeststrs

    		#---Skip if this request station has been already performed due to restart processing etc.---#
    		s_str = string(DateTime(starttimelist[startid]))[1:19]
	        # compute end time:
	        et = DateTime(starttimelist[startid]) + Second(download_time_unit)
	        e_str=string(et)[1:19]
	        # split network, station, location and channel
	        net, sta, loc, cha = split(requeststr[1], ".")
	        varname = join([net, sta], ".")*"__"*s_str*"__"*e_str

	        # make hierarchical directory
	        dir1 = joinpath(InputDict["pgvpga_dir"], join([net, sta], ".")) #ex. BP.LCCB..BP1
	        ts_year = split(split(varname, "__")[2], "-")[1]# start year
	        dir2 = joinpath(dir1, ts_year) #ex. 2014
	        mkpath(dir2)
	        # dump file
	        if ispath(joinpath(dir2, varname*".dat"))
	        	println("skip $(varname*".dat").")
	        	continue;
	        end

	        #NOTE: method is currently fixed "FDSN", which downloads waveforms.
	        method = "FDSN"

	        # including download margin
	        starttime = string(DateTime(starttimelist[startid]) - Second(InputDict["download_margin"]))
	        dltime = download_time_unit + 2 * InputDict["download_margin"]

	        ex = :(get_data($(method), $(requeststr), s=$(starttime), t=$(dltime),
	         v=$(0), src=$(reqsrc), unscale=$(InputDict["get_data_opt"][1]),
	          demean=$(InputDict["get_data_opt"][2]), detrend=$(InputDict["get_data_opt"][3]),taper=$(InputDict["get_data_opt"][4]),
	          ungap=$(InputDict["get_data_opt"][5]), rr=true))

	        t_dl = @elapsed Stemp = check_and_get_data(ex)
	        # If download is failed, skip
	        Stemp == 1 && continue
	        # If Stemp does not have 3 components, skip
	        Stemp.n != 3 && continue

	        # manipulate download_margin
	        # NOTE: using SeisIO.sync()
	        stsync = DateTime(starttimelist[startid])
	        etsync = u2d(d2u(stsync) + download_time_unit)
	        SeisIO.sync!(Stemp, s=stsync, t=etsync, v=0)

	        #---downsampling---#
	        if InputDict["sampling_frequency"] isa Number
	            for j = 1:Stemp.n
	                if Stemp.misc[j]["dlerror"] == 0 && !isempty(Stemp[j].t)
	                    # downsample
	                    if Stemp.fs[j] > InputDict["sampling_frequency"]
	                        SeisIO.resample!(Stemp, chans=j, fs=float(InputDict["sampling_frequency"]))
	                    end
	                end
	            end
	        end

	        #---compute pgv, pga and d_rms---#
	        pgv, pga, drms_mean = compute_pgvpga_drms(Stemp, InputDict["bp_fmin"], InputDict["bp_fmax"],
	                                                    InputDict["drms_fmin"], InputDict["drms_fmax"])

	        #---output result into dat file---#
	        # dump file
	        open(joinpath(dir2, varname*".dat"), "w") do io
	            write(io, "id,starttime,endtime,PGV,PGA,drms_mean\n")
	            write(io, "$(join([net, sta], ".")), $(s_str), $(e_str), $(pgv), $(pga), $(drms_mean)")
	        end;

	        #---save seisdata---#
	        if InputDict["save_seisdata"]
	            # make hierarchical directory
	            dir1 = joinpath(InputDict["seisdata_dir"], join([net, sta], ".")) #ex. BP.LCCB..BP1
	            ts_year = split(split(varname, "__")[2], "-")[1]# start year
	            dir2 = joinpath(dir1, ts_year) #ex. 2014
	            mkpath(dir2)
	            # dump file
	            wseis(joinpath(dir2, varname*".seisio"), Stemp)
	        end
	    end
	end

end

"""
    map_readPGVPGA(fi_csv::String)

read pgv, pga and d_rms from dat file.
"""
function map_readPGVPGA(fi_csv::String)
    df = CSV.File(fi_csv,  header=1, delim=',') |> DataFrame
    df["midtime"] = u2d((d2u(df["starttime"][1]) + d2u(df["endtime"][1]))/2)
    return df
end


text:
#=
Functionalities to manage files in the context of code evaluation.

- script files (can be read and written)
- output files (written, = where stdout is redirected)
=#

"""
$(SIGNATURES)

Internal function to take a relative path and return a unix version of the path
(if it isn't already). Used in [`resolve_rpath`](@ref).
"""
function unixify(rpath::AS)::AS
    # if empty, return "/"
    isempty(rpath) && return "/"
    # if windows, replace "\\" by "/"
    Sys.isunix() || (rpath = replace(rpath, "\\" => "/"))
    # if it has an extension e.g.: /blah.txt, return
    isempty(splitext(rpath)[2]) || return rpath
    # if it doesn't have an extension, check if it ends with `/` e.g. : /blah/
    # if it doesn't end with "/", add one and return
    endswith(rpath, "/") || return rpath * "/"
    return rpath
end

"""
$(SIGNATURES)

Internal function to take a unix path, split it along `/` and re-join it
(which will lead to the same path on unix but not on Windows).
This undoes what [`unixify`](@ref) does.
"""
join_rpath(rpath::AS) = joinpath(split(rpath, '/')...)

"""
$(SIGNATURES)

Internal function to parse a relative path and return either a path relative
to the website root folder or the website assets folder or a full system path.
This doesn't check whether the  path points to  an existing asset (for that
see [`resolve_rpath`](@ref).

* `/[path]`:  this is a path relative to the website root folder, e.g.:
              `/foo/figs/img1.png`.
* `[path]`:   this is a path relative to the site assets folder, e.g.:
              `/assets/foo`.
* `./[path]`: this is a path relative to the assets folder with the same path
              as the calling page. For instance, `./im1.png` in
              `src/pages/pg1.md` will point to `/assets/pages/pg1/im1.png`.

Otherwise, if in `code` mode, recurse with `./code/path`.

## Argument

1. `rpath`: the relative path (always in Unix format)

## Keywords

* `canonical=false`: whether to return a full path on the system (`true`) or
                     a unix-path relative to the website root folder (`false`).
* `code=false`:      whether we're in a code context or not (in which case an
                     additional shortcut form is allowed).
"""
function parse_rpath(rpath::AS; canonical::Bool=false, code::Bool=false)::AS

    FS2 = FD_ENV[:STRUCTURE] >= v"0.2"

    # path from the website root folder
    # examples with FS2:
    # /foo.bar          ==> { /foo.bar        ; [__site]/foo.bar        }
    # /_assets/foo.bar  ==> { /assets/foo.bar ; [__site]/assets/foo.bar }
    #
    if startswith(rpath, "/")
        if length(rpath) <= 1
            throw(RelativePathError("Relative path `$rpath` doesn't " *
                                    "look right."))
        end
        if FS2
            rpath = replace(rpath, "/_assets" => "/assets")
        end
        canonical || return rpath
        if FS2
            full_path = joinpath(PATHS[:site], join_rpath(rpath[2:end]))
        else
            full_path = joinpath(PATHS[:folder], join_rpath(rpath[2:end]))
        end
        return normpath(full_path)

    # path relative to the assets folder with the same path as the parent file
    # example with FS2 assuming we're on page /page/
    # ./foo.bar ==> { /assets/page/foo.bar  ; [__site]/assets/page/foo.bar}
    elseif startswith(rpath, "./")
        if length(rpath) <= 2
            throw(RelativePathError("Relative path `$rpath` doesn't " *
                                    "look right."))
        end
        # if the locvar ends with `index.md` but is not strictly
        # index.md then remove it
        loc_rpath = locvar("fd_rpath")
        if endswith(loc_rpath, "index.md") && loc_rpath != "index.md"
            loc_rpath = replace(loc_rpath, r"index\.md$" => "")
        end
        if !canonical
            # here we want to remain unix-style, so we don't use joinpath
            # note that rpath never starts with "/" so there's
            # no doubling of "//"
            rpath = unixify(splitext(loc_rpath)[1]) * rpath[3:end]
            return "/assets/" * rpath
        else
            if FS2
                full_path = joinpath(PATHS[:site], "assets",
                                     splitext(loc_rpath)[1],
                                     join_rpath(rpath[3:end]))
            else
                full_path = joinpath(PATHS[:assets],
                                     splitext(locvar("fd_rpath"))[1],
                                     join_rpath(rpath[3:end]))
            end
            return normpath(full_path)
        end
    end

    # if the rpath didn't start with "/" or "./", check if it's not empty
    # and if we're in a code context
    if iszero(length(rpath))
        throw(RelativePathError("Relative path `$rpath` (empty) isn't right."))
    end
    if code # recurse adding "./"
        return parse_rpath("./code/" * rpath; canonical=canonical)
    end
    # otherwise, we consider the path as a full path relative to the
    # /assets/ folder. For instance `blah/img1.png` => `/assets/blah/img1.png`
    canonical || return "/assets/" * rpath # no doubling bc rpath[1] != '/'
    if FS2
        full_path = joinpath(PATHS[:site], "assets", join_rpath(rpath))
    else
        full_path = joinpath(PATHS[:assets], join_rpath(rpath))
    end
    return normpath(full_path)
end

"""
$(SIGNATURES)

Internal function to resolve an input (usually a script in a given language)
given a relative path (in unix format) and return a tuple with:

* the full system path to the asset
* the full system path to the directory containing the asset
* the name of the asset without its extension

If an extension is not provided, Franklin will try to find something with an
adequate extension or, failing that, will try to find anything that matches
the file name.
"""
function resolve_rpath(rpath::AS, lang::AS="")::NTuple{3,String}
    # parse the relative path
    fpath = parse_rpath(rpath; canonical=true, code=!isempty(lang))

    # check if an extension is given, if not, consider it's
    # `.xx` with language `nothing`
    fp, ext = splitext(fpath)
    if isempty(ext)
        ext, = get(CODE_LANG, lang) do
            (".xx", nothing)
        end
        fpath *= ext
    end
    dir, fname = splitdir(fp)

    # If there's an extension, check the file exists and return
    if ext != ".xx"
        isfile(fpath) && return fpath, dir, fname
        throw(FileNotFoundError("Couldn't find a file when trying to " *
                "resolve an input request with relative path: `$rpath`."))
    end

    # Otherwise there's no extension, try to find a file with
    # a matching name
    if !isdir(dir)
        throw(FileNotFoundError("Couldn't find a file when trying to " *
                "resolve an input request with relative path: `$rpath`."))
    end
    files = readdir(dir)
    k     = findfirst(e -> splitext(e)[1] == fname, files)

    # if nothing is found, error
    if isnothing(k)
        throw(FileNotFoundError("Couldn't find a file when trying to " *
                "resolve an input request with relative path: `$rpath`."))
    end

    # otherwise resolve and return
    fname_ext = files[k]
    fpath     = joinpath(dir, fname_ext)
    fname     = splitext(fname_ext)[1]

    return fpath, dir, fname
end

"""
$SIGNATURES

A convenience internal function to get the system paths related to a code
script these paths are given raw, without checking whether things exist.
"""
function form_codepaths(rpath::AS)::NamedTuple
    path        = parse_rpath(rpath; canonical=true, code=true)
    # extension-less names are allowed
    script_path = ifelse(endswith(path, ".jl"), path, path * ".jl")
    dir, fname  = splitdir(script_path)
    fname_noext = splitext(fname)[1]
    out_dir     = joinpath(dir, "output")
    out_path    = joinpath(out_dir, fname_noext * ".out")
    res_path    = joinpath(out_dir, fname_noext * ".res")
    # return everything
    return (script_path = script_path,
            script_dir  = dir,
            script_name = fname,
            out_dir     = out_dir,
            out_path    = out_path,
            res_path    = res_path)
end


text:
using CUDArt, Knet, Base.Test
include("isapprox.jl")
@time include(Pkg.dir("Knet/test/mnist.jl"))

ytype(X::DataType)=
    (X <: KUsparse{Array} ? KUdense{Array} :
     X <: KUsparse{CudaArray} ? KUdense{CudaArray} :
     X <: SparseMatrixCSC ? Array : X)

net=nothing
xtrn=xtst=nothing
ytrn=ytst=nothing
w0 = w1 = b0 = b1 = nothing

for X in (
          SparseMatrixCSC{Float32,Int32},
          KUsparse{Array},
          KUsparse{CudaArray},
          Array,
          CudaArray,
          KUdense{Array},
          KUdense{CudaArray},
          )
    Y = ytype(X)
    @show (X,Y)
    net = [Mmul(10), Bias(), PercLoss()]
    setparam!(net, average=true, init=initzero)
    # xtrn = convert(X, copy(MNIST.xtrn))
    xtst = convert(X, copy(MNIST.xtst))
    # ytrn = convert(Y, copy(MNIST.ytrn))
    ytst = convert(Y, copy(MNIST.ytst))
    # @show map(summary, (xtrn, ytrn, xtst, ytst))
    for i=1:3
        # train(net, xtrn, ytrn)
        train(net, xtst, ytst)
        println((i, accuracy(ytst, predict(net, xtst)),
                 )) # accuracy(ytrn, predict(net, xtrn))))
    end
    if w0 == nothing
        w0 = convert(Array, net[1].w.arr)
        w1 = convert(Array, net[1].w.avg)
        b0 = convert(Array, net[2].b.arr)
        b1 = convert(Array, net[2].b.avg)
    else
        @test isapprox(w0, net[1].w.arr)
        @test isapprox(w1, net[1].w.avg)
        @test isapprox(b0, net[2].b.arr)
        @test isapprox(b1, net[2].b.avg)
    end
end # for X

# DEAD CODE:

# using Knet
# using Knet: accuracy
# @time require(Pkg.dir("Knet/test/mnist.jl"))

# Knet.gpu(false)
# xtrn = float32(255*MNIST.xtrn)
# xtst = float32(255*MNIST.xtst)
# ytrn = float32(MNIST.ytrn)
# ytst = float32(MNIST.ytst)
# strn = sparse(xtrn)
# stst = sparse(xtst)

# if true
# info("Dense input")
# xnet = Op[Perceptron(10)]
# @time for i=1:5
#     train(xnet, xtrn, ytrn)
#     println((i,
#              accuracy(ytst, predict(xnet, xtst)),
#              accuracy(ytrn, predict(xnet, xtrn))))
# end
# end # if false

# if true
# info("Sparse input")
# snet = Op[Perceptron(10)]
# @time for i=1:5
#     train(snet, strn, ytrn)
#     println((i,
#              accuracy(ytst, predict(snet, stst)),
#              accuracy(ytrn, predict(snet, strn))))
# end
# end # if false

# # This is the stuff that needs to be defined for CudaArrays:
# # We probably need gpu kernels for forw, update
# using CUDArt
# import Base: .+, getindex, ./, -

# if false
# .+(a::CudaArray{Float32,2}, b::CudaArray{Float32,2})=a
# getindex(::CudaArray{Float32,2}, ::Int64, ::Int64)=0
# ./(a::CudaArray{Float32,2}, ::Int64)=a
# -(a::CudaArray{Float32,2}, b::CudaArray{Float32,2})=a

# # Knet.atype(CudaArray)
# Knet.gpu(true)
# cnet = Op[Perceptron(10)]
# @time for i=1:5
#     train(cnet, xtrn, ytrn; iters=100)
#     println((i,
#              accuracy(ytst, predict(cnet, xtst)),
#              accuracy(ytrn, predict(cnet, xtrn))))
# end
# # Knet.atype(Array)
# end # if false

# # Similar stuff needs to be defined for sparse cuda arrays.


text:
export Loss, Context, getContext, setContext, resetContext, hasgpu, is_on_gpu, KorA

const ϵ = 10e-8

"""
A meter is reponsible for presenting metric values. This can be
printing it to the console output, showing it on TensorBoard of storing it
in a database.
"""
abstract type Meter end

abstract type MetricStore end
abstract type Layer end
abstract type Metric end
abstract type Optimizer end

"""
Abstract type for the loss functions. However Photon accepts
any function as a loss function as long as it is callable and
returns the loss value as a scalar type.
"""
abstract type Loss end

const Shape = Tuple{Vararg{Int}}

const Tensor{T,N} = Union{
						AbstractArray{T,N},
						Knet.KnetArray{T,N},
						Knet.Value{Array{T,N}},
						Knet.Value{Knet.KnetArray{T,N}}}

const Tensors = Union{Tensor, Tuple{Tensor}}

"""
Context is used by various parts of Photon to determine what the device and
datatype should be for data. It allows for quickly switches between GPU and CPU
based environments.

# Attributes
- device::Symbol the type of device. For now supported :cpu and :gpu
- deviceId::Int the id of the device, useful for example if you multiple GPU's
- dtype::Type the type of data you want to use. Most common are Float32, Float16 or Float64.
"""
mutable struct Context
	device::Symbol
	deviceId::Int
	dtype::Type

	function Context()
		device = Knet.gpu() >= 0 ? :gpu : :cpu
		new(device, 0, Float32)
	end
end

global ctx = Context()

function is_on_gpu()::Bool
	ctx.device == :gpu
end

hasgpu() = Knet.gpu() >= 0

function getContext()::Context
  ctx
end

function setContext(;device=ctx.device, deviceId=ctx.deviceId, dtype=ctx.dtype)
  ctx.device = device
  ctx.deviceId = deviceId
  ctx.dtype= dtype
  getContext()
end

function resetContext()
	global ctx
	ctx = Context()
end

"""
KorA makes it easy to move an array to the GPU or the other way around
"""
KorA(arr::Array) = (ctx.device == :gpu) ? Knet.KnetArray(arr) : arr
KorA(arr::Knet.KnetArray)= (ctx.device == :cpu) ? Array(arr) : arr
KorA(arr::Tuple)= (KorA(elem) for elem in arr)

# TODO: Make more generic
toFloat32(arr::Array) = convert(Array{Float32}, arr)
toFloat32(arr::Knet.KnetArray) = convert(Knet.KnetArray{Float32}, arr)
toFloat32(arr::Tuple) = (toFloat32(elem) for elem in arr)

addlast(x) = reshape(x, (size(x)...,1))
droplast(x) = reshape(x, (size(x)[1:end-1]...))

"""
autoConvertor converts data to the right format for a model.
It uses the context to determine the device (cpu or gpu) and datatype
that the data needs to be.

It supports Tuples, Arrays and KnetArrays and a combination of those.
"""
function autoConvertor(arr::Array)
	arr = convert(Array{ctx.dtype},arr)
	ctx.device == :gpu ? Knet.KnetArray(arr) : arr
end

function autoConvertor(arr::Knet.KnetArray)
	arr = convert(Knet.KnetArray{ctx.dtype}, arr)
	ctx.device == :gpu ? arr : Array(arr)
end

autoConvertor(arr::Tuple)= (autoConvertor(elem) for elem in arr)

# small util
makeArray(x::AbstractArray) = x
makeArray(x) = Vector(x)


text:
@connector function HeatPort(; name, T_start=273.15 + 20.0, Q_flow_start=0.0)
    @variables T(t)=T_start
    @variables Q_flow(t)=Q_flow_start [connect = Flow]
    ODESystem(Equation[], t, [T, Q_flow], [], name=name)
end
Base.@doc """
    HeatPort(; name, T_start=273.15 + 20.0, Q_flow_start=0.0)

Port for a thermal system.

# Parameters:
- `T_start`: [K] Temperature of the port
- `Q_flow_start`: [W] Heat flow rate at the port

# States:
- `T`: [K] Temperature of the port
- `Q_flow`: [W] Heat flow rate at the port
""" HeatPort

"""
    Element1D(;name, dT0=0.0, Q_flow0=0.0)

This partial model contains the basic connectors and variables to allow heat transfer models to be created that do not
store energy. This model defines and includes equations for the temperature drop across the element, `dT`, and the heat
flow rate through the element from `port_a` to `port_b`, `Q_flow`.

# Parameters:
- `dT_start`:  [K] Initial temperature difference across the component a.T - b.T
- `Q_flow_start`: [W] Initial heat flow rate from port a -> port b

# States:
- `dT`:  [K] Temperature difference across the component a.T - b.T
- `Q_flow`: [W] Heat flow rate from port a -> port b
"""
function Element1D(;name, dT_start=0.0, Q_flow_start=0.0)
    @named port_a = HeatPort()
    @named port_b = HeatPort()
    sts = @variables begin
        dT(t)=dT_start
        Q_flow(t)=Q_flow_start
    end
    eqs = [
        dT ~ port_a.T - port_b.T
        port_a.Q_flow ~ Q_flow
        port_a.Q_flow + port_b.Q_flow ~ 0
    ]

    return compose(ODESystem(eqs, t, sts, []; name=name), port_a, port_b)
end


text:
sz = (4,5)
filename = tempname()

# Create a simple dataset with some attributes to the variables
ds = NCDatasets.Dataset(filename,"c")

nlon = 120;
nlat = 50;
# Define the dimension "lon" and "lat" with the size 100 and 110 resp.
NCDatasets.defDim(ds,"lon",nlon)
NCDatasets.defDim(ds,"lat",nlat)

# Define a variable
lon = NCDatasets.defVar(ds,"lon",Float32,("lon",))
lat= NCDatasets.defVar(ds,"lat",Float32,("lat",))
# create two variables with the same standard name
vort = NCDatasets.defVar(ds,"vort",Float32,("lon","lat"))
vort2 = NCDatasets.defVar(ds,"vort2",Float32,("lon","lat"))

# write attributes
lon.attrib["standard_name"] = "longitude"
lat.attrib["standard_name"] = "latitude"
vort.attrib["units"] = "s^-1"
vort.attrib["standard_name"] = "ocean_relative_vorticity"
vort2.attrib["standard_name"] = "ocean_relative_vorticity"

close(ds)

ds = NCDatasets.Dataset(filename,"r")

# Single variables by standard name or units
@test NCDatasets.name.(NCDatasets.varbyattrib(ds, standard_name = "longitude")) == ["lon"]
@test NCDatasets.name.(NCDatasets.varbyattrib(ds, standard_name = "latitude")) == ["lat"]
@test NCDatasets.name.(NCDatasets.varbyattrib(ds, units = "s^-1")) == ["vort"]

# Two variables with the same attribute
@test length(NCDatasets.varbyattrib(ds, standard_name = "ocean_relative_vorticity")) == 2
@test NCDatasets.name(NCDatasets.varbyattrib(ds, standard_name = "ocean_relative_vorticity")[1]) == "vort"

# Empty result
@test NCDatasets.varbyattrib(ds, standard_name = "time") == []

# Bad input types (not Strings)
#@test_throws MethodError NCDatasets.varbyattrib("test.nc", 1, "lon")
#@test_throws MethodError NCDatasets.varbyattrib("test.nc", "standard_name", 12.2)

NCDatasets.close(ds)
rm(filename)


text:
"""
$(TYPEDEF)
"""
struct HermiteInterpolation{T1,T2,T3} <: AbstractDiffEqInterpolation
  t::T1
  u::T2
  du::T3
end

"""
$(TYPEDEF)
"""
struct LinearInterpolation{T1,T2} <: AbstractDiffEqInterpolation
  t::T1
  u::T2
end

"""
$(TYPEDEF)
"""
struct ConstantInterpolation{T1,T2} <: AbstractDiffEqInterpolation
  t::T1
  u::T2
end

interp_summary(::AbstractDiffEqInterpolation) = "Unknown"
interp_summary(::HermiteInterpolation) = "3rd order Hermite"
interp_summary(::LinearInterpolation) = "1st order linear"
interp_summary(::ConstantInterpolation) = "Piecewise constant interpolation"
interp_summary(::Nothing) = "No interpolation"
interp_summary(sol::DESolution) = interp_summary(sol.interp)

(id::HermiteInterpolation)(tvals,idxs,deriv,p,continuity::Symbol=:left) = interpolation(tvals,id,idxs,deriv,p,continuity)
(id::HermiteInterpolation)(val,tvals,idxs,deriv,p,continuity::Symbol=:left) = interpolation!(val,tvals,id,idxs,deriv,p,continuity)
(id::LinearInterpolation)(tvals,idxs,deriv,p,continuity::Symbol=:left) = interpolation(tvals,id,idxs,deriv,p,continuity)
(id::LinearInterpolation)(val,tvals,idxs,deriv,p,continuity::Symbol=:left) = interpolation!(val,tvals,id,idxs,deriv,p,continuity)
(id::ConstantInterpolation)(tvals,idxs,deriv,p,continuity::Symbol=:left) = interpolation(tvals,id,idxs,deriv,p,continuity)
(id::ConstantInterpolation)(val,tvals,idxs,deriv,p,continuity::Symbol=:left) = interpolation!(val,tvals,id,idxs,deriv,p,continuity)

@inline function interpolation(tvals,id,idxs,deriv,p,continuity::Symbol=:left)
  t = id.t; u = id.u
  typeof(id) <: HermiteInterpolation && (du = id.du)
  tdir = sign(t[end]-t[1])
  idx = sortperm(tvals,rev=tdir<0)
  i = 2 # Start the search thinking it's between t[1] and t[2]
  tdir*tvals[idx[end]] > tdir*t[end] && error("Solution interpolation cannot extrapolate past the final timepoint. Either solve on a longer timespan or use the local extrapolation from the integrator interface.")
  tdir*tvals[idx[1]] < tdir*t[1] && error("Solution interpolation cannot extrapolate before the first timepoint. Either start solving earlier or use the local extrapolation from the integrator interface.")
  if typeof(idxs) <: Number
    vals = Vector{eltype(first(u))}(undef, length(tvals))
  elseif typeof(idxs) <: AbstractVector
     vals = Vector{Vector{eltype(first(u))}}(undef, length(tvals))
  else
    vals = Vector{eltype(u)}(undef, length(tvals))
  end
  @inbounds for j in idx
    tval = tvals[j]
    i = searchsortedfirst(@view(t[i:end]),tval,rev=tdir<0)+i-1 # It's in the interval t[i-1] to t[i]
    avoid_constant_ends = deriv != Val{0} #|| typeof(tval) <: ForwardDiff.Dual
    avoid_constant_ends && i==1 && (i+=1)
    if !avoid_constant_ends && t[i] == tval
      lasti = lastindex(t)
      k = continuity == :right && i+1 <= lasti && t[i+1] == tval ? i+1 : i
      if idxs === nothing
        vals[j] = u[k]
      else
        vals[j] = u[k][idxs]
      end
    elseif !avoid_constant_ends && t[i-1] == tval # Can happen if it's the first value!
      if idxs === nothing
        vals[j] = u[i-1]
      else
        vals[j] = u[i-1][idxs]
      end
    else
      dt = t[i] - t[i-1]
      Θ = (tval-t[i-1])/dt
      idxs_internal = idxs
      if typeof(id) <: HermiteInterpolation
        vals[j] = interpolant(Θ,id,dt,u[i-1],u[i],du[i-1],du[i],idxs_internal,deriv)
      else
        vals[j] = interpolant(Θ,id,dt,u[i-1],u[i],idxs_internal,deriv)
      end
    end
  end
  DiffEqArray(vals,tvals)
end

"""
$(SIGNATURES)

Get the value at tvals where the solution is known at the
times t (sorted), with values u and derivatives ks
"""
@inline function interpolation!(vals,tvals,id,idxs,deriv,p,continuity::Symbol=:left)
  t = id.t; u = id.u
  typeof(id) <: HermiteInterpolation && (du = id.du)
  tdir = sign(t[end]-t[1])
  idx = sortperm(tvals,rev=tdir<0)
  i = 2 # Start the search thinking it's between t[1] and t[2]
  tdir*tvals[idx[end]] > tdir*t[end] && error("Solution interpolation cannot extrapolate past the final timepoint. Either solve on a longer timespan or use the local extrapolation from the integrator interface.")
  tdir*tvals[idx[1]] < tdir*t[1] && error("Solution interpolation cannot extrapolate before the first timepoint. Either start solving earlier or use the local extrapolation from the integrator interface.")
  @inbounds for j in idx
    tval = tvals[j]
    i = searchsortedfirst(@view(t[i:end]),tval,rev=tdir<0)+i-1 # It's in the interval t[i-1] to t[i]
    avoid_constant_ends = deriv != Val{0} #|| typeof(tval) <: ForwardDiff.Dual
    avoid_constant_ends && i==1 && (i+=1)
    if !avoid_constant_ends && t[i] == tval
      lasti = lastindex(t)
      k = continuity == :right && i+1 <= lasti && t[i+1] == tval ? i+1 : i
      if idxs === nothing
        vals[j] = u[k]
      else
        vals[j] = u[k][idxs]
      end
    elseif !avoid_constant_ends && t[i-1] == tval # Can happen if it's the first value!
      if idxs === nothing
        vals[j] = u[i-1]
      else
        vals[j] = u[i-1][idxs]
      end
    else
      dt = t[i] - t[i-1]
      Θ = (tval-t[i-1])/dt
      idxs_internal = idxs
      if eltype(u) <: Union{AbstractArray,ArrayPartition}
        if typeof(id) <: HermiteInterpolation
          interpolant!(vals[j],Θ,id,dt,u[i-1],u[i],du[i-1],du[i],idxs_internal,deriv)
        else
          interpolant!(vals[j],Θ,id,dt,u[i-1],u[i],idxs_internal,deriv)
        end
      else
        if typeof(id) <: HermiteInterpolation
          vals[j] = interpolant(Θ,id,dt,u[i-1],u[i],du[i-1],du[i],idxs_internal,deriv)
        else
          vals[j] = interpolant(Θ,id,dt,u[i-1],u[i],idxs_internal,deriv)
        end
      end
    end
  end
end

"""
$(SIGNATURES)

Get the value at tval where the solution is known at the
times t (sorted), with values u and derivatives ks
"""
@inline function interpolation(tval::Number,id,idxs,deriv,p,continuity::Symbol=:left)
  t = id.t; u = id.u
  typeof(id) <: HermiteInterpolation && (du = id.du)
  tdir = sign(t[end]-t[1])
  tdir*tval > tdir*t[end] && error("Solution interpolation cannot extrapolate past the final timepoint. Either solve on a longer timespan or use the local extrapolation from the integrator interface.")
  tdir*tval < tdir*t[1] && error("Solution interpolation cannot extrapolate before the first timepoint. Either start solving earlier or use the local extrapolation from the integrator interface.")
  @inbounds i = searchsortedfirst(t,tval,rev=tdir<0) # It's in the interval t[i-1] to t[i]
  avoid_constant_ends = deriv != Val{0} #|| typeof(tval) <: ForwardDiff.Dual
  avoid_constant_ends && i==1 && (i+=1)
  @inbounds if !avoid_constant_ends && t[i] == tval
    lasti = lastindex(t)
    k = continuity == :right && i+1 <= lasti && t[i+1] == tval ? i+1 : i
    if idxs === nothing
      val = u[k]
    else
      val = u[k][idxs]
    end
  elseif !avoid_constant_ends && t[i-1] == tval # Can happen if it's the first value!
    if idxs === nothing
      val = u[i-1]
    else
      val = u[i-1][idxs]
    end
  else
    dt = t[i] - t[i-1]
    Θ = (tval-t[i-1])/dt
    idxs_internal = idxs
    if typeof(id) <: HermiteInterpolation
      val = interpolant(Θ,id,dt,u[i-1],u[i],du[i-1],du[i],idxs_internal,deriv)
    else
      val = interpolant(Θ,id,dt,u[i-1],u[i],idxs_internal,deriv)
    end
  end
  val
end

"""
$(SIGNATURES)

Get the value at tval where the solution is known at the
times t (sorted), with values u and derivatives ks
"""
@inline function interpolation!(out,tval::Number,id,idxs,deriv,p,continuity::Symbol=:left)
  t = id.t; u = id.u
  typeof(id) <: HermiteInterpolation && (du = id.du)
  tdir = sign(t[end]-t[1])
  tdir*tval > tdir*t[end] && error("Solution interpolation cannot extrapolate past the final timepoint. Either solve on a longer timespan or use the local extrapolation from the integrator interface.")
  tdir*tval < tdir*t[1] && error("Solution interpolation cannot extrapolate before the first timepoint. Either start solving earlier or use the local extrapolation from the integrator interface.")
  @inbounds i = searchsortedfirst(t,tval,rev=tdir<0) # It's in the interval t[i-1] to t[i]
  avoid_constant_ends = deriv != Val{0} #|| typeof(tval) <: ForwardDiff.Dual
  avoid_constant_ends && i==1 && (i+=1)
  @inbounds if !avoid_constant_ends && t[i] == tval
    lasti = lastindex(t)
    k = continuity == :right && i+1 <= lasti && t[i+1] == tval ? i+1 : i
    if idxs === nothing
      copy!(out,u[k])
    else
      copy!(out,u[k][idxs])
    end
  elseif !avoid_constant_ends && t[i-1] == tval # Can happen if it's the first value!
    if idxs === nothing
      copy!(out,u[i-1])
    else
      copy!(out,u[i-1][idxs])
    end
  else
    dt = t[i] - t[i-1]
    Θ = (tval-t[i-1])/dt
    idxs_internal = idxs
    if typeof(id) <: HermiteInterpolation
      interpolant!(out,Θ,id,dt,u[i-1],u[i],du[i-1],du[i],idxs_internal,deriv)
    else
      interpolant!(out,Θ,id,dt,u[i-1],u[i],idxs_internal,deriv)
    end
  end
end

##################### Hermite Interpolants

"""
Hairer Norsett Wanner Solving Ordinary Differential Equations I - Nonstiff Problems Page 190

Hermite Interpolation
"""
@inline function interpolant(Θ,id::HermiteInterpolation,dt,y₀,y₁,dy₀,dy₁,idxs,T::Type{Val{0}})
  if idxs === nothing
    out = @. (1-Θ)*y₀+Θ*y₁+Θ*(Θ-1)*((1-2Θ)*(y₁-y₀)+(Θ-1)*dt*dy₀ + Θ*dt*dy₁)
  else
    out = similar(y₀,axes(idxs))
    @views @. out = (1-Θ)*y₀[idxs]+Θ*y₁[idxs]+
                    Θ*(Θ-1)*((1-2Θ)*(y₁[idxs]-y₀[idxs])+
                    (Θ-1)*dt*dy₀[idxs] + Θ*dt*dy₁[idxs])
  end
  out
end

"""
Hermite Interpolation
"""
@inline function interpolant(Θ,id::HermiteInterpolation,dt,y₀,y₁,dy₀,dy₁,idxs,T::Type{Val{1}})
  if idxs === nothing
    out = @. dy₀ + Θ*(-4*dt*dy₀ - 2*dt*dy₁ - 6*y₀ + Θ*(3*dt*dy₀ + 3*dt*dy₁ + 6*y₀ - 6*y₁) + 6*y₁)/dt
  else
    out = similar(y₀,axes(idxs))
    @views @. out = dy₀[idxs] + Θ*(-4*dt*dy₀[idxs] -
                    2*dt*dy₁[idxs] - 6*y₀[idxs] +
                    Θ*(3*dt*dy₀[idxs] + 3*dt*dy₁[idxs] +
                    6*y₀[idxs] - 6*y₁[idxs]) + 6*y₁[idxs])/dt
  end
  out
end

"""
Hermite Interpolation
"""
@inline function interpolant(Θ,id::HermiteInterpolation,dt,y₀,y₁,dy₀,dy₁,idxs,T::Type{Val{2}})
  if idxs === nothing
    out = @. (-4*dt*dy₀ - 2*dt*dy₁ - 6*y₀ + Θ*(6*dt*dy₀ + 6*dt*dy₁ + 12*y₀ - 12*y₁) + 6*y₁)/(dt*dt)
  else
    out = similar(y₀,axes(idxs))
    @views @. out = (-4*dt*dy₀[idxs] - 2*dt*dy₁[idxs] - 6*y₀[idxs] +
                    Θ*(6*dt*dy₀[idxs] + 6*dt*dy₁[idxs] + 12*y₀[idxs] -
                    12*y₁[idxs]) + 6*y₁[idxs])/(dt*dt)
  end
  out
end

"""
Hermite Interpolation
"""
@inline function interpolant(Θ,id::HermiteInterpolation,dt,y₀,y₁,dy₀,dy₁,idxs,T::Type{Val{3}})
  if idxs === nothing
    out = @. (6*dt*dy₀ + 6*dt*dy₁ + 12*y₀ - 12*y₁)/(dt*dt*dt)
  else
    out = similar(y₀,axes(idxs))
    @views @. out = (6*dt*dy₀[idxs] + 6*dt*dy₁[idxs] +
                    12*y₀[idxs] - 12*y₁[idxs])/(dt*dt*dt)
  end
  out
end

"""
Hairer Norsett Wanner Solving Ordinary Differential Euations I - Nonstiff Problems Page 190

Hermite Interpolation
"""
@inline function interpolant!(out,Θ,id::HermiteInterpolation,dt,y₀,y₁,dy₀,dy₁,idxs,T::Type{Val{0}})
  if out === nothing
    return (1-Θ)*y₀[idxs]+Θ*y₁[idxs]+Θ*(Θ-1)*((1-2Θ)*(y₁[idxs]-y₀[idxs])+(Θ-1)*dt*dy₀[idxs] + Θ*dt*dy₁[idxs])
  elseif idxs === nothing
    @. out = (1-Θ)*y₀+Θ*y₁+Θ*(Θ-1)*((1-2Θ)*(y₁-y₀)+(Θ-1)*dt*dy₀ + Θ*dt*dy₁)
  else
    @views @. out = (1-Θ)*y₀[idxs]+Θ*y₁[idxs]+Θ*(Θ-1)*((1-2Θ)*(y₁[idxs]-y₀[idxs])+(Θ-1)*dt*dy₀[idxs] + Θ*dt*dy₁[idxs])
  end
end

"""
Hermite Interpolation
"""
@inline function interpolant!(out,Θ,id::HermiteInterpolation,dt,y₀,y₁,dy₀,dy₁,idxs,T::Type{Val{1}})
  if out === nothing
    return dy₀[idxs] + Θ*(-4*dt*dy₀[idxs] - 2*dt*dy₁[idxs] - 6*y₀[idxs] + Θ*(3*dt*dy₀[idxs] + 3*dt*dy₁[idxs] + 6*y₀[idxs] - 6*y₁[idxs]) + 6*y₁[idxs])/dt
  elseif idxs === nothing
    @. out = dy₀ + Θ*(-4*dt*dy₀ - 2*dt*dy₁ - 6*y₀ + Θ*(3*dt*dy₀ + 3*dt*dy₁ + 6*y₀ - 6*y₁) + 6*y₁)/dt
  else
    @views @. out = dy₀[idxs] + Θ*(-4*dt*dy₀[idxs] - 2*dt*dy₁[idxs] - 6*y₀[idxs] + Θ*(3*dt*dy₀[idxs] + 3*dt*dy₁[idxs] + 6*y₀[idxs] - 6*y₁[idxs]) + 6*y₁[idxs])/dt
  end
end

"""
Hermite Interpolation
"""
@inline function interpolant!(out,Θ,id::HermiteInterpolation,dt,y₀,y₁,dy₀,dy₁,idxs,T::Type{Val{2}})
  if out === nothing
    return (-4*dt*dy₀[idxs] - 2*dt*dy₁[idxs] - 6*y₀[idxs] + Θ*(6*dt*dy₀[idxs] + 6*dt*dy₁[idxs] + 12*y₀[idxs] - 12*y₁[idxs]) + 6*y₁[idxs])/(dt*dt)
  elseif idxs === nothing
    @. out = (-4*dt*dy₀ - 2*dt*dy₁ - 6*y₀ + Θ*(6*dt*dy₀ + 6*dt*dy₁ + 12*y₀ - 12*y₁) + 6*y₁)/(dt*dt)
  else
    @views @. out = (-4*dt*dy₀[idxs] - 2*dt*dy₁[idxs] - 6*y₀[idxs] + Θ*(6*dt*dy₀[idxs] + 6*dt*dy₁[idxs] + 12*y₀[idxs] - 12*y₁[idxs]) + 6*y₁[idxs])/(dt*dt)
  end
end

"""
Hermite Interpolation
"""
@inline function interpolant!(out,Θ,id::HermiteInterpolation,dt,y₀,y₁,dy₀,dy₁,idxs,T::Type{Val{3}})
  if out === nothing
    return (6*dt*dy₀[idxs] + 6*dt*dy₁[idxs] + 12*y₀[idxs] - 12*y₁[idxs])/(dt*dt*dt)
  elseif idxs === nothing
    @. out = (6*dt*dy₀ + 6*dt*dy₁ + 12*y₀ - 12*y₁)/(dt*dt*dt)
  else
    @views @. out = (6*dt*dy₀[idxs] + 6*dt*dy₁[idxs] + 12*y₀[idxs] - 12*y₁[idxs])/(dt*dt*dt)
  end
end

############################### Linear Interpolants

"""
Linear Interpolation
"""
@inline function interpolant(Θ,id::LinearInterpolation,dt,y₀,y₁,idxs,T::Type{Val{0}})
  if idxs === nothing
    out = @. (1-Θ)*y₀ + Θ*y₁
  else
    out = similar(y₀,axes(idxs))
    Θm1 = (1-Θ)
    @views @. out = Θm1*y₀[idxs] + Θ*y₁[idxs]
  end
  out
end

@inline function interpolant(Θ,id::LinearInterpolation,dt,y₀,y₁,idxs,T::Type{Val{1}})
  if idxs === nothing
    out = @. (y₁ - y₀)/dt
  else
    out = similar(y₀,axes(idxs))
    @views @. out = (y₁[idxs] - y₀[idxs])/dt
  end
  out
end

"""
Linear Interpolation
"""
@inline function interpolant!(out,Θ,id::LinearInterpolation,dt,y₀,y₁,idxs,T::Type{Val{0}})
  Θm1 = (1-Θ)
  if out === nothing
    return Θm1*y₀[idxs] + Θ*y₁[idxs]
  elseif idxs === nothing
    @. out = Θm1*y₀ + Θ*y₁
  else
    @views @. out = Θm1*y₀[idxs] + Θ*y₁[idxs]
  end
end

"""
Linear Interpolation
"""
@inline function interpolant!(out,Θ,id::LinearInterpolation,dt,y₀,y₁,idxs,T::Type{Val{1}})
  if out === nothing
    return (y₁[idxs] - y₀[idxs])/dt
  elseif idxs === nothing
    @. out = (y₁ - y₀)/dt
  else
    @views @. out = (y₁[idxs] - y₀[idxs])/dt
  end
end

############################### Linear Interpolants

"""
Constant Interpolation
"""
@inline function interpolant(Θ,id::ConstantInterpolation,dt,y₀,y₁,idxs,T::Type{Val{0}})
  if idxs === nothing
    out = @. y₀
  else
    out = similar(y₀,axes(idxs))
    @views @. out = y₀[idxs]
  end
  out
end

@inline function interpolant(Θ,id::ConstantInterpolation,dt,y₀,y₁,idxs,T::Type{Val{1}})
  if idxs === nothing
    out = zeros(eltype(y₀),length(y₀))
  else
    out = similar(y₀,axes(idxs))
    @views @. out = 0
  end
  out
end

"""
Constant Interpolation
"""
@inline function interpolant!(out,Θ,id::ConstantInterpolation,dt,y₀,y₁,idxs,T::Type{Val{0}})
  if out === nothing
    return y₀[idxs]
  elseif idxs === nothing
    @. out = y₀
  else
    @views @. out = y₀[idxs]
  end
end

"""
Constant Interpolation
"""
@inline function interpolant!(out,Θ,id::ConstantInterpolation,dt,y₀,y₁,idxs,T::Type{Val{1}})
  if out === nothing
    return zeros(eltype(y₀),length(idxs))
  else
    @. out = 0
  end
end


text:
"""
Parameter state in Gibbs sampler for feature allocation model
with feature selection. This augments `State` with parameters:
`r`, `W_star`, and `p`-or-`omega`.

Usage example:
==============

```julia
julia> s = StateFS{Float64}()
julia> s.r = rand(Bool, 3, 10)
```

NOTE:
=====
The update order should be:
Z -> v -> alpha -> p -> r -> lam -> W* -> gamma -> eta -> delta -> sig^2 -> y*

Particularly, `lam` must be updated before W* and preferrably after p & r.
And `gamma` should be updated after Z and before delta.
"""
mutable struct StateFS{F <: AbstractFloat}
  r::Matrix{Bool}  # dim: IxK
  W_star::Matrix{F}  # dim: IxK
  omega::Vector{F}  # dim: number of covariates (`length(x_i)`) + 1
  theta::State{F}  # State

  # Primary constructor returns an uninitialized stateFS.
  StateFS{F}() where {F <: AbstractFloat} = new()
end

function StateFS{F}(theta::State{F}, dfs::DataFS) where {F <: AbstractFloat}
  sfs = StateFS{F}()
  I, K = size(theta.W)

  sfs.theta = theta
  # sfs.r = rand(Bool, I, K)
  sfs.r = ones(Bool, I, K)
  # sfs.W_star = ones(I, K)
  sfs.W_star = theta.W * 5
  sfs.omega = zeros(dfs.P)

  return sfs
end


text:
using TKTDbayes
using Test

max_stars_repo_name: virgile-baudrot/TKTDbayes

max_issues_repo_name: virgile-baudrot/TKTDbayes

max_forks_repo_name: virgile-baudrot/TKTDbayes