text:
############################################
#
#    ufi Operator
#
# https://fenicsproject.org/docs/ufl/1.6.0/ufl.html
# < Operator < Expr
############################################

# no need to export
# not exported in `feincs` and `ufl`
@pyclass ufl ListTensor Expression
@pyclass ufl ComponentTensor Expression
@pyclass ufl Indexed Expression

# export
@pyfunc ufl variable

export variable

############################################
#
#    ufi : Tensor algebra operators
#
# https://fenicsproject.org/docs/ufl/1.6.0/ufl.html
# < CompoundTensorOperator < Operator < Expr
############################################

# no need to export
# not exported in `dolfin` or `ufl`
@pyclass ufl Outer Expression
@pyclass ufl Inner Expression
@pyclass ufl Conj Expression
@pyclass ufl Dot Expression
@pyclass ufl Cross Expression
@pyclass ufl Determinant Expression
@pyclass ufl Inverse Expression
@pyclass ufl Cofactor Expression
@pyclass ufl Transposed Expression
@pyclass ufl Trace Expression
@pyclass ufl Deviatoric Expression
@pyclass ufl Skew Expression
@pyclass ufl Sym Expression
@pyclass ufl Product Expression
@pyclass ufl Sum Expression
@pyclass ufl Division Expression
@pyclass ufl PositiveRestricted Expression
@pyclass ufl NegativeRestricted Expression

# export
@pyfunc ufl outer
@pyfunc ufl inner
@pyfunc ufl dot
@pyfunc ufl cross
@pyfunc ufl perp
@pyfunc dolfin det
inv(u::OpType) = Inverse(ufl.inv(u.pyobject)) # Base.inv
@pyfunc ufl cofac
transpose(u::OpType) = Transposed(ufl.transpose(u.pyobject)) # Base.transpose
@pyfunc ufl tr
@pyfunc ufl diag
@pyfunc ufl diag_vector
@pyfunc ufl dev
@pyfunc ufl skew
@pyfunc ufl sym
@pyfunc ufl avg
@pyfunc ufl jump

export outer, inner, dot, cross, perp, det, cofac, tr, diag, diag_vector, dev, skew, sym, avg, jump

############################################
# Differential operators
#
# https://fenicsproject.org/docs/ufl/1.6.0/ufl.html
# < CompoundDerivative < Derivative < Operator < Expr
############################################

# no need to export
# not exported in `dolfin` or `ufl`
@pyclass ufl VariableDerivative Expression
@pyclass ufl Grad Expression
@pyclass ufl Div Expression
@pyclass ufl NablaGrad Expression
@pyclass ufl NablaDiv Expression
@pyclass ufl Curl Expression

# export
diff(f::OpType,v::Variable) = VariableDerivative(ufl.diff(f.pyobject, v.pyobject))  # Base.diff
@pyfunc ufl grad
div(u::OpType) = Div(ufl.div(u.pyobject)) # Base.div
@pyfunc ufl nabla_grad
@pyfunc ufl nabla_div
@pyfunc ufl Dx
@pyfunc ufl Dn
@pyfunc ufl curl
@pyfunc ufl rot

export diff, grad, nabla_div, nabla_grad, Dx, Dn, curl, rot

############################################
# Nonlinear functions
#
# https://fenicsproject.org/docs/ufl/1.6.0/ufl.html
############################################

# no need to export
# not exported in `dolfin` or `ufl`
@pyclass ufl Ln Expression
@pyclass ufl Erf Expression

@pyfunc ufl max_value
@pyfunc ufl min_value
abs(u::OpType) = Expression(ufl.abs(u.pyobject) * u.pyobject ) # Base.abs
sign(u::OpType) = Expression(ufl.sign(u.pyobject)) # Base.sign
sqrt(u::OpType) = Expression(ufl.sqrt(u.pyobject)) # Base.sqrt
exp(u::OpType) = Expression(ufl.exp(u.pyobject)) # Base.exp
@pyfunc ufl ln
@pyfunc ufl erf
cos(u::OpType) = Expression(ufl.cos(u.pyobject)) # Base.cos
sin(u::OpType) = Expression(ufl.sin(u.pyobject)) # Base.sin
tan(u::OpType) = Expression(ufl.tan(u.pyobject)) # Base.tan
acos(u::OpType) = Expression(ufl.acos(u.pyobject)) # Base.acos
asin(u::OpType) = Expression(ufl.asin(u.pyobject)) # Base.asin
atan(u::OpType) = Expression(ufl.atan(u.pyobject)) # Base.atan
@pyfunc ufl atan_2
cosh(u::OpType) = Expression(ufl.cosh(u.pyobject)) # Base.cosh
sinh(u::OpType) = Expression(ufl.sinh(u.pyobject)) # Base.sinh
tanh(u::OpType) = Expression(ufl.tanh(u.pyobject)) # Base.tanh
#@pyfunc dolfin bessel_J
#@pyfunc dolfin bessel_Y
#@pyfunc dolfin bessel_I
#@pyfunc dolfin bessel_K

export max_value, min_value, ln, erf, atan_2, bessel_J, bessel_Y, bessel_I, bessel_K

############################################
# Form transformations
#
# https://fenicsproject.org/docs/ufl/1.6.0/ufl.html
############################################

@pyfunc ufl lhs
@pyfunc ufl rhs
@pyfunc ufl system
@pyfunc ufl functional
replace(u::Expression) = to_fetype(ufl.replace(u.pyobject))  # Base.replace
@pyfunc ufl replace_integral_domains
adjoint(u::Expression) = to_fetype(dolfin.adjoint(u.pyobject)) # Base.adjoint
@pyfunc ufl action
@pyfunc ufl energy_norm
@pyfunc ufl sensitivity_rhs
@pyfunc dolfin derivative

export lhs, rhs, system, functional, replace, replace_integral_domains, adjoint, action,
       energy_norm, sensitivity_rhs, derivative

############################################
#    Operator
############################################

#import Base: +, -, *, /, ^, ==, <<   # must be explicitly imported to be extended

*(expr1::OpType , expr2::Measure) = Expression(expr2.pyobject.__rmul__(expr1.pyobject) )

*(expr1::OpType , expr2::OpType ) = Expression(expr1.pyobject * expr2.pyobject)
*(expr1::Real, expr2::OpType) = Expression(expr1 * expr2.pyobject)
*(expr1::OpType , expr2::Real) = Expression(expr2 * expr1.pyobject)

+(expr1::OpType, expr2::Real) = Expression(expr1.pyobject + expr2)
+(expr1::Real, expr2::OpType) = Expression(expr2.pyobject + expr1)
+(expr1::OpType, expr2::OpType) = Expression(expr1.pyobject + expr2.pyobject)

-(expr::OpType) = (-1)*expr
-(expr1::OpType, expr2::Real) = Expression(expr1.pyobject - expr2)
-(expr1::Real, expr2::OpType) = (-1)*(Expression(expr2.pyobject - expr1) )
-(expr1::OpType, expr2::OpType) = Expression(expr1.pyobject - expr2.pyobject)

/(expr1::OpType, expr2::Real) = Expression(expr1.pyobject / expr2)
/(expr1::OpType, expr2::OpType) = Expression(expr1.pyobject / expr2.pyobject)
/(expr1::Real,expr2::OpType) = Constant(expr1) / expr2

^(expr1::OpType, expr2::Real) = Expression(expr1.pyobject.__pow__(expr2) )
^(expr1::OpType, expr2::OpType) = Expression(expr1.pyobject.__pow__(expr2.pyobject) )

==(expr1::OpType, expr2::OpType) = Expression(expr1.pyobject == expr2.pyobject)
==(expr1::OpType, expr2::Real) = Expression(expr1.pyobject == expr2)
==(expr1::Real, expr2::OpType) = Expression(expr2.pyobject == expr1)

# Base: split
function split(fun::OpType)
     vec = dolfin.split(fun.pyobject)
     expr_vec = [to_fetype(spl) for spl in vec]
     return expr_vec
end


text:
@testset "1488.avoid-flood-in-the-city.jl" begin
    @test avoid_flood([1, 2, 3, 4]) == [-1, -1, -1, -1]
    @test avoid_flood([1, 2, 0, 0, 2, 1]) == [-1, -1, 2, 1, -1, -1]
    @test avoid_flood([1, 2, 0, 1, 2]) == []
    @test avoid_flood([69, 0, 0, 0, 69]) == [-1, 69, 1, 1, -1]
    @test avoid_flood([10, 20, 20]) == []
    @test avoid_flood([0, 1, 1]) == []
    @test avoid_flood([1, 0, 2, 0, 2, 1]) == [-1, 1, -1, 2, -1, -1]
    @test avoid_flood([1, 0, 2, 3, 0, 1, 2]) == [-1, 1, -1, -1, 2, -1, -1]
    @test avoid_flood([1, 0, 2, 0, 3, 0, 2, 0, 0, 0, 1, 2, 3]) ==
          [-1, 1, -1, 2, -1, 3, -1, 2, 1, 1, -1, -1, -1]
end

ext: jl


text:
using Base.Test
using OptimTools

if !isdefined(:probs) || isempty(probs)
	println("getting test problems")
	include("getTestFunctions.jl")
end

his = []; flag = []; x0 = [];

secondOrderMethods = (dampedNewton,newton,modnewton,newtoncg)
# try stopping based on atol
for method in secondOrderMethods
	println("testing $(string(method))")
	for k=1:length(probs)
		prob = probs[k]

		println("\t\t$(prob[1])")

		fk   = prob[2]
		dfk  = prob[3]
		d2fk = prob[4]
		xsol = prob[5]
		x0   = float(xsol) + 1*randn(size(xsol))*norm(xsol)
		atol = 1e-10

		xnt,flag,his = method(fk,dfk,d2fk, x0,maxIter=200,atol=atol,out=1*(k==2))
		@test (flag==-2) | (norm(dfk(xnt)) .< atol)
		if flag==-2
			@test minimum(eig(d2fk(xnt))[1]) < 0
		end
	end
end


text:
include("utils.jl")
#just a holder for the z partitions.
#to allow split_model to work correctly

abstract type SAFTgammaMieModel <: GCSAFTModel end

struct γMieZ <: EoSModel
    components::Vector{String}
    z::SingleParam{Vector{Float64}}
end

struct SAFTgammaMieParam <: EoSParam
    segment::SingleParam{Int}
    mixedsegment::SingleParam{Vector{Float64}}
    shapefactor::SingleParam{Float64}
    lambda_a::PairParam{Float64}
    lambda_r::PairParam{Float64}
    sigma::PairParam{Float64}
    epsilon::PairParam{Float64}
    epsilon_assoc::AssocParam{Float64}
    bondvol::AssocParam{Float64}
end

@registermodel γMieZ

struct SAFTgammaMie{I,VR} <: SAFTgammaMieModel
    components::Vector{String}
    groups::GroupParam
    sites::SiteParam
    params::SAFTgammaMieParam
    idealmodel::I
    vrmodel::VR
    mie_zfractions::γMieZ
    assoc_options::AssocOptions
    references::Array{String,1}
end

function split_model(model::SAFTgammaMie,subset = nothing)
    if !(subset === nothing)
        error("SAFTgammaMie does not support custom subsets")
    end
    pures = auto_split_model(model)
    for (i,pure) in pairs(pures)
        filter!(!iszero,pure.mie_zfractions.z.values[1])
        mixm = pure.params.mixedsegment.values
        for mi in mixm
            xi = mi[i]
            resize!(mi, 1)
            mi[1] = xi
        end
    end
    return pures
end

function SAFTgammaMie(components;
    idealmodel=BasicIdeal,
    userlocations=String[],
    ideal_userlocations=String[],
    verbose=false,
    assoc_options = AssocOptions())

    groups = GroupParam(components, ["SAFT/SAFTgammaMie/SAFTgammaMie_groups.csv"]; verbose=verbose)
    params,sites = getparams(groups, ["SAFT/SAFTgammaMie","properties/molarmass_groups.csv"]; userlocations=userlocations, verbose=verbose)
    _sites = sites.i_sites
    components = groups.components
    gc = groups.i_flattenedgroups
    comps = 1:length(components)
    gc_segment = params["vst"]
    shapefactor = params["S"]
    S = shapefactor.values
    vst = gc_segment.values
    v  = groups.n_flattenedgroups
    comp_segment = zeros(Float64,length(comps))
    for i ∈ comps
        res_i = 0.0
        vi = v[i]
        groups_i = groups.i_groups[i]
        for idx ∈ 1:length(groups_i)
            k = groups_i[idx]
            res_i += vi[k]*S[k]*vst[k]
        end
        comp_segment[i] = res_i
    end
    segment = SingleParam("segment",components,comp_segment)

    #used in x_S:
    #x_S(group i) = dot(z,mixsegment[i])/dot(z,m_vr)
    mixsegment =  [[v[i][k]*vst[k]*S[k] for i ∈ comps] for k ∈ gc]
    gc_mixedsegment = SingleParam("mixed segment",groups.flattenedgroups,mixsegment,[false for i ∈ gc],String[],String[])
    gc_sigma = sigma_LorentzBerthelot(params["sigma"])
    gc_sigma.values .*= 1E-10
    gc_epsilon = epsilon_HudsenMcCoubrey(params["epsilon"], gc_sigma)
    ϵ = gc_epsilon.values
    σ = gc_sigma.values
     #helper functions
     function ẑ(i, k)
        return v[i][k]*vst[k]*S[k] / ∑(v[i][l]*vst[l]*S[l] for l ∈ gc)
    end

    zz = [[0.0 for k ∈ 1:length(v[i])] for i ∈ comps]
    zzparam = SingleParam("z fraction",components,zz,[false for i ∈ 1:length(components)],String[],String[])
    for i ∈ 1:length(zz)
        zzi = zz[i]
        for k ∈ 1:length(zzi)
            zzi[k] = ẑ(i, k)
        end
    end

    function σ̄(i)
        return cbrt(∑(∑(ẑ(i,k)*ẑ(i,l)*σ[k,l]^3 for l ∈ gc) for k ∈ gc))
    end

    function σ̄(i, j)
        return (σ̄(i) + σ̄(j))/2
    end

    function ϵ̄(i)
        return ∑(∑(ẑ(i,k)*ẑ(i,l)*ϵ[k,l] for l ∈ gc) for k ∈ gc)
    end

    function ϵ̄(i, j)
        if i == j
            return ϵ̄(i)
        else
            return sqrt(σ̄(i)*σ̄(j))/σ̄(i,j) * sqrt(ϵ̄(i)*ϵ̄(j))
        end
    end

    comp_ϵ = [ϵ̄(i, j) for (i,j) ∈ Iterators.product(comps,comps)]
    epsilon = PairParam("epsilon",components,comp_ϵ)

    comp_σ = [σ̄(i, j) for (i,j) ∈ Iterators.product(comps,comps)]
    sigma = PairParam("sigma",components,comp_σ)

    gc_lambda_a = lambda_LorentzBerthelot(params["lambda_a"])
    gc_lambda_r = lambda_LorentzBerthelot(params["lambda_r"])
    #@show gc_lambda_a
    gc_λa = gc_lambda_a.values
    gc_λr = gc_lambda_r.values

    function λi(ll,i)
        return ∑(∑(ẑ(i,k)*ẑ(i,l)*ll[k,l] for l ∈ gc) for k ∈ gc)
    end
    comp_lambda_a = [λi(gc_λa,i) for i ∈ comps]
    comp_lambda_r = [λi(gc_λr,i) for i ∈ comps]
    lambda_a = lambda_LorentzBerthelot(SingleParam("lambda_a",components,comp_lambda_a))
    lambda_r = lambda_LorentzBerthelot(SingleParam("lambda_r",components,comp_lambda_r))

    #assoc
    gc_epsilon_assoc = params["epsilon_assoc"]
    gc_bondvol = params["bondvol"]
    comp_sites,idx_dict = gc_to_comp_sites(sites,groups)
    assoc_idx = gc_to_comp_assoc_idx(gc_bondvol,comp_sites,idx_dict)
    assoc_idxs,outer,inner,outer_size,inner_size = assoc_idx.values,assoc_idx.outer_indices,assoc_idx.inner_indices,assoc_idx.outer_size,assoc_idx.inner_size
    _comp_bondvol = [gc_bondvol.values.values[i] for i ∈ assoc_idxs]
    _comp_epsilon_assoc = [gc_epsilon_assoc.values.values[i] for i ∈ assoc_idxs]
    compval_bondvol = CompressedAssocMatrix(_comp_bondvol,outer,inner,outer_size,inner_size)
    compval_epsilon_assoc = CompressedAssocMatrix(_comp_epsilon_assoc,outer,inner,outer_size,inner_size)
    comp_bondvol = AssocParam{Float64}("epsilon assoc",components,compval_bondvol,comp_sites.sites,String[],String[])
    comp_epsilon_assoc = AssocParam{Float64}("bondvol",components,compval_epsilon_assoc,comp_sites.sites,String[],String[])
    #mixing of Mw
    comp_mw = zeros(Float64,length(components))
    gc_mw = params["Mw"].values
    for i ∈ 1:length(components)
        vi = v[i]
        gi = groups.i_groups[i]
        mwi = 0.0
        for idx ∈ 1:length(gi)
            j = gi[idx]
            mwi += gc_mw[j]*vi[j]
        end
        comp_mw[i] =mwi
    end
    _mw = SingleParam("molecular_weight",components,comp_mw)
    gcparams = SAFTgammaMieParam(gc_segment,gc_mixedsegment, shapefactor,gc_lambda_a,gc_lambda_r,gc_sigma,gc_epsilon,gc_epsilon_assoc,gc_bondvol)
    vrparams = SAFTVRMieParam(segment,sigma,lambda_a,lambda_r,epsilon,comp_epsilon_assoc,comp_bondvol,_mw)
    idmodel = init_model(idealmodel,components,ideal_userlocations,verbose)
    vr = SAFTVRMie(vrparams, comp_sites, idmodel; ideal_userlocations, verbose, assoc_options)
    mie_z = γMieZ(components,zzparam)
    γmierefs = ["10.1063/1.4851455", "10.1021/je500248h"]
    gmie = SAFTgammaMie(components,groups,sites,gcparams,idmodel,vr,mie_z,assoc_options,γmierefs)
    return gmie
end
@registermodel SAFTgammaMie

const SAFTγMie = SAFTgammaMie
export SAFTgammaMie,SAFTγMie

function a_res(model::SAFTgammaMieModel, V, T, z)
    _data = @f(data)
    return @f(a_hs,_data) + @f(a_disp,_data) + @f(a_chain,_data) + @f(a_assoc,_data)
end

function a_mono(model::SAFTgammaMieModel, V, T, z,_data = @f(data))
    return @f(a_hs,_data) + @f(a_disp,_data)
end

function data(model::SAFTgammaMieModel, V, T, z)
    _d_gc = @f(d_gc)
    _d = @f(d,_d_gc)
    ζi = @f(ζ0123,_d_gc)
    _ζ_X,σ3x = @f(ζ_X_σ3,_d_gc)
    _ρ_S = @f(ρ_S)
    _ζst = _ζst = σ3x*_ρ_S*π/6
    return (_d_gc,(_d,_ρ_S,ζi,_ζ_X,_ζst,σ3x))
end

function d_gc(model::SAFTgammaMieModel, V, T, z)
    _ϵ = model.params.epsilon.diagvalues
    _σ = model.params.sigma.diagvalues
    _λa = model.params.lambda_a.diagvalues
    _λr = model.params.lambda_r.diagvalues
    u = SAFTγMieconsts.u
    w = SAFTγMieconsts.w
    _0 = zero(T+first(z))
    n = length(_ϵ)
    _d_gc = zeros(typeof(_0),n)
    for k ∈ 1:n
        ϵ = _ϵ[k]
        σ = _σ[k]
        λa = _λa[k]
        λr = _λr[k]
        θ = Cλ(model.vrmodel,V,T,z,λa,λr)*ϵ/T
        di = _0
        λrinv = 1/λr
        λaλr = λa/λr
        for j ∈ 1:5
            di += w[j]*(θ/(θ+u[j]))^λrinv*(exp(θ*(1/(θ/(θ+u[j]))^λaλr-1))/(u[j]+θ)/λr)
        end
        _d_gc[k] = σ*(1-di)
    end
    return _d_gc
end

function d(model::SAFTgammaMieModel,V,T,z,_d_gc = @f(d_gc))
    _d = zeros(eltype(_d_gc),length(z))
    _0 = zero(eltype(_d))
    _zz = model.mie_zfractions.z.values
    for i ∈ @comps
        _z = _zz[i]
        m = length(_z)
        di = _0
        for k ∈ 1:m
            zk = _z[k]
            iszero(zk) && continue
            dk = _d_gc[k]
            di += zk*zk*dk^3
            for l ∈ 1:k - 1
                di += 0.25*zk*_z[l]*(dk + _d_gc[l])^3
            end
        end
        _d[i] = cbrt(di)
    end
    #d̂[i] = cbrt(∑(∑(@f(ẑ,i,k)*@f(ẑ,i,l)*@f(d,k,l)^3 for l ∈ @groups) for k ∈ @groups))
    return _d
end

function ζ0123(model::SAFTgammaMieModel, V, T, z,_d_gc=@f(d_gc))
    vrmodel = model.vrmodel
    m = vrmodel.params.segment.values
    m̄ = dot(z, m)
    m̄inv = 1/m̄
    _0 = zero(V+T+first(z))
    ζ0,ζ1,ζ2,ζ3 = _0,_0,_0,_0
    _ϵ = model.params.epsilon.diagvalues
    σ = model.params.sigma.diagvalues
    _λr = model.params.lambda_r.diagvalues
    _λa = model.params.lambda_a.diagvalues
    u = SAFTVRMieconsts.u
    _mi = model.params.mixedsegment.values
    w = SAFTVRMieconsts.w
    for i ∈ @groups
        λa = _λa[i]
        λr = _λr[i]
        ϵ = _ϵ[i]
        di =_d_gc[i]
        mi = _mi[i]
        xS = dot(mi,z)*m̄inv
        ζ0 += xS
        ζ1 += xS*di
        ζ2 += xS*di*di
        ζ3 += xS*di*di*di
    end
    c = π/6*N_A*m̄/V
    ζ0,ζ1,ζ2,ζ3 = c*ζ0,c*ζ1,c*ζ2,c*ζ3
    return ζ0,ζ1,ζ2,ζ3
end

#return π/6*@f(ρ_S)*∑(@f(x_S,k)*@f(d,k)^m for k ∈ @groups)

function ρ_S(model::SAFTgammaMieModel, V, T, z)
    return ρ_S(model.vrmodel,V,T,z)
end

function ζ_X_σ3(model::SAFTgammaMieModel, V, T, z,_d = @f(d_gc))
    vrmodel = model.vrmodel
    m = vrmodel.params.segment.values
    m̄ = dot(z, m)
    m̄inv = 1/m̄
    σ = model.params.sigma.values
    _mi = model.params.mixedsegment.values
    ρS = N_A/V*m̄
    _ζ_X = zero(first(_d))
    kρS = ρS* π/6
    σ3_x = _ζ_X
    r1 =_ζ_X
    for i ∈ @groups
        mi = _mi[i]
        x_Si = dot(mi,z)*m̄inv
        σ3_x += x_Si*x_Si*(σ[i,i]^3)
        _ζ_X += x_Si*x_Si*(_d[i])^3
        for j ∈ 1:i-1
            mj = _mi[j]
            x_Sj = dot(mj,z)*m̄inv
            σ3_x += 2*x_Si*x_Sj*(σ[i,j]^3)
            dij = 0.5*(_d[i] + _d[j])
            r1 = 2*x_Si*x_Sj*dij^3
            _ζ_X += r1
        end
    end

    #return π/6*@f(ρ_S)*∑(@f(x_S,i)*@f(x_S,j)*(@f(d,i)+@f(d,j))^3/8 for i ∈ comps for j ∈ comps)
    return kρS*_ζ_X,σ3_x
end

function σ3x(model::SAFTgammaMieModel, V, T, z)
    vrmodel = model.vrmodel
    m = vrmodel.params.segment.values
    m̄ = dot(z, m)
    m̄inv = 1/m̄
    σ = model.params.sigma.values
    _mi = model.params.mixedsegment.values
    σ3_x =  zero(first(z))
    for i ∈ @groups
        mi = _mi[i]
        x_Si = dot(mi,z)*m̄inv
        σ3_x += x_Si*x_Si*(σ[i,i]^3)
        for j ∈ 1:i-1
            mj = _mi[j]
            x_Sj = dot(mj,z)*m̄inv
            σ3_x += 2*x_Si*x_Sj*(σ[i,j]^3)
        end
    end
    return σ3_x
end

function  a_hs(model::SAFTgammaMieModel, V, T, z,_data = @f(data))
    _,vrdata = _data
    return a_hs(model.vrmodel,V,T,z,vrdata)
end

function a_chain(model::SAFTgammaMieModel, V, T, z,_data = @f(data))
    _,vrdata = _data
    return a_chain(model.vrmodel,V,T,z,vrdata)
end

function a_assoc(model::SAFTgammaMieModel, V, T, z,_data = @f(data))
    _,vrdata = _data
    return a_assoc(model.vrmodel,V,T,z,vrdata)
end

function a_disp(model::SAFTgammaMieModel, V, T, z,_data = @f(data))
    groups = @groups
    gcmodel = model
    model = gcmodel.vrmodel
    _d,vrdata = _data
    _,ρS,ζi,_ζ_X,_ζst,_ = vrdata
    comps = @comps
    l = length(comps)
    ∑z = ∑(z)
    m = model.params.segment.values
    _mi = gcmodel.params.mixedsegment.values
    _ϵ = gcmodel.params.epsilon.values
    _λr = gcmodel.params.lambda_r.values
    _λa = gcmodel.params.lambda_a.values
    _σ = gcmodel.params.sigma.values
    m̄ = dot(z, m)
    m̄inv = 1/m̄
    a₁ = zero(V+T+first(z))
    a₂ = a₁
    a₃ = a₁
    _ζst5 = _ζst^5
    _ζst8 = _ζst^8
    _KHS = @f(KHS,_ζ_X,ρS)
    for i ∈ groups
        j = i
        mi = _mi[i]
        x_Si = dot(mi,z)*m̄inv
        ϵ = _ϵ[i,j]
        λa = _λa[i,j]
        λr = _λr[i,j]
        σ = _σ[i,j]
        _C = @f(Cλ,λa,λr)
        dij = _d[i]
        x_0ij = σ/dij
        dij3 = dij^3
        x_0ij = σ/dij
        #calculations for a1 - diagonal
        aS_1_a = @f(aS_1,λa,_ζ_X)
        aS_1_r = @f(aS_1,λr,_ζ_X)
        B_a = @f(B,λa,x_0ij,_ζ_X)
        B_r = @f(B,λr,x_0ij,_ζ_X)
        a1_ij = (2*π*ϵ*dij3)*_C*ρS*
        (x_0ij^λa*(aS_1_a+B_a) - x_0ij^λr*(aS_1_r+B_r))

        #calculations for a2 - diagonal
        aS_1_2a = @f(aS_1,2*λa,_ζ_X)
        aS_1_2r = @f(aS_1,2*λr,_ζ_X)
        aS_1_ar = @f(aS_1,λa+λr,_ζ_X)
        B_2a = @f(B,2*λa,x_0ij,_ζ_X)
        B_2r = @f(B,2*λr,x_0ij,_ζ_X)
        B_ar = @f(B,λr+λa,x_0ij,_ζ_X)
        α = _C*(1/(λa-3)-1/(λr-3))
        f1,f2,f3,f4,f5,f6 = @f(f123456,α)
        _χ = f1*_ζst+f2*_ζst5+f3*_ζst8
        a2_ij = π*_KHS*(1+_χ)*ρS*ϵ^2*dij3*_C^2 *
        (x_0ij^(2*λa)*(aS_1_2a+B_2a)
        - 2*x_0ij^(λa+λr)*(aS_1_ar+B_ar)
        + x_0ij^(2*λr)*(aS_1_2r+B_2r))

        #calculations for a3 - diagonal
        a3_ij = -ϵ^3*f4*_ζst * exp(f5*_ζst+f6*_ζst^2)
        #adding - diagonal
        a₁ += a1_ij*x_Si*x_Si
        a₂ += a2_ij*x_Si*x_Si
        a₃ += a3_ij*x_Si*x_Si
        for j ∈ 1:i-1
            mj = _mi[j]
            x_Sj = dot(mj,z)*m̄inv
            ϵ = _ϵ[i,j]
            λa = _λa[i,j]
            λr = _λr[i,j]
            σ = _σ[i,j]
            _C = @f(Cλ,λa,λr)
            dij = 0.5*(_d[i]+_d[j])
            x_0ij = σ/dij
            dij3 = dij^3
            x_0ij = σ/dij
            #calculations for a1
            a1_ij = (2*π*ϵ*dij3)*_C*ρS*
            (x_0ij^λa*(@f(aS_1,λa,_ζ_X)+@f(B,λa,x_0ij,_ζ_X)) - x_0ij^λr*(@f(aS_1,λr,_ζ_X)+@f(B,λr,x_0ij,_ζ_X)))

            #calculations for a2
            α = _C*(1/(λa-3)-1/(λr-3))
            f1,f2,f3,f4,f5,f6 = @f(f123456,α)
            _χ = f1*_ζst+f2*_ζst5+f3*_ζst8
            a2_ij = π*_KHS*(1+_χ)*ρS*ϵ^2*dij3*_C^2 *
            (x_0ij^(2*λa)*(@f(aS_1,2*λa,_ζ_X)+@f(B,2*λa,x_0ij,_ζ_X))
            - 2*x_0ij^(λa+λr)*(@f(aS_1,λa+λr,_ζ_X)+@f(B,λa+λr,x_0ij,_ζ_X))
            + x_0ij^(2*λr)*(@f(aS_1,2λr,_ζ_X)+@f(B,2*λr,x_0ij,_ζ_X)))

            #calculations for a3
            a3_ij = -ϵ^3*f4*_ζst * exp(f5*_ζst+f6*_ζst^2)
            #adding
            a₁ += 2*a1_ij*x_Si*x_Sj
            a₂ += 2*a2_ij*x_Si*x_Sj
            a₃ += 2*a3_ij*x_Si*x_Sj
        end
    end
    a₁ = a₁*m̄/T/∑z
    a₂ = a₂*m̄/(T*T)/∑z
    a₃ = a₃*m̄/(T*T*T)/∑z
    #@show (a₁,a₂,a₃)
    adisp =  a₁ + a₂ + a₃
    return adisp
end

const SAFTγMieconsts =(
    A = SA[ 0.81096    1.7888   -37.578   92.284;
          1.02050  -19.341    151.26  -463.50 ;
         -1.90570   22.845   -228.14   973.92 ;
          1.08850   -6.1962   106.98  -677.64  ],

    ϕ = [[7.5365557, -359.440,  1550.9, -1.199320, -1911.2800,  9236.9  ],
        [-37.604630,  1825.60, -5070.1,  9.063632,  21390.175, -129430 ],
        [71.745953, -3168.00,  6534.6, -17.94820, -51320.700,  357230  ],
        [-46.835520,  1884.20, -3288.7,  11.34027,  37064.540, -315530 ],
        [-2.4679820,- 0.82376, -2.7171,  20.52142,  1103.7420,  1390.2 ],
        [-0.5027200, -3.19350,  2.0883, -56.63770, -3264.6100, -4518.2 ],
        [8.0956883,  3.70900,  0.0000,  40.53683,  2556.1810,  4241.6  ]],

    u = [0.26356031971814109102031,1.41340305910651679221800,3.59642577104072208122300,7.08581000585883755692200,12.6408008442757826594300],
    w = [0.5217556105828086524759,0.3986668110831759274500,7.5942449681707595390e-2,3.6117586799220484545e-3,2.3369972385776227891e-5],

         c  = [0.0756425183020431	-0.128667137050961	 0.128350632316055	-0.0725321780970292	   0.0257782547511452  -0.00601170055221687	  0.000933363147191978  -9.55607377143667e-05  6.19576039900837e-06 -2.30466608213628e-07 3.74605718435540e-09
                0.134228218276565	    -0.182682168504886 	 0.0771662412959262	-0.000717458641164565 -0.00872427344283170	0.00297971836051287	 -0.000484863997651451	 4.35262491516424e-05 -2.07789181640066e-06	4.13749349344802e-08 0
               -0.565116428942893	     1.00930692226792   -0.660166945915607	 0.214492212294301	  -0.0388462990166792	0.00406016982985030	 -0.000239515566373142	 7.25488368831468e-06 -8.58904640281928e-08	0	                 0
               -0.387336382687019	    -0.211614570109503	 0.450442894490509	-0.176931752538907	   0.0317171522104923  -0.00291368915845693	  0.000130193710011706  -2.14505500786531e-06  0	                0	                 0
                2.13713180911797	    -2.02798460133021 	 0.336709255682693	 0.00118106507393722  -0.00600058423301506	0.000626343952584415 -2.03636395699819e-05	 0	                   0	                0	                 0
               -0.300527494795524	     2.89920714512243   -0.567134839686498	 0.0518085125423494	  -0.00239326776760414	4.15107362643844e-05  0	                     0	                   0	                0                    0
               -6.21028065719194	    -1.92883360342573	 0.284109761066570	-0.0157606767372364	   0.000368599073256615	0 	                  0	                     0	                   0	                0	                 0
                11.6083532818029	     0.742215544511197  -0.0823976531246117	 0.00186167650098254   0	                0	                  0	                     0	                   0	                0	                 0
               -10.2632535542427	    -0.125035689035085	 0.0114299144831867	 0	                   0	                0	                  0	                     0	                   0	                0	                 0
                4.65297446837297	    -0.00192518067137033 0	                 0	                   0	                0	                  0	                     0	                   0	                0	                 0
               -0.867296219639940	     0	                 0	                 0	                   0	                0	                  0	                     0	                   0	                0	                 0],
)


text:
using Distributed
addprocs(2)

using Test
using MLJTuning
using MLJBase

include("test_utilities.jl")

print("Loading some models for testing...")
# load `Models` module containing models implementations for testing:
@everywhere include("models.jl")
print("\r                                           \r")

@testset "utilities" begin
    @test include("utilities.jl")
end

@testset "tuned_models.jl" begin
    @test include("tuned_models.jl")
end

@testset "ranges" begin
    @test include("ranges.jl")
end

@testset "grid" begin
    @test include("strategies/grid.jl")
end

@testset "learning curves" begin
        @test include("learning_curves.jl")
end

# @testset "julia bug" begin
#     @test include("julia_bug.jl")
# end


text:
function save_sparse_matrix(fid, A::SparseMatrixCSC, name::String)

    n, m = size(A)
    Is, Js, Vs = findnz(A)

    fid["$(name)_Is"] = Int32[Is...]
    fid["$(name)_Js"] = Int32[Js...]
    fid["$(name)_Vs"] = Float64[Vs...]
    fid["$(name)_size"] = Int32[n, m]

    @info "Sparse Matrix '$(name)' saved to disk."
end

function load_sparse_matrix(fid, name::String)

    if isempty(fid["$(name)_Is"])
        Is = Int32[1]
        Js = Int32[1]
        Vs = Float64[0]
    else
        Is = fid["$(name)_Is"][:]
        Js = fid["$(name)_Js"][:]
        Vs = fid["$(name)_Vs"][:]
    end
    n, m = fid["$(name)_size"][:]

    @info "Loaded Sparse Matrix '$(name)'."

    return sparse(Is, Js, Vs, n, m)
end

function save_grid(fid, grd::Grid, name::String)

    # create group for grid
    create_group(fid, name)

    # save projection
    save_sparse_matrix(fid, grd.Proj, "$(name)/Proj")
    fid["$(name)/r0"] = grd.r0[:]
    fid["$(name)/r1"] = grd.r1[:]
    fid["$(name)/points"] = grd.points[:,:]
    fid["$(name)/spacing"] = Float64[grd.dx, grd.dy]
    fid["$(name)/_size"] = Int64[grd.Nk, grd.Nz, grd._Nx, grd._Ny, grd._Nbuffer]
    fid["$(name)/_NaNs"] = grd._nan_outside_boundaries[:,:]

    @info "Grid '$(name)' saved to disk."

end

function load_grid(fid, name::String)

    Proj = load_sparse_matrix(fid, "$(name)/Proj")
    r0 = fid["$(name)/r0"][:]
    r1 = fid["$(name)/r1"][:]
    points = fid["$(name)/points"][:,:]
    dx, dy = fid["$(name)/spacing"][:]
    Nk, Nz, Nx, Ny, Nbuffer = fid["$(name)/_size"][:]
    NaNs = fid["$(name)/_NaNs"][:,:]

    @info "Loaded Grid '$(name)'."

    return Grid(r0, r1, points, Proj, dx, dy, Nk, Nz, Nx, Ny, Nbuffer, NaNs)
end

function save_ghost_conditions(fid, gc::GhostConditions, name::String)

    # create group
    create_group(fid, name)

    save_sparse_matrix(fid, gc.Proj, "$(name)/Proj")
    save_sparse_matrix(fid, gc.Mirror, "$(name)/Mirror")
    fid["$(name)/swaps"] = gc.swaps[:,:]
end

function load_ghost_conditions(fid, name::String)

    Proj = load_sparse_matrix(fid, "$(name)/Proj")
    Mirror = load_sparse_matrix(fid, "$(name)/Mirror")
    swaps = fid["$(name)/swaps"][:,:]

    return GhostConditions(Proj, Mirror, swaps)
end

function example_geometry_setup(path::String, Nx, Ny)

    # define the magnetic field
    psi, bx, by, bz = solovev_flux_function(0.1, 0.1, 0.3, 1.7, 10.0, downsep=[1, -0.561])

    # values for flux surfaces
    psi_in = psi(1 - 0.28, 0)
    psi_out = psi(1 - 0.34, 0)
    psi_priv = psi(1, -(0.561 + 0.01))

    # define barriers for flux surfaces
    inner_flux = Barrier() do
        func(x,y) = (y > -0.561) ? psi(x,y) : abs(psi(x,y))
        rmap(x,y) = trace_reflection(x, y, psi, psi_in, 0.001)
        func, psi_in, rmap, 1
    end

    outer_flux = Barrier() do
        func(x,y) = (norm([x-1,y+0.6]) > 0.05) ? psi(x,y) : psi_out + psi_priv - psi(x,y)
        rmap(x,y) = trace_reflection(x, y, func, psi_out, 0.001)
        func, psi_out, rmap, -1
    end

    # define barrier for divertor target plates
    target = Barrier() do
        func(x,y) = y
        val = -(0.561 + 0.04)
        rmap(x,y) = x, 2*val - y
        func, val, rmap, 1
    end

    # compute step function widths
    dr = 0.01
    dp_in = abs(ForwardDiff.derivative(u -> psi(u,0), 0.7)) * dr
    dp_out = abs(ForwardDiff.derivative(u -> psi(u,0), 0.6)) * dr
    dp_priv = abs(ForwardDiff.derivative(u -> psi(1,u), -(0.561 + 0.05))) * dr

    # create Grid
    deltas = Float64[dp_in, 0.5*(dp_out + dp_priv), dr]
    barrs = [inner_flux, outer_flux, target]
    grd = Grid(Nx, Ny, 1) do
        inside(x,y) = check_if_inside(x, y, 10*deltas, barrs)
        r0 = Float64[0.5, -1]
        r1 = Float64[1.5, 1]
        inside, r0, r1
    end

    @info grd.Nk

    # compute penalization values as step functions
    pen = zeros(grd.Nk, 3)
    for i=1:3
        pen[:,i] = f_to_grid(grd) do x,y
            smoothstep(x, y, deltas[i], barrs[i])
        end
    end

    # K1 and K2 are operators for penalized forward-euler timesteps
    K1 = sparse(Diagonal(pen[:,1] .* pen[:,2] .* pen[:,3]))
    K2 = 100 * (I - K1)

    # H1, H2, and H3 are for defining boundary values.
    h1 = ones(grd.Nk)
    h2 = ones(grd.Nk)
    h3 = ones(grd.Nk)
    h1[pen[:,1] .!= 0.0] .= 0.0
    h2[pen[:,2] .!= 0.0] .= 0.0
    h3[pen[:,3] .!= 0.0] .= 0.0
    H1 = sparse(Diagonal(h1))
    H2 = sparse(Diagonal(h2))
    H3 = sparse(Diagonal(h3))

    # compute ghost-conditions
    GC_nmann = GhostConditions(barrs, grd) #XXX this is backwards now I think.
    GC_dchlt = swap_sign(GC_nmann, 1, 2, 3)
    GC_u = swap_sign(GC_nmann, 3)

    # define vector for function that is +1 on left plate and -1 on right plate
    trgt_sgn = f_to_grid((x,y) -> sign(1 - x), grd)

    # lcfs average
    lcfs_avg = flux_surface_average(psi, 0.0, 150, Float64[1, 0], 4, 4, stencil2d(4, 4), grd)

    # trace field-line images
    img = fieldline_images(bx, by, bz, 0.01, grd)

    fid = h5open(path, "w")
    save_grid(fid, grd, "Grid")
    save_sparse_matrix(fid, K1, "K1")
    save_sparse_matrix(fid, K2, "K2")
    save_sparse_matrix(fid, H1, "H1")
    save_sparse_matrix(fid, H2, "H2")
    save_sparse_matrix(fid, H3, "H3")
    save_ghost_conditions(fid, GC_dchlt, "GC_dchlt")
    save_ghost_conditions(fid, GC_nmann, "GC_nmann")
    save_ghost_conditions(fid, GC_u, "GC_u")
    fid["trgt_sgn"] = trgt_sgn[:]
    save_sparse_matrix(fid, lcfs_avg, "lcfs_avg")
    fid["FL_img"] = img[:,:]
    close(fid)
end

function example_physics_setup(R0, n0, T0, B0, init_path, geo_path)

    am, ad, ki, ke, er, eg, ev, de2, eta = dimensionless_parameters(0.3*R0, R0, n0, T0, B0)
    psi, bx, by, bz = solovev_flux_function(0.1, 0.1, 0.3, 1.7, 10.0, downsep=[1, -0.561])

    geofid = h5open(geo_path, "r")
    grd = load_grid(geofid, "Grid")
    GC_dchlt = load_ghost_conditions(geofid, "GC_dchlt")
    lcfs_avg = load_sparse_matrix(geofid, "lcfs_avg")
    H1 = load_sparse_matrix(geofid, "H1")
    H2 = load_sparse_matrix(geofid, "H2")
    H3 = load_sparse_matrix(geofid, "H3")
    FL_img = geofid["FL_img"][:,:]
    close(geofid)

    Minv = stencil1d(5)
    MinvT = stencil2d(4, 4)
    Dx = derivative_matrix(1, 0, Minv, Minv, grd) * 0.3
    Dy = derivative_matrix(0, 1, Minv, Minv, grd) * 0.3
    Dxy = derivative_matrix(1, 1, Minv, Minv, grd) * (0.3)^2
    Dxx = derivative_matrix(2, 0, Minv, Minv, grd) * (0.3)^2
    Dyy = derivative_matrix(0, 2, Minv, Minv, grd) * (0.3)^2
    Dxxx = derivative_matrix(3, 0, Minv, Minv, grd) * (0.3)^3
    Dyyy = derivative_matrix(0, 3, Minv, Minv, grd) * (0.3)^3
    Dxxy = derivative_matrix(2, 1, Minv, Minv, grd) * (0.3)^3
    Dxyy = derivative_matrix(1, 2, Minv, Minv, grd) * (0.3)^3
    Ds, Dss = fieldline_derivatives(FL_img, 4, 4, MinvT, grd)

    psi_in = psi(1 - 0.28, 0)
    dp = abs(psi_in)

    n = Te = Ti = f_to_grid(grd) do x,y
        psi_xy = psi(x,y)
        (y > -0.561) ? (3/4)*(tanh((psi_in/2 - psi_xy)/dp)+1) + 0.5 : 0.5
    end
    uvec = f_to_grid(grd) do x,y
        0.5*(tanh((-0.4 - y)/0.05) + 1)*tanh((1-x)/0.05)
    end
    u = hcat(uvec, uvec, uvec)
    Sp = f_to_grid(grd) do x,y
        psi_xy = psi(x,y)
        (y > -0.561) ? 1/2 * (tanh((psi_in/2 - psi_xy)/dp) + 1) : 0.0
    end

    Pe = n .* Te
    Pi = n .* Ti

    lnn = lnTe = lnTi = begin v = log.(n); hcat(v, v, v) end
    j = jn = psi = zeros(grd.Nk)
    u = w = A = zeros(grd.Nk, 3)

    lnn_x = Dx * lnn[:,2]
    lnn_y = Dy * lnn[:,2]
    Pi_xx = Dxx * Pi
    Pi_yy = Dyy * Pi
    phi_b = bval_phi(w[:,2], Te, lcfs_avg, H1, H2, H3)
    ad = 1.0 # XXX
    phi = solve_vorticity_eqn(phi_b, w[:,2], n, lnn_x, lnn_y, Pi_xx, Pi_yy, ad, Dx, Dy, Dxx, Dyy, GC_dchlt)

    fid = h5open(init_path, "w")
    save_sparse_matrix(fid, Dx, "Dx")
    save_sparse_matrix(fid, Dy, "Dy")
    save_sparse_matrix(fid, Dxy, "Dxy")
    save_sparse_matrix(fid, Dxx, "Dxx")
    save_sparse_matrix(fid, Dyy, "Dyy")
    save_sparse_matrix(fid, Dxxx, "Dxxx")
    save_sparse_matrix(fid, Dyyy, "Dyyy")
    save_sparse_matrix(fid, Dxxy, "Dxxy")
    save_sparse_matrix(fid, Dxyy, "Dxyy")
    save_sparse_matrix(fid, Ds, "Ds")
    save_sparse_matrix(fid, Dss, "Dss")
    fid["lnn"] = lnn[:,:]
    fid["lnTe"] = lnTe[:,:]
    fid["lnTi"] = lnTi[:,:]
    fid["u"] = u[:,:]
    fid["w"] = w[:,:]
    fid["A"] = A[:,:]
    fid["phi"] = phi[:]
    fid["psi"] = psi[:]
    fid["n"] = n[:]
    fid["Te"] = Te[:]
    fid["Ti"] = Ti[:]
    fid["Pe"] = Pe[:]
    fid["Pi"] = Pi[:]
    fid["j"] = j[:]
    fid["jn"] = jn[:]
    fid["Sp"] = Sp[:]
    fid["params"] = ones(Float64, 9)
    # fid["params"] = Float64[am, ad, ki, ke, er, eg, ev, de2, eta]
    close(fid)
end

function test_simulate(Nt, sn, ste, sti, kdiff, dt, output_path, init_path, geo_path)

    geofid = h5open(geo_path, "r")
    grd = load_grid(geofid, "Grid")
    K1 = load_sparse_matrix(geofid, "K1")
    K2 = load_sparse_matrix(geofid, "K2")
    H1 = load_sparse_matrix(geofid, "H1")
    H2 = load_sparse_matrix(geofid, "H2")
    H3 = load_sparse_matrix(geofid, "H3")
    GC_dchlt = load_ghost_conditions(geofid, "GC_dchlt")
    GC_nmann = load_ghost_conditions(geofid, "GC_nmann")
    GC_u = load_ghost_conditions(geofid, "GC_u")
    trgt_sgn = geofid["trgt_sgn"][:]
    lcfs_avg = load_sparse_matrix(geofid, "lcfs_avg")
    FL_img = geofid["FL_img"][:,:]
    close(geofid)

    fid = h5open(init_path, "r")
    Dx = load_sparse_matrix(fid, "Dx")
    Dy = load_sparse_matrix(fid, "Dy")
    Dxy = load_sparse_matrix(fid, "Dxy")
    Dxx = load_sparse_matrix(fid, "Dxx")
    Dyy = load_sparse_matrix(fid, "Dyy")
    Dxxx = load_sparse_matrix(fid, "Dxxx")
    Dyyy = load_sparse_matrix(fid, "Dyyy")
    Dxxy = load_sparse_matrix(fid, "Dxxy")
    Dxyy = load_sparse_matrix(fid, "Dxyy")
    Ds = load_sparse_matrix(fid, "Ds")
    Dss = load_sparse_matrix(fid, "Dss")
    lnn = fid["lnn"][:,:]
    lnTe = fid["lnTe"][:,:]
    lnTi = fid["lnTi"][:,:]
    u = fid["u"][:,:]
    w = fid["w"][:,:]
    A = fid["A"][:,:]
    phi = fid["phi"][:]
    psi = fid["psi"][:]
    n = fid["n"][:]
    Te = fid["Te"][:]
    Ti = fid["Ti"][:]
    Pe = fid["Pe"][:]
    Pi = fid["Pi"][:]
    j = fid["j"][:]
    jn = fid["jn"][:]
    Sp = fid["Sp"][:]
    am, ad, ki, ke, er, eg, ev, de2, eta = fid["params"][:]
    close(fid)
    am, ad, ki, ke, er, eg, ev, de2, eta = ones(Float64, 9)

    fid = h5open(output_path, "w")
    create_dataset(fid, "n", datatype(Float64), dataspace(grd.Nk*grd.Nz,Nt), chunk=(grd.Nk*grd.Nz,1))
    create_dataset(fid, "Te", datatype(Float64), dataspace(grd.Nk*grd.Nz,Nt), chunk=(grd.Nk*grd.Nz,1))
    create_dataset(fid, "Ti", datatype(Float64), dataspace(grd.Nk*grd.Nz,Nt), chunk=(grd.Nk*grd.Nz,1))
    create_dataset(fid, "u", datatype(Float64), dataspace(grd.Nk*grd.Nz,Nt), chunk=(grd.Nk*grd.Nz,1))
    create_dataset(fid, "phi", datatype(Float64), dataspace(grd.Nk*grd.Nz,Nt), chunk=(grd.Nk*grd.Nz,1))
    create_dataset(fid, "psi", datatype(Float64), dataspace(grd.Nk*grd.Nz,Nt), chunk=(grd.Nk*grd.Nz,1))
    fid["n"][:,1] = n[:]
    fid["Te"][:,1] = Te[:]
    fid["Ti"][:,1] = Ti[:]
    fid["u"][:,1] = u[:,2]
    fid["phi"][:,1] = phi[:]
    fid["psi"][:,1] = psi[:]
    close(fid)

    Sn = sn * Sp
    STe = ste * Sp
    STi = sti * Sp
    kdiff_lnn, kdiff_lnTe, kdiff_lnTi, kdiff_u, kdiff_w, kdiff_A = kdiff[:]

    Nskip = 10
    p = Progress(Nskip*(Nt-1), desc="Running GDB... ")
    for t=2:Nt
        for _=1:10
            leapfrog!(lnn, lnTe, lnTi, u, w, A, phi, psi, n, Te, Ti, Pe, Pi, j, jn, Sn, STe, STi, Dx, Dy, Dxy, Dxx, Dyy, Dxxx, Dyyy,
                      Dxxy, Dxyy, Ds, Dss, dt, 10, K1, K2, H1, H2, H3, trgt_sgn, lcfs_avg, GC_nmann, GC_dchlt, GC_u, kdiff_lnn,
                      kdiff_lnTe, kdiff_lnTi, kdiff_u, kdiff_w, kdiff_A, am, ad, ki, ke, er, eg, ev, de2, eta)
            next!(p)
        end

        fid = h5open(output_path, "r+")
        fid["n"][:,t] = n[:]
        fid["Te"][:,t] = Te[:]
        fid["Ti"][:,t] = Ti[:]
        fid["u"][:,t] = u[:,2]
        fid["phi"][:,t] = phi[:]
        fid["psi"][:,t] = psi[:]
        close(fid)

    end

end


text:
#=
P28 (**) Sorting a list of lists according to length of sublists
a) We suppose that a list contains elements that are lists themselves. The objective is to sort the elements of this list according to their length. E.g. short lists first, longer lists later, or vice versa.

Example:
* (lsort '((a b c) (d e) (f g h) (d e) (i j k l) (m n) (o)))
((O) (D E) (D E) (M N) (A B C) (F G H) (I J K L))

b) Again, we suppose that a list contains elements that are lists themselves. But this time the objective is to sort the elements of this list according to their length frequency; i.e., in the default, where sorting is done ascendingly, lists with rare lengths are placed first, others with a more frequent length come later.

Example:
* (lfsort '((a b c) (d e) (f g h) (d e) (i j k l) (m n) (o)))
((i j k l) (o) (a b c) (f g h) (d e) (d e) (m n))

Note that in the above example, the first two lists in the result have length 4 and 1, both lengths appear just once. The third and forth list have length 3 which appears twice (there are two list of this length). And finally, the last three lists have length 2. This is the most frequent length.
=#

function lsort(colls)
    sort(colls, by=length)
end

colls = Vector{}[[:a, :b, :c], [:d, :e], [:f, :g, :h], [:d, :e], [:i, :j, :k, :l], [:m, :n], [:o]]
sorted_by_length = lsort(colls)
@assert sorted_by_length == Vector{}[[:o], [:d, :e], [:d, :e], [:m, :n], [:a, :b, :c], [:f, :g, :h], [:i, :j, :k, :l]]

println("Tests passed: JL-28.jl")


text:
export simplify_expi

SymEngine.free_symbols(syms::Union{Real, Complex}) = Basic[]
SymEngine.free_symbols(syms::AbstractArray{T}) where {T<:Union{Real, Complex}} = Basic[]
SymEngine.free_symbols(syms::AbstractArray{T}) where {T<:Union{Basic, SymEngine.BasicType}} = SymEngine.free_symbols(Matrix(syms))
function rand_assign(syms...)
    fs = union(free_symbols.(syms)...)
    Dict(zip(fs, randn(length(fs))))
end

function _basic_approx(x, y; atol=1e-8)
    diff = x-y
    assign = rand_assign(x, y)
    length(assign) > 0 && (diff = subs.(diff, Ref.((assign...,))...))
    nres = ComplexF64.(diff)
    all(isapprox.(nres, 0; atol=atol))
end

Base.:≈(x::AbstractArray{<:Basic}, y::AbstractArray; atol=1e-8) = _basic_approx(x, y, atol=atol)
Base.:≈(x::AbstractArray, y::AbstractArray{<:Basic}; atol=1e-8) = _basic_approx(x, y, atol=atol)
Base.:≈(x::AbstractArray{<:Basic}, y::AbstractArray{<:Basic}; atol=1e-8) = _basic_approx(x, y, atol=atol)
function Base.Complex{T}(a::Basic) where T
    a = SymEngine.evalf(a)
    T(real(a)) + im*T(imag(a))
end

# simplification rules
function simplify_expi(expr::Basic)
    subs(expr, exp(im*Basic(π)/2)=>Basic(im), exp(im*Basic(π))=>-1, exp(im*Basic(π)*3/2)=>-Basic(im))
end


text:
using POMDPs
using ARDESPOT
using POMDPToolbox
using POMDPModels
using ProgressMeter

T = 50
N = 50

pomdp = BabyPOMDP()

bounds = IndependentBounds(DefaultPolicyLB(FeedWhenCrying()), 0.0)
# bounds = IndependentBounds(reward(pomdp, false, true)/(1-discount(pomdp)), 0.0)

solver = DESPOTSolver(epsilon_0=0.1,
                      K=100,
                      D=50,
                      lambda=0.01,
                      bounds=bounds,
                      T_max=Inf,
                      max_trials=500,
                      rng=MersenneTwister(4),
                      # random_source=SimpleMersenneSource(50),
                      random_source=FastMersenneSource(100, 10),
                      # random_source=MersenneSource(50, 10, MersenneTwister(10))
                     )

@show solver.lambda

rsum = 0.0
fwc_rsum = 0.0
@showprogress for i in 1:N
    planner = solve(solver, pomdp)
    sim = RolloutSimulator(max_steps=T, rng=MersenneTwister(i))
    fwc_sim = deepcopy(sim)
    rsum += simulate(sim, pomdp, planner)
    fwc_rsum += simulate(fwc_sim, pomdp, FeedWhenCrying())
end

@show rsum/N
@show fwc_rsum/N


text:
include("ConstantSpace.jl")
include("SumSpace.jl")
include("ArraySpace.jl")
include("ProductSpaceOperators.jl")
include("SubSpace.jl")
include("QuotientSpace.jl")

⊕(A::Space,B::Space)::Any = domainscompatible(A,B) ? SumSpace(A,B) : PiecewiseSpace(A,B)
⊕(f::Fun,g::Fun) = Fun(space(f) ⊕ space(g), interlace(coefficients(f),coefficients(g)))

⊕(f::Fun,g::Fun,h::Fun...) = ⊕((f ⊕ g), h...)

+(A::Space,B::Space) = A ⊕ B

#split the cfs into component spaces
function coefficients(cfs::AbstractVector,A::SumSpace,B::SumSpace)
    if spacescompatible(A,B)
        cfs
    else
        mapreduce(f->Fun(f,B),+,components(Fun(A,cfs))).coefficients
    end
end
function coefficients(cfs::AbstractVector,A::PiecewiseSpace,B::PiecewiseSpace)
    if spacescompatible(A,B)
        cfs
    else
        mapreduce(f->Fun(f,B),+,components(Fun(A,cfs))).coefficients
    end
end

# spread a single space into a sum space by placing
# its coefficients depending on k
function interlacewithzeros(cfs::AbstractVector,k,it)
    n = length(cfs)

    ret = Array{eltype(cfs)}(undef, 0)
    n == 0 && return ret

    for (K,j) in it
        if K == k
            push!(ret,cfs[j])
            n == j && break
        else
            push!(ret,0)
        end
    end

    ret
end

interlacewithzeros(cfs::AbstractVector,k,B::DirectSumSpace) = interlacewithzeros(cfs,k,interlacer(B))

function sumspacecoefficients(cfsin::AbstractVector,A::Space,B::SumSpace)
    m=length(components(B))

    for k=1:m
        if isconvertible(A,component(B,k))
            cfs = coefficients(cfsin,A,component(B,k))
            return interlacewithzeros(cfs,k,B)
        end
    end

    defaultcoefficients(cfsin,A,B)
end

function sumspacecoefficients(cfsin::AbstractVector,A::Space,B::PiecewiseSpace)
    m=length(components(B))

    for k=1:m
        if domain(component(B,k)) == domain(A) && isconvertible(A,component(B,k))
            cfs = coefficients(cfsin,A,component(B,k))
            return interlacewithzeros(cfs,k,B)
        end
    end

    defaultcoefficients(cfsin,A,B)
end

for TYP in (:SumSpace,:PiecewiseSpace), ATYP in (:ConstantSpace,:(ConstantSpace{<:Domain{<:Number}}),:Space)
    @eval coefficients(cfsin::AbstractVector,A::$ATYP,B::$TYP) = sumspacecoefficients(cfsin,A,B)
end

## LowRank Constructors

# convert a vector of functionals and an operator to a LowRnakPertOperator
# the rangespace is a DirectSumSpace specified by ST of the input rangespaces
# the default is a  ArraySpace, but support is there for PiecewiseSpace
# for bcs
for TYP in (:PiecewiseSpace,:VectorSpace)
    @eval function LowRankPertOperator(A::AbstractVector{OT},::Type{$TYP}) where OT<:Operator
        A=promotedomainspace(A)
        for k=1:length(A)-1
            @assert isafunctional(A[k])
        end
        @assert isbanded(A[end])
        L=LowRankOperator(A[1:end-1],$TYP)
        # add zero functionals to shift down
        BB=[fill(ZeroOperator(domainspace(BB),ConstantSpace()),length(A)-1);A[end]]
        S=InterlaceOperator(BB,domainspace(BB),$TYP(map(rangespace,A)))
        L+S
    end
end

LowRankPertOperator(A::AbstractVector{OT}) where {OT<:Operator}=LowRankPertOperator(A,VectorSpace)

function LowRankOperator(Bin::AbstractVector{FT},::Type{PiecewiseSpace}) where FT<:Operator
    B=promotedomainspace(Bin)
    rsp=PiecewiseSpace(map(rangespace,B))
    LowRankOperator(
        VFun{typeof(rsp),Float64}[Fun(rsp,[zeros(k-1);1]) for k=1:length(B)],
        B)
end

function LowRankOperator(Bin::AbstractVector{FT},::Type{VectorSpace}) where FT<:Operator
    B=promotedomainspace(Bin)
    rsp=Space([map(rangespace,B);ZeroSpace()])  #TODO: Why the hack?
    LowRankOperator(
        VFun{typeof(rsp),Float64}[Fun(rsp,[zeros(k-1);1]) for k=1:length(B)],
        B)
end

LowRankOperator(Bin::AbstractVector{FT}) where {FT<:Operator} = LowRankOperator(Bin,VectorSpace)

"""
 WeightSpace represents a space that weights another space.
 Overload weight(S,x).
"""
abstract type WeightSpace{S,DD,RR} <: Space{DD,RR} end

domain(S::WeightSpace) = domain(S.space)
points(sp::WeightSpace,n) = points(sp.space,n)

weight(S::WeightSpace, x...) = error("Override `weight(::$S, $x)`")

struct WeightSpacePlan{S,P,T,V}
    space::S
    plan::P
    points::Vector{T}
    weights::Vector{V}
end

struct IWeightSpacePlan{S,P,T,V}
    space::S
    plan::P
    points::Vector{T}
    weights::Vector{V}
end

function plan_transform(S::WeightSpace,vals::AbstractVector)
    pts=points(S,length(vals))
    WeightSpacePlan(S,plan_transform(S.space,vals),pts,weight.(Ref(S),pts))
end

function plan_itransform(S::WeightSpace,vals::AbstractVector)
    pts=points(S,length(vals))
    IWeightSpacePlan(S,plan_itransform(S.space,vals),pts,weight.(Ref(S),pts))
end

*(P::WeightSpacePlan,vals::AbstractVector) = P.plan*(vals./P.weights)
*(P::IWeightSpacePlan,cfs::AbstractVector) = P.weights.*(P.plan*cfs)

# used for ProductFun
transform(sp::WeightSpace,vals::AbstractVector,plan::WeightSpacePlan) =
    transform(sp.space,vals./(sp==plan.space ? plan.weights : weight.(sp,plan.points)),plan.plan)
itransform(sp::WeightSpace,cfs::AbstractVector,plan::WeightSpacePlan) =
    itransform(sp.space,cfs,plan.plan).*(sp==plan.space ? plan.weights : weight.(sp,plan.points))

evaluate(f::AbstractVector,S::WeightSpace,x) = weight(S,x)*evaluate(f,S.space,x)

function evaluate(f::AbstractVector,S::WeightSpace,x...)
    fv=evaluate(f,S.space,x...)
    weight(S,x...).*fv
end

# recurrence is inhereted
for FUNC in (:recα,:recβ,:recγ)
    @eval $FUNC(T,ws::WeightSpace,k) = $FUNC(T,ws.space,k)
end

include("HeavisideSpace.jl")
include("DiracSpace.jl")


text:
"""
get_particle_string(PART_TYPE, p::Vector{<:Real})

Description of get_particle_string
--------------------------------
Returns one-line string for particle in form "TYPE pₓ p𝔂 p𝐳 0"

example:

"electron px = 1.0 py = 1.0 pz = 1.0 dt = 0"

"3 1.0 1.0 1.0 0"
"""
function get_particle_string(PART_TYPE, p::Vector{<:Real})

    return string(PART_TYPE," ",p[2]," ",p[3]," ",p[4]," 0")
end

function get_particle_string(PART_TYPE, p::DataFrameRow)

    return string(PART_TYPE," ",p.px," ",p.py," ",p.pz," 0")
end

"""
    get_header_string(ID)

Description of get_header_string
------------------------------
Returns one-line string for header in form "ID 0 2" (2 stands for number of particles)

example:

"ID = 0 t = 0 particles = 2"

"0 0 2"
"""
function get_header_string(ID)
    return string(ID," 0 2")
end

"""
```get_event_string(ID, PART_TYPE, p1::Vector{<:Real}, p2::Vector{<:Real})```
```get_event_string(ID, PART_TYPE, p1::DataFrameRow, p2::DataFrameRow)```

Description of get_event_string
------------------------------
Returns one event, separated at the end by #.
"""
function get_event_string(ID, PART_TYPE, p1::Vector{<:Real}, p2::Vector{<:Real})
    sHeader= get_header_string(ID)
    sPart1 = get_particle_string(PART_TYPE, p1)
    sPart2 = get_particle_string(PART_TYPE, p2)

    return sHeader*"\n"*sPart1*"\n"*sPart2*"\n#\n"
end

function get_event_string(ID, PART_TYPE, p1::DataFrameRow, p2::DataFrameRow)
    sHeader= get_header_string(ID)
    sPart1 = get_particle_string(PART_TYPE, p1)
    sPart2 = get_particle_string(PART_TYPE, p2)

    return sHeader*"\n"*sPart1*"\n"*sPart2*"\n#\n"
end


text:
using MathProgBase, GLPKMathProgInterface

include("misc.jl")

function leastAbsolutes(X, y)
  (n, d) = size(X)
  X = [ones(n, 1) X]
  d += 1
  y = y[:,1]
  c = [ones(n); zeros(d)]
  A = [eye(n) -X; eye(n) X]
  b = [-y; y]

  display(size(y))
  display(size(c))
  display(size(A))
  display(size(b))

  l = [zeros(n); ones(d) * -Inf]

  sol = linprog(c, A, '>', b, l, Inf, GLPKSolverLP())
  w = sol.sol[n+1:n+d]
  if sol.status == :Optimal
      println("Optimal objective value is $(sol.objval)")
      println("Optimal solution vector is: $(w)")
  else
      println("Error: solution status $(sol.status)")
  end
  predict(Xtilde) = [ones(size(Xtilde,1),1) Xtilde]*w

  # Return model
  return LinearModel(predict, w)
end


text:
using ProgressMeter:Progress, next!

export stop_at_step, stop_at_episode, stop_when_done, steps_per_episode, rewards_of_each_episode

"""
    stop_at_step(n::Int, is_show_progress::Bool=true)

Return a function, which will return false after been called `n` times.
`is_show_progress` will control whether print the progress meter or not.
"""
function stop_at_step(n::Int, is_show_progress::Bool=true)
    i, p = 0, Progress(n)
    function f(env, agent)
        i += 1
        is_show_progress && next!(p; showvalues = [(:step, i)])
        i ≥ n
    end
    f() = i
end

"""
    stop_at_episode(n::Int, is_show_progress::Bool=true)

Return a function, which will return false after `n` episodes.
`isend(env)` is used to check if it is the end of an episode.
`is_show_progress` will control whether print the progress meter or not.
"""
function stop_at_episode(n::Int, is_show_progress::Bool=true)
    i, p = 0, Progress(n)
    function f(env, agent)
        if isend(env)
            i += 1
            is_show_progress && next!(p; showvalues = [(:episode, i)])
            i ≥ n
        else
            false
        end
    end
    f() = i
end

"Return false when encountered an end of an episode"
stop_when_done(env, agent) = isend(env)

"A callback(closure) which will record the length of each episode"
function steps_per_episode()
    steps = []
    count = 0
    function acc(env, agent)
        count += 1
        if isend(env)
            push!(steps, count)
            count = 0
        end
    end
    function acc()
        steps
    end
end

"""
A callback(closure) which will record the total reward of each episode.
Only support single agent yet.
"""
function rewards_of_each_episode()
    rewards = []
    r = 0.
    function acc(env, agent)
        r += buffer(agent).reward[end]
        if isend(env)
            push!(rewards, r)
            r = 0.
        end
    end
    function acc()
        rewards
    end
end


text:
using Feather
using Arrow  # TODO: review whether Arrow needs to be a dependency
# using Missings

function create_folder(p::AbstractString)
    if isdir(p)
        return joinpath(pwd(),p)
    else
        mkdir(p)
        return joinpath(pwd(),p)
    end
end

function add_info!(dtr::DieterModel)
    r = dtr.results

    for (k,df) in r
        for s in [:min_res, :ev, :heat, :h2, :timestep, :scen]
            if !(dtr.settings[s] |> ismissing)
                df[!,s] .= dtr.settings[s]
            else
                @debug "Setting missing column $(s) to have entries with value -1"
                df[!,s] .= -1
            end
        end
        if !isempty(df)
            df[!,:Value] = df[!,:Value] ./ 1000
        end
    end

    for k  in [:N, :G]
        r[k][!,:fuel] = [dtr.parameters[:Fuel][x] for x in r[k][!,:Technologies]]
        r[k][!,:renewable] = [dtr.parameters[:Renewable][x] == 1 for x in r[k][!,:Technologies]]
        r[k][!,:dispatchable] = [dtr.parameters[:Dispatchable][x] == 1 for x in r[k][!,:Technologies]]
    end

    return dtr
end

"""
Save the results of the solved DieterModel `dtr` to `<rdir>/<singledir>/`
"""
function save_results(dtr::DieterModel, rdir::AbstractString; singledir::AbstractString="single_results")

    add_info!(dtr)

    scen = dtr.settings[:scen]
    min_res = dtr.settings[:min_res]
    ev = coalesce(dtr.settings[:ev],"NUL")
    heat = coalesce(dtr.settings[:heat],"NUL")
    h2 = coalesce(dtr.settings[:h2],"NUL")
    timestep = dtr.settings[:timestep]

    fn  = "-scen_$(scen)-res_$(min_res)%-ev_$(ev)-heat_$(heat)%-h2_$(h2)twh-timestep_$(timestep).csv"
    fn2 = "-scen_$(scen)-res_$(min_res)%-ev_$(ev)-heat_$(heat)%-h2_$(h2)twh-timestep_$(timestep).feather"

    create_folder(rdir)
    srdir = create_folder(joinpath(rdir,singledir))

    for (k,v) in dtr.results
        if !isempty(v)
            path = create_folder(joinpath(srdir,string(k)))

            # fp = joinpath(path, string(k)*fn)
            # CSV.write(fp, v)

            fp2 = joinpath(path, string(k)*fn2)
            # fp2 = joinpath(srdir, string(k)*fn2)
            Feather.write(fp2, v)
        end
    end

    # CSV.write(joinpath(rdir, "load.csv"), DataFrame(Load=dtr.parameters[:Load]))
    fpath = joinpath(rdir, "load.feather")
    isfile(fpath) || Feather.write(fpath, DataFrame(Load=dtr.parameters[:Load]))

    return nothing
end


text:
# This file is part of InvertedFiles.jl

export WeightedInvertedFile

"""
    struct WeightedInvertedFile <: AbstractInvertedFile

An inverted index is a sparse matrix representation of with floating point weights, it supports only positive non-zero values.
This index is optimized to efficiently solve `k` nearest neighbors (cosine distance, using previously normalized vectors).

# Parameters

- `lists`: posting lists (non-zero id-elements in rows)
- `weights`: non-zero weights (in rows)
- `sizes`: number of non-zero values in each element (non-zero values in columns)
- `locks`: per-row locks for multithreaded construction
"""
struct WeightedInvertedFile <: AbstractInvertedFile
    lists::Vector{Vector{UInt32}}  ## posting lists
    weights::Vector{Vector{Float32}}  ## associated weights
    sizes::Vector{Int32}  ## number of non zero elements per vector
    locks::Vector{SpinLock}
end

"""
    WeightedInvertedFile(vocsize::Integer)

Convenient function to create an empty `WeightedInvertedFile` with the given vocabulary size.
"""
function WeightedInvertedFile(vocsize::Integer)
    vocsize > 0 || throw(ArgumentError("voc must not be empty"))
    WeightedInvertedFile(
        [Vector{Vector{UInt32}}(undef, 0) for i in 1:vocsize],
        [Vector{Float32}(undef, 0) for i in 1:vocsize],
        Vector{UInt32}(undef, 0),
        [SpinLock() for i in 1:vocsize]
    )
end

function _internal_push!(idx::WeightedInvertedFile, tokenID, objID, weight, sort)
    push!(idx.lists[tokenID], objID)
    push!(idx.weights[tokenID], weight)
    sort && sortlastpush!(idx.lists[tokenID], idx.weights[tokenID])
end


text:
#!/usr/bin/env julia

function main()
    if length(ARGS) > 0
        n = parse(Int64, ARGS[1])
        save_n = n
        factors = Int64[]
        d = 2
        while n > 1
            while n % d == 0
                push!(factors, d)
                n = div(n,d)
            end
            d += 1
            if d * d > n
                if n > 1
                    push!(factors, n)
                end
                break
            end
        end

        @printf("The largest prime factor of %d is %d\n",
                save_n, factors[length(factors)])
    else
        error("this program needs a single numeric argument")
    end
end

main()


text:
using DashBootstrapComponents

dropdown = html_div([
    dbc_dropdownmenu(
        [
            dbc_dropdownmenuitem("A button", id = "dropdown-button", n_clicks = 0),
            dbc_dropdownmenuitem("Internal link", href = "/docs/components/dropdown_menu"),
            dbc_dropdownmenuitem("External Link", href = "https://github.com"),
            dbc_dropdownmenuitem(
                "External relative",
                href = "/docs/components/dropdown_menu",
                external_link = true,
            ),
        ],
        label = "Menu",
    ),
    html_p(id = "item-clicks", className = "mt-3"),
]);

callback!(app, Output("item-clicks", "children"), Input("dropdown-button", "n_clicks")) do n
    return n > 0 ? "Button clicked $n times." : "Button not clicked yet."
end;

max_issues_repo_name: sthagen/dash-bootstrap-components

max_forks_repo_name: sthagen/dash-bootstrap-components


text:
using RGlasso
using Base.Test

# none yet


text:
# # Modal analysis of a free-floating steel circle: three-dimensional solid model

# Source code: [`circle_3d_modal_conv_tut.jl`](circle_3d_modal_conv_tut.jl)

# ## Description

# Vibration analysis of a free-floating steel ring. This is a
# benchmark from the NAFEMS Selected Benchmarks for Natural Frequency Analysis,
# publication: Test VM09: Circular Ring --  In-plane and Out-of-plane
# Vibration.

# The results can be compared with analytical expressions, but the main purpose
# is to compute data for extrapolation to the limit to predict the true natural
# frequencies.

# ## Reference frequencies

# There will be 6 rigid body modes (zero natural frequencies).

# The numerical results are due to the publication:
# NAFEMS Finite Element Methods & Standards, Abbassian, F., Dawswell, D. J., and
# Knowles, N. C. Selected Benchmarks for Natural Frequency Analysis, Test No.
# 6. Glasgow: NAFEMS, Nov., 1987.

# The reference values were analytically determined (Blevins, FORMULAS FOR
# DYNAMICS, ACOUSTICS AND VIBRATION, Table 4.16). Note that shear flexibility
# was neglected when computing the reference values.

# | Mode       |         Reference Value (Hz)  |  NAFEMS Target Value (Hz) |
# | -------   |     -------  |  ------- |
# | 7, 8 | (out of plane)   |        51.85          |         52.29  |
# | 9, 10 |  (in plane)       |       53.38         |          53.97  |
# | 11, 12 |  (out of plane)   |     148.8          |         149.7  |
# | 13, 14 |  (in plane)       |     151.0          |         152.4  |
# | 15, 16 |  (out of plane)   |     287.0          |         288.3  |
# |  17, 18 |  (in plane)      |      289.5         |          288.3  |

# ## Goals

# - Show convergence relative to reference values.
# - Compute data for extrapolation to the limit to predict the true natural
#   frequencies.

# ## Definition of the basic inputs

# The finite element code realize on the basic functionality implemented in this
# package.
using FinEtools
using FinEtoolsDeforLinear
using SymRCM
using Arpack
using LinearAlgebra

# The material parameters may be defined with the specification of the units.
# The elastic properties are:
E = 200.0 * phun("GPa")
nu = 0.3;

# The mass density is
rho = 8000 * phun("kg/m^3")

# The geometry consists of the radius of the circle, and the diameter of the
# circular cross-section.
radius = 1.0 * phun("m"); diameter = 0.1 * phun("m");

# We shall calculate these eigenvalues, but we are mostly interested in the
# first three  natural frequencies.
neigvs = 18;

# The mass shift needs to be applied since the structure is free-floating.
oshift = (2*pi*15)^2

# We will generate this many elements per radius of the cross-section, and along
# the length of the circular ring.
results = let
    results = []
    for i in 1:3
        nperradius, nL = 2^i, 20*2^i
        tolerance = diameter/nperradius/100
        # Generate the mesh in a straight cylinder.
        fens, fes = H8cylindern(diameter/2, 2*pi, nperradius, nL)
        # Select the nodes in the bases of the cylinder for future merging of the nodes.
        z0l = selectnode(fens, plane = [0, 0, 1, 0], inflate = tolerance)
        z2l = selectnode(fens, plane = [0, 0, 1, 2*pi], inflate = tolerance)
        # Twist the straight cylinder into a ring.
        for i in 1:count(fens)
            a = fens.xyz[i, 3]
            x = fens.xyz[i, 1]
            y = fens.xyz[i, 2] + radius
            fens.xyz[i, :] .= (x, y*cos(a), y*sin(a))
        end
        # Merge the nodes of the bases, which involves renumbering the connectivity.
        fens, fes = mergenodes(fens, fes, tolerance, vcat(z0l, z2l))

        # If desired, visualize the mesh.
        # File = "ring.vtk"
        # vtkexportmesh(File, fens, fes)
        # @async run(`"paraview.exe" $File`)

        # Generate the material and the FEM machine.
        MR = DeforModelRed3D
        material = MatDeforElastIso(MR, rho, E, nu, 0.0)
        femm = FEMMDeforLinearESNICEH8(MR, IntegDomain(fes, NodalTensorProductRule(3)), material)
        # Set up the nodal fields for the geometry in the displacements.
        geom = NodalField(fens.xyz)
        u = NodalField(zeros(size(fens.xyz,1),3)) # displacement field
        # Renumber the nodes to produce quicker linear equation solves.
        p = let
            C = connectionmatrix(femm, count(fens))
            p = symrcm(C)
        end
        numberdofs!(u, p)
        # Now set up the discrete model.
        associategeometry!(femm,  geom)
        K  = stiffness(femm, geom, u)
        M = mass(femm, geom, u)

        # Solve the free vibration problem.
        evals, evecs, nconv = eigs(Symmetric(K+oshift*M), Symmetric(M); nev=neigvs, which=:SM, explicittransform = :none)
        # Correct for the mass shift.
        evals = evals .- oshift;
        sigdig(n) = round(n * 10000) / 10000
        fs = real(sqrt.(complex(evals)))/(2*pi)
        println("Eigenvalues: $(sigdig.(fs)) [Hz]")
        push!(results, (fs[6+1], fs[6+3], fs[6+5]))
    end
    results # return these values from the block
end

@show results

# ## Richardson extrapolation

# Here we will use Richardson extrapolation from the three sets of data. This
# will allow us to predict the convergence rate and the true solution for each
# of the three frequencies (or rather the pairs of frequencies, 7 and 8, 9 and
# 10, and 11 and 12).

# We will immediately set up the convergence plots. We will extrapolate and then
# compute from that the normalized error to be plotted with respect to the
# refinement factors (which in this case are 4, 2, and 1). We use the
# refinement factor as a convenience: we will calculate the element size by
# dividing the circumference of the ring with a number of elements generated
# circumferentially.

using FinEtools.AlgoBaseModule: richextrapol

using Gnuplot
@gp  "set terminal windows 0 "  :-

# Modes 7 and 8
sols = [r[1] for r in results]
resextrap = richextrapol(sols, [4.0, 2.0, 1.0])
print("Predicted frequency 7 and 8: $(resextrap[1])\n")
errs = abs.(sols .- resextrap[1])./resextrap[1]
@gp  :- 2*pi*radius./[80, 160, 320] errs " lw 2 lc rgb 'red' with lp title 'Mode 7, 8' "  :-

# Modes 9 and 10
sols = [r[2] for r in results]
resextrap = richextrapol(sols, [4.0, 2.0, 1.0])
print("Predicted frequency 9 and 10: $(resextrap[1])\n")
errs = abs.(sols .- resextrap[1])./resextrap[1]
@gp  :- 2*pi*radius./[80, 160, 320] errs " lw 2 lc rgb 'green' with lp title 'Mode 9, 10' "  :-

# Modes 11 and 12
sols = [r[3] for r in results]
resextrap = richextrapol(sols, [4.0, 2.0, 1.0])
print("Predicted frequency 11 and 12: $(resextrap[1])\n")
errs = abs.(sols .- resextrap[1])./resextrap[1]
@gp  :- 2*pi*radius./[80, 160, 320] errs " lw 2 lc rgb 'blue' with lp title 'Mode 11, 12' "  :-

@gp  :- "set xrange [0.01:0.1]" "set logscale x" :-
@gp  :- "set logscale y" :-
@gp  :- "set xlabel 'Element size'" :-
@gp  :- "set ylabel 'Normalized error [ND]'" :-
@gp  :- "set title '3D: Convergence of modes 7, ..., 12'"


text:
"""
ДАНО: Робот находится в произвольной клетке ограниченного прямоугольного поля (присутствуют перегородки)
РЕЗУЛЬТАТ: Робот — в исходном положении в центре прямого креста из маркеров, расставленных вплоть до внешней рамки.
"""
function new_cross(r::Robot)
    side=Nord
    for a ∈ 1:4
        x=1
        while x!=0
            x=find_decision(r,side)
            putmarker!(r)
        end
        while ismarker(r)
            find_decision(r,inverse(side))
        end
        side=inverse(next(side))
    end
    putmarker!(r)
end

function find_decision(r::Robot, side::HorizonSide) #Обход локальных преград
    num_steps=0
    while isborder(r,side)==true&&isborder(r,next(side))==false
        move!(r,next(side))
        num_steps+=1
    end
    x=0 #Счётчик
    if isborder(r,side)==false
        move!(r,side)
        x+=1
    end
    if num_steps !=0
        while isborder(r,inverse(next(side)))==true
            move!(r,side)
            x+=1
        end
        for a ∈ 1:num_steps
            move!(r,inverse(next(side)))
        end
    end
    return x #Возвращение значения счётчика
end

function next(side::HorizonSide)
    (HorizonSide(mod(Int(side)+1,4)))
end
function inverse(side::HorizonSide)
    (HorizonSide(mod(Int(side)+2,4)))
end

max_stars_repo_name: alexgruzberg/mirea-progs

max_issues_repo_name: alexgruzberg/mirea-progs

max_forks_repo_name: alexgruzberg/mirea-progs


text:
#================================================================================================
    functions common to both lists and primitives
================================================================================================#
Base.values(A::Primitive) = A
Base.values(A::ArrowVector) = A.values
bitmask(A::ArrowVector) = A.bitmask
offsets(A::ArrowVector) = A.offsets

Base.copy(A::T) where T<:ArrowVector = A[:]

"""
    valuespointer(A::ArrowVector)

Returns a pointer to the start of the values buffer for `A`.
"""
valuespointer(A::ArrowVector) = valuespointer(A.values)

"""
    bitmaskpointer(A::ArrowVector)

Returns a pointer to the start of the bitmask buffer for `A`.
"""
bitmaskpointer(A::ArrowVector{Union{J,Missing}}) where J = valuespointer(A.bitmask)

"""
    offsetspointer(A::AbstractList)

Returns a pointer to the start of the offsets buffer for `A`.
"""
offsetspointer(A::AbstractList) = valuespointer(A.offsets)

"""
    rawvalues(A::ArrowVector, idx)

Gets a `Vector{UInt8}` of the raw values associated with the indices `idx`.
"""
rawvalues(A::ArrowVector, i) = rawvalues(A.values, i)
rawvalues(A::ArrowVector) = rawvalues(A.values)

# TODO this actually gets fucked up if bitmask has trailing ones (that's also why this is backwards)
"""
    nullcount(A::ArrowVector)

Return the number of nulls (`missing`s) in `A`.
"""
nullcount(A::ArrowVector) = 0
function nullcount(A::ArrowVector{Union{T,Missing}}) where T
    s = 0
    for i ∈ 1:minbitmaskbytes(A)
        s += count_ones(A.bitmask[i])
    end
    length(A) - s
end

function unsafe_getbit(A::Primitive{J}, i::Integer) where J
    a, b = divrem(i-1, 8sizeof(J)) .+ (1,1)
    getbit(unsafe_getvalue(A, a), b)
end
unsafe_getbit(A::Primitive, idx::AbstractVector{<:Integer}) = Bool[unsafe_getbit(A, i) for i ∈ idx]
function unsafe_getbit(A::Primitive, idx::AbstractVector{Bool})
    Bool[unsafe_getbit(A,i) for i ∈ 1:length(A) if idx[i]]
end

"""
    unsafe_isnull(A::ArrowVector, idx)

Check whether element(s) `idx` of `A` are null.
"""
unsafe_isnull(A::ArrowVector, i::Integer) = false
unsafe_isnull(A::ArrowVector{Union{J,Missing}}, i::Integer) where J = !unsafe_getbit(bitmask(A), i)
function unsafe_isnull(A::ArrowVector{Union{J,Missing}}, idx::AbstractVector{<:Integer}) where J
    Bool[unsafe_isnull(A, i) for i ∈ idx]
end
function unsafe_isnull(A::AbstractVector{Union{J,Missing}}, idx::AbstractVector{Bool}) where J
    Bool[unsafe_isnull(A, i) for i ∈ 1:length(A) if idx[i]]
end

"""
    getbit(A::Primitive, i)

Get the `i`th bit of the primitive.  For example, for a `Primitive{UInt8]` the 10th bit is the second
bit of the second value.
"""
function getbit(A::Primitive{J}, i::Integer) where J
    a, b = divrem(i-1, 8sizeof(J)) .+ (1,1)
    getbit(getvalue(A, a), b)
end
getbit(A::Primitive, idx::AbstractVector{<:Integer}) = Bool[getbit(A, i) for i ∈ idx]
getbit(A::Primitive, idx::AbstractVector{Bool}) = Bool[getbit(A, i) for i ∈ 1:lenght(A) if idx[i]]

function setbit!(A::Primitive{J}, x::Bool, i::Integer) where J
    a, b = divrem(i-1, 8sizeof(J)) .+ (1,1)
    byte = setbit(getvalue(A, a), x, b)
    setvalue!(A, byte, a)
end
function setbit!(A::Primitive, x::AbstractVector{Bool}, idx::AbstractVector{<:Integer})
    j = 1
    for i ∈ idx
        setbit!(A, x[j], i)
        j += 1
    end
    x
end
function setbit!(A::Primitive, x::AbstractVector{Bool}, idx::AbstractVector{Bool})
    j = 1
    for i ∈ 1:length(A)
        if idx[i]
            setbit!(A, x[j], i)
            j += 1
        end
    end
    x
end

"""
    isnull(A::ArrowVector, idx)

Check whether element(s) `idx` of `A` are null.
"""
isnull(A::ArrowVector, i) = false
isnull(A::ArrowVector{Union{J,Missing}}, i::Integer) where J = !getbit(bitmask(A), i)
isnull(A::ArrowVector, idx::AbstractVector{<:Integer}) = Bool[isnull(A, i) for i ∈ idx]
isnull(A::ArrowVector, idx::AbstractVector{Bool}) = Bool[isnull(A, i) for i ∈ 1:length(A) if idx[i]]

"""
    unsafe_rawbitmask(p::ArrowVector{Union{J,Missing}})

Retrieve the raw value of the null bit mask for `p`.
"""
function unsafe_rawbitmask(p::ArrowVector{Union{J,Missing}}) where J
    unsafe_rawpadded(bitmaskpointer(A), bitmaskbytes(p))
end

setnull!(A::ArrowVector{Union{J,Missing}}, x::Bool, i::Integer) where J = setbit!(bitmask(A), !x, i)

# TODO these are flipped relative to each other which is super confusing
function setnulls!(A::ArrowVector{Union{J,Missing}}, bytes::Vector{UInt8}) where J
    setvalue!(A.bitmask, bytes, 1:length(bytes))
end
function setnulls!(A::ArrowVector{Union{J,Missing}}, nulls::AbstractVector{Bool}) where J
    setnulls!(A, bitpack(.!nulls))
end

"""
    unsafe_setnull!(A::ArrowVector{Union{J,Missing}}, x::Bool, i::Integer)

Set element `i` of `A` to be null. This involves no bounds checking and a call to `unsafe_store!`.
"""
function unsafe_setnull!(A::ArrowVector{Union{J,Missing}}, x::Bool, i::Integer) where J
    a, b = divrem(i-1, 8) .+ (0,1)
    ptr = bitmaskpointer(A) + a
    byte = setbit(unsafe_load(ptr), !x, b)
    unsafe_store!(ptr, byte)
end

"""
    unsafe_setnulls!(A::ArrowVector, nulls::AbstractVector{Bool})

Set *all* the nulls for the `ArrowVector`. This does not check bounds and contains a call to
`unsafe_copy!` (but does not copy directly from `nulls`).
"""
function unsafe_setnulls!(A::ArrowVector{Union{J,Missing}}, bytes::Vector{UInt8}) where J
    unsafe_copy!(bitmaskpointer(A), pointer(bytes), length(bytes))
end
function unsafe_setnulls!(A::ArrowVector{Union{J,Missing}}, nulls::AbstractVector{Bool}) where J
    unsafe_setnulls!(A, bitpack(.!nulls))
end

macro _make_fillmissings_funcs(name::Symbol, func::Symbol)
esc(quote
    function $name(v::AbstractVector{Union{J,Missing}}, A::ArrowVector{Union{J,Missing}},
                   idx::AbstractVector{<:Integer}) where J
        for (i, j) ∈ enumerate(idx)
            $func(A, j) && (v[i] = missing)
        end
    end
    function $name(v::AbstractVector{Union{J,Missing}}, A::ArrowVector{Union{J,Missing}},
                   idx::AbstractVector{Bool}) where J
        j = 1
        for i ∈ 1:length(A)
            if idx[i]
                $func(A, i) && (v[j] = missing)
                j += 1
            end
        end
    end
    function $name(v::AbstractVector{Union{J,Missing}}, A::ArrowVector{Union{J,Missing}}) where J
        $name(v, A, 1:length(A))
    end
end)
end
@_make_fillmissings_funcs(unsafe_fillmissings!, unsafe_isnull)
@_make_fillmissings_funcs(fillmissings!, isnull)

function setnonmissing!(A::ArrowVector{J}, v::AbstractVector{T}) where {J,T<:Union{J,Union{J,Missing}}}
    @boundscheck length(A) == length(v) || throw(ArgumentError("trying to set from wrong sized array"))
    for i ∈ 1:length(A)
        !ismissing(v[i]) && (A[i] = v[i])
    end
    A
end

Base.length(A::ArrowVector) = A.length
Base.size(A::ArrowVector) = (length(A),)
function Base.size(A::ArrowVector, i::Integer)
    if i == 1
        return length(A)
    else
        return 1
    end
    throw(ArgumentError("arraysize: dimension $i out of range"))
end

Base.convert(::Type{Array{T}}, A::ArrowVector{T}) where T = A[:]
Base.convert(::Type{Vector{T}}, A::ArrowVector{T}) where T = A[:]

Base.IndexStyle(::Type{<:ArrowVector}) = IndexLinear()

# macro for creating arrowformat functions
macro _formats(constructor, argtype, w...)
esc(quote
    arrowformat(x::$argtype) where {$(w...)} = $constructor(x)
    function arrowformat(::Type{<:Array}, x::$argtype) where {$(w...)}
        $constructor(Array, x)
    end
    function arrowformat(data::Vector{UInt8}, i::Integer, x::$argtype) where {$(w...)}
        $constructor(data, i, x)
    end
end)
end

"""
    arrowformat(v::AbstractVector)
    arrowformat(Array, v::AbstractVector)

Convert a vector to the appropriate arrowformat.  If `Array` is passed, a contiguous array will
be used for the data buffer.
"""
function arrowformat end
@_formats Primitive AbstractVector{J} J
@_formats NullablePrimitive AbstractVector{Union{J,Missing}} J
@_formats List AbstractVector{J} J<:AbstractString
@_formats NullableList AbstractVector{Union{J,Missing}} J<:AbstractString
@_formats Primitive{Datestamp} AbstractVector{T} T<:Date
@_formats Primitive{Timestamp{Dates.Millisecond}} AbstractVector{T} T<:DateTime
@_formats Primitive{TimeOfDay{Dates.Nanosecond,Int64}} AbstractVector{T} T<:Dates.Time
@_formats NullablePrimitive{Datestamp} AbstractVector{Union{T,Missing}} T<:Date
@_formats NullablePrimitive{Timestamp{Dates.Millisecond}} AbstractVector{Union{T,Missing}} T<:DateTime
@_formats NullablePrimitive{TimeOfDay{Dates.Nanosecond,Int64}} AbstractVector{Union{T,Missing}} T<:Dates.Time
@_formats DictEncoding CategoricalArray{T,1,U} T<:Any U
@_formats BitPrimitive AbstractVector{J} J<:Bool
@_formats NullableBitPrimitive AbstractVector{J} J<:Union{Bool,Missing}
arrowformat(A::T) where {T<:ArrowVector} = T(A)

# TODO clean up bounds checking

@inline function Base.getindex(l::ArrowVector{J}, i::Union{<:Integer,AbstractVector{<:Integer}}) where J
    @boundscheck checkbounds(l, i)
    @inbounds o = unsafe_getvalue(l, i)
    o
end
@inline function Base.getindex(l::ArrowVector{Union{J,Missing}}, i::Integer)::Union{J,Missing} where J
    @boundscheck checkbounds(l, i)
    @inbounds o = unsafe_isnull(l, i) ? missing : unsafe_getvalue(l, i)
    o
end
@inline function Base.getindex(l::ArrowVector{Union{J,Missing}}, idx::AbstractVector{<:Integer}) where J
    @boundscheck checkbounds(l, idx)
    @inbounds v = convert(Vector{Union{J,Missing}}, unsafe_getvalue(l, idx))
    @inbounds unsafe_fillmissings!(v, l, idx)
    v
end
@inline Base.getindex(l::ArrowVector, ::Colon) = l[1:end]

"""
    writepadded(io::IO, A::Primitive)

Write the values of `A` to `io` with Arrow padding (8 bytes).
"""
function writepadded(io::IO, A::Primitive, idx::Union{<:Integer,AbstractVector{<:Integer}})
    vals = rawvalues(A, idx)
    write(io, vals)
    pad = zeros(UInt8, padding(length(vals)) - length(vals))
    write(io, pad)
    padding(length(vals))
end
writepadded(io::IO, A::Primitive) = writepadded(io, A, 1:length(A))
function writepadded(io::IO, A::ArrowVector, subbuffs::Function...)
    s = 0
    for sb ∈ subbuffs
        s += writepadded(io, sb(A))
    end
    s
end

"""
    rawpadded(A::Primitive)

Get a `Vector{UInt8}` describing the data of `A` with Arrow padding (8 bytes).
"""
function rawpadded(A::Primitive, idx::Union{<:Integer,AbstractVector{<:Integer}})
    vals = rawvalues(A, idx)
    vcat(vals, zeros(UInt8, padding(length(vals)) - length(vals)))
end
rawpadded(A::Primitive) = rawpadded(A, 1:length(A))
rawpadded(A::ArrowVector, subbuffs::Function...) = reduce(vcat, (rawpadded(sb(A)) for sb ∈ subbuffs))


text:
using PihVfi890
using Test

function euler_test()
    @testset "Euler" begin
        m = init_test_model();
        t = 2;
        k = 5.0;
        kPrimeTomorrow(kPrime) = 0.5 .* kPrime;

        kPrime = solve_one_point(m, t, k, kPrimeTomorrow);
        dev = euler_dev_one_point(m, t, k, kPrime, kPrimeTomorrow);
        @test kPrime > kprime_min(m, t);
        @test kPrime < kprime_max(m, t);
        @test abs(dev) < 1e-4
        dev2 = euler_dev_one_point2(m, t, k, kPrime, kPrimeTomorrow);
        @test abs(dev2) < 1e-4
    end
end

function solve_test()
    @testset "Solve" begin
        m = init_test_model();
        kPrime_tV = solve(m);
    end
end

@testset "VFI" begin
    euler_test();
    solve_test();
end

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


text:
using Test
using TrigPolys

@testset "Arithmetic" begin
    p1 = random_trig_poly(102)
    p2 = random_trig_poly(201)
    x = range(0; stop=2*pi, length=2000)

    @test isapprox((-p1).(x), -(p1.(x)))
    @test isapprox((p1 + p2).(x), p1.(x) .+ p2.(x))
    @test isapprox((p1 - p2).(x), p1.(x) .- p2.(x))
    @test isapprox((p1 * p2).(x), p1.(x) .* p2.(x))

    @test TrigPolys.a0(p1 + 3) == TrigPolys.a0(p1) + 3
    @test TrigPoly(pi) == TrigPoly(pi)
    @test p1 + p2 == p2 + p1
    @test p1 * p2 == p2 * p1
    @test p1 + pi == pi + p1
    @test p1 * pi == pi * p1
    @test p1 / pi == p1 * (1/pi)
end

@testset "Type Conversion and Promotion" begin
    F1 = Float32
    F2 = Float64
    VF1 = Vector{F1}
    VF2 = Vector{F2}
    @test typeof(TrigPoly(F2(4))+F1(3)) == TrigPoly{F2, VF2}
    @test typeof(convert(TrigPoly{F2, VF2}, F1(.3))) == TrigPoly{F2, VF2}
    @test typeof(convert(TrigPoly{F1, VF1}, F2(.3))) == TrigPoly{F1, VF1}
    @test typeof(TrigPoly(F1(2))+TrigPoly(F2(3))) == TrigPoly{F2, VF2}
    @test typeof(TrigPoly(F1(2))+F1(3)) == TrigPoly{F1, VF1}
end

@testset "Evaluate and interpolate" begin
    n = 1000
    p = random_trig_poly(n)
    s = p.n*2+1
    samples = [2*pi*(i-1)/s for i in 1:s]
    B = hcat([TrigPolys.basis(p.n, xi) for xi in samples]...)

    u = randn(s)
    @test isapprox(evaluate(u), B'*u)
    @test isapprox(evaluate(p), p.(samples))
    @test isapprox(evaluateT(vec(u)), B*u)
    @test isapprox(interpolatev(evaluate(u)), u)
end


text:
using DataStructures
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
using Plasm

###	Data:  1D complex embedded in 2D
#####################################################################

V = [732.725 1123.87 1123.87 1124.49 1124.49 776.005 776.005 731.819 731.819 732.725 732.725 732.725 1284.06 1592.17 1592.17 1592.32 1592.32 1281.66 1281.66 1284.06 1284.06 1284.06 789.106 789.236 789.236 1030.28 1030.28 1031.02 1031.02 789.106 789.106 789.106 734.892 1087.13 1087.13 1087.63 1087.63 1227.73 1227.73 1226.84 1226.84 1648.35 1648.35 1649.66 1649.66 1591.08 1591.08 1590.19 1590.19 1469.27 1469.27 1470.21 1470.21 1414.42 1414.42 1412.64 1412.64 1178.75 1178.75 1179.1 1179.1 1123.36 1123.36 1123.49 1123.49 677.788 677.788 677.027 677.027 736.012 736.012 734.892 734.892 734.892 692.591 692.119 692.119 710.549 710.549 748.971 748.971 748.398 748.398 774.963 774.963 773.044 773.044 1046.13 1046.13 1046.3 1046.3 1073.16 1073.16 1075.13 1075.13 1240.87 1240.87 1239.87 1239.87 1266.97 1266.97 1266.83 1266.83 1606.23 1606.23 1605.59 1605.59 1632.89 1632.89 1634.74 1634.74 1608.8 1608.8 1605.94 1605.94 1455.92 1455.92 1456.28 1456.28 1431.23 1431.23 1428.3 1428.3 1163.49 1163.49 1164.92 1164.92 1138.7 1138.7 1136.96 1136.96 934.882 934.882 775.749 775.749 776.12 776.12 761.45 761.45 718.072 718.072 719.533 719.533 1137.14 1137.14 1137.32 1137.32 719.091 719.091 719.221 719.221 692.591 692.591 692.591 776.005 1047.47 1047.47 1046.64 1046.64 776.368 776.368 776.005 776.005 776.005 1269.2 1268.89 1268.89 1606.74 1606.74 1606.72 1606.72 1269.2 1269.2 1269.2; 781.1 782.673 782.673 1048.94 1048.94 1045.53 1045.53 1005.26 1005.26 781.1 781.1 781.1 1108.07 1106.08 1106.08 1174.09 1174.09 1172.27 1172.27 1108.07 1108.07 1108.07 1106.08 1288.98 1288.98 1289.83 1289.83 1106.9 1106.9 1106.08 1106.08 1106.08 1349.71 1349.71 1349.71 1107.42 1107.42 1105.77 1105.77 1230.35 1230.35 1230.86 1230.86 615.61 615.61 615.61 615.61 1051.75 1051.75 1052.69 1052.69 615.61 615.61 615.61 615.61 1049.7 1049.7 1050.6 1050.6 615.61 615.61 615.61 615.61 726.053 726.053 726.777 726.777 1024.77 1024.77 1078.88 1078.88 1349.71 1349.71 1349.71 730.096 1024.88 1024.88 1025.99 1025.99 1067.18 1067.18 1345.41 1345.41 1345.84 1345.84 1333.11 1333.11 1333.92 1333.92 1349.71 1349.71 1347.57 1347.57 1094.47 1094.47 1093.37 1093.37 1226.19 1226.19 1226.85 1226.85 1214.78 1214.78 1215.55 1215.55 1227.84 1227.84 1227.51 1227.51 615.946 615.946 615.61 615.61 1066.73 1066.73 1067.74 1067.74 616.213 616.213 616.301 616.301 1066.68 1066.68 1065.78 1065.78 616.712 616.712 617.245 617.245 1067.13 1067.13 1064.62 1064.62 1065.23 1065.23 1057.26 1057.26 1057.24 1057.24 1011.09 1011.09 768.914 768.914 770.027 770.027 742.676 742.676 742.388 742.388 729.39 729.39 730.096 730.096 730.096 1090.67 1092.32 1092.32 1306.35 1306.35 1304.4 1304.4 1090.67 1090.67 1090.67 1093.67 1187.53 1187.53 1187.41 1187.41 1093.67 1093.67 1093.67 1093.67 1093.67]

EV = Array{Int64,1}[[1, 2], [3, 4], [5, 6], [7, 8], [9, 10], [11, 12], [13, 14], [15, 16], [17, 18], [19, 20], [21, 22], [23, 24], [25, 26], [27, 28], [29, 30], [31, 32], [33, 34], [35, 36], [37, 38], [39, 40], [41, 42], [43, 44], [45, 46], [47, 48], [49, 50], [51, 52], [53, 54], [55, 56], [57, 58], [59, 60], [61, 62], [63, 64], [65, 66], [67, 68], [69, 70], [71, 72], [73, 74], [75, 76], [77, 78], [79, 80], [81, 82], [83, 84], [85, 86], [87, 88], [89, 90], [91, 92], [93, 94], [95, 96], [97, 98], [99, 100], [101, 102], [103, 104], [105, 106], [107, 108], [109, 110], [111, 112], [113, 114], [115, 116], [117, 118], [119, 120], [121, 122], [123, 124], [125, 126], [127, 128], [129, 130], [131, 132], [133, 134], [135, 136], [137, 138], [139, 140], [141, 142], [143, 144], [145, 146], [147, 148], [149, 150], [151, 152], [153, 154], [155, 156], [157, 158], [159, 160], [161, 162], [163, 164], [165, 166], [167, 168], [169, 170], [171, 172], [173, 174]]

###	execution
#####################################################################

using Plasm

V,bases,coboundaries = Lar.chaincomplex(V,EV)
objs = Lar.lar2obj2D(V, [coboundaries...])
open("./villa.obj", "w") do f
	write(f, objs)
end

V = convert(Lar.Points,V)
Plasm.view(V,bases[1])
Plasm.view(V,bases[2])

ev,fv = bases
VV = [[k] for k=1:size(V,2)]
model = (V, [VV,ev,fv])
Plasm.view(Plasm.numbering(80.)(model))

objs = Lar.lar2obj2D(convert(Lar.Points,V'), [coboundaries...])
open("./villa.obj", "w") do f
	write(f, objs)
end

V,(EV,FV) = Lar.obj2lar2D("./villa.obj")

V = convert(Lar.Points,V')
Plasm.view(V,EV)
Plasm.view(V,FV)