text:
using Test
using SequentialRaytrace
using LsqFit

Air2(t) = ConstantIndex(convert(t, 1.0))

function mylens(imageDist)
    Lens(
        "",
        Object(Air, 10.0),
        [OpticalComponent(
            :base,
            Nothing,
            nothing,
            [
                EvenAsphere(0.0, 0.0, 0.0, 0.0, 0.0, ClearDiameter(100.0), Silica, 1.5, :S1),
                Plano(nothing, Silica, 1.5),
                Sphere(-1.0 / 98.3, nothing, Air, imageDist, :S2),
            ],
        )],
    )
end

function model(x, p)
    imageDist = p[1]
    t = typeof(imageDist)
    testray = ray_from_NA(0.0, 0.125)
    lens = mylens(imageDist)
    result = gen_result(Vector{Ray{typeof(imageDist)}}, lens)
    r = trace!(result, lens, testray, 1.0)
    r[end].y
end

@testset "LsqFit" begin
    @testset "optimize without automatic derivatives" begin
        xdata = [0.0]
        ydata = [0.0]

        vect_model(x, p) = model.(x, p)
        r = curve_fit(vect_model, xdata, ydata, [10.0])
        # r = curve_fit(vect_model, xdata, ydata, [10.0], autodiff = :forwarddiff)
        @test r.param[1] ≈ -12.7875519 atol = 1e-7
    end
end


text:
using Makie
using HDF5
using PETScBinaryIO
using BALUtils
using ArgParse
using LinearAlgebra
using SparseArrays
using ColorSchemes
using Colors
using LightGraphs
using NetworkLayout
using SimpleWeightedGraphs
using RCall
using DataFrames
using Statistics

function poses2centers(poses)
    cameras = mapslices(x -> BALUtils.Camera(x), poses, dims=1)
    hcat(Vector.(BALUtils.center.(cameras))...)
end

function draw_diff!(scene :: Scene, initial :: Array, change :: Array, scale :: Float64; kwargs...)
    moved = initial .+ change * scale
    paired = reshape(vcat(initial, moved), size(initial,1), :)
    # linesegments!(scene, paired[1,:], paired[2,:], paired[3,:], transparency=true, show_axis=false; kwargs...)
    c = change .* scale
    arrows!(scene, initial[1,:], initial[2,:], initial[3,:], c[1,:], c[2,:], c[3,:], transparency=true, show_axis=false; kwargs...)
end

function setup_scene!(s :: Scene, points)
    plot_center = mapslices(median, points, dims=2)[:]
    old_center = mapslices(mean, points, dims=2)[:]
    plot_center = old_center
    function setcam()
        cam3d_cad!(s)
        update_cam!(s)
        off = s.data_limits[].widths[2]
        cam = cameracontrols(s)
        cam.near[] = 0.0001
        update_cam!(s, plot_center .+ [0., 170, 0.], plot_center)
    end
    lift(a -> setcam(), s.events.window_area)
    lift(a -> setcam(), s.events.window_open)
end

function gplot(A :: SparseMatrixCSC; vertexcolors=nothing, vertexcolors2=nothing)
    net = NetworkLayout.Circular.layout(A)
    s = Scene()
    is, js, ks = findnz(A)
    ls = map((i,j) -> net[i] => net[j], is, js)
    (min, max) = extrema(ks)
    linecolors = RGBA.(0.0, 0.0, 0.0, (ks .- min) ./ (max - min) .* 0.9 .+ 0.1)
    linesegments!(s, ls, color=linecolors, transparecy=true, linewidth=3)
    scatter!(s, net, color=vertexcolors, transparency=true, colormap=:pu_or, markersize=0.2)
    #scatter!(s, net, transparency=true, colormap=:blues, markersize=0.1)
    s
end

function gplot3d(A :: SparseMatrixCSC, poses :: Matrix; vertexcolors=nothing, vertexcolors2=nothing)
    centers = mapslices(x -> BALUtils.center(BALUtils.Camera(x)), poses, dims=1)
    s = Scene()
    is, js, ks = findnz(A)
    ls = mapreduce((i,j) -> hcat(centers[:,i], centers[:,j]), hcat, is, js)
    (min, max) = extrema(ks)
    linecolors = RGBA.(0.0, 0.0, 0.0, (ks .- min) ./ (max - min) .* 0.9 .+ 0.1)
    #linesegments!(s, ls[1,:], ls[2,:], ls[3,:], transparency=true, color=linecolors, linewidth=3)
    scatter!(s, centers[1,:], centers[2,:], centers[3,:], color=vertexcolors, transparency=true, colormap=:pu_or, markersize=0.02)
    #scatter!(s, net, transparency=true, colormap=:blues, markersize=0.1)
    s
end

function plot_eigen_distribution(d, index)
    bal_file = joinpath(d, "problem_noised.bal")
    bbal_file = joinpath(d, "problem_noised.bbal")
    ba = if isfile(bal_file)
        readbal(bal_file)
    else
        readbal(bbal_file)
    end
    dscale = map(x -> x[1:num_cameras(ba)*9], readpetsc(joinpath(d, "dscale.petsc")))

    r = 0:300

    x = map(r) do i
        eig = h5read(joinpath(d, "eig_$(index)_smallest.h5"), "/eigenvector$i") .* reshape(dscale[index], 9, :)
        sort(mapslices(norm, eig[1:3,:], dims=1)[:]),
        sort(mapslices(norm, eig[4:6,:], dims=1)[:]),
        sort(mapslices(norm, eig[7:9,:], dims=1)[:])
    end
    rots = getindex.(x, 1)
    trans = getindex.(x, 2)
    intrin = getindex.(x, 3)

   #  df2 = vcat(
   #             DataFrame(norm=map(norm, rots), count=map(x -> count(x .> 1e-4)/length(x), rots),
   #                       eig=1:length(rots), group="rotation"),
   #             DataFrame(norm=map(norm, trans), count=map(x -> count(x .> 1e-2)/length(x), trans),
   #                       eig=1:length(trans), group="translation"),
   #             DataFrame(norm=map(norm, intrin), count=map(x -> count(x .> 1e-3)/length(x), intrin),
   #                       eig=1:length(intrin), group="intrinsic")
   # )
   #  ggplot(df2, aes(x=:eig, y=:count, color=:group)) + geom_point()
    df = vcat(
         mapreduce(x -> DataFrame(magnitude=x[2], x=1:length(x[2]), group="rotation", eig=x[1]), vcat, enumerate(rots)),
         mapreduce(x -> DataFrame(magnitude=x[2], x=1:length(x[2]), group="translation", eig=x[1]), vcat, enumerate(trans)),
         mapreduce(x -> DataFrame(magnitude=x[2], x=1:length(x[2]), group="intrinsic", eig=x[1]), vcat, enumerate(intrin)))
    ggplot(df, aes(x=:x, y=:magnitude, color=:eig)) + geom_point(size=0.1) + scale_y_log10()
end

function visualize(d, opts, index, eigindex, scale)
    solution = readpetsc(joinpath(d, "solution.petsc"))
    bal_file = joinpath(d, "problem_noised.bal")
    bbal_file = joinpath(d, "problem_noised.bbal")
    ba = if isfile(bal_file)
        readbal(bal_file)
    else
        readbal(bbal_file)
    end
    dscale = map(x -> x[1:num_cameras(ba)*9], readpetsc(joinpath(d, "dscale.petsc")))
    poses = map(x -> reshape(x[1:num_cameras(ba)*9], 9, :), solution)[index]
    centers = mapslices(x -> BALUtils.center(BALUtils.Camera(x)), poses, dims=1)

    S = readpetsc(joinpath(d, "bamg_$(index)_$(opts)_S.petsc"))[1]
    # network = NetworkLayout.SFDP.layout(S, Point3f0)
    # centers[1:3,:] = hcat(map(x -> Vector(x), network)...)

    eig_ = reshape(h5read(joinpath(d, "cr_eigs_$(index)_$(opts).h5"), "/eigenvector$(eigindex-1)"), 9, :)
    eig = eig_ .* reshape(dscale[index][1:num_cameras(ba)*9], 9, :)
    # To visualize modes of main matrix
    # eig_ = h5read(joinpath(d, "eig_$(index)_smallest.h5"), "/eigenvector$(eigindex-1)")
    # eig = eig_ .* reshape(dscale[index], 9, :)

    # eig = reshape(h5read(joinpath(d, "cr_eigs_$(index)_$(opts).h5"), "/eigenvector$eigindex") .* dscale[index], 9, :)
    displaced = poses .+ eig
    displaced_centers = mapslices(x -> BALUtils.center(BALUtils.Camera(x)), displaced, dims=1)
    eignorms = mapslices(x -> norm(x, Inf), eig, dims=1)[:]

    # get aggregate ids for each vertex
    P = readpetsc(joinpath(d, "bamg_$(index)_$(opts)_P.petsc"))[1]
    aggs = fill(-1, num_cameras(ba))
    is, js, _ = findnz(P)
    for (i, j) in zip(is, js)
        if (i-1) % 9 == 0 && (j - 1) % 9 == 0
            aggs[(i-1)÷9+1] = (j-1)÷9+1
        end
    end
    aggs_ = copy(aggs)
    aggs_[findall(aggs_ .== -1)] = maximum(aggs_)+1:maximum(aggs_)+count(aggs_ .== -1)
    colors = distinguishable_colors(maximum(aggs_))[aggs_]

    # verify that the eigenvector is orthogonal to prolongation on the aggregate
    aggids = unique(aggs)
    for i in aggids
        if i == -1
            continue
        end
        inds = findall(aggs .== i)
        binds = mapreduce(x -> (x-1)*9+1:x*9, vcat, inds)
        n = norm(P[binds, :]' * eig_[:, inds][:], Inf)
        if length(inds) > 1 && n > 1e-10
            @warn "Aggregate is in coarse space (inner product: $n)"
        end
    end

    aggids = unique(aggs)
    diams = map(aggids) do agg
        if agg == -1
            0
        else
            inds = aggs .== agg
            diameter(SimpleWeightedGraph(1 ./ S[inds, inds]))
        end
    end
    daggs = Dict(zip(aggids, diams))

    println("Largest parts of eigenvector")
    r = sortperm(eignorms, rev=true)[1:10]
    display(eig[:, r])
    println("")
    # println("A eigenvector")
    # display(A_eig[:, r])
    # println("")
    println("Associated poses")
    display(poses[:, r])

    # i = eignorms .> 0.1
    # @show S[i, i]
    # @show findall(i)
    # @show aggs[i]
    # j = aggs .== aggs[i][1]
    # @show findall(j)
    # @show eignorms[j]
    # @show S[j, j]
    # @show sum(S[j,:], dims=2)
    # @show sum(S[j,j], dims=2)
    # @show sum(S[j,j], dims=2) ./ sum(S[j,:], dims=2)
    # @show sum(S[j,j], dims=2) ./ sum(S[j,j])
    # @show eignorms[j]
    # @show diameter(SimpleWeightedGraph(1 ./ S[j, j]))
    # net = NetworkLayout.Spring.layout(S[j,j], Point2f0, iterations=100)
    # s2 = Scene()
    # is, js, ks = findnz(S[j,j])
    # inds = vcat(is, js)[:]
    # (min, max) = extrema(ks)
    # linecolors = RGBA.(0.0, 0.0, 0.0, (ks .- min) ./ (max - min))
    # linesegments!(s2, net[inds], color = linecolors)
    # scatter!(s2, net, color=eignorms[j],transparency=true,colormap=:blues)
    # # display(s2)
    # is, js, ks = findnz(S)
    # S_ = sparse(is, js, ones(size(ks)), size(S)...)
    # k = (S_ * S_* j) .!= 0.0
    # # display(gplot(S[k,k], vertexcolors=colors[k], vertexcolors2=eignorms[k]))

    # colors = eignorms

    s = Scene()
    scatter!(s, centers[1, :], centers[2, :], centers[3, :], color=colors
            , transparency=true, markersize=1.1, colormap=:blues)
    draw_diff!(s, centers, displaced_centers - centers, scale, transparency=true, linewidth=2.0)

    # draw strength of connection
    is, js, ks = findnz(S)
    (min, max) = extrema(ks)
    linecolors = RGBA.(0.0, 0.0, 0.0, (ks .- min) ./ (max - min))
    inds = vcat(is, js)[:]
    # linesegments!(s, centers[1, inds], centers[2, inds], centers[3, inds], color=linecolors, transparency=true)

    s
end

if PROGRAM_FILE == @__FILE__
    s = ArgParseSettings(autofix_names = true)
    @add_arg_table s begin
        "input_dir"
        help = "Bundle adjustment problem input directory"
        required = true
        "eigs"
        required = true
        "index"
        required = true
        arg_type = Int
    end
    args = parse_args(ARGS, s)

    visualize(args["input_dir"], args["eigs"], args["index"])
end

ext: jl


text:
### accessing libc functions (XOPEN_SOURCE >= 700, POSIX_C_SOURCE >= 200809L glibc>=2.24)

if Sys.isunix() && get(Base.ENV, "NO_CLOCALE", "") != "1"

    function newlocale_c(mask::Cint, clocale::AbstractString, base::CLocaleType)
        ccall(:newlocale, CLocaleType, (Cint, Cstring, CLocaleType), mask, clocale, base)
    end

    function duplocale(ploc::CLocaleType)
        ccall(:duplocale, CLocaleType, (CLocaleType,), ploc)
    end

    function freelocale(ploc::CLocaleType)
        ccall(:freelocale, Nothing, (CLocaleType,), ploc)
    end

    function uselocale(newloc::CLocaleType)
        ccall(:uselocale, CLocaleType, (CLocaleType,), newloc)
    end

    function strcoll_c(s1::AbstractString, s2::AbstractString, ploc::CLocaleType)
        res = 0
        if ploc == CL0
            ploc = current_clocale()
            res = ccall(:strcoll_l, Cint, (Cstring, Cstring, CLocaleType), s1, s2, ploc)
            freelocale(ploc)
        else
            res = ccall(:strcoll_l, Cint, (Cstring, Cstring, CLocaleType), s1, s2, ploc)
        end
        Int(res)
    end

    function nl_langinfo_c(nlitem::Cint, ploc::CLocaleType)
        ccall(:nl_langinfo_l, Ptr{UInt8}, (Cint, CLocaleType), nlitem, ploc)
    end

else # for the non case provide dummy methods

    newlocale_c(mask::Cint, clocale::AbstractString, base::CLocaleType) = CL0
    duplocale(ploc::CLocaleType) = CL0
    freelocale(ploc::CLocaleType) = nothing
    uselocale(newloc::CLocaleType) = CL0
    function strcoll_c(s1::AbstractString, s2::AbstractString, ploc::CLocaleType)
        s1 == s2 ? 0 : s1 < s2 ? -1 : 1
    end
    nl_langinfo_c(nlitem::Cint, ploc::CLocaleType) = Ptr{UInt8}(0)

end
#################################


text:
@windows? (
    begin
        if haskey(ENV, "VS120COMNTOOLS")
            vs_cmd_prompt = string(ENV["VS120COMNTOOLS"], "..\\..\\VC\\vcvarsall.bat")
        elseif haskey(ENV, "VS110COMNTOOLS")
            vs_cmd_prompt = string(ENV["VS110COMNTOOLS"], "..\\..\\VC\\vcvarsall.bat")
        elseif haskey(ENV, "VS100COMNTOOLS")
            vs_cmd_prompt = string(ENV["VS100COMNTOOLS"], "..\\..\\VC\\vcvarsall.bat")
        else
            error("Cannot find proper Visual Studio installation. VS 2013, 2012, or 2010 is required.")
        end

        # check whether 32 or 64 bit archtecture
        # NOTE: Actually, nvcc in x86 visual studio command prompt doesn't make 32-bit binary
        #       It depends on whether CUDA toolkit is 32bit or 64bit
        if Int == Int64
            arch = "amd64"
        else
            arch = "x86"
        end

        # Run nmake -f Windows.mk under visual studio command prompt
        run(`cmd /C "$vs_cmd_prompt" $arch & nmake -f Windows.mk`)
        cd("../test") do
            run(`cmd /C "$vs_cmd_prompt" $arch & nmake -f Windows.mk`)
        end
    end
    : # for linux or mac
    begin
        run(`make`)
        cd("../test") do
            run(`make`)
        end
    end)


text:
# Contains various membership function types
# ------------------------------------------

"""	 Membership function type
"""
abstract type MF end

"""	 Triangular membership function type

	 TriangularMF(l_vertex, center, r_vertex)

	 Properties
	 ----------
	 `l_vertex`, `center` and `r_vertex` are the vertices of the triangle, in order

	 `eval` function returns membership value at a point
	 `mean_at` function returns mean value at line clipped by given firing strength
"""
mutable struct TriangularMF<:MF

	l_vertex::Real
	center::Real
	r_vertex::Real

	eval::Function
	mean_at::Function
	get_n::Function
	cut::Function

	function TriangularMF(l_vertex::Real, center::Real, r_vertex::Real)

		if l_vertex <= center <= r_vertex

			this = new()

			this.l_vertex = l_vertex
			this.center = center
			this.r_vertex = r_vertex

			function inv(α)
				x1 = (this.center - this.l_vertex) * α + this.l_vertex
				x2 = (this.center - this.r_vertex) * α + this.r_vertex
				return x1, x2
			end
			this.eval = function eval(x)
				maximum([minimum([((x - this.l_vertex) / (this.center - this.l_vertex)), ((this.r_vertex - x) / (this.r_vertex - this.center))]), 0])
			end

			this.mean_at = function mean_at(firing_strength)

				if firing_strength != 1
					p1 = (this.center - this.l_vertex) * firing_strength + this.l_vertex
					p2 = (this.center - this.r_vertex) * firing_strength + this.r_vertex
					(p1 + p2) / 2
				elseif firing_strength == 1
					return this.center
				end

			end

			this.n_opt = function get_n(tol)
				return 0
			end

			this.cut = function cut(α)
				x1, x2 = inv(α)
				return CutMF(x1, x2, α, this)
			end

			this

		else

			error("invalid vertices")

		end

	end

end

"""
 Gaussian membership function type

	GaussianMF(center, sigma)

	 Properties
	 ----------
	 `center` is the center of the distribution
	 `sigma` determines width of the distribution

	 `eval` function returns membership value at a point
	 `mean_at` function returns mean value at line clipped by given firing strength

"""
mutable struct GaussianMF<:MF

	center::Real
	sigma::Real

	eval::Function
	mean_at::Function
	get_n::Function
	cut::Function

	function GaussianMF(center::Real, sigma::Real)

		this = new()

		this.center = center
		this.sigma = sigma

		function inv(α)
			x1 = this.center - this.sigma*sqrt(-2*log(α))
			x2 = this.center + this.sigma*sqrt(-2*log(α))
			return x1, x2
		end

		this.eval = function eval(x)
			exp( - 0.5 * ((x - this.center) / this.sigma) ^ 2)
		end

		this.mean_at = function mean_at(firing_strength)
			this.center
		end

		this.get_n = function get_n(tol)
			l, u = inv(0.01)
			c = this.center
			s = this.sigma
			sdiff(x) = this.eval(x)*(c^2 - 2*c*x - s^2 + x^2)/s^4
			xs = collect(l:(u-l)/100:u)
			ys = sdiff.(xs)
			ym = maximum(abs.(ys))
			sqrt((u - l)^3*ym/(12*tol))
		end

		this.cut = function cut(α)
			x1, x2 = inv(α)
			return CutMF(x1, x2, α, this)
		end

		this

	end

end

"""
	 Generalised Bell membership function type

	BellMF(a, b, c)

	 Properties
	 ----------
	 `a`, `b` and `c` the usual bell parameters with `c` being the center

	 `eval` function returns membership value at a point
	 `mean_at` function returns mean value at line clipped by given firing strength

"""
mutable struct BellMF<:MF

	a::Real
	b::Real
	c::Real

	eval::Function
	mean_at::Function
	get_n::Function
	cut::Function

	function BellMF(a::Real, b::Real, c::Real)

		this = new()

		this.a = a
		this.b = b
		this.c = c

		function inv(α)
			x1 = this.c - this.a*((1-α)/α)^(1/(2*b))
			x2 = this.c + this.a*((1-α)/α)^(1/(2*b))
			return x1, x2
		end
		this.eval = function eval(x)
			1 / (1 + abs((x - this.c) / this.a) ^ (2 * this.b))
		end

		this.mean_at = function mean_at(firing_strength)
			this.c
		end

		this.get_n = function get_n(tol)
			l, u = inv(α)
			return (u-l)^3/(12*tol)
		end
		this.cut = function cut(α)
			x1, x2 = inv(α)
			return CutMF(x1, x2, α, this)
		end

		this

	end

end

"""
	 Trapezoidal membership function type

	TrapezoidalMF(l_bottom_vertex, l_top_vertex, r_top_vertex, r_bottom_vertex)

	 Properties
	 ----------
	 `l_bottom_vertex`, `l_top_vertex`, `r_top_vertex` and `r_bottom_vertex` are the vertices of the trapezoid, in order

	 `eval` function returns membership value at a point
	 `mean_at` function returns mean value at line clipped by given firing strength

"""
mutable struct TrapezoidalMF<:MF

	l_bottom_vertex::Real
	l_top_vertex::Real
	r_top_vertex::Real
	r_bottom_vertex::Real

	eval::Function
	mean_at::Function
	get_n::Function
	cut::Function

	function TrapezoidalMF(l_bottom_vertex::Real, l_top_vertex::Real, r_top_vertex::Real, r_bottom_vertex::Real)

		if l_bottom_vertex <= l_top_vertex <= r_top_vertex <= r_bottom_vertex

			this = new()

			this.l_bottom_vertex = l_bottom_vertex
			this.l_top_vertex = l_top_vertex
			this.r_top_vertex = r_top_vertex
			this.r_bottom_vertex = r_bottom_vertex

			function inv(α)
				x1 = (this.l_top_vertex - this.l_bottom_vertex) * α + this.l_bottom_vertex
				x2 = (this.r_top_vertex - this.r_bottom_vertex) * α + this.r_bottom_vertex
				return x1, x2
			end

			this.eval = function eval(x)
				maximum([minimum([((x - this.l_bottom_vertex) / (this.l_top_vertex - this.l_bottom_vertex)), 1, ((this.r_bottom_vertex - x) / (this.r_bottom_vertex - this.r_top_vertex))]), 0])
			end

			this.mean_at = function mean_at(firing_strength)
				p1 = (this.l_top_vertex - this.l_bottom_vertex) * firing_strength + this.l_bottom_vertex
				p2 = (this.r_top_vertex - this.r_bottom_vertex) * firing_strength + this.r_bottom_vertex
				(p1 + p2) / 2
			end

			this.get_n = function get_n(tol)
				return 0
			end

			this.cut = function cut(α)
				x1, x2 = inv(α)
				return CutMF(x1, x2, α, this)
			end

			this

		else

			error("invalid vertices")

		end

	end

end

"""
	 Sigmoid membership function type

	SigmoidMF(a, c, limit)

	 Properties
	 ----------
	 `a` controls slope
	 `c` is the crossover point
	 `limit` sets the extreme limit

	 `eval` function returns membership value at a point
	 `mean_at` function returns mean value at line clipped by given firing strength

"""
mutable struct SigmoidMF<:MF

	a::Real
	c::Real
	limit::Real

	eval::Function
	mean_at::Function
	get_n::Function
	cut::Function

	function SigmoidMF(a::Real, c::Real, limit::Real)

		if (a > 0 && limit > c) || (a < 0 && limit < c)

			this = new()

			this.a = a
			this.c = c
			this.limit = limit

			function inv(α)
				return this.c - (1/this.a)*log((1-α)/α)
			end

			this.eval = function eval(x)
				1 / (1 + exp(-this.a * (x - this.c)))
			end

			this.mean_at = function mean_at(firing_strength)

				if firing_strength == 1
					p_firing_strength = 0.999
				elseif firing_strength == 0
					p_firing_strength = 0.001
				else
					p_firing_strength = firing_strength
				end

				p1 = -log((1 / p_firing_strength) - 1) / this.a + this.c
				p2 = this.limit
				(p1 + p2) / 2

			end

			this.get_n = function get_n(tol)
				a = this.a
				c = this.c
				if a > 0
					u = this.limit
					l = inv(0.01)
				else
					l = this.limit
					u = inv(0.01)
				end
				aux1(x) = exp(-a*(x-c))
				aux2(x) = 2*a^2*(aux1(x))^2/(aux1(x) + 1)^3
				sdiff(x) = aux2(x) - a^2*aux1(x)/(aux1(x) + 1)^2
				xs = collect(l:(u-l)/100:u)
				ys = sdiff.(xs)
				ym = maximum(abs.(ys))
				sqrt((u - l)^3*ym/(12*tol))
			end

			this.cut = function cut(α)
				x = inv(α)
				if this.a > 0
					x1 = x
					x2 = this.limit
				else
					x1 = this.limit
					x2 = x
				end
				return CutMF(x1, x2, α, this)
			end

			this

		else

			error("invalid parameters")

		end

	end
end

"""
	Cut MF

	CutMF(a, b, c)

	 Properties
	 ----------
	 `a`, `b` and `c` the usual bell parameters with `c` being the center

	 `eval` function returns membership value at a point
	 `mean_at` function returns mean value at line clipped by given firing strength

"""
mutable struct CutMF<:MF

	x1::Real
	x2::Real
	α::Real
	toCutMF::MF

	eval::Function
	mean_at::Function
	get_n::Function
	cut::Function

	function CutMF(x1::Real, x2::Real, α::Real, toCutMF::MF)

		this = new()

		this.x1 = x1
		this.x2 = x2
		this.α = α
		this.toCutMF = toCutMF

		this.eval = function eval(x)
			if x1 <= x <= x2
				return this.α
			else
				return this.toCutMF.eval(x)
			end
		end

		this.mean_at = function mean_at(firing_strength)
			if firing_strength >= α
				return this.toCutMF.mean_at(α)
			else
				return this.toCutMF.mean_at(firing_strength)
			end
		end

		this.get_n = function get_n(tol)
			return this.toCutMF.get_n(tol)
		end

		this.cut = function cut(α)
			if α >= this.α
				return this
			else
				return this.toCutMF.cut(α)
			end
		end

		this

	end

end


text:
# test for symmetry on distances

using Manifolds
using Test

import Rotations as _Rot

##
@testset "test symmetry of Manifolds.distance" begin
##

M = TranslationGroup(2)
a,b = randn(2), randn(2)
@test isapprox( distance(M, a, b), distance(M, b, a), atol = 1e-5)

M = SpecialOrthogonal(2)
a,b = _Rot.RotMatrix(randn()), _Rot.RotMatrix(randn())
@test isapprox( distance(M, a, b), distance(M, b, a), atol = 1e-5)

M = SpecialEuclidean(2)
a = ProductRepr(randn(2),_Rot.RotMatrix(randn()))
b = ProductRepr(randn(2),_Rot.RotMatrix(randn()))
@test isapprox( distance(M, a, b), distance(M, b, a), atol = 1e-5)

M = SpecialOrthogonal(3)
a = _Rot.RotZ(randn())*_Rot.RotY(randn())*_Rot.RotX(randn())
b = _Rot.RotZ(randn())*_Rot.RotY(randn())*_Rot.RotX(randn())
@test isapprox( distance(M, a, b), distance(M, b, a), atol = 1e-5)

##
end


text:
const REGISTRY_PATH = joinpath(@__DIR__, "Registry.toml")

using Pkg
using UUIDs
using REPL.TerminalMenus

function yesno(query)
    menu = RadioMenu(["Yes", "No"])
    return request(query, menu) == 1
end

function get_pkgs_in_general()
    registries = Pkg.Types.collect_registries()
    general = registries[findfirst(reg -> reg.uuid == UUID("23338594-aafe-5451-b93e-139f81909106"), registries)]

    return Pkg.Types.read_registry(joinpath(general.path, "Registry.toml"))["packages"]
end

function valid_location(method, location)
    isempty(location) && return false

    if method == "vendored"
        if occursin(r"^(/|\\)", location)
            printstyled("\n`location` for `vendored` docs must be path relative to the package's root directory.\n", color=:red)
            return false
        end
    else
        if !occursin(r"^(http|https|git)\://", location)
            printstyled("\n`location` for `hosted` or `git-repo` docs must be a valid URL.\n", color=:red)
            return false
        end
    end
    return true
end

function get_config(pkg; novendor = false)
    methods = novendor ? ["hosted", "git-dir"] : ["vendored", "hosted", "git-dir"]
    menu = RadioMenu(methods)
    method = methods[request("Choose a `method` for $(pkg):", menu)]
    println()

    location = ""
    while !valid_location(method, location)
        println("Please choose a `location` for $(pkg):")
        location = strip(readline())
    end

    println()
    return Dict(
        "method" => method,
        "location" => location
    )
end

function add_by_uuid(uuid, pkg, config)
    if config == nothing
        config = get_config(pkg)
    end
    config["name"] = pkg

    toml = Pkg.TOML.parsefile(REGISTRY_PATH)
    if haskey(toml, uuid)
        println("UUID `$(uuid)` for package `$(pkg)` already in registry.")
        overwrite = yesno("Do you want to overwrite the existing entry?")
        println()
        if !overwrite
            return false
        end
    end
    toml[uuid] = config
    open(REGISTRY_PATH, "w") do io
        Pkg.TOML.print(io, toml)
    end
    return true
end

function add_to_registry(pkg, general_pkgs, config = config)
    pkgs_with_name = String[]
    for (uuid, _pkg) in general_pkgs
        _pkg["name"] == pkg && push!(pkgs_with_name, uuid)
    end

    if length(pkgs_with_name) == 0
        printstyled("No package with name `$(pkg)` found in registry.", color=:red)
        return false
    elseif length(pkgs_with_name) > 1
        menu = MultiSelectMenu(pkgs_with_name)
        indices = request("Choose one or multiple UUIDs for which you want to add a documentation hosting method:", menu)
        pkgs_with_name = pkgs_with_name[indices]
    end

    for uuid in pkgs_with_name
        success = add_by_uuid(uuid, pkg, config)
    end
end

function main(pkgnames)
    general_pkgs = get_pkgs_in_general()

    config = nothing
    if length(pkgnames) > 1
        one_for_all = yesno("Will all packages have the same configuration?")
        println()

        if one_for_all
            config = get_config("all provided packages", novendor = one_for_all)
        end
    end

    for pkg in pkgnames
        success = add_to_registry(pkg, general_pkgs, config)
    end
end

if length(ARGS) == 0 || any(isempty.(ARGS))
    println("Please provide one or multiple package names as arguments.")
else
    try
        main(ARGS)
    catch err
        println("\n")
        printstyled("An error occured. Please revert any changes and try again.", color=:red)
    end
    println("Please commit your changes and open a PR.")
end

max_forks_repo_name: pfitzseb/DocumentationGeneratorRegistry


text:
#=
max_percent_diff.jl
Copyright 2000-9-16, Jeff Fessler, University of Michigan
=#

export max_percent_diff

using Test: @test

"""
    d = max_percent_diff(s1, s2, [options])

compute the "maximum percent difference" between two signals: `s1, s2`
default is to normalize by `maximum(abs.(s1))`

options
- `maxboth::Bool` use max of both arguments to normalize; default `false`
- `normalize::Bool` normalize each before comparing; default `false`

"""
function max_percent_diff(s1, s2 ;
        maxboth::Bool = false, normalize::Bool = false)

    #check to make sure that we have comparable signals
    size(s1) != size(s2) &&
        throw(DomainError("Dimension mismatch $(size(s1)) $(size(s2))"))

    any(isnan.(s1)) && @warn("NaN values in s1 input to max_percent_diff")
    any(isnan.(s2)) && @warn("NaN values in s2 input to max_percent_diff")

    if normalize
        normFactor = sum(abs.(s2)) / sum(abs.(s1))
        print("Normalization factor (matrix 2 / matrix 1): $normFactor \n")
        s1 /= sum(abs.(s1))
        s2 /= sum(abs.(s2))
    end

    if maxboth
        denom = max(maximum(abs.(s1)), maximum(abs.(s2)))
        d = denom == 0 ? 0 : maximum(abs.(s1 - s2)) / denom
    else
        denom = maximum(abs.(s1))
        denom = denom == 0 ? denom : maximum(abs.(s2))
        d = denom == 0 ? 0 : maximum(abs.(s1 - s2)) / denom
    end

    return 100 * d
end

function max_percent_diff_test()
    @test max_percent_diff([200], [202]) ≈ 100/101
    @test max_percent_diff([0 200], [0 202]) ≈ 100/101
    @test max_percent_diff([0 200], [0 198], maxboth=true) ≈ 1
    #≈ does not work with 0 values (it's proportional to the second arg)
    @test max_percent_diff([0 100], [0 200], normalize=true) < .01
    true
end

"""
    max_percent_diff(:test)
self test
"""
function max_percent_diff(test::Symbol)
    test != :test && throw(DomainError("Invalid signal"));
    return max_percent_diff_test()
end

# max_percent_diff(:test)


text:
using AffineArithmetic, IntervalArithmetic, Polynomials, StaticArrays
using Test

using AffineArithmetic: Aff

@testset "Construction from intervals" begin
    X = 1..3
    X_a = Affine(X)
    X_a isa Affine

    Y = 2..4
    Y_a = Affine(Y)
    @test Y_a.range == 2..4
    @test Y_a.affine == Aff(3.0, SVector{1, Float64}(1.0), 0..0)

    sum_a = X_a + Y_a
    @test sum_a.range == 3..7
    @test sum_a.affine == Aff(5.0, SVector{1, Float64}(2.0), 0..0)

    prod_a = X_a * Y_a
    @test prod_a.range == 2..12
    @test prod_a.affine == Aff(6.0, SVector{1, Float64}(5.0), 0..1)
end

@testset "Small powers and range" begin
    X = 1..3
    Y = Affine(X)

    f(x) = x^2 - x + 1

    @test range(f(X)) == -1..9
    @test range(f(Y)) == 0..7   # affine is a little better

    g(x) = (x - 1)^3

    @test range(g(X)) == 0..8
    @test range(g(Y)) == 0..8   # affine gives the same
end

reset_affine_index()

@testset "Matrix multiplication" begin
    A = 0.5 * [1 2; -1 1]

    X = Affine(-1..1)
    Y = Affine(-1..1)

    XX = [X, Y]
    @test XX == [Affine(Aff(0.0, SVector{1, Float64}(1.0), 0..0), -1..1),
                 Affine(Aff(0.0, SVector{1, Float64}(1.0), 0..0), -1..1)]

    @test A * XX == [Affine(Aff(0.0, SVector{1, Float64}(1.5), 0..0), -1.5 .. 1.5),
                     Affine(Aff(0.0, SVector{1, Float64}(0.0), 0..0), 0..0)]
end

reset_affine_index()

@testset "Range of polynomial" begin
    # example from Rump 2015
    p = Polynomial([-3, 1])   # x - 3
    p2 = p^8

    x = 4 ± 1e-4
    y = Affine(x)

    @test (-70..70) ⊆ range(p2(x))
    @test range(p2(y)) ⊆ 0.998..1.002   # affine is extremely much better!
end

@testset "Constructors with intervals" begin
    X = Affine(1..2)
    @test X isa Affine
    @test X - X == Affine(0..0)

    x, y = affine(1..2, 3..4)
    @test interval(x + y - x - y) == 0..0
end

max_stars_repo_name: dpsanders/AffineArithmetic

max_issues_repo_name: dpsanders/AffineArithmetic

max_forks_repo_name: dpsanders/AffineArithmetic


text:
using Gadfly, DataArrays

spy(randn((10, 10)), Coord.cartesian(fixed=true))


text:
#####################################
# Helper functions for Dual numbers #
#####################################

@inline realpart(r::Real)             = r
@inline realpart(d::ForwardDiff.Dual) = d.value
@inline realpart(ds::Union{Vector,SubArray}) = Float64[realpart(d) for d in ds] # NOTE: the function below is assumed to return a Vector now
@inline realpart!(arr::Union{Array,SubArray}, ds::Union{Array,SubArray}) = for i = 1:length(ds) arr[i] = realpart(ds[i]) end
@inline realpart(ds::Matrix{T}) where {T <: Real} = Float64[realpart(col) for col in ds]
@inline realpart(ds::Matrix{Any}) = [realpart(col) for col in ds]
@inline realpart(ds::Array)  = map(d -> realpart(d), ds)  # NOTE: this function is not optimized
# @inline realpart(ds::TArray) = realpart(Array(ds))    # TODO: is it disabled temporary
@inline realpart(ta::Tracker.TrackedReal) = ta.data

@inline dualpart(d::ForwardDiff.Dual)       = d.partials.values
@inline dualpart(ds::Union{Array,SubArray}) = map(d -> dualpart(d), ds)

# Base.promote_rule(D1::Type{Real}, D2::Type{ForwardDiff.Dual}) = D2
import Base: <=
<=(a::Tracker.TrackedReal, b::Tracker.TrackedReal) = a.data <= b.data

#####################################################
# Helper functions for vectorize/reconstruct values #
#####################################################

vectorize(d::UnivariateDistribution, r::Real) = Vector{Real}([r])
vectorize(d::MultivariateDistribution, r::AbstractVector{<:Real}) = Vector{Real}(r)
vectorize(d::MatrixDistribution, r::AbstractMatrix{<:Real}) = Vector{Real}(vec(r))

# NOTE:
# We cannot use reconstruct{T} because val is always Vector{Real} then T will be Real.
# However here we would like the result to be specifric type, e.g. Array{Dual{4,Float64}, 2},
# otherwise we will have error for MatrixDistribution.
# Note this is not the case for MultivariateDistribution so I guess this might be lack of
# support for some types related to matrices (like PDMat).
@inline reconstruct(d::Distribution, val::Union{Vector,SubArray}) = reconstruct(d, val, typeof(val[1]))
# NOTE: the implementation below was: `T(val[1])`; it is changed to `val[1]` due to the need of ReverseDiff.jl, no side effect found yet
@inline reconstruct(d::UnivariateDistribution,   val::Union{Vector,SubArray}, T::Type) = val[1]
@inline reconstruct(d::MultivariateDistribution, val::Union{Vector,SubArray}, T::Type) = Array{T, 1}(val)
@inline reconstruct(d::MatrixDistribution,       val::Union{Vector,SubArray}, T::Type) = Array{T, 2}(reshape(val, size(d)...))

@inline reconstruct!(r, d::Distribution, val::Union{Vector,SubArray}) = reconstruct!(r, d, val, typeof(val[1]))
@inline reconstruct!(r, d::MultivariateDistribution, val::Union{Vector,SubArray}, T::Type) = (r[:] = val; r)

@inline reconstruct(d::Distribution, val::Union{Vector,SubArray}, n::Int) = reconstruct(d, val, typeof(val[1]), n)
@inline reconstruct(d::UnivariateDistribution,   val::Union{Vector,SubArray}, T::Type, n::Int) = Array{T, 1}(val)
@inline reconstruct(d::MultivariateDistribution, val::Union{Vector,SubArray}, T::Type, n::Int) = Array{T, 2}(reshape(val, size(d)[1], n))

@inline reconstruct(d::MatrixDistribution,       val::Union{Vector,SubArray}, T::Type, n::Int) = begin
  orig = Vector{Matrix{T}}(undef, n)
  tmp = Array{T, 3}(reshape(val, size(d)[1], size(d)[2], n))
  for i = 1:n
    orig[i] = tmp[:,:,i]
  end
  orig
end

@inline reconstruct!(r, d::Distribution, val::Union{Vector,SubArray}, n::Int) = reconstruct!(r, d, val, typeof(val[1]), n)
@inline reconstruct!(r, d::MultivariateDistribution, val::Union{Vector,SubArray}, T::Type, n::Int) = (r[:] = val; r)


text:
## --- Test simple metropolis_min and metropolis_minmax functions

    mu, sigma = collect(100:0.1:101), 0.01*ones(11);
    @test Chron.check_dist_ll(MeltsVolcanicZirconDistribution, mu, sigma, 100,101) ≈ -3.6933357793356576

    tmindist, tmaxdist, lldist, acceptancedist = metropolis_minmax(2*10^5, MeltsVolcanicZirconDistribution, mu, sigma, burnin=10^5)
    @test isapprox(nanmean(tmindist), 99.92, atol=0.015)
    @test isapprox(nanmean(tmaxdist), 101.08, atol=0.015)
    @test isapprox(nanmean(metropolis_min(2*10^5, MeltsVolcanicZirconDistribution, mu, sigma, burnin=10^5)), 99.9228, atol=0.015)

## -- Test distribution bootstrapping

     dist = BootstrapCrystDistributionKDE(1:10)
     @test isapprox(dist, [0.465404, 0.503149, 0.540308, 0.576485, 0.611314, 0.644472, 0.675691, 0.704756, 0.731516, 0.75588, 0.777815, 0.797345, 0.81454, 0.82951, 0.842398, 0.85337, 0.862605, 0.87029, 0.876611, 0.881749, 0.885873, 0.889136, 0.891676, 0.893612, 0.895043, 0.896047, 0.896684, 0.896993, 0.896993, 0.896684, 0.896047, 0.895043, 0.893612, 0.891676, 0.889136, 0.885873, 0.881749, 0.876611, 0.87029, 0.862605, 0.85337, 0.842398, 0.82951, 0.81454, 0.797345, 0.777815, 0.75588, 0.731516, 0.704756, 0.675691, 0.644472, 0.611314, 0.576485, 0.540308, 0.503149, 0.465404, 0.427489, 0.389826, 0.352829, 0.316888, 0.282361, 0.249559, 0.218742, 0.19011, 0.163804, 0.139901, 0.118424, 0.099338, 0.082567, 0.067992, 0.055465, 0.044819, 0.03587, 0.028432, 0.022318, 0.017347, 0.01335, 0.010173, 0.007675, 0.005732, 0.004238, 0.003102, 0.002247, 0.001612, 0.001145, 0.000806, 0.000563, 0.000392, 0.000275, 0.000197, 0.00015, 0.000128], atol=0.00001)

     dist = BootstrapCrystDistributionKDE(1:10, ones(10))
     @test isapprox(dist, [0.500228, 0.535143, 0.569056, 0.601647, 0.632641, 0.661805, 0.688961, 0.713982, 0.736794, 0.757375, 0.775748, 0.79198, 0.806169, 0.818443, 0.828947, 0.837842, 0.845291, 0.851459, 0.856502, 0.86057, 0.863796, 0.866299, 0.868179, 0.869517, 0.870375, 0.870794, 0.870794, 0.870375, 0.869517, 0.868179, 0.866299, 0.863796, 0.86057, 0.856502, 0.851459, 0.845291, 0.837842, 0.828947, 0.818443, 0.806169, 0.79198, 0.775748, 0.757375, 0.736794, 0.713982, 0.688961, 0.661805, 0.632641, 0.601647, 0.569056, 0.535143, 0.500228, 0.464656, 0.428798, 0.39303, 0.357728, 0.323251, 0.289935, 0.258079, 0.227938, 0.199721, 0.173583, 0.149624, 0.127894, 0.108393, 0.091075, 0.075858, 0.062627, 0.051244, 0.041553, 0.03339, 0.026585, 0.020972, 0.016391, 0.012691, 0.009734, 0.007396, 0.005566, 0.004149, 0.003064, 0.002241, 0.001623, 0.001166, 0.00083, 0.000587, 0.000415, 0.000295, 0.000215, 0.000166, 0.000143], atol=0.00001)

## -- Test systematic uncertainty propagation

    dist = 100 .+ randn(100000)

    dist_upb = add_systematic_uncert_UPb(dist)
    @test std(dist_upb) > std(dist)
    @test isapprox(std(dist_upb) - std(dist), 0.0015, atol=0.0006)

    dist_ar = add_systematic_uncert_ArAr(dist)
    @test std(dist_ar) > std(dist)
    @test isapprox(std(dist_ar) - std(dist), 0.005, atol=0.002)

    dist2d = 100 .+ randn(1000,1000)

    dist_upb = add_systematic_uncert_UPb(dist2d)
    # Check that uncertainty is added in vertical bands
    diff = dist_upb .- dist2d
    @test all(std(diff,dims=2) .> std(diff,dims=1)')

    dist_ar = add_systematic_uncert_ArAr(dist2d)
    # Check that uncertainty is added in vertical bands
    diff = dist_ar .- dist2d
    @test all(std(diff,dims=2) .> std(diff,dims=1)')

## ---


text:
import Pkg; Pkg.activate(joinpath(@__DIR__, "../../../../../"))
using Distributions
using PyPlot

"""
Paper: Bayesian inference for finite mixtures of univariate and multivariate
       skew-normal and skew-t distributions, Biostatistics 2010.
skew (delta): a real number in (-1, 1)
"""
function rand_skewnormal(loc, scale, skew)
    z = rand(TruncatedNormal(0, 1, 0, Inf))
    return loc + scale * skew * z + scale * sqrt(1 - skew ^ 2) * randn()
end

function rand_skewnormal(loc, scale, skew, dims...)
    z = rand(TruncatedNormal(0, 1, 0, Inf), dims...)
    return loc .+ scale * skew * z + scale * sqrt(1 - skew ^ 2) * randn(dims...)
end

x = rand_skewnormal(1, .5, -.97, 100000)
plt.hist(x, bins=100);

mean(x .< 0)

mean(x), std(x)


text:
module GitTest

using ForwardDiff

include("extra_file01.jl")

export f, f2, dfdx

end # module


text:
using Parameters

##############################
# Parameters
##############################

# conversion factor
const day2sec = 24*60*60
const year2sec = 365*day2sec

# "Convert m^3/s in Gt/a"
# qms2Gta(flux, rho) = flux*rho/1e3/1e9*year2sec

# A ice flow factor lookup table (fixed):
# Temperature dependence of A (Table 3.4, p.75), for n=3 [Pa^-3 s^-1]
const As = (0.0:-5:-50,[ 2.4e-24
                         9.3e-25
                         3.5e-25
                         2.1e-25
                         1.2e-25
                         6.8e-26
                         3.7e-26
                         2.0e-26
                         1.0e-26
                         5.2e-27
                         2.6e-27])
const _A_itr = Interp.interpolate((reverse(As[1]),), reverse(As[2]), Interp.Gridded(Interp.Linear()))
"""
Ice flow rate factor A as function of temperature.  Only valid for n=3

Interpolated from a lookup table (C&P p. 75)

Input:
- temp: ice temp
- pp:: physical para

Out:
- A
"""
ice_flow_factor(temp) = _A_itr[temp]

"""
$(TYPEDEF)
$(FIELDS)

Fairly fixed physical parameters, not used in tuning.

Units:
- kg, m, s, Celsius
"""
@with_kw struct Phys{FN<:Function} <:APara @deftype F
    # general
    g=9.81 # [m/s^2]
    # ice (references are to Cuffey&Patterson 2010)
    rhow=1000      # water density
    rhoi=910       # ice density
    rho=870        # glacier density including firn layer (kg m-3)
    rhosw=1025     # sea water density (kg m-3)
    n=3            # Glen's n
    @assert n==3   # needs to be ==3! for the A-table to be valid
    A::FN = ice_flow_factor
    r=(n+1)/(n+2) # (average deformational speed) / (surface deformational speed) for shallow ice
    @assert r==(n+1)/(n+2)
end

"""
What error is fitted:
- sqrerr: corresponds to Normal distribution of errors
- abserr: corresponds to Laplac distribution of errors
"""
baremodule ErrorModel
using Base: @enum, getindex
@enum(_ErrorModel, sqrerr, abserr)
end

"""
$(TYPEDEF)
$(FIELDS)

Model parameters.  These are not particularly fixed, ranging from
educated to wild guesses.  Some may be used in the inversion, some define the
inversion.

TODO:
"""
@with_kw struct MPara<:APara @deftype F
    # Standard deviation of errors due forward model limitations
    # (assuming perfectly tuned parameters):
    sigma_h_model = 10 # of ice thickness
    sigma_iv_model = 10 # of surface velocity
    @assert sigma_h_model>=0
    @assert sigma_iv_model>=0
    # total sigma is: sqrt(pm.sigma_h^2 + gl.h.sigma^2)

    # whether to assume Normal-distributed errors (:sqrerr)
    # or Laplace-distributed errors (:abserr)
    error_model::ErrorModel._ErrorModel=[ErrorModel.sqrerr, ErrorModel.abserr][1]

    # poor man's avoidance of spatial correlation: only use every n'th datapoint:
    nthin_h::Int=1
    nthin_iv::Int=10

    ## 2D->1D
    #########

    # Smooth dem to get smooth alpha, smoother bands.  This is in
    # particular important when there is cross-flow bumpiness, such as
    # on Unteraar.  However, it can also be bad.  YMMV, check!
    # This is then also used for h and iv extrapolation. The window
    # should probably be chosen around the size of the offending surface
    # features, or zero.
    window_dem_smooth = 100

    # Allowed minimal/maximal slope of bands
    slope_min = deg2rad(0.4)
    slope_max = deg2rad(60.0)

    bandsize = 30 # (m) arguably this should go into Num

    # Smooth the width of the bands over +/- this elevation range (m)
    # This makes the 1D iv much smoother, however it can make the
    # malphas much noisier!  In particular it can introduce bad
    # artifacts at the edges!
    window_width_smooth = 50

    ## 1D
    #####
    ## Scalars
    # Distance d over which to average tau in units of local ice thickness.
    # Takes average over elevation bands within +/- mean_tau_dist*h; no distance weighting is done.
    mean_tau_dist = 3.0  # Kamb and Echelmeyer 1986 (Fig 2 & Conclusions) suggest a triangularly-weighted
                         # window of length ±2*(1.5 to 2.5 ice thicknesses)

    # use Nye's ice-flow shape factor, or not
    shapeF::Bool = true

    # Whether to use the magic flux corrector
    use_flux_corrector::Bool=false

    # The ice thickness at the top of the glacier.  If set to
    # negative, set ice thickness at top equal to ice thickness one
    # below (good for ice-caps).
    # If set to NaN it will use -1 if isicecap==true and 0 otherwise.
    h_at_top = NaN

    ## 1D->2D
    #########

    # smoothing bed-dem with boxcar filter, probably around mean ice thickness:
    window = 500 # [m]

    # threshold on slopes used in 2D extrapolation
    alpha_min = deg2rad(1.5)

    # max distance to land-margin to consider
    dist_land_maxdist = 3e3 # [m] in Matthias' code: crit_distnu
    dist_exp = 0.3 # exponent to use in the distance to land-margin.  <1 gives more U-shape
    @assert dist_exp>0

    # Cap ice thickness at flotation thickness.  (Setting to true will violate mass-cons)
    cap_at_floatation::Bool = false

    ## IV extrapolation
    #
    # Trial u is made with the formula:
    #
    # !iv_nye *(
    # taus[i]^iv_tau_exp * # surface slope term on band
    # hs2d[jj]^iv_h_exp *  # ice thickness term
    # d[jj]^iv_dist_exp1) + # distance to land term 1
    # iv_nye * (1 - ( (len_scale-d2)/len_scale )^iv_dist_exp2 ) # Nye 1965 inspired
    #
    # Shallow ice is: tau^n*h^(n+1)
    iv_tau_exp = [0, 3][2] # no-effect at 0
    @assert iv_tau_exp>=0
    iv_h_exp = [0, 4][2] # no-effect at 0
    @assert iv_h_exp>=0
    iv_dist_exp1 = [0, 0.3][2] # no-effect at 0
    @assert 0<=iv_dist_exp1<=1

    # Nye 1965 extrapolation:
    iv_nye::Bool = true # to switch-on Nye 1965-extrapolation (and switch off Huss-extrapolation)
    iv_dist_exp2 = [3+1, 1][1] # larger makes IV more even flow-field across glacier, n+1 according to theory.
                               # Ceases to have much effect when >10.
    @assert iv_dist_exp2>0

    # How the trial u is scaled.  Best :masscons
    iv_extrapolation_scaling::Symbol = [:huss, :masscons, :none][2]
    # :huss has no parameters

    # :masscons parameters
    # Over how many grid cells to average the flux direction.  Does not
    # seem to make much difference (at least when using a smoothed DEM).
    iv_flux_dir_window = 300 # [m]
    iv_window_frac = 0.7 # smoothing distance as fraction of max band width

    # :none parameters
    iv_scale_factor = 1 # a scale factor.  Only needed if :none is used.

    # Afl_all::Vector{F}=[0.0255,0.0175,0.0405,0.0485,0.05,0.025,0.014,0.0255,0.0255]               # flow rate factor

    # # cs parabelflörmig über den Gletscher verteilt
    # cs=0.              # minimaler Wert der Korrektur-Funktion
    # rancurv::Vector{F}=[4,6]      # Exponent

    # # Options for redistribution of ice thickness within elevation band
    # fact_ip=0.75
    # #fact_ip=0.5           # exponent regulating distribution of thickness values in the cross section
    # # (0<fact<1) low: U-valley, high: V-Valley
    # slp_crit=1.5             # Cut off threshold for slope influence [deg]

    # rho_glacier=870.     # kg m-3 (glacier density including firn layer)

    # # surface elevation changes
    # dh_dtelev_boundary=30    # % lowermost x %
    # dh_dtsurf::Vector{F}=[-1,-7.5]       # m/a (buttressed / not buttressed)
    # #dh_dtsurf::Vector{F}=[-0.25,-3]       # m/a (buttressed / not buttressed)
    # #dh_dtsurf::Vector{F}=[-2,-15]       # m/a (buttressed / not buttressed)

    # meanslope_lower=.4      # [deg] Lower boundary for slope in elevation band ...
    # coldice_temp=-8
    # fact_a=6.
    # crit_distnu=3000 # [m] critical distance from nunataks - no influence when further away!

end

"""
$(TYPEDEF)
$(FIELDS)

Numerical parameters for deterministic model.
"""
@with_kw mutable struct Num<:APara @deftype F
#    bstep=10       # elevation band discretisation
    reltol=0.05    # relative tolerance for iterative solver
    minthick=1.0    # minimal ice thickness, acts as a abstol
    niter::Int=100 # max number of iterations in Picard loop
    relaxfac=0.1   # relax the Picard iteration: 0 -> no relaxation, 0.1 recommended
    relaxmin=0.1  # relax the Picard iteration: 0 -> no relaxation, 0.05  recommended (TODO?)
    @assert relaxmin>=2*reltol
    # performance
    calc_boxcar_M::Bool=true # set to false when only doing a few forward runs.
    # output
    plotyes::Bool = false
    verbose::Bool = false
    # misc flags
    test::Bool=true # activate some tests
    error_on_neg_qtot::Bool = false # if false return h=0 where qtot<0, otherwise error
end

##############################
# Data
##############################

const SomeDataValType = Union{F,AGridded{F},Traj{F},Void}
"""
$(TYPEDEF)

Holds measurements or synthetic guesses (scalar, 1D-trajectory,
1D-elevation band data, 2D or missing) and their errors expressed as
standard-deviation and ranges (once collapsed into 1D elevation
bands).

See function `logprior_scalar` and `logprior_1dfields` for implementations.

Fields:

    vals::Union{F,AGridded{F},Traj{F},Void}}
    sigma::F # Standard deviation of error
    range_rel::Tuple{F,F} # Range of allowed values relative to local value
    range_rel_cutoff::F # only check range_rel for local values `lv` with `abs(lv)>range_rel_cutoff`
    range_abs::Tuple{F,F} # Range of allowed values in absolute values
    smooth_window::F # elevation range to smooth its 1D version over
                       (not applicable to scalars).  This essentially gives
                       a smoothness constraint on the data
    vals_prior::Union{F,AGridded{F},Traj{F},Void}} # use this value if the vals produce FILL or NaN

TODO:
- add a field which give the time/date range of the measurements.
- think about how to transform a, say, sigma which applies to 2D data,
  into a sigma which applies to the corresponding elevation band data.
- allow non-scalar sigma, etc.
"""
@with_kw struct SomeData{V<:SomeDataValType, VV<:SomeDataValType}
    vals::V  # mean/expectation value
    sigma::F=Inf # Standard deviation of vals distribution (Gaussian)
    range_rel::Tuple{F,F}=(-Inf,Inf) # Range of allowed values relative to local expectation value
    range_rel_cutoff::F=Inf # only check range_rel for local values `lv` with `abs(lv)>range_rel_cutoff`
    range_abs::Tuple{F,F}=(-Inf,Inf) # Range of allowed values in absolute values
    smooth_window::F=0 # elevation range to smooth its 1D version over
    vals_prior::VV=F(NaN) # use this value if the vals produce FILL or NaN
end
SomeData(val::V,sigma,range_rel,range_rel_cutoff,range_abs,smooth_window, vals_prior::VV) where {V,VV} =
    SomeData{V,VV}(val,sigma,range_rel,range_rel_cutoff,range_abs,smooth_window, vals_prior)
SomeData(val::V) where {V} = SomeData{V,F}(val, Inf, (-Inf,Inf), Inf, (-Inf,Inf), 0.0, NaN)
SomeData(val::V,sigma,range_rel,range_rel_cutoff,range_abs) where {V<:Union{F,Void}} =
    SomeData{V,F}(val,sigma,range_rel,range_rel_cutoff,range_abs,0.0,NaN)
hasdata(::SomeData{Void}) = false
hasdata(::SomeData{F}) = true
hasdata(sd::SomeData) = length(sd.vals.v)>0

"""
Regions of the world according to RGI and some extras (with numbers >99)

See also: symbol2region, region2symbol
"""
baremodule Region
using Base: @enum, getindex
@enum(_Region,
      Alaska=1,
      WesternCanadaUSA,
      ArcticCanadaNorth,
      ArcticCanadaSouth,
      Greenland,
      Iceland,
      Svalbard, # and Jan Mayen
      Scandinavia,
      RussianArctic,
      NorthAsia,
      CentralEurope,
      CaucasusMiddleEast,
      CentralAsia,
      SouthAsiaWest,
      SouthAsiaEast,
      LowLatitudes,
      SouthernAndes,
      NewZealand,
      SubAntarctic,
      # Misc
      NotKnown=100,
      Synthetic=101,
      # Other regions of the world
      AntarcticPeninsula=110,
      # ITMIX
      ITMIX=201,
      ITMIX2=202,
      )
end
let
    sym = Symbol.(instances(Region._Region))
    int = Int.(instances(Region._Region))
    global symbol2region(sym) = Region._Region(int[findfirst(x->x==sym, Region.sym)])
    global region2symbol(reg) = Symbol(reg)
end

# """
# The different kinds of surfaces we need to consider:

# - `this_glacier` -- the glacier under study
# - `other_glacier` -- where it borders on other glaciers
# - `land_margin` -- where it borders on solid land
# - `calving_margin` -- where it borders on something it can calf into
#                       (water for tide-water glaciers, air for hanging glaciers)
# """
# @enum SurfaceTypes this_glacier other_glacier_margin land_margin calving_margin

"""
$(TYPEDEF)

Glacier ID, type to dispatch on a particular glacier.  As
type-parameter use a symbol of the RGI, GLIMS id or some other ID.
"""
abstract type GlacierID end
Base.:(==)(a::GlacierID, b::GlacierID) = getid(a)==getid(b)

"Return glacier ID, usually a symbol or Int"
getid(gid::GlacierID) = gid.id
"Return glacier RGI number, a string"
getrgi(::GlacierID)::String = error("to be defined")
"Return the name of a glacier, if it has one (string)"
getname(gid::GlacierID)::String = gid.name
"Return the region (an instance of the `Region` enumeration)"
getregion(gid::GlacierID)::Region = gid.region

"""
$(TYPEDEF)

Glacier outline and outcrops.

The outline gives the glacier outline, usually including outcrops.
However, outcrops can also be specified with separate polygons.

The polygons follow the "bigpoly" convention as in
https://github.com/mauro3/VAWTools.jl/blob/67db8fe62daabfc10d313a46e827349644738ca1/src/VAWTools.jl#L607

    gid::GID
    outline::Matrix{Float64} # size==(2,n)
    splits::Vector{Int} # i-th poly has indices splits[i]:splits[i+1]-1
    outcrops::Matrix{Float64}
    outcrops_splits::Vector{Int}

"""
struct Outline{GID<:GlacierID}
    gid::GID
    outline::Matrix{Float64} # size==(2,n)
    splits::Vector{Int}
    outcrops::Matrix{Float64}
    outcrops_splits::Vector{Int}
end
Outline(gid::GID, o, s) where {GID} =
    Outline{GID}(gid, o, s, Matrix{Float64}(0,0), Int[])
Outline(gid::GID, o) where {GID} =
    Outline{GID}(gid, o, [1,size(o,2)+1] , Matrix{Float64}(0,0), Int[])
hasoutcrops(outl::Outline) = length(outl.outcrops)>0

"Approximate area of an outline"
function area(outline::Outline, gridpts=100)
    ol = outline.outline
    dx = (maximum(ol[1,:])-minimum(ol[1,:]))/gridpts
    dy = (maximum(ol[2,:])-minimum(ol[2,:]))/gridpts
    dx = min(dx,dy)
    x = minimum(ol[1,:]):dx:maximum(ol[1,:])
    y = minimum(ol[2,:]):dx:maximum(ol[2,:])
    pts = 0
    for xx in x, yy in y
        if inpoly((xx,yy), ol)
            pts +=1
        end
    end
    return pts*dx^2
    # TODO: remove outcrops
end

# ## refactor Glacier to use this instead of iscalving at some point
# """
# $(TYPEDEF)

# Holds glacier attributes, roughly following RGI.
# """
# @with_kw_noshow struct GlacierAttrs
#     form::GlacierAttrs.Form
#     terminus::GlacierAttrs.Terminus
#     surging::GlacierAttrs.surging = GlacierAttrs.not_surging
# end

# module GlacierAttrs
# @enum(Form,
#       glacier
#       glacier_stystem
#       icecap
#       icecap_catchment
#       )
# @enum(Terminus,
#       land,
#       marine,
#       lake,
#       dry_calving,
#       regenerated,
#       shelf)
# @enum(Surging,
#       not_surging,
#       possible,
#       probable,
#       observed)
# end

"""
Enumeration for the different types of teminus-flux specification available:

- abs: absolute value [m^3/s]
- rel_bdot_teminus: as ratio in terms of flux from integrating
                    expectation of prior of bdot
- rel_btilde_teminus: as ratio in terms of flux from integrating
                      expectation of prior of btilde
- rel_bdot_max: as ratio in terms of maximum flux from integrating
                expectation of prior of bdot
- none: do not look at the flux
"""
baremodule FluxKind
using Base
@enum _FluxKind abs rel_bdot_teminus rel_btilde_teminus rel_bdot_max none
end

"""
$(TYPEDEF)

Holds all input for a particular glacier, i.e. all our prior
knowledge of a glacier. This data stays fixed during fitting.

Note that no fields should contain NaNs
"""
@with_kw_noshow struct Glacier
    gid::GlacierID # RGI 4.0 number except for Peninsula
    iscalving::Bool # true for calving glacier (tide-water, fresh-water, hanging glacier)
    isicecap::Bool # true for ice-cap like glacier with several outlets
    proj4::String # map projection, a Proj4 string

    # DEM (assumed error free)
    dem::Gridded{F}         # surface DEM (cropped) [m].  Note that
                            # the dem.x,dem.y is also the
                            # computational 2D grid.
    outline::Outline # glacier outline and possibly additional outcrops
    cropbox::Tuple{Float64,Float64,Float64,Float64} # crop-box
    # "True on the glacier"
    glaciermask::Matrix{Bool}
    @assert size(dem)==size(glaciermask)
    # "True on solid ground"
    landmask::Matrix{Bool}
    @assert size(dem)==size(landmask)
    # "distance to closest land point (SurfaceTypes: land_margin)"
    dist_land::Matrix{F}
    # "dist_land was calculated up to this distance"
    dist_land_maxdist_calc::F
    # DEM products, all calculated
    alpha::Matrix{F}       # surface slope [radians] (prod)
    elemean::F
    elemax::F
    elemin::F
    elemedian::F

    # Distributed priors
    ###############
    # Measurements/estimates of input fields (1D, 2D or nothing)
    bdot::SomeData        # mass balance  [m ice /a]
    dhdt::SomeData        # surface elevation change [m ice/a]
    temp::SomeData        # ice temperature [C]
    fsl::SomeData         # basal sliding fraction []

    # Scalar priors
    ###############
    # Flux at terminus can be set to a few different settings:
    # - (:abs, SomeData(flux...)) : absolute flux at terminus [m^3/s].
    # - (:rel_bdot_teminus, SomeData(flux...)) : as fraction of the
    #    integrated bdot at the terminus --> calving
    # - (:rel_btilde_teminus, SomeData(flux...)) : as fraction of the
    #    integrated btilde at the terminus --> calving
    # - (:rel_bdot_max, SomeData(flux...)) : as fraction of the
    #    maximum of the integrated bdot --> calving or land terminating
    #
    # Probably (:abs, SomeData(0.0)) for land-terminating glaciers,
    # and (:rel_bdot_terminus, SomeData(1.0)), or similar, for
    # tide-water glaciers.  Sigma is small for land-terminating and
    # big for tide water.  Range: can restrict it to be non-negative
    # (thus h>0).
    #
    # Note that for all rel_* kinds the SomeData fields vals, sigma and range_abs,
    # are in fact relative to the "*".  E.g. sigma for rel_bdot_max, corresponds to
    # an absolute sigma of bdot_max*sigma.
    terminus_flux::Tuple{FluxKind._FluxKind, SomeData} # [m^3/s] or []

    ## Measurements of forward model output fields
    #####
    h::SomeData   # Traj{F}: measured ice thickness at radar lines, abuses .vals.err field to hold surface elevation [m]
    iv::SomeData    # iv velocity: measured surface flow speed: on different grid to DEM [m/a], possibly at points
    ivmask::Matrix{Bool} # true on the glacier, matches grid of .iv

    ###
    misc::Dict{Symbol,Any}=Dict{Symbol,Any}() # for stuff
end
## This does not work with @with_kw:
# function Glacier(args...)
#     mask = args[8]
#     for (i,arg) in enumerate(args)
#         if _hasnan(arg)
#             error("Argument #$i of Glacier constructor call contains NaNs! Not allowed to enable certain performance tricks.")
#         end
#         if _notmasked(arg, mask)
#             error("Argument #$i of Glacier constructor call contains unmasked FILL=$FILL values!")
#         end
#     end
#     Glacier(args...)
# end

getname(gl::Glacier) = getname(gl.gid)
getregion(gl::Glacier) = getregion(gl.gid)

Base.show(io::IO, gl::Glacier) = print(io, """Glacier $(getid(gl.gid)) with name "$(getname(gl))" in region "$(getregion(gl))" """)
area(gl::Glacier) = sum(gl.glaciermask) * step(gl.dem.x)^2

volume_area_mean_h(gl::Glacier) = volume_area_mean_h(area(gl), isicecap=gl.isicecap)

"Make sure there are no NaN lurking"
_hasnan(arg::AbstractArray{T}) where {T<:AbstractFloat} = any(isnan.(arg))
function _hasnan(arg::AGridded{T}) where {T<:AbstractFloat}
    VAWTools.haserror(arg) ? (_hasnan(arg.v) && _hasnan(arg.v)) : _hasnan(arg.v)
end
_hasnan(arg::SomeData) = hasdata(arg) ? _hasnan(arg.vals) : false
_hasnan(arg::AbstractFloat) = isnan(arg)
_hasnan(arg::Tuple) = mapreduce(_hasnan, |, false, arg)
_hasnan(arg::Any) = false

function _notmasked(arg, mask)
    if isa(arg,AbstractArray) && eltype(arg)<:AbstractFloat
        return any(arg[mask].==FILL)
    elseif isa(arg,AGridded) && eltype(arg.v)<:AbstractFloat
        return any(arg.v[mask]==FILL)
    elseif isa(arg,SomeData)
        hasdata(arg) ? _notmasked(arg.vals) : false
    elseif isa(arg, Tuple)
        return mapreduce(x->_notmasked(x,mask), |, false, arg)
    else
        return false
    end
end

# """
# Container which holds all needed parameters for a model run.
# """
# mutable struct Model
#     ph::Phys
#     pm:MPara
#     dirs::Dirs
#     gl::Glacier
#     num::Num
#     # res::Result
#     # Model(ph,pa,dirs,gl,num,res) = new(ph,pa,dirs,gl,num,res)
#     Model(ph,pa,dirs,gl,num)     = new(ph,pa,dirs,gl,num)
# end


text:
module toJuMP

using JuMP
using JSON

include("const.jl")
include("types.jl")
include("writer.jl")
include("utility.jl")

include("gms2jump.jl")
include("mod2jump.jl")

!isdir("$(Pkg.dir())/toJuMP/.jls") && mkdir("$(Pkg.dir())/toJuMP/.jls")

m_tester = Model()

export gms2jump, mod2jump

end # module


text:
# LibTOAST.jl: interface to the TOAST++ library
# Copyright (C) 2019 Samuel Powell

# Imports

# Export types
export Regul, RegulTK0, RegulTK1, RegulTV, RegulPM

# Export methods
export val, grad, hess

"""
    Regul

Abstract supertype of regulariers capable of calculating the value, gradient,
and Hessian of of a particular regularisation functional, with parameters
defined in a SolutionCoeff raster basis.

See RegulTK0, RegulTK1, RegulPM, RegulTV
"""
@compat abstract type Regul end

"""
  RegulTK0(x₀)

Construct a zeroth-order Tikhonov regularisation functional with baseline
parameters x₀ (defined in a SolutionCoeff basis).

``f(x) = ‖x-x₀‖²``

"""
mutable struct RegulTK0 <: Regul

  ptr::Cxx.CxxCore.CppPtr
  rast::Raster

  function RegulTK0(x0::SolutionCoeff)

    rast = x0.rast
    xptr = pointer(x0.data)
    xlen = length(x0)

    regulptr = icxx"""
      RVector x0($(xlen), $(xptr), SHALLOW_COPY);
      RVector xs($(xlen), 1);
      Regularisation *reg = new Tikhonov0 (1.0, &x0, &xs);
      return reg;
    """

    regul = new(regulptr, rast)
    finalizer(_regul_delete, regul)
    return regul
  end

end

"""
  RegulTK1(x₀)

Construct a first-order Tikhonov regularisation functional with baseline
parameters x₀ (defined in a SolutionCoeff basis).

``f(x) = ‖∇(x-x₀)‖²``
"""
mutable struct RegulTK1 <: Regul

  ptr::Cxx.CxxCore.CppPtr
  rast::Raster

  function RegulTK1(x0::SolutionCoeff)

    rast = x0.rast
    rptr = rast.ptr
    xptr = pointer(x0.data)
    xlen = length(x0)

    regulptr = icxx"""
      void *kapref = 0;
      bool istensor = false;
      RVector x0($(xlen), $(xptr), SHALLOW_COPY);
      Regularisation *reg = new TK1 (1.0, &x0, $(rptr), kapref, istensor);
      return reg;
    """

    regul = new(regulptr, rast)
    finalizer(_regul_delete, regul)
    return regul
  end

end

"""
    RegulTV(x₀, β=1.0)

Construct a soft Total-Variation regularisation functional with parameter
β and baseline parameters x₀ (defined in a SolutionCoeff basis).
"""
mutable struct RegulTV <: Regul

  ptr::Cxx.CxxCore.CppPtr
  rast::Raster

  function RegulTV(x0::SolutionCoeff, β::Float64)

    rast = x0.rast
    rptr = rast.ptr
    xptr = pointer(x0.data)
    xlen = length(x0)

    regulptr = icxx"""
      RVector x0($(xlen), $(xptr), SHALLOW_COPY);
      Regularisation *reg = new TV(1.0, $(β), &x0, $(rptr));
      return reg;
    """

    regul = new(regulptr, rast)
    finalizer(_regul_delete, regul)
    return regul

  end

end

function RegulTV(x0::SolutionCoeff)
  @info "Using default TV β = 1.0"
  return RegulTV(x0, 1.0)
end

"""
    RegulPM(x₀, T=1.0)

Construct a Perona-Malik regularisation functional with parameter T, and
baseline parameters x₀ (defined in a SolutionCoeff basis).
"""
mutable struct RegulPM <: Regul

  ptr::Cxx.CxxCore.CppPtr
  rast::Raster

  function RegulPM(x0::SolutionCoeff, T::Float64)

    rast = x0.rast
    rptr = rast.ptr
    xptr = pointer(x0.data)
    xlen = length(x0)

    regulptr = icxx"""
      RVector x0($(xlen), $(xptr), SHALLOW_COPY);
      Regularisation *reg = new PM(1.0, $(T), &x0, $(rptr));
      return reg;
    """

    regul = new(regulptr, rast)
    finalizer(_regul_delete, regul)
    return regul

  end

end

function RegulPM(x0::SolutionCoeff)
  @info "Using default PM T = 1.0"
  return RegulPM(x0, 1.0)
end

# Delete a raster map
_regul_delete(regul::T) where {T<:Regul} = finalize(regul.ptr)

# Regul val
"""
  val(regul, x)

Return the value of the regularisation functional defined by `regul` evaluated
with parameters `x` (defined in a SolutionCoeff basis).
"""
function val(regul::T, x::SolutionCoeff) where {T<:Regul}

  @assert regul.rast == x.rast

  xptr = pointer(x.data)
  xlen = length(x)

  icxx"""
    RVector x($xlen, $(xptr), SHALLOW_COPY);
    return $(regul.ptr)->GetValue(x);
  """

end

# grad of regul
"""
  grad(regul, x)

Return the gradient of the regularisation functional defined by `regul` evaluated
with parameters `x` (defined in a SolutionCoeff basis).
"""
function grad(regul::T, x::SolutionCoeff) where {T<:Regul}

  @assert regul.rast == x.rast

  g = SolutionCoeff(regul.rast)

  gptr = pointer(g.data)
  xptr = pointer(x.data)
  xlen = length(x)

  icxx"""
    RVector x($xlen, $xptr, SHALLOW_COPY);
    RVector g($xlen, $gptr, SHALLOW_COPY);
    g = $(regul.ptr)->GetGradient(x);
  """

  return g

end

# Hessian of regul
"""
  val(regul, x)

Return the Hessian of the regularisation functional defined by `regul` evaluated
with parameters `x` (defined in a SolutionCoeff basis).
"""
function hess(regul::T, x::SolutionCoeff) where {T<:Regul}
  #
  # assert(regul.rast == x.rast)
  #
  # xptr = pointer(x)
  # xlen = length(x)

  # icxx"""
  #
  #   RVector x(xlen, xvec, SHALLOW_COPY);
  #   RCompRowMatrix *H = new RCompRowMatrix;
  #
  #   int nprm = regul->GetNParam();
  #   int n = x.Dim();
  #   int n0 = n/nprm;
  #   int i, j;
  #
  #   RCompRowMatrix Hi, Hij;
  #   for (i = 0; i < nprm; i++) {
  #       for (j = 0; j < nprm; j++) {
  #           if (!j) {
  #               Hi.New(n0,n0);
  #               if (j==i) regul->SetHess1 (Hi, x, j);
  #           } else {
  #               Hij.New(n0,n0);
  #               if (j==i) regul->SetHess1 (Hij, x, j);
  #               Hi = cath (Hi, Hij);
  #           }
  #       }
  #       if (!i) *H = Hi;
  #       else    *H = catv (*H, Hi);
  #   }
  #
  #   int nval = H->nVal();
  #
  #   lengthptrs[0] = nval;
  #   lengthptrs[1] = n0+1;
  #
  #   return H;
  # """
  #
  # icxx"""
  #   memcpy(colidx, H->colidx, H->nVal()*sizeof(int));
  #   memcpy(rowptr, H->rowptr, (H->nRows()+1)*sizeof(int));
  #   memcpy(val, H->ValPtr(), H->nVal()*sizeof(double));
  #
  #   delete H;
  # """

end


text:
using Distributions

if VERSION < v"0.5.0-"
    function Markdown.rstinline(io::IO, md::Markdown.Link)
        if ismatch(r":(func|obj|ref|exc|class|const|data):`\.*", md.url)
            Markdown.rstinline(io, md.url)
        else
            Markdown.rstinline(io, "`", md.text, " <", md.url, ">`_")
        end
    end
end

function printrst(io,md)
    mdd = md.content[1]
    sigs = shift!(mdd.content)

    decl = ".. function:: "*replace(sigs.code, "\n","\n              ")
    body = Markdown.rst(mdd)
    println(io, decl)
    println(io)
    for l in split(body, "\n")
        ismatch(r"^\s*$", l) ? println(io) : println(io, "   ", l)
    end
end

cd(joinpath(dirname(@__FILE__),"source")) do
    for (name,T) in [("Univariate Discrete",DiscreteUnivariateDistribution),
                     ("Univariate Continuous",ContinuousUnivariateDistribution),
                     ]

        fname = replace(lowercase(name),' ','-')
        open("$fname.rst","w") do f
            println(f,".. _$fname:")
            println(f)
            println(f,"$name Distributions")
            println(f,"----------------------------------------------------")
            println(f)
            println(f, ".. DO NOT EDIT: this file is generated from Julia source.")
            println(f)

            for D in subtypes(T)
                if isleaftype(D)
                    md = Base.doc(D)
                    if isa(md,Markdown.MD)
                        isa(md.content[1].content[1],Markdown.Code) || error("Incorrect docstring format: $D")

                        printrst(f,md)
                    else
                        warn("$D is not documented.")
                    end
                end
            end
        end
    end
end


text:
#############################################################################
#############################################################################
#
# This file implements polynomial addition
#
#############################################################################
#############################################################################

"""
Add two polynomials.
"""
function +(p1::Polynomial, p2::Polynomial)::Polynomial
    p1, p2 = deepcopy(p1), deepcopy(p2)
    p3 = Polynomial()
    while !iszero(p1) && !iszero(p2)
        t1, t2 = trailing(p1), trailing(p2)
        if t1.degree == t2.degree
            push!(p3, popfirst!(p1)+popfirst!(p2))
        elseif t1.degree < t2.degree
            push!(p3,popfirst!(p1))
        else
            push!(p3,popfirst!(p2))
        end
    end
    while !iszero(p1)
        push!(p3,popfirst!(p1))
    end
    while !iszero(p2)
        push!(p3,popfirst!(p2))
    end
    return p3
end

"""
Add a polynomial and a term.
"""
+(p::Polynomial, t::Term) = p + Polynomial(t)
+(t::Term, p::Polynomial) = p + t

"""
Add a polynomial and an integer.
"""
+(p::Polynomial, n::Int) = p + Term(n,0)
+(n::Int, p::Polynomial) = p + Term(n,0)


text:
using Roots
using LinearAlgebra
using SparseArrays
import FiniteDiff.finite_difference_jacobian
using Plots
pyplot()
# %%
using ReiterJulia
# %%
n_p = 100
n_d = 200
n_ξ = 41
# %%
d = DiscretizedModel()
params = Parameters()
initParameters!(params)
P = Prices()
A = Aggregates()
z = 0.
@save A z
# %%
initDiscretizedModel!(d, n_p, n_d, n_ξ, params)
# %%
s = 0.8*d.x  # houdeholds' policy function for capital taken at points in the earnings grid x
# %%
Kss = find_zero(x->excessCapital!(d, A, P, params, s, x, verbose=false), (0.95, 1.05), verbose=true, xatol = eps())
# %%
A.K = Kss
# %%
Π = setTransition!(d, P, params, s)
M = Π - I(d.n_d+1)
M[:,end] .= 1
b = zeros(d.n_d+1)
b[end] = 1
p = M'\b
# %%
println("Check steady-state residual: ", A.K - capitalSupply(p,d))
# %%
plot(d.x_fine, p, color="black", title = "Steady-state capital distribution")
xgrid!(:off)
ygrid!(:off)
# %%
x = vcat(p[2:end], [A.K, A.z])
y = s
# %%
H_x = finite_difference_jacobian(u->H(u,x,y,y,d,params), x)
# %%
H_xp = finite_difference_jacobian(u->H(x,u,y,y,d,params), x)
# %%
H_y = finite_difference_jacobian(u->H(x,x,u,y,d,params), y)
# %%
H_yp = finite_difference_jacobian(u->H(x,x,y,u,d,params), y)
# %%
AA = hcat(H_xp, H_yp)
BB = -hcat(H_x, H_y)
# %%
# g_x is the jacobian of the decision function (controls) wrt states
# h_x is the jacobian of the transition function (states) wrt states
g_x, h_x, eu = solve_eig(AA, BB, d.n_d+2)
# %%
T = 60
dx = zeros(n_d+2,T)
dx[end,1] = 0.1
# %%
for t in 2:T
    dx[:,t] .= h_x*dx[:,t-1]
end
# %%
plot(d.k[2:end], dx[1:n_d,1]+x[1:n_d], label="t=0")
plot!(d.k[2:end], dx[1:n_d,20]+x[1:n_d], label="t=20")
plot!(d.k[2:end], dx[1:n_d,40]+x[1:n_d], label="t=40")
xgrid!(:off)
ygrid!(:off)
xlabel!("Capital distribution")
# %%
plot(1:T, x[end] .+ dx[end-1,:])
xlabel!("Time")
ylabel!("Aggregate capital")
xgrid!(:off)
ygrid!(:off)
# %%

max_stars_repo_name: NormannR/ReiterJulia

max_issues_repo_name: NormannR/ReiterJulia

max_forks_repo_name: NormannR/ReiterJulia


text:
module DpForPfc
  using Plots
  pyplot()

  include(joinpath(split(@__FILE__, "src")[1], "src/decision_making/markov_decision_process/dynamic_programming.jl"))
  include(joinpath(split(@__FILE__, "src")[1], "src/model/goal/goal.jl"))
  include(joinpath(split(@__FILE__, "src")[1], "src/model/puddle/puddle.jl"))

  function main(is_test=false)
    puddles = []

    dp = DynamicProgramming([0.2, 0.2, pi/18], Goal(-1.5, -1.5), puddles, 0.1, 10)

    counter = 0
    delta = 1e100

    if is_test == false
      while delta > 0.01
        delta = value_iteration_sweep(dp)
        counter += 1
        println("$(counter), $(delta)")
      end

      txt_path = joinpath(split(@__FILE__, "src")[1], "src/decision_making/partially_observable_mdp/policy.txt")
      fp = open(txt_path, "w")
      for i in dp.indexes
        p = dp.policy[i[1]+1, i[2]+1, i[3]+1, :]
        write(fp, "$(i[1]) $(i[2]) $(i[3]) $(p[1]) $(p[2])\n")
      end
      close(fp)

      txt_path = joinpath(split(@__FILE__, "src")[1], "src/decision_making/partially_observable_mdp/value.txt")
      fp = open(txt_path, "w")
      for i in dp.indexes
        v = dp.value_function[i[1]+1, i[2]+1, i[3]+1]
        write(fp, "$(i[1]) $(i[2]) $(i[3]) $(v)\n")
      end
      close(fp)
    else
      for i in 1:10
        delta = value_iteration_sweep(dp)
        counter += 1
      end
    end
  end
end

max_stars_repo_name: ShisatoYano/JuliaAutonomy

max_issues_repo_name: ShisatoYano/JuliaAutonomy

max_forks_repo_name: ShisatoYano/JuliaAutonomy


text:
"""
``ProjectedRotationAngle`` is an abstract type to be used
as the first argument for functions.
It represents the rotation angle
with respect to a reference axis.
"""
abstract ProjectedRotationAngle


text:
# GA seclections
# ==============

# Rank-based fitness assignment
# sp - selective linear presure in [1.0, 2.0]
function ranklinear(sp::Float64)
    @assert 1.0 <= sp <= 2.0 "Selective pressure has to be in range [1.0, 2.0]."
    function rank(fitness::Vector{<:Real}, N::Int)
        λ = length(fitness)
        idx = sortperm(fitness)
        ranks = zeros(λ)
        for i in 1:λ
            ranks[i] = ( 2 - sp + 2*(sp - 1)*(idx[i] - 1) / (λ - 1) ) / λ
        end
        return pselection(ranks, N)
    end
    return rank
end

# (μ, λ)-uniform ranking selection
function uniformranking(μ::Int)
    function uniformrank(fitness::Vector{<:Real}, N::Int)
        λ = length(fitness)
        idx = sortperm(fitness, rev=true)
        @assert μ < λ "μ should be less then $(λ)"
        ranks = similar(fitness, Float64)
        for i in 1:μ
            ranks[idx[i]] = 1/μ
        end
        return pselection(ranks, N)
    end
    return uniformrank
end

# Roulette wheel (proportionate selection) selection
function roulette(fitness::Vector{<:Real}, N::Int)
    prob = fitness./sum(fitness)
    return pselection(prob, N)
end

# Stochastic universal sampling (SUS)
function sus(fitness::Vector{<:Real}, N::Int)
    F = sum(fitness)
    P = F/N
    start = P*rand()
    pointers = [start+P*i for i = 0:(N-1)]
    selected = Array{Int}(undef, N)
    i = c = 1
    for P in pointers
        while sum(fitness[1:i]) < P
            i += 1
        end
        selected[c] = i
        c += 1
    end
    return selected
end

# Truncation selection
function truncation(fitness::Vector{<:Real}, N::Int)
    λ = length(fitness)
    @assert λ >= N "Cannot select more then $(λ) elements"
    idx = sortperm(fitness, rev=true)
    return idx[1:N]
end

# Tournament selection
function tournament(groupSize :: Int)
    @assert groupSize > 0 "Group size must be positive"
    function tournamentN(fitness::Vector{<:Real}, N::Int)
        selection = fill(0,N)

        nFitness = length(fitness)

        for i in 1:N
            contender = unique(rand(1:nFitness, groupSize))
            while length(contender) < groupSize
                contender = unique(vcat(contender, rand(1:nFitness, groupSize - length(contender))))
            end

            winner = first(contender)
            winnerFitness = fitness[winner]
            for idx = 2:groupSize
                c = contender[idx]
                if winnerFitness < fitness[c]
                    winner = c
                    winnerFitness = fitness[c]
                end
            end

            selection[i] = winner
        end
        return selection
    end
    return tournamentN
end

# Utils: selection
function pselection(prob::Vector{<:Real}, N::Int)
    cp = cumsum(prob)
    selected = Array{Int}(undef, N)
    for i in 1:N
        j = 1
        r = rand()
        while cp[j] < r
            j += 1
        end
        selected[i] = j
    end
    return selected
end


text:
using Downloads: download
using OrdinaryDiffEq
using Trixi

###############################################################################
# semidiscretization of the acoustic perturbation equations

equations = AcousticPerturbationEquations2D(v_mean_global=(0.0, -0.5), c_mean_global=1.0,
                                            rho_mean_global=1.0)

# Create DG solver with polynomial degree = 4 and (local) Lax-Friedrichs/Rusanov flux
solver = DGSEM(polydeg=4, surface_flux=flux_lax_friedrichs)

# Create unstructured quadrilateral mesh from a file
default_mesh_file = joinpath(@__DIR__, "mesh_five_circles_in_circle.mesh")
isfile(default_mesh_file) || download("https://gist.githubusercontent.com/andrewwinters5000/3c79baad6b4d73bb26ec6420b5d16f45/raw/22aefc4ec2107cf0bffc40e81dfbc52240c625b1/mesh_five_circles_in_circle.mesh",
                                       default_mesh_file)
mesh_file = default_mesh_file

mesh = UnstructuredMesh2D(mesh_file)

"""
    initial_condition_gauss_wall(x, t, equations::AcousticPerturbationEquations2D)

A Gaussian pulse, used in the `gauss_wall` example elixir in combination with
[`boundary_condition_wall`](@ref). Uses the global mean values from `equations`.
"""
function initial_condition_gauss_wall(x, t, equations::AcousticPerturbationEquations2D)
  v1_prime = 0.0
  v2_prime = 0.0
  p_prime = exp(-log(2) * (x[1]^2 + (x[2] - 25)^2) / 25)

  prim = SVector(v1_prime, v2_prime, p_prime, global_mean_vars(equations)...)

  return prim2cons(prim, equations)
end
initial_condition = initial_condition_gauss_wall

boundary_conditions = Dict( :OuterCircle  => boundary_condition_slip_wall,
                            :InnerCircle1 => boundary_condition_slip_wall,
                            :InnerCircle2 => boundary_condition_slip_wall,
                            :InnerCircle3 => boundary_condition_slip_wall,
                            :InnerCircle4 => boundary_condition_slip_wall,
                            :InnerCircle5 => boundary_condition_slip_wall )

# A semidiscretization collects data structures and functions for the spatial discretization
semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
                                    boundary_conditions=boundary_conditions)

###############################################################################
# ODE solvers, callbacks etc.

# Create ODE problem with time span from 0.0 to 300.0
tspan = (0.0, 300.0)
ode = semidiscretize(semi, tspan)

# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup
# and resets the timers
summary_callback = SummaryCallback()

# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results
analysis_callback = AnalysisCallback(semi, interval=100)

# The SaveSolutionCallback allows to save the solution to a file in regular intervals
save_solution = SaveSolutionCallback(interval=50, solution_variables=cons2state)

# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver
callbacks = CallbackSet(summary_callback, analysis_callback, save_solution)

###############################################################################
# run the simulation

# use a Runge-Kutta method with automatic (error based) time step size control
sol = solve(ode, RDPK3SpFSAL49(), abstol=1.0e-6, reltol=1.0e-6,
            save_everystep=false, callback=callbacks);
# Print the timer summary
summary_callback()


text:
module TestBase

using CUDA
using Test
using UnionArrays
using UnionArrays: unionof

const ==′ = isequal

function test_unionvector()
    xs = UnionVector(undef, CuVector, unionof(Float32, Missing), 3)
    fill!(xs, 1)
    CUDA.@allowscalar xs[2] = missing
    @test (CUDA.@allowscalar xs[1]) === 1.0f0
    @test (CUDA.@allowscalar xs[2]) === missing
    @test (CUDA.@allowscalar xs[3]) === 1.0f0
    @test collect(xs) ==′ [1.0f0, missing, 1.0f0]
end

function test_unionmatrix()
    A = UnionArray(undef, CuVector, unionof(Float32, Missing), (3, 2))
    fill!(A, 1)
    CUDA.@allowscalar A[2, 2] = missing
    @test (CUDA.@allowscalar A[1, 1]) === 1.0f0
    @test (CUDA.@allowscalar A[2, 2]) === missing
    @test (CUDA.@allowscalar A[3, 2]) === 1.0f0
    @test collect(A) ==′ [
        1.0f0 1.0f0
        1.0f0 missing
        1.0f0 1.0f0
    ]
end

end  # module


text:
module DebugPrinting

export setdebug!, withdebug, getdebug

const DEBUG_LEVEL = Dict{Symbol, Int}()

getdebug(mod::Symbol) = DEBUG_LEVEL[mod]
getdebug(mod::Module) = getdebug(nameof(mod))

function setdebug!(mod::Symbol, level::Integer)
    DEBUG_LEVEL[mod] = Int(level)
    return
end

setdebug!(mod::Module, level::Integer) = setdebug!(nameof(mod), level)

function withdebug(f, mod, level)
    lvl = DEBUG_LEVEL[]
    try
        setdebug!(mod, level)
        f()
    finally
        setdebug!(mod, lvl)
    end
end

macro debug(level, msg)
    esc(quote
        if getdebug(__module__) >= $level
            println(string("DEBUG: ", $(QuoteNode(__source__.file)), ":", $(QuoteNode(__source__.line)), " ", $msg))
        end
    end)
end

end # module


text:
using StatsBase

xs = sort(rand(1000))
ys = diff(rand(1001));
T = eltype(ys)

g = -0.3:0.025:1.2
nbins = length(g) - 1

ybinned = bin(median, g, xs, ys)
@test length(ybinned) == nbins
@test ybinned isa Vector{T}

ybinned = bin(quantile, g, xs, ys, [0.5])
@test length(ybinned) == nbins
@test ybinned isa Vector{T}

ybinned = bin_mean(g, xs, ys)
@test ybinned isa Vector
@test length(ybinned) == nbins


text:
using ..CUDASupport
using ..TopOpt: @init_cuda, CPU, GPU
@init_cuda()
import ..TopOpt: whichdevice
using ..GPUUtils, ForwardDiff

abstract type AbstractGPUPenalty{T} <: AbstractPenalty{T} end

whichdevice(::AbstractCPUPenalty) = CPU()
whichdevice(::AbstractGPUPenalty) = GPU()

function CUDAnative.pow(d::TD, p::AbstractFloat) where {T,TV,TD<:ForwardDiff.Dual{T,TV,1}}
    return ForwardDiff.Dual{T}(
        CUDAnative.pow(d.value, p), p * d.partials[1] * CUDAnative.pow(d.value, p - 1)
    )
end

struct GPUPowerPenalty{T} <: AbstractGPUPenalty{T}
    p::T
end
@inline (P::GPUPowerPenalty)(x) = CUDAnative.pow(x, P.p)

struct GPURationalPenalty{T} <: AbstractGPUPenalty{T}
    p::T
end
@inline (P::GPURationalPenalty)(x) = x / (1 + R.p * (1 - x))

struct GPUSinhPenalty{T} <: AbstractGPUPenalty{T}
    p::T
end
@inline (P::GPUSinhPenalty)(x) = CUDAnative.sinh(R.p * x) / CUDAnative.sinh(R.p)

for T in (:PowerPenalty, :RationalPenalty)
    fname = Symbol(:GPU, T)
    @eval @inline CuArrays.cu(p::$T) = $fname(p.p)
end
CuArrays.cu(p::AbstractGPUPenalty) = p

@define_cu(IterativeSolvers.CGStateVariables, :u, :r, :c)
whichdevice(ra::RaggedArray) = whichdevice(ra.offsets)
@define_cu(RaggedArray, :offsets, :values)


text:
"""
    interwine{T}(v1::Vector{T}, v2::Vector{T})

Interwine the elements of two vectors. The resulting
vector has double the length of the original ones.

**Example**
```
interwine([1,1,1], [2,2,2]) == [1,2,1,2,1,2]
```
"""
function interwine(v1::Vector{T}, v2::Vector{T}) where T
    @assert size(v1) == size(v2)
    v = Vector{T}(2length(v1))
    v[1:2:end] = v1
    v[2:2:end] = v2
    return v
end

"""
    interwine{T}(m1::Matrix{T}, m2::Matrix{T})

Interwine the columns of two matrices.
"""
function interwine(m1::Matrix{T}, m2::Matrix{T}) where T
    @assert size(m1) == size(m2)
    v = Vector{Vector{T}}(2size(m1,2))
    for j=1:size(m1,2)
        v[2j-1] = m1[:,j]
        v[2j] = m2[:,j]
    end
    return hcat(v...)
end

max_stars_repo_name: JuliaTagBot/LittleScienceTools

max_issues_repo_name: JuliaTagBot/LittleScienceTools

max_forks_repo_name: JuliaTagBot/LittleScienceTools


text:
using LinearAlgebra
using Surrogates
using Test
#1D
lb = 0.0
ub = 10.0
f = x -> log(x)*exp(x)
x = sample(5,lb,ub,SobolSample())
y = f.(x)
p = 1.9
my_k = Kriging(x,y,p)
add_point!(my_k,4.0,75.68)
add_point!(my_k,[5.0,6.0],[238.86,722.84])
pred = my_k(5.5)

#ND
lb = [0.0,0.0,1.0]
ub = [5.0,7.5,10.0]
x = sample(5,lb,ub,SobolSample())
f = x -> x[1]+x[2]*x[3]
y = f.(x)
theta = [2.0,2.0,2.0]
p = [1.9,1.9,1.9]
my_k = Kriging(x,y,p,theta)
add_point!(my_k,(4.0,3.2,9.5),34.4)
add_point!(my_k,[(1.0,4.65,6.4),(2.3,5.4,6.7)],[30.76,38.48])
pred = my_k((3.5,5.5,6.5))

#test sets
#WITHOUT ADD POINT
x = [(1.0,2.0,3.0),(4.0,5.0,6.0),(7.0,8.0,9.0)]
y = [1.0,2.0,3.0]
p = [1.0,1.0,1.0]
theta = [2.0,2.0,2.0]
my_k = Kriging(x,y,p,theta)
est = my_k((1.0,2.0,3.0))
std_err = std_error_at_point(my_k,(1.0,2.0,3.0))

#WITH ADD POINT adding singleton
x = [(1.0,2.0,3.0),(4.0,5.0,6.0),(7.0,8.0,9.0)]
y = [1.0,2.0,3.0]
p = [1.0,1.0,1.0]
theta = [2.0,2.0,2.0]
my_k = Kriging(x,y,p,theta)
add_point!(my_k,(10.0,11.0,12.0),4.0)
est = my_k((10.0,11.0,12.0))
std_err = std_error_at_point(my_k,(10.0,11.0,12.0))
@test std_err < 10^(-6)

#WITH ADD POINT ADDING MORE
x = [(1.0, 2.0, 3.0),(4.0,5.0,6.0),(7.0, 8.0, 9.0)]
y = [1.0,2.0,3.0]
p = [1.0,1.0,1.0]
theta = [2.0,2.0,2.0]
my_k = Kriging(x,y,p,theta)
add_point!(my_k,[(10.0, 11.0, 12.0),(13.0,14.0,15.0)],[4.0,5.0])
est = my_k((10.0,11.0,12.0))
std_err = std_error_at_point(my_k,(10.0,11.0,12.0))
@test std_err < 10^(-6)

#WITHOUT ADD POINT
x = [1.0,2.0,3.0]
y = [4.0,5.0,6.0]
p = 1.3
my_k = Kriging(x,y,p)
est = my_k(1.0)
std_err = std_error_at_point(my_k,1.0)
@test std_err < 10^(-6)

#WITH ADD POINT adding singleton
x = [1.0,2.0,3.0]
y = [4.0,5.0,6.0]
p = 1.3
my_k = Kriging(x,y,p)
add_point!(my_k,4.0,9.0)
est = my_k(4.0)
std_err = std_error_at_point(my_k,4.0)
@test std_err < 10^(-6)

#WITH ADD POINT adding more
x = [1.0,2.0,3.0]
y = [4.0,5.0,6.0]
p = 1.3
my_k = Kriging(x,y,p)
add_point!(my_k,[4.0,5.0,6.0],[9.0,13.0,15.0])
est = my_k(4.0)
std_err = std_error_at_point(my_k,4.0)
@test std_err < 10^(-6)


text:
"""
onesite infinite dmrg
"""
@with_kw struct Idmrg1{} <: Algorithm
    tol_galerkin::Float64 = Defaults.tol
    tol_gauge::Float64 = Defaults.tolgauge
    maxiter::Int = Defaults.maxiter
    verbose::Bool = Defaults.verbose
end

function find_groundstate(ost::InfiniteMPS, ham::Hamiltonian,alg::Idmrg1,oenvs=environments(ost,ham))
    st = copy(ost);
    envs = IDMRGEnv(ost,oenvs);

    delta::Float64 = 2*alg.tol_galerkin;

    for topit in 1:alg.maxiter
        delta = 0.0;

        curc = st.CR[0];

        for pos = 1:length(st)
            (eigvals,vecs) = @closure eigsolve(st.AC[pos],1,:SR,Arnoldi()) do x
                ac_prime(x,pos,st,envs)
            end

            st.AC[pos] = vecs[1]
            (st.AL[pos],st.CR[pos]) = leftorth(vecs[1]);

            setleftenv!(envs,pos+1,transfer_left(leftenv(envs,pos),ham[pos],st.AL[pos],st.AL[pos]));
        end

        for pos = length(st):-1:1

            (eigvals,vecs) = @closure eigsolve(st.AC[pos],1,:SR,Arnoldi()) do x
                ac_prime(x,pos,st,envs)
            end
            st.AC[pos] = vecs[1]
            (st.CR[pos-1],temp) = rightorth(_transpose_tail(vecs[1]));
            st.AR[pos] = _transpose_front(temp);

            setrightenv!(envs,pos-1,transfer_right(rightenv(envs,pos),ham[pos],st.AR[pos],st.AR[pos]));
        end

        delta = norm(curc-st.CR[0]);
        delta<alg.tol_galerkin && break;
        alg.verbose && @info "idmrg iter $(topit) err $(delta)"
    end

    nst = InfiniteMPS(st.AR[1:end],tol=alg.tol_gauge);
    nenvs = environments(nst, ham, solver=oenvs.solver)
    return nst,nenvs,delta;
end

"""
twosite infinite dmrg
"""
@with_kw struct Idmrg2{} <: Algorithm
    tol_galerkin::Float64 = Defaults.tol
    tol_gauge::Float64 = Defaults.tolgauge
    maxiter::Int = Defaults.maxiter
    verbose::Bool = Defaults.verbose
    trscheme = truncerr(1e-6);
end

function find_groundstate(ost::InfiniteMPS, ham::Hamiltonian,alg::Idmrg2,oenvs=environments(ost,ham))
    length(ost) < 2 && throw(ArgumentError("unit cell should be >= 2"))

    st = copy(ost);
    envs = IDMRGEnv(ost,oenvs);

    delta::Float64 = 2*alg.tol_galerkin;

    for topit in 1:alg.maxiter
        delta = 0.0;

        curc = st.CR[0];

        #sweep from left to right
        for pos = 1:length(st)-1
            ac2 = st.AC[pos]*_transpose_tail(st.AR[pos+1]);

            (eigvals,vecs) = @closure eigsolve(ac2,1,:SR,Arnoldi()) do x
                ac2_prime(x,pos,st,envs)
            end

            (al,c,ar,ϵ) = tsvd(vecs[1],trunc=alg.trscheme,alg=TensorKit.SVD())
            normalize!(c);

            st.AL[pos] = al
            st.CR[pos] = complex(c);
            st.AR[pos+1] = _transpose_front(ar);
            st.AC[pos+1] = _transpose_front(c*ar);

            setleftenv!(envs,pos+1,transfer_left(leftenv(envs,pos),ham[pos],st.AL[pos],st.AL[pos]));
            setrightenv!(envs,pos,transfer_right(rightenv(envs,pos+1),ham[pos+1],st.AR[pos+1],st.AR[pos+1]))
        end

        #update the edge
        @plansor ac2[-1 -2;-3 -4] := st.AC[end][-1 -2;1]*inv(st.CR[0])[1;2]*st.AL[1][2 -4;3]*st.CR[1][3;-3]
        (eigvals,vecs) = @closure eigsolve(ac2,1,:SR,Arnoldi()) do x
            ac2_prime(x,0,st,envs)
        end
        (al,c,ar,ϵ) = tsvd(vecs[1],trunc=alg.trscheme,alg=TensorKit.SVD())
        normalize!(c);

        st.AC[end] = al*c;
        st.AL[end] = al;
        st.CR[end] = complex(c);
        st.AR[1] = _transpose_front(ar);
        st.AC[1] = _transpose_front(c*ar);
        st.AL[1] = _transpose_front(c*ar)*inv(st.CR[1]);

        curc = complex(c);

        #update environments
        setleftenv!(envs,1,transfer_left(leftenv(envs,0),ham[0],st.AL[0],st.AL[0]));
        setrightenv!(envs,0,transfer_right(rightenv(envs,1),ham[1],st.AR[1],st.AR[1]));

        #sweep from right to left
        for pos = length(st)-1:-1:1
            ac2 = st.AL[pos]*_transpose_tail(st.AC[pos+1]);

            (eigvals,vecs) = @closure eigsolve(ac2,1,:SR,Arnoldi()) do x
                ac2_prime(x,pos,st,envs)
            end

            (al,c,ar,ϵ) = tsvd(vecs[1],trunc=alg.trscheme,alg=TensorKit.SVD())
            normalize!(c);

            st.AL[pos] = al
            st.AC[pos] = al*c
            st.CR[pos] = complex(c);
            st.AR[pos+1] = _transpose_front(ar);
            st.AC[pos+1] = _transpose_front(c*ar);

            setrightenv!(envs,pos,transfer_right(rightenv(envs,pos+1),ham[pos+1],st.AR[pos+1],st.AR[pos+1]))
            setleftenv!(envs,pos+1,transfer_left(leftenv(envs,pos),ham[pos],st.AL[pos],st.AL[pos]));
        end

        #update the edge
        @plansor ac2[-1 -2;-3 -4] :=  st.CR[end-1][-1;1]*st.AR[end][1 -2;2]*inv(st.CR[end])[2;3]*st.AC[1][3 -4;-3]
        (eigvals,vecs) = @closure eigsolve(ac2,1,:SR,Arnoldi()) do x
            ac2_prime(x,0,st,envs)
        end
        (al,c,ar,ϵ) = tsvd(vecs[1],trunc=alg.trscheme,alg=TensorKit.SVD())
        normalize!(c);

        st.AR[end] = _transpose_front(inv(st.CR[end-1])*_transpose_tail(al*c))
        st.AL[end] = al;
        st.CR[end] = complex(c);
        st.AR[1] = _transpose_front(ar);
        st.AC[1] = _transpose_front(c*ar);

        setleftenv!(envs,1,transfer_left(leftenv(envs,0),ham[0],st.AL[0],st.AL[0]));
        setrightenv!(envs,0,transfer_right(rightenv(envs,1),ham[1],st.AR[1],st.AR[1]));

        #update error
        smallest = infimum(_firstspace(curc),_firstspace(c));
        e1 = isometry(_firstspace(curc),smallest);
        e2 = isometry(_firstspace(c),smallest);
        delta = norm(e2'*c*e2-e1'*curc*e1)
        alg.verbose && @info "idmrg iter $(topit) err $(delta)"

        delta<alg.tol_galerkin && break;

    end

    nst = InfiniteMPS(st.AL[1:end],tol=alg.tol_gauge);
    nenvs = environments(nst, ham, solver=oenvs.solver)
    return nst,nenvs,delta;
end