text:
#=
Copyright (c) 2015, Intel Corporation
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
- Redistributions of source code must retain the above copyright notice,
  this list of conditions and the following disclaimer.
- Redistributions in binary form must reproduce the above copyright notice,
  this list of conditions and the following disclaimer in the documentation
  and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
THE POSSIBILITY OF SUCH DAMAGE.
=#

using ParallelAccelerator
using DocOpt
using Images

@acc function blur(img, iterations)
    buf = Array(Float32, size(img)...)
    runStencil(buf, img, iterations, :oob_skip) do b, a
       b[0,0] =
            (a[-2,-2] * 0.003  + a[-1,-2] * 0.0133 + a[0,-2] * 0.0219 + a[1,-2] * 0.0133 + a[2,-2] * 0.0030 +
             a[-2,-1] * 0.0133 + a[-1,-1] * 0.0596 + a[0,-1] * 0.0983 + a[1,-1] * 0.0596 + a[2,-1] * 0.0133 +
             a[-2, 0] * 0.0219 + a[-1, 0] * 0.0983 + a[0, 0] * 0.1621 + a[1, 0] * 0.0983 + a[2, 0] * 0.0219 +
             a[-2, 1] * 0.0133 + a[-1, 1] * 0.0596 + a[0, 1] * 0.0983 + a[1, 1] * 0.0596 + a[2, 1] * 0.0133 +
             a[-2, 2] * 0.003  + a[-1, 2] * 0.0133 + a[0, 2] * 0.0219 + a[1, 2] * 0.0133 + a[2, 2] * 0.0030)
       return a, b
    end
    return img
end

function main()
    doc = """gaussian-blur.jl

Gaussian blur image processing.

Usage:
  gaussian-blur.jl -h | --help
  gaussian-blur.jl [--img-file=<img-file>] [--iterations=<iterations>]

Options:
  -h --help                  Show this screen.
  --img-file=<img-file>      Specify a path to an input image file [default: ../example.jpg].
  --iterations=<iterations>  Specify a number of iterations [default: 100].
"""
    arguments = docopt(doc)

    dir = dirname(@__FILE__)
    img_file = joinpath(dir, arguments["--img-file"])
    iterations = parse(Int, arguments["--iterations"])

    (fname, ext) = splitext(img_file)
    out_file = string(fname, "-blur", ".jpg")

    println("input file = ", img_file)
    println("iterations = ", iterations)
    println("output file = ", out_file)

    function gaussian_blur(input_fname, output_fname, iterations)
        local img = convert(Array{Float32, 2}, load(input_fname))
        tic()
        blur(img, 0)
        println("SELFPRIMED ", toq())

        tic()
        img = blur(img, iterations)
        println("SELFTIMED ", toq())
        save(output_fname, img)
    end

    gaussian_blur(img_file, out_file, iterations)

end

main()


text:
using OffsetArrays
using FFTW

bwstoyan(pp) = 0.15 / sqrt(npoints(pp)/area(window(pp)))
"""
    pcf(xy, window, r, bandwidth=bwstoyan)

Return pair-correlation function estimate for point pattern xy using translation edge correction.

The kernel estimator only uses the r values provided plus padding.

Bandwidth can be specified as a function of xy and window or as a number.

Stoyan & Stoyan, 1994, p 284
"""
pcf(pp, r, bandwidth=bwstoyan) = pcf(pp, r, bandwidth(pp))
function pcf(pp, r, bandwidth::Number)
    xy = pp.data
    window = pp.window
    if length(xy) == 0
        return fill(0.0/0.0, length(r))
    end
    #xy = getxy.(xy)
    λ = length(xy)/area(window)
    a = bandwidth

    buf = 4*ceil(Int, a/step(r))
    A = pi * 2 * step(r) * r * length(xy) * λ
    h = disthisttrans(xy, step(r), length(r) + buf, window, true)
    p1 = vcat(fill(0.0, buf), h)
    rk = range(-step(r)*(div(length(r),2) + buf - 1), step=step(r), length=length(r) + 2buf)
    kernel = (x -> ifelse(-a < x < a, 3/4 * (1 - (x/a)^2)/a, 0.0)).(rk)
    #println(OffsetArray(rk, -div(length(kernel), 2))[0])
    kernel = OffsetArray(kernel, -div(length(kernel), 2))
    kernel ./= sum(kernel)
    irfft(rfft(p1) .* rfft(kernel), length(kernel))[(1+buf:end-buf) .- 1] ./ A
end

ext: jl


text:
#===============
    STRUCTS
===============#

mutable struct sphere <: RayTracer.RTObject
    id::UUID
    radius::Float64
    position::rt_tuple
    transform::Array{Float64,2}
    material::RayTracer.RTMaterial
    function sphere()
        return new(uuid4(),1.0,RayTracer.point(0.,0.,0.),Matrix(I,4,4),RayTracer.material())
    end
end

#===============
    FUNCTIONS
===============#

function intersect(s::sphere,r::ray)::Array{intersection,1}
    r2 = transform(r,inv(s.transform))
    x0 =  r2.origin - s.position
    a = RayTracer.dot(r2.direction,r2.direction)
    b = 2*RayTracer.dot(r2.direction,x0)
    c = RayTracer.dot(x0,x0) - s.radius^2

    discriminant = b^2 - 4*a*c
    if discriminant<0
        return []
    end
    t1 = -(b+sqrt(discriminant))/(2*a)
    t2 = -(b-sqrt(discriminant))/(2*a)
    return [intersection(t1,s),intersection(t2,s)]
end

function normal_at(s::sphere,p::RayTracer.rt_tuple)::RayTracer.rt_tuple
    object_normal = inv(s.transform)*p - s.position
    return_transform = zeros(4,4)
    return_transform[1:3,1:3] = inv(RayTracer.submatrix(s.transform,4,4))'
    return_transform[4,4] = 1
    return RayTracer.normalize(return_transform*object_normal)
end


text:
"""
$(TYPEDEF)

Abstract type for transient solution
"""
abstract type AbstractTransientSolution{T,N,A,B}  <: AbstractDiffEqArray{T, N, A} end

"""
$(TYPEDEF)

Transient solution structure

## Fields
$(TYPEDFIELDS)

## Interface

Object of this type adhere to the `AbstractDiffEqArray`  interface.
For indexing and interpolation, see [https://diffeq.sciml.ai/stable/basics/solution/](https://diffeq.sciml.ai/stable/basics/solution/).

In particular, a TransientSolution `sol` can be accessed as follows:
- `sol[it]` contains the solution for timestep `i`
- `sol[ispec,:,it]` contains the solution for component `ispec` at timestep `i`
- `sol(t)` returns a (linearly) interpolated solution value for `t`.
- `sol.t[it]` is the corresponding time
- `sol[ispec,ix,it]` refers to solution of component `ispec` at node `ix` at moment `it`

"""
mutable struct TransientSolution{T,N,A,B} <: AbstractTransientSolution{T,N,A,B}
    """
    Vector of solutions
    """
    u::A
    """
    Vector of times
    """
    t::B
end

TransientSolution(vec::AbstractVector{T}, ts, ::NTuple{N}) where {T, N} = TransientSolution{eltype(T), N, typeof(vec), typeof(ts)}(vec, ts)

TransientSolution(vec::AbstractVector,ts::AbstractVector) = TransientSolution(vec, ts, (size(vec[1])..., length(vec)))

Base.append!(s::TransientSolution,t::Real,sol::AbstractArray)=push!(s.t,t), push!(s.u,copy(sol))

(sol::TransientSolution)(t) = _interpolate(sol,t)

function _interpolate(sol,t)
    if isapprox(t,sol.t[1],atol=1.0e-10*abs(sol.t[2]-sol.t[1]))
        return sol[1]
    end
    idx=searchsortedfirst(sol.t,t)
    if idx==1 || idx > length(sol)
        return nothing
    end
    if t==sol.t[idx-1]
        return sol[idx-1]
    else
        retval=similar(sol[idx])
        dt=sol.t[idx]-sol.t[idx-1]
        a=(sol.t[idx]-t)/dt
        b=(t-sol.t[idx-1])/dt
        retval.=a*sol[idx-1] + b*sol[idx]
    end
end

mutable struct VectorOfDiskArrays{T} <: AbstractVector{T}
    fname::String
    file::Union{JLD2.JLDFile,Nothing}
    n::Int64
end

function Base.push!(v::VectorOfDiskArrays, obj)
    v.n+=1
    if isnothing(v.file)
        jldopen(v.fname,"a+") do file
            file[string(v.n)]=obj
        end
    else
        v.file[string(v.n)]=obj
    end
end

Base.size(v::VectorOfDiskArrays)=(v.n,)
Base.length(v::VectorOfDiskArrays)=v.n
Base.eltype(v::VectorOfDiskArrays{T}) where T = T
function Base.getindex(v::VectorOfDiskArrays,i)
    if isnothing(v.file)
        jldopen(v.fname,"r") do file
            file[string(i)]
        end
    else
        v.file[string(i)]
    end
end

_tempname()=Base.VERSION<v"1.4" ? tempname() : tempname(pwd())*".jld2"

"""
````
VectorOfDiskArrays(firstobj:AbstractArray;
                   keep_open=true,
                   fname= fname=tempname(pwd())*".jld2")
````
Constructor of vector of arrays stored on disk (via JLD2).

- `keep_open`: if true, disk file is not closed during the existence of the object
- `fname`: file name for the disk file

The disk file is automatically removed if the object is garbage collected.
"""
function VectorOfDiskArrays(obj::AbstractArray{T}; keep_open=true, fname=_tempname()) where T
    file=jldopen(fname,"a+")
    file[string(1)]=obj
    if !keep_open
        close(file)
        file=nothing
    end
    v=VectorOfDiskArrays{T}(fname,file,1)
    finalizer(v-> (isnothing(v.file) ? nothing : close(v.file); rm(v.fname,force=true)),v)
    v
end

Base.append!(s::TransientSolution{T,N,VectorOfDiskArrays{T},B}, t::Real,sol::AbstractArray) where {T,N,B} = push!(s.t,t), push!(s.u,sol)

"""
````
TransientSolution(t0,inival;
                  in_memory=true,
                  keep_open=true,
                  fname=tempname(pwd())*".jld2"
````
Constructor of transient solution with initial value and inital time.

- `in_memory`: if true (default), data are kept in main memory, otherwise on disk (via JLD2)
- `keep_open`: if true, disk file is not closed during the existence of the object
- `fname`: file name for the disk file
"""
function TransientSolution(t0::Number,
                           inival::AbstractArray{T};
                           in_memory=true,
                           keep_open=true,
                           fname=_tempname()) where T
    if !in_memory && !isa(inival,SparseSolutionArray)
        TransientSolution(VectorOfDiskArrays(inival,keep_open=keep_open,fname=fname),[t0])
    else
        TransientSolution([copy(inival)], [t0])
    end
end


text:
module bankNote

#--- functions

function _5dollar(value=5)
    return value
end
function _10dollar(value=10)
    return value
end
function _20dollar(value=20)
return value
end

#TODO: time of start counting of Inflation -i.e. inflation, since When?
```
a generic dollar amount
has a value
Changing based on inflation

a generic dollar amount , has a value Changing based on inflation
TODO: simple math interpolation
```
function _Xdollar(value, inflation=10^-1)

    intrinsic_value =  (value - value * inflation)
    #if has
end

#The must-have factor #TODO: wehich library of differentiation to choose from?
```
velocity [1t derivative of distance ]
acceleration [2nd derivative of distance ]: based on netwon's renowned f = m * a , where a = f/m
#Edit: either : 1. Elaborate on the Context(newton), change line above to somewhere
(more)
a dollar has a physical value, that is changing by changin the
1. velocity: change of a dollar; how fast dollar rises (falls) in t days (months,years)
2. acceleration: how fast dollar velocity changes in t days (months,years)

```
function _Xdollar(value, velocity, acceleration,inflation=10^-1)

    intrinsic_value =  (value - value * inflation)
    #if has
end
#end

max_stars_repo_name: adamwillisXanax/accountingJulia

max_issues_repo_name: adamwillisXanax/accountingJulia

max_forks_repo_name: adamwillisXanax/accountingJulia


text:
@inline function indmin3(a,b,c,i,j)
    if a <= b
        if a <= c
            return 1,i-1,j-1
        else
            return 3,i,j-1
        end
    else
        if b <= c
            return 2,i-1,j
        else
            return 3,i,j-1
        end
    end
end


text:
if VERSION == v"1.3.1"
    using Pkg.Artifacts: do_artifact_str, find_artifacts_toml, load_artifacts_toml

    # A copy of `Pkg.Artifacts.@artifact_str` where `name` is properly escaped on
    # Julia 1.3.1
    # https://github.com/JuliaLang/Pkg.jl/issues/1912
    # https://github.com/JuliaLang/Pkg.jl/pull/1580
    macro artifact_str(name)
        # Load Artifacts.toml at compile time, so that we don't have to use `__source__.file`
        # at runtime, which gets stale if the `.ji` file is relocated.
        local artifacts_toml = find_artifacts_toml(string(__source__.file))
        if artifacts_toml === nothing
            error(string(
                "Cannot locate '(Julia)Artifacts.toml' file when attempting to use artifact '",
                name,
                "' in '",
                __module__,
                "'",
            ))
        end

        local artifact_dict = load_artifacts_toml(artifacts_toml)
        return quote
            # Invalidate .ji file if Artifacts.toml file changes
            Base.include_dependency($(artifacts_toml))

            # Use invokelatest() to introduce a compiler barrier, preventing many backedges from being added
            # and slowing down not only compile time, but also `.ji` load time.  This is critical here, as
            # artifact"" is used in other modules, so we don't want to be spreading backedges around everywhere.
            Base.invokelatest(do_artifact_str, $(esc(name)), $(artifact_dict), $(artifacts_toml), $__module__)
        end
    end
elseif VERSION >= v"1.3.0"
    # Issue does not exist in 1.3.0, and >= 1.4.0. Also, assume that the issue would also be
    # fixed on >= 1.3.2
    using Pkg.Artifacts: @artifact_str
end


text:
#Demo 6: Smith chart for publication (labels etc removed - will be in caption).
#-------------------------------------------------------------------------------
using InspectDR
using Colors

#==Input
===============================================================================#

#Constants
#-------------------------------------------------------------------------------
black = RGB24(0, 0, 0)
white = RGB24(1, 1, 1)
red = RGB24(1, 0, 0)
green = RGB24(0, 1, 0)
blue = RGB24(0, 0, 1)

#Input data
#-------------------------------------------------------------------------------
f = collect(50e6:50e6:10e9)
C1 = 5e-12; GC1 = 100
C2 = .2e-12; RC2 = 10
L = 1.5e-9; RL = 15

#==Equations
===============================================================================#
Γ(Z; Zref::Real=50.0) = (Z .- Zref) ./ (Z .+ Zref)

#Calculations
#-------------------------------------------------------------------------------
ω = 2π*f; j = im
YC1 = (j*ω*C1) .+ 1/GC1
ZC1 = 1 ./ YC1
ZC2 = 1 ./ (j*ω*C2) .+ RC2
ZL = j*ω*L .+ RL

#==Generate plot
===============================================================================#
mplot = InspectDR.Multiplot(title="Smith Chart (Publication Layout)")
smithext = InspectDR.PExtents1D(min=-1.1,max=1.1) #Padded a bit

plot = add(mplot, InspectDR.smithchart(:Z, ref=50, title="", xlabel="", ylabels=String[]))
	strip = plot.strips[1]
	plot.xext_full = smithext; strip.yext_full = smithext
	#plot.layout[:enable_legend] = true
	plot.layout[:halloc_legend] = 110
	lyt = plot.layout

#Zero-out title areas, etc:
	lyt[:halloc_left] = 0
	lyt[:halloc_right] = 0
	lyt[:valloc_top] = 0
	lyt[:valloc_bottom] = 0

#Make data area sqare when saved:
	lyt[:halloc_data] = 500
	lyt[:valloc_data] = 500

	wfrm = add(plot, f, Γ(ZC1), id="cap: $(C1/1e-12)pF")
	wfrm.line = line(color=blue, width=3)
	wfrm = add(plot, f, Γ(ZC2), id="cap: $(C2/1e-12)pF")
	wfrm.line = line(color=red, width=3)
	wfrm = add(plot, f, Γ(ZL), id="ind: $(L/1e-9)nH")
	wfrm.line = line(color=green, width=3)

gplot = display(InspectDR.GtkDisplay(), mplot)
InspectDR.write_png("export_pubsmith.png", plot)
#InspectDR.write_svg("export_pubsmith.svg", plot)
#InspectDR.write_eps("export_pubsmith.eps", plot)
#InspectDR.write_pdf("export_pubsmith.pdf", plot)

:DONE


text:
#-------------------------------------------------------------------------------
# This function creates an instance of SNEASY+Hector.
#-------------------------------------------------------------------------------

# Load required packages.
using Mimi
using MimiHector
using MimiSNEASY
using DataFrames
using CSVFiles

function create_sneasy_hectorch4(;rcp_scenario::String="RCP85", start_year::Int=1765, end_year::Int=2300, etminan_ch4_forcing::Bool=true)

    # ---------------------------------------------
    # Load and clean up necessary data.
    # ---------------------------------------------

    # Load RCP emissions and concentration scenario values (RCP options = "RCP26" and "RCP85").
    rcp_emissions      = DataFrame(load(joinpath(@__DIR__, "..", "..", "data", "model_data", rcp_scenario*"_emissions.csv"), skiplines_begin=36))
    rcp_concentrations = DataFrame(load(joinpath(@__DIR__, "..", "..", "data", "model_data", rcp_scenario*"_concentrations.csv"), skiplines_begin=37))
    rcp_forcing        = DataFrame(load(joinpath(@__DIR__, "..", "..", "data", "model_data", rcp_scenario*"_midyear_radforcings.csv"), skiplines_begin=58))

    # Find start and end year indices to crop RCP scenario data to correct model time horizon.
    rcp_indices = findall((in)(collect(start_year:end_year)), rcp_emissions.YEARS)

    # Set pre-industrial atmospheric CO₂, CH₄, and N₂O concentrations to RCP values in 1765.
    CO₂_0 = rcp_concentrations[rcp_concentrations.YEARS .== 1765, :CO2][1]
    CH₄_0 = rcp_concentrations[rcp_concentrations.YEARS .== 1765, :CH4][1]
    N₂O_0 = rcp_concentrations[rcp_concentrations.YEARS .== 1765, :N2O][1]

    # Calculate carbon dioxide emissions as well as aerosol and exogenous RCP radiative forcing scenarios.
    rcp_co2_emissions     = rcp_emissions.FossilCO2 .+ rcp_emissions.OtherCO2
    rcp_aerosol_forcing   = rcp_forcing.TOTAER_DIR_RF .+ rcp_forcing.CLOUD_TOT_RF
    rcp_exogenous_forcing = rcp_forcing.TOTAL_INCLVOLCANIC_RF .- rcp_forcing.CO2_RF .- rcp_forcing.CH4_RF .- rcp_forcing.TROPOZ_RF .- rcp_forcing.CH4OXSTRATH2O_RF .- rcp_aerosol_forcing

    # ------------------------------------------------------------
    # Initialize Mimi-SNEASY and add new CH₄ components to model.
    # ------------------------------------------------------------

    # Get an instance of Mimi-SNEASY.
    m = MimiSNEASY.get_model(start_year=start_year, end_year=end_year)

    # Remove old radiative forcing components.
    delete!(m, :rfco2)
    delete!(m, :radiativeforcing)

    # Add in new components.
    add_comp!(m, rf_total, before = :doeclim)
    add_comp!(m, rf_co2_etminan, before = :rf_total)
    add_comp!(m, MimiHector.rf_ch4h2o, before = :rf_co2_etminan)
    add_comp!(m, MimiHector.rf_o3, before = :rf_ch4h2o)
    add_comp!(m, MimiHector.ch4_cycle, before = :rf_o3)
    add_comp!(m, MimiHector.oh_cycle, before = :ch4_cycle)

    # Add in user-specified CH₄ radiative forcing component.
    # Note: If not using Etminan et al. equations, use original forcing equations from parent CH₄ model.
    if etminan_ch4_forcing == true
        add_comp!(m, rf_ch4_etminan, before = :rf_ch4h2o)
    else
        add_comp!(m, MimiHector.rf_ch4, before = :rf_ch4h2o)
    end

    # ---------------------------------------------
    # Set component parameters.
    # ---------------------------------------------

    # ---- Common parameters ----
    Mimi.set_external_param!(m, :CH₄_0, CH₄_0)
    Mimi.set_external_param!(m, :N₂O_0, N₂O_0)
    Mimi.set_external_param!(m, :N₂O, rcp_concentrations.N2O[rcp_indices], param_dims=[:time])
    Mimi.set_external_param!(m, :M0, CH₄_0)
    Mimi.set_external_param!(m, :CO_emissions, rcp_emissions.CO[rcp_indices], param_dims=[:time])
    Mimi.set_external_param!(m, :NMVOC_emissions, rcp_emissions.NMVOC[rcp_indices], param_dims=[:time])

    # ---- Carbon Cycle ---- #
    update_param!(m, :CO2_emissions, rcp_co2_emissions[rcp_indices])
    update_param!(m, :atmco20, CO₂_0)

    # ---- Tropospheric Sink (OH) Lifetime ---- #
    set_param!(m, :oh_cycle, :CNOX, 0.0042)
    set_param!(m, :oh_cycle, :CCO, -0.000105)
    set_param!(m, :oh_cycle, :CNMVOC, -0.000315)
    set_param!(m, :oh_cycle, :CCH4, -0.32)
    set_param!(m, :oh_cycle, :NOX_emissions, rcp_emissions.NOx[rcp_indices])
    connect_param!(m, :oh_cycle, :CO_emissions, :CO_emissions)
    connect_param!(m, :oh_cycle, :NMVOC_emissions, :NMVOC_emissions)
    set_param!(m, :oh_cycle, :TOH0, 6.586)
    connect_param!(m, :oh_cycle, :M0, :M0)

    # ---- Methane Cycle ---- #
    set_param!(m, :ch4_cycle, :UC_CH4, 2.78)
    set_param!(m, :ch4_cycle, :CH4N, 300.)
    set_param!(m, :ch4_cycle, :Tsoil, 160.)
    set_param!(m, :ch4_cycle, :Tstrat, 120.)
    connect_param!(m, :ch4_cycle, :M0, :M0)
    set_param!(m, :ch4_cycle, :CH4_emissions, rcp_emissions.CH4[rcp_indices])

    # ---- Methane Radiative Forcing ---- #
    if etminan_ch4_forcing == true
        connect_param!(m, :rf_ch4_etminan, :CH₄_0, :CH₄_0)
        connect_param!(m, :rf_ch4_etminan, :N₂O_0, :N₂O_0)
        set_param!(m, :rf_ch4_etminan, :scale_CH₄, 1.0)
        set_param!(m, :rf_ch4_etminan, :a₃, -1.3e-6)
        set_param!(m, :rf_ch4_etminan, :b₃, -8.2e-6)
        connect_param!(m, :rf_ch4_etminan, :N₂O, :N₂O)
    else
        connect_param!(m, :rf_ch4, :N₂O_0, :N₂O_0)
        connect_param!(m, :rf_ch4, :CH₄_0, :CH₄_0)
        set_param!(m, :rf_ch4, :scale_CH₄, 1.0)
    end

    # ---- Straospheric Water Vapor From Methane Radiative Forcing ---- #
    connect_param!(m, :rf_ch4h2o, :M0, :M0)
    set_param!(m, :rf_ch4h2o, :H₂O_share, 0.05)

    # ---- Tropospheric Ozone Radiative Forcing ---- #
    set_param!(m, :rf_o3, :O₃_0, 32.38)
    set_param!(m, :rf_o3, :NOx_emissions, rcp_emissions.NOx[rcp_indices])
    connect_param!(m, :rf_o3, :CO_emissions, :CO_emissions)
    connect_param!(m, :rf_o3, :NMVOC_emissions, :NMVOC_emissions)

    # ---- Carbon Dioxide Radiative Forcing ---- #
    set_param!(m, :rf_co2_etminan, :a₁, -2.4e-7)
    set_param!(m, :rf_co2_etminan, :b₁, 7.2e-4)
    set_param!(m, :rf_co2_etminan, :c₁, -2.1e-4)
    set_param!(m, :rf_co2_etminan, :CO₂_0, CO₂_0)
    connect_param!(m, :rf_co2_etminan, :N₂O_0, :N₂O_0)
    connect_param!(m, :rf_co2_etminan, :N₂O, :N₂O)
    set_param!(m, :rf_co2_etminan, :rf_scale_CO₂, co2_rf_scale(3.7, CO₂_0, N₂O_0))

    # ---- Total Radiative Forcing ---- #
    set_param!(m, :rf_total, :α, 1.0)
    # TODO It would be nice if `rf_aerosol` wouldn't be a model
    # parameter at this point
    update_param!(m, :rf_aerosol, rcp_aerosol_forcing[rcp_indices])
    connect_param!(m, :rf_total, :rf_aerosol, :rf_aerosol)
    set_param!(m, :rf_total, :rf_exogenous, rcp_exogenous_forcing[rcp_indices])

    # ----------------------------------------------------------
    # Create connections between Mimi SNEASY+Hector components.
    # ----------------------------------------------------------
    connect_param!(m, :doeclim,        :forcing, :rf_total,  :total_forcing)
    connect_param!(m, :ccm,            :temp,    :doeclim,   :temp)
    connect_param!(m, :rf_co2_etminan, :CO₂,     :ccm,       :atmco2)
    connect_param!(m, :oh_cycle,       :CH4,     :ch4_cycle, :CH4)
    connect_param!(m, :ch4_cycle,      :TOH,     :oh_cycle,  :TOH)
    connect_param!(m, :rf_o3,          :CH₄,     :ch4_cycle, :CH4)
    connect_param!(m, :rf_ch4h2o,      :CH4,     :ch4_cycle, :CH4)

    # Create different connections if using updated Etminan et al. CH₄ forcing equations.
    if etminan_ch4_forcing == true
        connect_param!(m, :rf_ch4_etminan, :CH₄,    :ch4_cycle,      :CH4)
        connect_param!(m, :rf_total,       :rf_CH₄, :rf_ch4_etminan, :rf_CH₄)
    else
        connect_param!(m, :rf_ch4,   :CH4,    :ch4_cycle, :CH4)
        connect_param!(m, :rf_total, :rf_CH₄, :rf_ch4,    :rf_CH4)
    end

    connect_param!(m, :rf_total, :rf_O₃,      :rf_o3,          :rf_O₃)
    connect_param!(m, :rf_total, :rf_CH₄_H₂O, :rf_ch4h2o,      :rf_ch4h2o)
    connect_param!(m, :rf_total, :rf_CO₂,     :rf_co2_etminan, :rf_CO₂)

    # Return constructed model.
    return m
end


text:
# --- Day 1: Sonar Sweep ---
# You're minding your own business on a ship at sea when the overboard alarm
# goes off! You rush to see if you can help. Apparently, one of the Elves
# tripped and accidentally sent the sleigh keys flying into the ocean!
#
# Before you know it, you're inside a submarine the Elves keep ready for
# situations like this. It's covered in Christmas lights (because of course it
# is), and it even has an experimental antenna that should be able to track the
# keys if you can boost its signal strength high enough; there's a little meter
# that indicates the antenna's signal strength by displaying 0-50 stars.
#
# Your instincts tell you that in order to save Christmas, you'll need to get
# all fifty stars by December 25th.
#
# Collect stars by solving puzzles. Two puzzles will be made available on each
# day in the Advent calendar; the second puzzle is unlocked when you complete
# the first. Each puzzle grants one star. Good luck!
#
# As the submarine drops below the surface of the ocean, it automatically
# performs a sonar sweep of the nearby sea floor. On a small screen, the sonar
# sweep report (your puzzle input) appears: each line is a measurement of the
# sea floor depth as the sweep looks further and further away from the
# submarine.
#
# For example, suppose you had the following report:
#
# 199
# 200
# 208
# 210
# 200
# 207
# 240
# 269
# 260
# 263
#
# This report indicates that, scanning outward from the submarine, the sonar
# sweep found depths of 199, 200, 208, 210, and so on.
#
# The first order of business is to figure out how quickly the depth increases,
# just so you know what you're dealing with - you never know if the keys will
# get carried into deeper water by an ocean current or a fish or something.
#
# To do this, count the number of times a depth measurement increases from the
# previous measurement. (There is no measurement before the first measurement.)
# In the example above, the changes are as follows:
#
# 199 (N/A - no previous measurement)
# 200 (increased)
# 208 (increased)
# 210 (increased)
# 200 (decreased)
# 207 (increased)
# 240 (increased)
# 269 (increased)
# 260 (decreased)
# 263 (increased)
#
# In this example, there are 7 measurements that are larger than the previous
# measurement.
#
# How many measurements are larger than the previous measurement?

readings = map(l -> parse(Int, l), eachline("input.txt"))
lead = @view readings[2:end]
println("Count: ", count(((a, b),) -> a < b, zip(readings, lead)))


text:
"""
    function disp_plots(args...; kwargs)

This function plots the output results from "plots_TCRDD.jl"

#Arguments
- `bids::Array{Float64, 1}`                     Bids
- `gen_profit::Array{Float64, 1}`               Generator Profit
- `gen_profit_tcrdd::Array{Float64, 1}`         Approximated Generator Profit
- `gen_profit_noCost::Array{Float64, 1}`        Generator Profit no Cost
- `gen_profit_tcrdd_noCost::Array{Float64, 1}`  Approximated Generator Profit no Cost
- `gen_cost_slack::Array{Float64, 1}`           Slack Generator Cost
- `Pg_slack::Array{Float64, 1}`                 Slack Generator Active Power Output
- `lmp_slack::Array{Float64, 1}`                Slack Bus LMP
- `lmp_slack_apx_tcrdd::Array{Float64, 1}`      Slack Bus Approximated LMP
- `residualD_slack_Pg::Array{Float64, 1}`       Slack Generator Residual Demand
- `tcrdd::Array{Float64, 1}`                    Transmission Constrained Residual Demand
- `nbind_lines::Array{Int64, 1}`                Number of binding lines
- `div::Int64`                                  Number of divisions for the samples

#Keywords
- `examplecase4::Bool = false`                  If Case 4 example, provide an extra plot
- `branch1_3::Array{Float64, 1}`                Active Power Flow for Branch (Bus1 to Bus3)

"""
function disp_plots(
    bids::Array{Float64, 1},
    gen_profit::Array{Float64, 1},
    gen_profit_tcrdd::Array{Float64, 1},
    gen_profit_noCost::Array{Float64, 1},
    gen_profit_tcrdd_noCost::Array{Float64, 1},
    gen_cost_slack::Array{Float64, 1},
    Pg_slack::Array{Float64, 1},
    lmp_slack::Array{Float64, 1},
    lmp_slack_apx_tcrdd::Array{Float64, 1},
    residualD_slack_Pg::Array{Float64, 1},
    tcrdd::Array{Float64, 1},
    nbind_lines::Array{Int64, 1},
    div::Int64;
    examplecase4::Bool = false,
    branch1_3::Array{Float64, 1}
    )
    # Plot using Plotly
    plotly()

    # Slack Generator Profit
    gen_profit_both = Array{Float64}(undef, (div + 1, 2));
    gen_profit_both[:, 1] = gen_profit;
    gen_profit_both[:, 2] = gen_profit_tcrdd;
    plot(bids, gen_profit_both , title = "Generator Profit",
        label = ["Profit Function Real" "Profit Function Aprox TCRDD"], legend = :topleft);
    xlabel!("Bid q(pu)");
    a = ylabel!("profit(\$/hr)")
    display(a)

    # Slack Generator Profit No Cost
    gen_profit_both_noCost = Array{Float64}(undef, (div + 1, 2));
    gen_profit_both_noCost[:, 1] = gen_profit_noCost;
    gen_profit_both_noCost[:, 2] = gen_profit_tcrdd_noCost;
    plot(bids, gen_profit_both_noCost , title = "Generator Profit no Cost",
        label = ["Profit Function Real_noCost" "Profit Function Aprox TCRDD_noCost"],
        legend = :bottomright);
    xlabel!("Bid q(pu)");
    a = ylabel!("profit(\$/hr)")
    display(a)

    # Slack Generator Cost
    plot(bids, gen_cost_slack, title = "Generator Cost",
        label = ["Cost Function Real" "Cost Function TCRDD"], legend = :bottomright);
    xlabel!("Bid q(pu)");
    a = ylabel!("Cost(\$/hr)")
    display(a)

    # Slack Generator Active Power Output
    plot(bids, Pg_slack, title = "Dispatched vs Bid on Slack",
        label = ["Dispatched vs Bid Slack"], legend = :bottomright);
    xlabel!("Bid q(pu)");
    a = ylabel!("P Gen (pu)")
    display(a)

    # Slack Bus LMP
    lmp_both = Array{Float64}(undef,(div + 1, 2));
    lmp_both[:, 1] = lmp_slack;
    lmp_both[:, 2] = lmp_slack_apx_tcrdd;
    plot(bids, lmp_both, title = "LMPs of slack vs bids",
        label = ["LMP Real" "LMP approx TCRDD"], legend = :topright);
    xlabel!("Bid q(pu)");
    a = ylabel!("LMP (\$/hrpu)")
    display(a)

    # Residual Demand Curve
    plot(bids, residualD_slack_Pg, title = " Residual Demand Curve vs Bids",
        label = ["RDC"], legend = :topright);
    xlabel!("Bid q(pu)");
    a = ylabel!("RDC (pu)")
    display(a)

    # TCRDD
    plot(bids, tcrdd ,title = "TCRDD vs bids",
        label = ["TCRDD vs Bid Slack"], legend = :bottomright);
    xlabel!("Bid q(pu)");
    a = ylabel!("TCRDD (pu)")
    display(a)

    # RDC & TCRDD vs Bid
    rdc_tcrdd = Array{Float64}(undef,(div + 1, 2));
    rdc_tcrdd[:, 1] = residualD_slack_Pg;
    rdc_tcrdd[:, 2] = tcrdd;
    plot(bids, rdc_tcrdd, title = "RDC & TCRDD vs Bid",
        label = ["RDC" "TCRDD"], legend = :topright);
    xlabel!("Bid q(pu)");
    a = ylabel!("RDC & TCRDD")
    display(a)

    # Residual Demand Curve vs LMP
    plot(lmp_slack, residualD_slack_Pg, title = " Residual Demand Curve vs LMP",
        label = ["RDC"], legend = :topright);
    xlabel!("LMP (\$/hrpu)");
    a = ylabel!("RDC (pu)")
    display(a)

    # TCRDD vs LMP Slack
    plot(lmp_slack, tcrdd, title = "TCRDD vs LMP", label = ["TCRDD vs LMP Slack"],
        legend = :bottomright);
    xlabel!("LMP (\$/hrpu)");
    a = ylabel!("TCRDD (pu)")
    display(a)

    # RDC & TCRDD vs LMP Slack
    rdc_tcrdd = Array{Float64}(undef,(div + 1, 2));
    rdc_tcrdd[:, 1] = residualD_slack_Pg;
    rdc_tcrdd[:, 2] = tcrdd;
    plot(lmp_slack, rdc_tcrdd, title = "RDC & TCRDD vs LMP",
        label = ["RDC" "TCRDD"], legend = :topright);
    xlabel!("LMP (\$/hrpu)");
    a = ylabel!("RDC & TCRDD")
    display(a)

    # Number of Binding lines
    plot(bids, nbind_lines, title = "Number of Binding lines", label = ["Binding Lines"],
        legend = :bottomright);
    xlabel!("Bid q(pu)");
    a = ylabel!("Number binding lines")
    display(a)

    if examplecase4
        # Only for example case4 Pf branch Bus1 to Bus3
        plot(bids, branch1_3, title = "Pf branch Bus1 to Bus3",
            label = ["Pf line 1 to 3"], legend = :bottomright);
        xlabel!("Bid q(pu)");
        a = ylabel!("Active Power Flow Pf (pu)")
        display(a)
    end
end

max_stars_repo_name: invenia/MaximumGeneratorProfitTCRDD

max_issues_repo_name: invenia/MaximumGeneratorProfitTCRDD

max_forks_repo_name: invenia/MaximumGeneratorProfitTCRDD


text:
export Sync

struct Sync
    xs::AbstractArray
end

function islocal(x)
    !any(startswith.(x, ("//", "https://", "http://", "ftp://")))
end

function path2url(path::AbstractString)
    if startswith(path, "/pkg/")
        @warn("/pkg/ URLs are deprecated, load files with their absolute path in Scope")
        return path
    elseif isfile(abspath(path))
        path = abspath(path)
        # first lookup to see if any of the file itself or any of the parent
        # directories are registered.
        AssetRegistry.isregistered(path) && AssetRegistry.getkey(path)
        cur_path = path
        while true
            if AssetRegistry.isregistered(cur_path) && isdir(cur_path)
                key = AssetRegistry.getkey(cur_path)
                url = key * "/" * replace(path, cur_path => "")
                return url
            end
            cur_path1 = dirname(cur_path)
            if cur_path1 == cur_path
                # this means we have reached root directory,
                # and none of the parents are in registry
                # register the original path uniquely
                return AssetRegistry.register(path)
            end
            cur_path = cur_path1
        end
    else
        error("File $path not found")
    end
end

function dep2url(dep::AbstractString)
    # if is an url, we are done :)
    islocal(dep) || return dep
    query_parts = split(dep, "?") # remove anything after ?
    file_path = first(query_parts)
    query_part = length(query_parts) == 2 ? query_parts[2] : ""
    url = path2url(file_path)
    return string(baseurl[], url, query_part)
end

function lowerdeps(name, imp)
    url = dep2url(imp)
    extension = split(url, '.')[end]
    if !(extension in ("js", "css", "html"))
        error("WebIO can't load dependency of unknown type $url")
    end
    return Dict{String,Any}(
        "type" => extension,
        "name" => name,
        "url" => url
    )
end

lowerdeps(x::String) = lowerdeps(nothing, x)
lowerdeps(x::Pair) = lowerdeps(x[1], x[2])
function lowerdeps(x::Dict)
    if haskey(x, "url")
        key = "url"
    elseif haskey(x, :url)
        key = :url
    else
        error("Dict import specification doesn't have key :url")
    end

    url = x[key]
    dict = lowerdeps(nothing, url)
    x1 = copy(x)
    x1[key] = dict["url"]
    x1
end

lowerdeps(xs::AbstractArray) = Dict(
    "type" => "async_block",
    "data" => map(lowerdeps, xs),
)

lowerdeps(x::Sync) = Dict("type"=>"sync_block",
                          "data" => map(lowerdeps, x.xs))

function import!(scope, x)
    push!(scope.imports, x)
end

Base.@deprecate adddeps!(scope, x) import!(scope, x)


text:
using Plots
include("unrate_btg.jl")
include("../../computation/tdist_mle.jl")

path = "../../datasets/"
data, target = unrate(path)

btg0 = load_unrate_btg()
(pdf0_raw, cdf0_raw, dpdf0_raw, quantInfo0_raw) = solve(btg0);
f = theta -> tdist_mle(btg0, theta, 0.1)
g(theta, lambda) = tdist_mle(btg0, theta, lambda)
h = lambda -> tdist_mle(btg0, 130.0, lambda)
####################################
########## MLE Estimation ##########
####################################

####################################
######## Optional Plotting #########
####################################
if true
    Plots.plot()
    plt(y -> fzero(x -> cdf0_raw(reshape([y], 1, 1), reshape([1.0], 1, 1), x) - 0.5, 0, 2), 0, 1)
    plt!(y -> fzero(x -> cdf0_raw(reshape([y], 1, 1), reshape([1.0], 1, 1), x) - 0.75, 0, 2), 0, 1)
    plt!(y -> fzero(x -> cdf0_raw(reshape([y], 1, 1), reshape([1.0], 1, 1), x) - 0.25, 0, 2), 0, 1)
end

####################################
############### Test ###############
####################################
if true
    @info "Start Test"
    before = Dates.now()
    count_test = 0
    error_abs = 0.
    error_sq = 0.
    nlpd = 0.
    id_fail = []
    id_nonproper = []
    CI_set = zeros(n_test, 2)
    median_set = zeros(n_test)
    for i in 1:n_test
        global error_abs, error_sq, nlpd, count_test
        # mod(i, 20) == 0 ? (@info i) : nothing
        # @info "i" i
        x_test_i = reshape(x_test[i, :], 1, d)
        Fx_test_i = reshape(Fx_test[i, :], 1, p)
        try
            #pdf_test_i, cdf_test_i, dpdf_test_i, quantbound_test_i, support_test_i = pre_process(x_test_i, Fx_test_i, pdf0_raw, cdf0_raw, dpdf0_raw, quantInfo0_raw)
            y_test_i = y_test[i]

            #median_test_i = max_train * quantile(cdf_test_i, quantbound_test_i, support_test_i)[1]
            median_test_i = fzero(x -> cdf0_raw(x_test_i, reshape([1.0], 1, 1), x) - 0.5, 0, 10)
            # @info "True, median " y_test_i_true, median_test_i
            median_set[i] = median_test_i
            try
                #CI_test_i = max_train .* credible_interval(cdf_test_i, quantbound_test_i, support_test_i; mode=:equal, wp=.95)[1]
                aa = fzero(x -> cdf0_raw(x_test_i, reshape([1.0], 1, 1), x) - 0.975, 0, 20)
                bb = fzero(x -> cdf0_raw(x_test_i, reshape([1.0], 1, 1), x) - 0.025, 0, 20)
                CI_test_i = [bb, aa]
                count_test += (y_test_i >= CI_test_i[1])&&(y_test_i <= CI_test_i[2]) ? 1 : 0
                # @info "CI" CI_test_i
                #@info "95% confidence interval:", CI_test_i
                CI_set[i, :] = CI_test_i
            catch err
                append!(id_fail, i)
            end
            error_abs += abs(y_test_i - median_test_i)
            error_sq += (y_test_i - median_test_i)^2
            nlpd += log(pdf0_raw(y_test_i))
        # @info "Count, id_fail" count_test, id_fail
        catch err
            append!(id_nonproper, i)
        end
    end
        count_test /= n_test
        error_abs  /= n_test
        error_sq   /= n_test
        nlpd       /= -n_test
        after = Dates.now()
        elapsedmin = round(((after - before) / Millisecond(1000))/60, digits=5)
        println("error_sq")
        display(error_sq)
        println("count_test")
        display(count_test)
        println("elapsedmin")
        display(elapsedmin)
        # Plot
        PyPlot.close("all") #close existing windows
        PyPlot.plot(x_test, median_set, label = "BTG median")
        PyPlot.plot(xs, label, label = "true")
        PyPlot.fill_between(dropdims(x_test; dims = 2), CI_set[:, 1], CI_set[:, 2], alpha = 0.3, label = "95% confidence interval")
        PyPlot.scatter(x_train, y_train, s = 10, c = "k", marker = "*")
        #PyPlot.scatter(x_test, y_test, s = 10, c = "k", marker = "*")
        PyPlot.legend(fontsize=8)
        PyPlot.grid()
        PyPlot.title("BTG $myquadtype", fontsize=10)
        PyPlot.savefig("exp_synthetic_mle_btg4.pdf")
    end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
# Raytracer.jl
# Raytracing for the generation of photorealistic images in Julia
# Copyright (c) 2021 Samuele Colombo, Paolo Galli

# Main package file

"""
    Raytracer

Raytracing package for the generation of photorealistic images in Julia.
"""
module Raytracer

##########
# Imports

import Base:
    (+), (-), (*), (≈), (==),
    Matrix, OneTo,
    axes, clamp, convert, eltype, fill!, firstindex, getindex, iterate,
    lastindex, length, one, print_matrix, rand, read, readline, setindex!, show,
    size, write, zero,
    union, intersect, setdiff,
    isless, isequal

import Base.Broadcast:
    BroadcastStyle, Broadcasted, Style,
    broadcastable, combine_eltypes, copy, similar

using ColorTypes:
    RGB, Fractional

using FileIO:
    File, @format_str, query
import FileIO:
    save, load

import ImagePFM:
    _read

import StaticArrays:
    @SMatrix,
    FieldVector, MMatrix, SMatrix, SVector, Size,
    similar_type

import LinearAlgebra:
    (⋅), (×),
    Diagonal, I,
    inv, norm, normalize

using Random
import Random:
    Random.CloseOpen01, Sampler, SamplerTrivial

using Intervals, ProgressMeter

##########
# Exports

export # Rendering
    RGB,
        luminosity, clamp, γ_correction,
    HdrImage,
        luminosity, clamp, γ_correction,
    Camera,
        PerspectiveCamera, OrthogonalCamera,
        aperture_deg,
        fire_ray,
    ImageTracer,
        fire_all_rays!,
    Ray,
    HitRecord,
        HitOrMiss,
    Pigment,
        UniformPigment,
        CheckeredPigment,
        ImagePigment,
    BRDF,
        DiffuseBRDF, SpecularBRDF,
        at,
    Material,
    Renderer,
        OnOffRenderer, FlatRenderer, PathTracer, PointLightRenderer

export # Scene
    Normal, Vec,
        normalize,
        norm,
        norm²,
        create_onb_from_z,
        NORMAL_X, NORMAL_Y, NORMAL_Z,
        VEC_X, VEC_Y, VEC_Z,
        NORMAL_X_false, NORMAL_Y_false, NORMAL_Z_false,
    Point,
        convert,
        ORIGIN,
    Vec2D,
    Transformation,
        isconsistent,
        rotationX, rotationY, rotationZ,
        scaling, translation,
    Shape,
        Sphere, Plane, Cube, Cylinder, Cone,
        CSG,
            UnionCSG,
            IntersectionCSG,
            DiffCSG,
            FusionCSG,
        fuse,
        Rule,
            UniteRule,
            IntersectRule,
            DiffRule,
            FuseRule,
        get_t, get_all_ts, get_uv, get_normal,
        ray_intersection, all_ray_intersections, quick_ray_intersection,
    World,
        is_point_visible,
    PointLight,
    Lights

export # Random number generation
    PCG

export # High level API
    tonemapping, render, render_from_script

export # image tools
    clamp_image, normalize_image, γ_correction, load, save

export # Colors
    BLACK, WHITE,
    RED, GREEN, BLUE,
    CYAN, MAGENTA, YELLOW

###########
# Includes

include("colors.jl")
include("hdrimage.jl")

include("pcg.jl")

include("geometry.jl")
include("transformations.jl")
include("ray.jl")
include("cameras.jl")
include("materials.jl")
include("lights.jl")
include("hitrecord.jl")
include("shapes.jl")
include("world.jl")
include("renderers.jl")
include("imagetracer.jl")

include("interpreter.jl")

using .Interpreter:
    print_subsequent_lexer_exceptions,
    Scene,
    open_stream,
    parse_variables_from_string,
    InterpreterException,
    parse_scene,
    SourceLocation,
    UndefinedSetting

include("user_utils.jl")

end # module Raytracer


text:
function new_cdatanode(xdoc::XMLDocument, txt::AbstractString)
    p = ccall((:xmlNewCDataBlock,libxml2), Xptr,
              (Xptr, Xstr, Cint),
              xdoc.ptr, txt, length(txt) + 1)
    XMLNode(p)
end

add_cdata(xdoc::XMLDocument, x::XMLElement, txt::AbstractString) =
    add_child(x, new_cdatanode(xdoc, txt))


text:
function rilt(t,y,smin,smax,N,α=missing;kwargs...)
    s = smin * (smax / smin).^range(0, 1, length=N + 1)
   invert(t,y,s,(t,s) -> exp(-t*s),Tikhonov(α),kwargs...)
end
function invert(t,y,s,kernel::Function,aamethod::AutoRegMethod;kwargs...)
    reg = auto_reg(t,y,s,kernel::Function,aamethod;kwargs...)
    invert(t,y,s,kernel,reg;kwargs...)
end

"""
    invert(s,y,t,k::Function,reg;yoffset=true,tdomain=:realplus)

Compute the discretized form a(s) in y(t) = ∫a(s)k(t,s)ds. `s` and `y` represents
the data wich we want to invert on discrete points `t`.

# Examples
```julia-repl
julia> bar([1, 2], [1, 2])
1
```
"""

function invert(s,y,t,k::Function,reg;yoffset=true,kwargs...)
    AR,yr = build_ar(s,y,t,k,reg,yoffset)
    yt=get_yt(AR,yr;kwargs...)
    reshape
    return (t,yt[1:length(t),1],s,(AR*yt)[1:length(s),1])

end

function get_yt(AR,yr;tdomain=:realplus)

    if tdomain ==:real
        sy = (AR \ yr)
        return  sy
    end
    return nonneg_lsq(AR,yr;alg=:fnnls)

end

function build_ar(t,y,s,kernel,reg,yoffset)
    #TODO: whole function does allocate like hell when AutoRegMethod is present
    #   cacheing must be implemented
    # central points

    sc=  (s[1:end - 1] + s[2:end]) / 2
    N=length(sc)
    ds = diff(s)
    dd=ones(length(ds))
    L = get_regularization_matrix(reg,N,dd)

    #Build Matrix approximation of kernel
    sType = eltype(s)
    lType = eltype(L)
    aType =promote_type(sType,lType,eltype(t),eltype(ds))
    A = aType[kernel(i, sc[j])*ds[j] for i in t, j in eachindex(sc)]

    if reg.α isa AbstractArray
        L= L*(Diagonal(reg.α))
    else
        L= L*(Diagonal(reg.α*ones(aType,length(sc))))
    end

    if yoffset
        #add a ofset column
        A= hcat(A,ones(size(A,1)))
        #Do not regularize the y-offset
        L= hcat(L,zeros(size(L,1)))
       # push!(ds, )
        push!(sc, 0)
    end

    AR = vcat(A, L)

    # add  entry to y to store the regularization
    yr = vcat(y, zeros(size(L,1 )))
    yr=convert.((eltype(AR),),yr) #allows for automatic dfferentation  NNLS

    return (AR,yr)
end


text:
using Pkg.Artifacts
using Pkg.GitTools
using Pkg.PlatformEngines
using Pkg.PlatformEngines: sha256
using Pkg

Pkg.PlatformEngines.probe_platform_engines!()

function bind!(toml, item; force=false)
    h = artifact_hash(item.name, toml)
    if force && (h != nothing) && ispath(artifact_path(h))
        rm(path, recursive=true)
    end
    h_tar = nothing
    if h == nothing
        h = create_artifact() do artifact_dir
            @info "Creating artifact $item"
            tarball = download(item.url, joinpath(tempdir(), "$(item.name).tar.bz2"))
            try
                h_tar = open(tarball) do io
                    bytes2hex(sha256(io))
                end
                unpack(tarball, artifact_dir)
            finally
                rm(tarball)
            end
        end
        bind_artifact!(toml, item.name, h;
                   download_info=[(item.url, h_tar)],
                   lazy=false,
                   force=true)
        @info "At path $(artifact_path(h))"
    end
    @info "Skipping existing artifact $item"
    h
end

function delete!(toml, item)
    h = artifact_hash(item.name, toml)
    path = artifact_path(h)
    @info "removing $path"
    @info "ispath(path) = $(ispath(path))"
    rm(artifact_path(h), recursive=true, force=true)
    @info "ispath(path) = $(ispath(path))"
end

toml = joinpath(@__DIR__, "..", "Artifacts.toml")
for item in [
             (name="dhall", url="https://github.com/dhall-lang/dhall-haskell/releases/download/1.30.0/dhall-1.30.0-x86_64-linux.tar.bz2"),
             (name="dhall-json", url="https://github.com/dhall-lang/dhall-haskell/releases/download/1.30.0/dhall-json-1.6.2-x86_64-linux.tar.bz2"),
            # (name="dhall-bash", url="https://github.com/dhall-lang/dhall-haskell/releases/download/1.30.0/dhall-bash-1.0.28-x86_64-linux.tar.bz2"),
             (name="dhall-yaml", url="https://github.com/dhall-lang/dhall-haskell/releases/download/1.30.0/dhall-yaml-1.0.2-x86_64-linux.tar.bz2")
            ]
    bind!(toml, item)
    delete!(toml, item)
end


text:
### A Pluto.jl notebook ###
# v0.12.4

using Markdown
using InteractiveUtils

# ╔═╡ b80a0fc0-0675-11eb-1fcd-277d5ec1c289
using Pkg, DrWatson

# ╔═╡ b80a447c-0675-11eb-319c-01dc82d3fb7c
begin
	@quickactivate "SR2TuringPluto"
	using Turing
	using StatisticalRethinking
end

# ╔═╡ b221c3e6-0675-11eb-1ea1-4307128c6de3
md"## Clip-04-07-25t.jl"

# ╔═╡ b80ac79e-0675-11eb-2e70-17f8aa994dc0
md"## Clip-04-97-25t.jl"

# ╔═╡ b81d124e-0675-11eb-16e7-91a1a32f0369
md"### snippets 4.7 - 4.11"

# ╔═╡ b81df8d2-0675-11eb-33ad-819faf448b8a
begin
	df = CSV.read(sr_datadir("Howell1.csv"), DataFrame)
	df = df[df.age .>= 18, :]
end;

# ╔═╡ b82e9b88-0675-11eb-27b5-69985a392610
md"### snippets 4.12, 4.13"

# ╔═╡ b82f4f42-0675-11eb-3f77-d9e53da11c70
begin
	fig1 = plot(p -> pdf(Normal(178, 20), p), 100:250, leg=false, title="prior mu")
	fig2 = plot(p -> pdf(Uniform(0, 50), p), -10:60, leg=false, title="prior sigma")
	plot(fig1, fig2, layout=(1, 2))
end

# ╔═╡ b83efd16-0675-11eb-2e8d-fd33be380806
md"### snippet 4.14"

# ╔═╡ b84a549a-0675-11eb-08c9-b908f9639415
begin
	N = 1000
	d_μ = Normal(178, 20)
	sample_mu = rand(d_μ, N)
	d_σ = Uniform(0, 50)
	sample_sigma = rand(d_σ, N)
	prior_h = rand.(Normal.(sample_mu, sample_sigma))
	density(prior_h)
end

# ╔═╡ b85bf01a-0675-11eb-17e9-67c60998b96b
md"### snippet 4.15"

# ╔═╡ b8664768-0675-11eb-07de-555dbc1de52b
begin
	d_μ_100 = Normal(178, 100)
	sample_μ_100 = rand(d_μ_100, N)
end;

# ╔═╡ b86fd274-0675-11eb-2564-dfb336f2f5e6
begin
	sample_height_100 = rand.(Normal.(sample_μ_100, sample_sigma))
	density(sample_height_100)
end

# ╔═╡ 97fbe892-0676-11eb-0a85-6d793b380f4e
md"### snippets 4.16, 4.17"

# ╔═╡ 9b8b89fe-0676-11eb-1eaa-e986dcba4794
begin
	μ_grid_1 = range(150, 160, length = 50)
	σ_grid_1 = range(7, 9, length = 50)
	post_1 = [(; :μ => μ, :σ => σ) for μ in μ_grid_1, σ in σ_grid_1]  # this is `post`
	ll_1 = [loglikelihood(Normal(μ, σ), df.height) for (μ, σ) in post_1]
	prior_1 = [logpdf(d_μ, μ) + logpdf(d_σ, σ) for (μ, σ) in post_1]
	lposterior_1 = ll_1 .+ prior_1
	lposterior_1 .-= maximum(lposterior_1)
	posterior_1 = exp.(lposterior_1)
	wireframe(μ_grid_1, σ_grid_1, posterior_1)
end

# ╔═╡ b897b90a-0676-11eb-1fe9-cbd663edb4f0
md"### snippets 4.19 - 4.22"

# ╔═╡ c3a4dba2-0676-11eb-07fc-bfb72963914b
begin
	samples = sample(post_1, pweights(posterior_1), N) |> DataFrame
	scatter(samples.μ, samples.σ, alpha = 0.05, ms = 5, xlabel = "μ", ylabel = "σ")
end

# ╔═╡ fdacb630-0676-11eb-0217-f1e07c127aa2
begin
	fig3 = density(samples.μ, leg=false, title="samples.mu")
	fig4 = density(samples.σ, leg=false, title="samples.sigma")
	plot(fig3, fig4, layout=(1,2))
end

# ╔═╡ 643e622c-0677-11eb-187b-033fe474252d
quantile(samples.μ, (0.1, 0.9))

# ╔═╡ 6eed8892-0677-11eb-3c04-e9cb30ea2b66
quantile(samples.σ, (0.1, 0.9))

# ╔═╡ 776436b0-0677-11eb-19bf-bd953d980210
md"### snippets 4.23 - 4.25"

# ╔═╡ 9ba86ca0-0677-11eb-3304-118f06fda2f0
begin
	d3_height = df.height[1:20]
	μ_grid = range(140, 160, length = 50)
	σ_grid = range(4, 20, length = 50)
	post = [(μ = μ, σ = σ) for μ in μ_grid, σ in σ_grid]
	ll = [loglikelihood(Normal(μ, σ), d3_height) for (μ, σ) in post]
	prior = [logpdf(d_μ, μ) + logpdf(d_σ, σ) for (μ, σ) in post]
	lposterior = ll .+ prior
	lposterior .-= maximum(lposterior)
	posterior = exp.(lposterior)
	wireframe(μ_grid, σ_grid, posterior)
end

# ╔═╡ 8613457c-0675-11eb-10bd-217d514940ef
begin
	samples2 = sample(post, pweights(posterior), N) |> DataFrame
	scatter(samples2.μ, samples2.σ, alpha = 0.05, ms = 5)
end

# ╔═╡ d8ab7426-067a-11eb-0bac-1bb4747061b1
begin
	fig5 = density(samples2.μ, leg=false, title="samples2.mu")
	fig6 = density(samples2.σ, leg=false, title="samples2.sigma")
	plot(fig1, fig2, fig3, fig4, fig5, fig6, layout=(3,2))
end

# ╔═╡ d8d3d31c-067a-11eb-2a92-ef47acab4194
md"## End of clip-04-07-25t.jl"

# ╔═╡ Cell order:
# ╠═b221c3e6-0675-11eb-1ea1-4307128c6de3
# ╠═b80a0fc0-0675-11eb-1fcd-277d5ec1c289
# ╠═b80a447c-0675-11eb-319c-01dc82d3fb7c
# ╟─b80ac79e-0675-11eb-2e70-17f8aa994dc0
# ╟─b81d124e-0675-11eb-16e7-91a1a32f0369
# ╠═b81df8d2-0675-11eb-33ad-819faf448b8a
# ╟─b82e9b88-0675-11eb-27b5-69985a392610
# ╠═b82f4f42-0675-11eb-3f77-d9e53da11c70
# ╟─b83efd16-0675-11eb-2e8d-fd33be380806
# ╠═b84a549a-0675-11eb-08c9-b908f9639415
# ╟─b85bf01a-0675-11eb-17e9-67c60998b96b
# ╠═b8664768-0675-11eb-07de-555dbc1de52b
# ╠═b86fd274-0675-11eb-2564-dfb336f2f5e6
# ╟─97fbe892-0676-11eb-0a85-6d793b380f4e
# ╠═9b8b89fe-0676-11eb-1eaa-e986dcba4794
# ╟─b897b90a-0676-11eb-1fe9-cbd663edb4f0
# ╠═c3a4dba2-0676-11eb-07fc-bfb72963914b
# ╠═fdacb630-0676-11eb-0217-f1e07c127aa2
# ╠═643e622c-0677-11eb-187b-033fe474252d
# ╠═6eed8892-0677-11eb-3c04-e9cb30ea2b66
# ╟─776436b0-0677-11eb-19bf-bd953d980210
# ╠═9ba86ca0-0677-11eb-3304-118f06fda2f0
# ╠═8613457c-0675-11eb-10bd-217d514940ef
# ╠═d8ab7426-067a-11eb-0bac-1bb4747061b1
# ╟─d8d3d31c-067a-11eb-2a92-ef47acab4194


text:
# TODO: better names for the functions
using ForwardDiff

function updatecoeffs!{T<:Real}(coeff0, coeff1, coeff2, rhs, model, v, t::T, x,
                                a::Vector{T}, Δτ::T, Δx::T)
    # Updates
    # coeffn = values in linear system
    # rhs    = value function at previous timestep + f at current timestep

    # Input
    # model  = HJBOneDim object
    # v      = value function at previous timestep
    # t      = value of (forward) time
    # x      = vector of x-values
    # a      = policy-values on interior
    # Δτ     = time-step size
    # Δx     = spacial step length
    taux = Δτ/Δx
    htaux2 = 0.5*Δτ/Δx^2

    n = length(coeff0)
    for j = 2:n-1
        @inbounds begin
            bval = model.b(t,x[j],a[j-1])
            sval2 = model.σ(t,x[j],a[j-1])^2
            coeff1[j] = -(sval2*htaux2 + max(bval,0.)*taux)
            coeff2[j-1] = -(sval2*htaux2 - min(bval,0.)*taux)
            coeff0[j] = 1.-coeff1[j]-coeff2[j-1]
            rhs[j] = v[j] + Δτ*model.f(t,x[j],a[j-1])
        end
    end
end

function updatepol!(pol, v, model::HJBOneDim, t, x, Δx)
    # Loops over each x[i] value and optimises the control
    # TODO: Should we instead optimize the whole control vector
    # by considering the sum of the individual objectives
    # (the gradient should be diagonal?)

    idx = 1./Δx
    hdx2 = 0.5/Δx^2

    let
        # Search for control on left boundary
        function objective(a)
            @inbounds begin
                bval = idx*model.b(t,x[1],a[1])
                sval2 = hdx2*model.σ(t,x[1],a[1])^2
                coeff1 = -(sval2-bval)
                coeff2 = -(bval-2*sval2)
                coeff3 = -sval2
                return coeff1*v[1] + coeff2*v[2] + coeff3*v[3] - model.f(t,x[1],a[1])
            end
        end
        res = optimize(objective, model.amin, model.amax)
        @inbounds pol[1] = Optim.minimizer(res)
    end

    let
        # Search for control on right boundary
        function objective(a)
            @inbounds begin
                bval = idx*model.b(t,x[end],a[1])
                sval2 = hdx2*model.σ(t,x[end],a[1])^2
                coeff1 = -(sval2+bval)
                coeff2 = bval+2*sval2
                coeff3 = -sval2
                return coeff3*v[end-2] + coeff2*v[end-1] + coeff1*v[end] - model.f(t,x[end],a[1])
            end
        end
        res = optimize(objective, model.amin, model.amax)
        @inbounds pol[end] = Optim.minimizer(res)
    end

    function hamiltonian(a, j::Int)
        @inbounds begin
            bval = model.b(t,x[j],a[1])
            sval2 = model.σ(t,x[j],a[1])^2
            coeff1 = sval2*hdx2 + max(bval,0.)*idx
            coeff2 = sval2*hdx2 - min(bval,0.)*idx
            return coeff1*(v[j]-v[j+1]) + coeff2*(v[j]-v[j-1]) - model.f(t,x[j],a[1])
        end
    end

    for j = 2:length(pol)-1
        objective(a) = hamiltonian(a, j)
        # g!(x, out) = ForwardDiff.gradient!(out, objective, x)
        # diffobj = DifferentiableFunction(objective, g!)

        # res = optimize(diffobj, [pol[j]], [model.amin], [model.amax],
        #                Fminbox(), optimizer=LBFGS,
        #                optimizer_o = OptimizationOptions(g_tol=1e-3,f_tol=1e-3,x_tol=1e-3))#,
        #@inbounds pol[j] = Optim.minimizer(res)[1]
        # TODO: Give the choice of Brent vs Fminbox()?
        res = optimize(objective, model.amin, model.amax)
        @inbounds pol[j] = Optim.minimizer(res)
    end
end

function policynewtonupdate{T<:Real}(model::HJBOneDim{T},
                                     v, a, x,
                                     Δx, Δτ, ti::Int;
                                     tol = 1e-3,
                                     scale = 1.0,
                                     maxpolicyiter::Int = 10)
    # v = value function at previous time-step
    # a = policy function at previous time-step / initial guess for update
    taux = Δτ/Δx
    htaux2 = 0.5*Δτ/Δx^2

    t = (ti-1)*Δτ
    n = length(x)
    @assert length(a) == n

    coeff0 = ones(x)   # v_i
    coeff1 = zeros(n-1) # v_{i+1} # TODO: type stability
    coeff2 = zeros(n-1) # v_{i-1} # TODO: type stability
    rhs = zeros(x)
    # Dirichlet conditions
    if model.bcond[1] == true
        @inbounds rhs[1] = model.Dfun[1](t)
    end
    if model.bcond[2] == true
        @inbounds rhs[end] = model.Dfun[2](t)
    end

    vnew = v
    pol = copy(a) # TODO: just update a instead?

    for k in 0:maxpolicyiter
        updatepol!(pol, vnew, model, t, x, Δx)
        updatecoeffs!(coeff0, coeff1, coeff2, rhs, model, v, t, x, pol, Δτ, Δx)
        Mat = spdiagm((coeff2, coeff0, coeff1), -1:1, n, n)

        @inbounds begin
            # Move to sparse(I,J,K) instead?
            if model.bcond[1] == false
                bval = taux*model.b(t,x[1],pol[1])
                sval2 = htaux2*model.σ(t,x[1],pol[1])^2
                Mat[1,1:3] = [1.-(sval2-bval), -(bval-2*sval2), -sval2]
                rhs[1] = v[1] + Δτ*model.f(t,x[1],pol[1])
            end

            if model.bcond[2] == false
                bval = taux*model.b(t,x[end],pol[end])
                sval2 = htaux2*model.σ(t,x[end],pol[end])^2
                Mat[end,end-2:end] = [-sval2, bval+2*sval2, 1.-(sval2+bval)]
                rhs[end] = v[end] + Δτ*model.f(t,x[end],pol[end])
            end
        end

        # TODO: Use Krylov solver for high-dimensional PDEs
        vold = vnew
        vnew = Mat\rhs

        vchange = maximum(abs(vnew-vold)./max(scale,abs(vnew)))
        if vchange < tol && k > 0
            break
        end
    end
    updatepol!(pol, vnew, model, t, x, Δx)

    return vnew, pol
end

function timeloopiteration(model::HJBOneDim, N::Int,
                           Δτ, vinit, x, Δx)
    # Pass v and pol by reference?
    K = length(x)
    v = zeros(K, N+1)
    pol = zeros(K, N) # No policy at t = T

    @inbounds v[:,N+1] = vinit # We use forward time t instead of backward time τ
    polinit = fill(0.5*(model.amax+model.amin), K) # initial guess for control
    @inbounds v[:,N], pol[:,N] = policynewtonupdate(model, v[:,N+1],
                                                    polinit, x, Δx, Δτ, N)

    for j = N-1:-1:1
        @inbounds begin
            # t = (N-j)*Δτ
            # TODO: pass v-column, pol-column by reference?
            v[:,j], pol[:,j] = policynewtonupdate(model, v[:,j+1],
                                                  pol[:,j+1], x, Δx, Δτ, j)
        end
    end

    return v, pol
end

max_stars_repo_name: anriseth/HJBSolver

max_issues_repo_name: anriseth/HJBSolver

max_forks_repo_name: anriseth/HJBSolver


text:
"""
    new_linear_constraint(model::Optimizer, func::SAF{T}, set) where {T<:Real}

Create a new linear constraint and return a `LinearConstraint` with already a correct index
such that it can be simply added with [`add_constraint!`](@ref)
"""
function new_linear_constraint(model::Optimizer, func::SAF{T}, set) where {T<:Real}
    lc = new_linear_constraint(func, set)
    lc.idx = length(model.inner.constraints) + 1
    return lc
end

function new_linear_constraint(func::SAF{T}, set) where {T<:Real}
    func = remove_zero_coeff(func)

    indices = [v.variable_index.value for v in func.terms]

    return LinearConstraint(0, func, set, indices)
end

function remove_zero_coeff(func::MOI.ScalarAffineFunction)
    terms = [term for term in func.terms if term.coefficient != 0]
    return MOI.ScalarAffineFunction(terms, func.constant)
end

"""
    get_indices(func::VAF{T}) where {T}

Get indices from the VectorAffineFunction
"""
function get_indices(func::VAF{T}) where {T}
    return [v.scalar_term.variable_index.value for v in func.terms]
end

"""
    get_inner_constraint

Create the inner constraint of a reified or indicator constraint
"""
function get_inner_constraint(func::VAF{T}, set::Union{ReifiedSet, IndicatorSet}, inner_set::MOI.AbstractVectorSet) where {T<:Real}
    f = MOIU.eachscalar(func)
    inner_internals = ConstraintInternals(
        0,
        f[2:end],
        set.set,
        get_indices(f[2:end]),
    )
    return init_constraint_struct(set.set, inner_internals)
end

function get_inner_constraint(func::VAF{T}, set::Union{ReifiedSet, IndicatorSet, MOI.IndicatorSet}, inner_set::MOI.AbstractScalarSet) where {T<:Real}
    inner_terms = [v.scalar_term for v in func.terms if v.output_index == 2]
    inner_constant = func.constants[2]
    inner_set = set.set

    inner_func = MOI.ScalarAffineFunction{T}(inner_terms, inner_constant)

    return new_linear_constraint(inner_func, inner_set)
end

function get_inner_constraint(vars::MOI.VectorOfVariables, set::Union{ReifiedSet, IndicatorSet}, inner_set)
    inner_internals = ConstraintInternals(
        0,
        MOI.VectorOfVariables(vars.variables[2:end]),
        set.set,
        Int[v.value for v in vars.variables[2:end]],
    )
    return init_constraint_struct(set.set, inner_internals)
end

function get_activator_internals(A, indices)
    ActivatorConstraintInternals(A, indices[1] in indices[2:end], false, 0)
end


text:
using BinaryBuilder, Pkg

name = "RDKit"
version = v"2022.03.1"

sources = [
    GitSource("https://github.com/rdkit/rdkit.git", "7e205e0d93a3046c1eaab37120c9f6971194ddf2"),
    DirectorySource("./bundled"),
]

script = raw"""
cd ${WORKSPACE}/srcdir/rdkit

# Fix name of static libraries dependencies of `librdkitcffi` when building for Windows.
atomic_patch -p1 ../patches/static-libraries-windows.patch
# To check whether to optimise popcnt you must check the _*TARGET*_ system, not
# the host one.
atomic_patch -p1 ../patches/popcnt-target-system.patch
# Windows build fails to link a test, despite the fact we don't want tests.
atomic_patch -p1 ../patches/do-not-build-cffi-test.patch

FLAGS=()
if [[ "${target}" == *-mingw* ]]; then
    FLAGS+=(-DRDK_BUILD_THREADSAFE_SSS=OFF)
fi

mkdir build
cd build
cmake \
    -DCMAKE_BUILD_TYPE=Release \
    -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
    -DCMAKE_INSTALL_PREFIX=${prefix} \
    -DRDK_INSTALL_INTREE=OFF \
    -DRDK_BUILD_INCHI_SUPPORT=ON \
    -DRDK_BUILD_PYTHON_WRAPPERS=OFF \
    -DRDK_BUILD_CFFI_LIB=ON \
    -DRDK_BUILD_FREETYPE_SUPPORT=ON \
    -DRDK_BUILD_CPP_TESTS=OFF \
    -DRDK_BUILD_SLN_SUPPORT=OFF \
    -DRDK_TEST_MULTITHREADED=OFF \
    "${FLAGS[@]}" \
    ..
make -j${nproc}
make install
"""

platforms = [
    Platform("x86_64", "linux"; libc="glibc"),
    Platform("i686", "linux"; libc="glibc"),
    Platform("aarch64", "linux"; libc="glibc"),
    Platform("x86_64", "macos"),
    Platform("aarch64", "macos"),
    Platform("x86_64", "windows"),
]

platforms = expand_cxxstring_abis(platforms)

products = [
    LibraryProduct("librdkitcffi", :librdkitcffi),
]

dependencies = [
    Dependency("FreeType2_jll"),
    Dependency("boost_jll"; compat="=1.76.0"),
    BuildDependency("Eigen_jll"),
    Dependency("Zlib_jll"),
]

build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies;
               preferred_gcc_version=v"7", julia_compat="1.6")

max_stars_repo_name: MichelJuillard/Yggdrasil

max_issues_repo_name: MichelJuillard/Yggdrasil

max_forks_repo_name: MichelJuillard/Yggdrasil


text:
@testset "test for voltage deviation" begin
    @testset "case_vd: cigre mv" begin
        result = run_pandamodels_vd(case_vd)

        pm = _PdM.load_pm_from_json(case_vd)
        params = _PdM.extract_params!(pm)

        @test string(result["termination_status"]) == "LOCALLY_SOLVED"
        @test string(result["dual_status"]) == "FEASIBLE_POINT"
        @test string(result["primal_status"]) == "FEASIBLE_POINT"

        bus = result["solution"]["bus"]
        for (idx, bus) in result["solution"]["bus"]
            if idx in keys(params[:setpoint_v])
                @test isapprox(bus["vm"], params[:setpoint_v][idx]["value"], atol=1e-1)
            end
        end

        @test isapprox(result["objective_lb"], -Inf)
        @test isapprox(result["objective"], 0.000453688; atol = 1e-8)

        @test result["solve_time"] > 0.0
    end

end


text:
# This problem finds a least cost shipping schedule that meets
# requirements at markets and supplies at factories.
#
# Dantzig, G B, Chapter 3.3. In Linear Programming and Extensions.
# Princeton University Press, Princeton, New Jersey, 1963.

using GAMS
using JuMP
using Test

function example_transport(; verbose = true)
   model = Model(GAMS.Optimizer)
   set_optimizer_attribute(model, MOI.Silent(), !verbose)

   a = [350, 600]
   b = [325, 300, 275]
   d = [2.5 1.7 1.8; 2.5 1.8 1.4]

   @variable(model, x[1:2,1:3] >= 0)

   @objective(model, Min, 0.09 * sum(d[i,j] * x[i,j] for i = 1:2, j = 1:3))

   @constraint(model, [i = 1:2], sum(x[i,j] for j = 1:3) <= a[i])
   @constraint(model, [j = 1:3], sum(x[i,j] for i = 1:2) >= b[j])

   if verbose
      print(model)
   end

   JuMP.optimize!(model)

   obj_opt = JuMP.objective_value(model)
   x_opt = JuMP.value.(x)

   if verbose
      println("Objective value: ", obj_opt)
      println("x = ", x_opt)
   end

   @test obj_opt ≈ 153.6750
   @test x_opt[1,1] ≈ 50.0
   @test x_opt[1,2] ≈ 300.0
   @test x_opt[1,3] ≈ 0.0
   @test x_opt[2,1] ≈ 275.0
   @test x_opt[2,2] ≈ 0.0
   @test x_opt[2,3] ≈ 275.0
end

example_transport(verbose = false)


text:
module MLJScientificTypes

# Dependencies
using ScientificTypes
using Tables, CategoricalArrays, ColorTypes, PrettyTables, Dates

# re-exports from ScientificTypes
export Scientific, Found, Unknown, Known, Finite, Infinite,
    OrderedFactor, Multiclass, Count, Continuous, Textual,
    Binary, ColorImage, GrayImage, Image, Table,
    ScientificTimeType, ScientificDate, ScientificDateTime,
    ScientificTime
export scitype, scitype_union, elscitype, nonmissing, trait

# exports
export coerce, coerce!, autotype, schema, info

# -------------------------------------------------------------
# Abbreviations

const ST   = ScientificTypes
const Arr  = AbstractArray
const CArr = CategoricalArray
const Cat  = CategoricalValue

# Indicate the convention, see init.jl where it is set.
struct MLJ <: Convention end

include("init.jl")

# -------------------------------------------------------------
# Includes

include("coerce.jl")
include("schema.jl")
include("autotype.jl")

include("convention/scitype.jl")
include("convention/coerce.jl")

end # module


text:
import YaoBase: measure!, measure
using LinearAlgebra: eigen!

"""
    eigenbasis(op::AbstractBlock{N})

Return the `eigenvalue` and `eigenvectors` of target operator.
By applying `eigenvector`' to target state,
one can swith the basis to the eigenbasis of this operator.
However, `eigenvalues` does not have a specific form.
"""
function eigenbasis(op::AbstractBlock{N}) where {N}
    m = mat(op)
    if m isa Diagonal || m isa IMatrix
        op, IdentityGate{N}()
    else
        E, V = eigen!(Matrix(m))
        matblock(Diagonal(E)), matblock(V)
    end
end

for BT in []
    @eval function eigenbasis(op::$BT)
        @warn "eigenbasis on blocktype `$($BT)` calls into the fallback implementation, which might be slow. Try using `kron`, `repeat` if items commute to each oher."
        invoke(eigenbasis, Tuple{AbstractBlock}, op)
    end
end

"""
Return true if operators commute to each other.
"""
function simple_commute_eachother(ops::Vector{<:AbstractBlock{N}}) where {N}
    occ = zeros(Bool, N)
    for op in ops
        for i in occupied_locs(op)
            if occ[i]
                return false
            else
                occ[i] = true
            end
        end
    end
    return true
end

function eigenbasis(op::ChainBlock{N}) where {N}
    # detect commute operators
    if simple_commute_eachother(subblocks(op))
        E = chain(N)
        blks = chain(N)
        for b in subblocks(op)
            Ei, Vi = eigenbasis(b)
            push!(E, Ei)
            push!(blks, Vi)
        end
        return E, blks
    else
        if N > 5
            @warn "eigenbasis on blocktype `ChainBlock` (size $N) calls into the fallback implementation, which might be slow. Try using `kron`, `repeat` if items commute to each oher. If this behavior is not what you expected, please file an issue here: https://github.com/QuantumBFS/Yao.jl/issues."
        end
        invoke(eigenbasis, Tuple{AbstractBlock}, op)
    end
end

for GT in [:PutBlock, :RepeatedBlock, :ControlBlock, :Daggered]
    @eval function eigenbasis(op::$GT)
        E, V = eigenbasis(content(op))
        chcontent(op, E), chcontent(op, V)
    end
end

for GT in [:RotationGate, :TimeEvolution, :Scale]
    @eval function eigenbasis(op::$GT)
        E, V = eigenbasis(content(op))
        chcontent(op, E), V
    end
end

for GT in [:CachedBlock]
    @eval function eigenbasis(op::$GT)
        eigenbasis(content(op))
    end
end

function eigenbasis(op::KronBlock{N}) where {N}
    E = []
    blks = []
    for (k, b) in op
        Ei, Vi = eigenbasis(b)
        push!(E, k => Ei)
        push!(blks, k => Vi)
    end
    kron(N, E...), kron(N, blks...)
end

function eigenbasis(op::XGate)
    Z, H
end

function eigenbasis(op::YGate)
    Z, ConstGate.S * H
end

function measure!(
    postprocess::NoPostProcess,  # operator measuring does not allow removing bits or resetting vaules.
    op::AbstractBlock,
    reg::AbstractRegister,
    locs::AllLocs;
    kwargs...,
)
    _check_msize(op, reg, locs)
    E, V = eigenbasis(op)
    res = measure!(postprocess, BlockedBasis(diag(mat(E))), reg |> V', locs; kwargs...)
    apply!(reg, V)
    return res
end

# `BlockedBasis` is for measuring operators on its eigen basis, where a `block` is a subspace with the same eigenvalue.
# * `perm` is the permutation that permute the basis by the ascending order of eigenvalues,
# * `values` are eigenvalues of target observable for each block,
# * `block_ptr` is the pointers for blocks, e.g. to index block `i`, one can use `block_ptr[i]:block_ptr[i-1]`, or `subblock(blockbasis, i)` for short.
struct BlockedBasis{VT}
    perm::Vector{Int}
    values::VT
    block_ptr::Vector{Int}
end

subblock(bb::BlockedBasis, i::Int) = bb.block_ptr[i]:bb.block_ptr[i+1]-1
nblocks(bb::BlockedBasis) = length(bb.block_ptr) - 1

function BlockedBasis(values::AbstractVector{T}) where {T}
    if length(values) == 1
        return BlockedBasis([1], values, [1, 2])
    elseif length(values) == 0
        return BlockedBasis([], values, [1])
    end
    order = sortperm(values; by = real)
    values = values[order]
    vpre = values[1]
    block_ptr = [1]
    unique_values = [vpre]
    k = 1
    for i in 2:length(values)
        v = values[i]
        if !isapprox(v, vpre)  # use approx in order to ignore the round off error
            k += 1
            push!(block_ptr, i)
            push!(unique_values, v)
        end
        vpre = v
    end
    push!(block_ptr, length(values) + 1)
    return BlockedBasis(order, unique_values, block_ptr)
end

function YaoBase.measure!(
    ::NoPostProcess,
    bb::BlockedBasis,
    reg::ArrayReg{B,T},
    ::AllLocs;
    rng::AbstractRNG = Random.GLOBAL_RNG,
) where {B,T}
    state = (reg|>rank3)[bb.perm, :, :]  # permute to make eigen values sorted
    pl = dropdims(sum(abs2, state, dims = 2), dims = 2)
    pl_block = zeros(eltype(pl), nblocks(bb), B)
    for ib in 1:B
        for i in 1:nblocks(bb)
            for k in subblock(bb, i)
                pl_block[i, ib] += pl[k, ib]
            end
        end
    end
    res = Vector{Int}(undef, B)
    @inbounds @views for ib in 1:B
        ires = sample(rng, 1:nblocks(bb), Weights(pl_block[:, ib]))
        # notice ires is `BitStr` type, can be use as indices directly.
        range = subblock(bb, ires)
        state[range, :, ib] ./= sqrt(pl_block[ires, ib])
        state[1:range.start-1, :, ib] .= zero(T)
        state[range.stop+1:size(state, 1), :, ib] .= zero(T)
        res[ib] = ires
    end
    # undo permute and assign back
    _state = reshape(state, 1 << nactive(reg), :)
    rstate = reshape(reg.state, 1 << nactive(reg), :)
    for j in 1:size(rstate, 2)
        for i in 1:size(rstate, 1)
            rstate[bb.perm[i], j] = _state[i, j]
        end
    end
    return B == 1 ? bb.values[res[]] : bb.values[res]
end

function measure(op::AbstractBlock, reg::AbstractRegister, locs::AllLocs; kwargs...)
    _check_msize(op, reg, locs)
    E, V = eigenbasis(op)
    res = measure(ComputationalBasis(), copy(reg) |> V', locs; kwargs...)
    diag(mat(E))[Int64.(res).+1]
end

render_mlocs(alllocs::AllLocs, locs) = locs
render_mlocs(alllocs, locs) = alllocs[locs]

function _check_msize(op, reg, locs)
    if (locs isa AllLocs ? nactive(reg) : length(locs)) != nqubits(op)
        throw(
            QubitMismatchError(
                "operator of size $(nqubits(op)) does not match register size $(nactive(reg))",
            ),
        )
    end
end

function measure(op::Scale, reg::AbstractRegister, locs::AllLocs; kwargs...)
    factor(op) .* measure(content(op), reg, locs; kwargs...)
end

function measure(op::CachedBlock, reg::AbstractRegister, locs::AllLocs; kwargs...)
    measure(content(op), reg, locs; kwargs...)
end

function measure(op::Daggered, reg::AbstractRegister, locs::AllLocs; kwargs...)
    conj(measure(content(op), reg, locs; kwargs...))
end

function measure(ab::Add, reg::AbstractRegister, locs::AllLocs; kwargs...)
    sum(subblocks(ab)) do op
        measure(op, reg, locs; kwargs...)
    end
end

function measure(op::PutBlock{N}, reg::AbstractRegister, locs; kwargs...) where {N}
    _check_msize(op, reg, locs)

    # get eigen basis
    E, V = eigenbasis(op)
    ai = AddressInfo(nactive(reg), locs)
    _E = map_address(E, ai)
    _V = map_address(V, ai)
    _reg = copy(reg) |> _V'

    # perform equivalent measure
    E = diag(mat(content(_E)))
    res = measure(ComputationalBasis(), _reg, _E.locs; kwargs...)
    map(ri -> E[Int64(ri)+1], res)
end

function measure(op::PutBlock{N}, reg::AbstractRegister, locs::AllLocs; kwargs...) where {N}
    invoke(measure, Tuple{PutBlock,AbstractRegister,Any}, op, reg, locs; kwargs...)
end


text:
export ups, downs

ups(::Type{ST}, dims::Int...) where ST = ups(ST, dims)
downs(::Type{ST}, dims::Int...) where ST = downs(ST, dims)

ups(::Type{ST}, dims::Dims) where ST = fill(up(ST), dims)
downs(::Type{ST}, dims::Dims) where ST = fill(down(ST), dims)

Random.rand(::Type{Bit{T}}, dims::Dims) where T = Bit{T}.(rand(Bool, dims))
Random.rand(::Type{Spin{T}}, dims::Dims) where T = Spin{T}.(2 * rand(Bool, dims) .- 1)
Random.rand(::Type{Half{T}}, dims::Dims) where T = Half{T}.(rand(Bool, dims) .- 0.5)

# static array
import Base: @_inline_meta

ups(::Type{SA}) where {SA <: StaticArray} = _ups(Size(SA), SA)

@generated function _ups(::Size{s}, ::Type{SA}) where {s, SA <: StaticArray}
    T = eltype(SA)
    if T == Any
        T = Bit{Float64} # default
    end
    v = [:(up($T)) for i = 1:prod(s)]
    return quote
        @_inline_meta
        $SA(tuple($(v...)))
    end
end

downs(::Type{SA}) where {SA <: StaticArray} = _downs(Size(SA), SA)

@generated function _downs(::Size{s}, ::Type{SA}) where {s, SA <: StaticArray}
    T = eltype(SA)
    if T == Any
        T = Bit{Float64} # default
    end
    v = [:(down($T)) for i = 1:prod(s)]
    return quote
        @_inline_meta
        $SA(tuple($(v...)))
    end
end


text:
# Autogenerated wrapper script for CUDA_jll for powerpc64le-linux-gnu
export compute_sanitizer, libcublas, libcublasLt, libcudadevrt, libcudart, libcufft, libcufftw, libcupti, libcurand, libcusolver, libcusolverMg, libcusparse, libdevice, libnppc, libnppial, libnppicc, libnppidei, libnppif, libnppig, libnppim, libnppist, libnppisu, libnppitc, libnpps, libnvblas, libnvtoolsext, libnvvm, nvdisasm

JLLWrappers.@generate_wrapper_header("CUDA")
JLLWrappers.@declare_executable_product(compute_sanitizer)
JLLWrappers.@declare_library_product(libcublas, "libcublas.so.11")
JLLWrappers.@declare_library_product(libcublasLt, "libcublasLt.so.11")
JLLWrappers.@declare_file_product(libcudadevrt)
JLLWrappers.@declare_library_product(libcudart, "libcudart.so.11.0")
JLLWrappers.@declare_library_product(libcufft, "libcufft.so.10")
JLLWrappers.@declare_library_product(libcufftw, "libcufftw.so.10")
JLLWrappers.@declare_library_product(libcupti, "libcupti.so.11.5")
JLLWrappers.@declare_library_product(libcurand, "libcurand.so.10")
JLLWrappers.@declare_library_product(libcusolver, "libcusolver.so.11")
JLLWrappers.@declare_library_product(libcusolverMg, "libcusolverMg.so.11")
JLLWrappers.@declare_library_product(libcusparse, "libcusparse.so.11")
JLLWrappers.@declare_file_product(libdevice)
JLLWrappers.@declare_library_product(libnppc, "libnppc.so.11")
JLLWrappers.@declare_library_product(libnppial, "libnppial.so.11")
JLLWrappers.@declare_library_product(libnppicc, "libnppicc.so.11")
JLLWrappers.@declare_library_product(libnppidei, "libnppidei.so.11")
JLLWrappers.@declare_library_product(libnppif, "libnppif.so.11")
JLLWrappers.@declare_library_product(libnppig, "libnppig.so.11")
JLLWrappers.@declare_library_product(libnppim, "libnppim.so.11")
JLLWrappers.@declare_library_product(libnppist, "libnppist.so.11")
JLLWrappers.@declare_library_product(libnppisu, "libnppisu.so.11")
JLLWrappers.@declare_library_product(libnppitc, "libnppitc.so.11")
JLLWrappers.@declare_library_product(libnpps, "libnpps.so.11")
JLLWrappers.@declare_library_product(libnvblas, "libnvblas.so.11")
JLLWrappers.@declare_library_product(libnvtoolsext, "libnvToolsExt.so.1")
JLLWrappers.@declare_library_product(libnvvm, "libnvvm.so.4")
JLLWrappers.@declare_executable_product(nvdisasm)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_executable_product(
        compute_sanitizer,
        "bin/compute-sanitizer",
    )

    JLLWrappers.@init_library_product(
        libcublas,
        "lib/libcublas.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libcublasLt,
        "lib/libcublasLt.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_file_product(
        libcudadevrt,
        "lib/libcudadevrt.a",
    )

    JLLWrappers.@init_library_product(
        libcudart,
        "lib/libcudart.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libcufft,
        "lib/libcufft.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libcufftw,
        "lib/libcufftw.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libcupti,
        "lib/libcupti.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libcurand,
        "lib/libcurand.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libcusolver,
        "lib/libcusolver.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libcusolverMg,
        "lib/libcusolverMg.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libcusparse,
        "lib/libcusparse.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_file_product(
        libdevice,
        "share/libdevice/libdevice.10.bc",
    )

    JLLWrappers.@init_library_product(
        libnppc,
        "lib/libnppc.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libnppial,
        "lib/libnppial.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libnppicc,
        "lib/libnppicc.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libnppidei,
        "lib/libnppidei.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libnppif,
        "lib/libnppif.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libnppig,
        "lib/libnppig.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libnppim,
        "lib/libnppim.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libnppist,
        "lib/libnppist.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libnppisu,
        "lib/libnppisu.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libnppitc,
        "lib/libnppitc.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libnpps,
        "lib/libnpps.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libnvblas,
        "lib/libnvblas.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libnvtoolsext,
        "lib/libnvToolsExt.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libnvvm,
        "lib/libnvvm.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        nvdisasm,
        "bin/nvdisasm",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using Images, StaticArrays, LinearAlgebra, Statistics
using FileIO, ProgressMeter, Printf, ImageSegmentation
#using Plots, ImageView

include("psf_fun.jl")
include("s_save_image.jl")
include("filter.jl")
include("psf_analyze.jl")
include("beadstest.jl")

THRESHOLD = 7

print("Disk label: ")
label = chomp(readline())

print("Date: ")
date = chomp(readline())
#date="20190315"
@show date

cd(@sprintf("/mnt/%s/jchang/%s/", label, date))
files = readdir()
filelist = filter(x -> occursin(".imagine", x), files)
@show filelist

print("add scale bar(Y as 1, N as 0): ")
r = chomp(readline())
r = parse(Int, r)

l = 0
if r == 1
    print("The length of scale bar(in px): ")
    l = chomp(readline())
    l = parse(Int, l)
end

print("Background Filename: ")
BG_filename = chomp(readline())

if BG_filename == ""
    BG_filename = filelist[1][1:end-8]
    BG_filename *= "_1"
end

print("Save every file during process (1 as Yes, 0 as No): ")
Save = chomp(readline())
Save = (Save == "1")

for k = 1:size(filelist, 1)
    filename = filelist[k][1:end-8]
    exp = load(@sprintf("%s.imagine", filename))

    transfertotif(exp, filename, r, l) # save .imagine into .tif
    files = readdir()
    framelist = filter(x -> occursin(".tif", x), files)
    framelist = filter(x -> occursin(@sprintf("%s", filename), x), framelist)

    for i = 1:size(framelist, 1)
        filename = framelist[i][1:end-4]
        exp = load(@sprintf("%s.tif", filename))

        #img = BG_subtraction(exp,BG_filename,filename,Save)# filter image with background signal
        img = convert(Array{N0f16}, exp)
        img = beads_segment(img, THRESHOLD, 10000, 3)
        #psf_analyze(img,filename)
        img3 = Gray.(convert(Array{N0f16}, OffsetArrays.no_offset_view(img)))
        img4 = Gray.(convert.(Normed{UInt16,16}, img3))
        img_save(img, "/home/jchang/image/result/", @sprintf("%s-fi.tif", filename))
    end
end


text:
using Test
using POMDPModels

let
    rng = MersenneTwister(41)
    N = 5
    p = LightDarkND{N, typeof(default_sigma)}()
    @test discount(p) == 0.95
    s0 = LightDarkNDState{N}(0,fill(0.0, N))
    s0, _, r = @gen(:sp, :o, :r)(p, s0, +1, rng)
    @test s0.loc[1] == 1.0
    @test r == 0
    s1, _, r = @gen(:sp, :o, :r)(p, s0, 0, rng)
    @test s1.status < 0
    @test r == -10.0
    s2 = LightDarkNDState{N}(0, fill(5.0, N))
    obs = rand(rng, observation(p, nothing, nothing, s2))
    @test sum(abs.(obs.-5.0)) <= 0.1

    sv = convert_s(Array{Float64}, s2, p)
    @time for i = 1:1000000 convert_s(Array{Float64}, s2, p) end
    @test sv == [0.0, fill(5.0, N)...]
    s = convert_s(LightDarkNDState{N}, sv, p)
    @time for i = 1:1000000 convert_s(LightDarkNDState{N}, sv, p) end
    @test s == s2
end


text:
struct ModelClosure{M,A} <: AbstractConditionalModel{M,A,NamedTuple{(), Tuple{}}}
    model::M
    argvals::A
end

function Base.show(io::IO, mc::ModelClosure)
    println(io, "ModelClosure given")
    println(io, "    arguments    ", keys(argvals(mc)))
    println(io, "    observations ", keys(observations(mc)))
    println(io, model(mc))
end

export argvals
argvals(c::ModelClosure) = c.argvals

export observations
observations(c::ModelClosure) = NamedTuple()

export observed
function observed(mc::ModelClosure{M,A}) where {M,A}
    NamedTuple()
end

model(c::ModelClosure) = c.model

ModelClosure(m::AbstractModel) = ModelClosure(m,NamedTuple())

(m::AbstractModel)(nt::NamedTuple) = ModelClosure(m, nt)

(mc::ModelClosure)(nt::NamedTuple) = ModelClosure(model(mc), merge(mc.argvals, nt))

argvalstype(mc::ModelClosure{M,A}) where {M,A} = A
argvalstype(::Type{MC}) where {M,A,MC<:ModelClosure{M,A}} = A

obstype(::ModelClosure) = NamedTuple{(), Tuple{}}
obstype(::Type{<:ModelClosure}) = NamedTuple{(), Tuple{}}

type2model(::Type{MC}) where {M,MC<:ModelClosure{M}} = type2model(M)

import Base

Base.:|(m::ModelClosure, nt::NamedTuple) = ModelPosterior(m, nt)


text: using JLD


text:
using Test
using MARY_fRG.Basics

@testset "长方形" begin
    rect1 = Rectangle(Point2D(0., 0.), 1, 2)
    println(rect1)
    @test area(rect1) == 2
end

@testset "正方形" begin
    squ1 = Square(Point2D(0., 0.), 2)
    println(squ1)
    @test area(squ1) == 4
end