text:
using JSON3
using Serialization
using TableScraper: scrape_tables

using CodecZlib: GzipDecompressor

export load_namesdb

"""
Save a namesdb file somehwere
"""
function load_namesdb(path="namesdb"; force=false)
    if force | !isfile(path)
        # download latest names database from ugo.net
        latest_row = scrape_tables("https://u-go.net/playerdb-archive/", identity)[1].rows[2]
        link = latest_row[3][1].attributes["href"]
        http = HTTP.request("GET", link)
        json = JSON3.read(transcode(GzipDecompressor, http.body))
        name_dict = create_names_db(json)
        serialize(path, name_dict)
        @info "The Names database is cached at path `$path`"
        return name_dict
    else
        deserialize(path)
    end
end

function create_names_db(json::JSON3.Array)
    # if isfile("c:/git/BadukGoWeiqiTools/src/namesdb") & !force
    #     deserialize("c:/git/BadukGoWeiqiTools/src/namesdb")
    # else
    simplenames = Dict{String,String}()
    for json1 in json
        key = json1.key_name
        for name in json1.names
            any_in_go4go = filter(jn -> "Go4go" in jn.databases, name.simplenames)
            if length(any_in_go4go) >= 1
                key = any_in_go4go[1].name
            end
        end

        for name in json1.names
            for sn in name.simplenames
                simplenames[sn.name] = key
            end
        end

    end

    # manuals
    manuals = Dict("金禹丞" => "Jin Yucheng",
        "王春晖" => "Wang Chunhui",
        "许一笛" => "Xu Yidi",
        "黄明宇" => "Huang Mingyu",
        "申旻俊" => "Shin Minjun",
        "陈一纯" => "Chen Yichun",
        "叶长欣" => "Ye Changxin",
        "李泽锐" => "Li Zerui",
        "吕立言" => "Lyu Liyan",
        "何語涵" => "He Yuhan",
        "王春暉" => "Wang Chunhui",
        "蒋其潤" => "Jiang Qirun",
        "張涛" => "Zhang Tao",
        "李維清" => "Li Weiqing",
        "屠暁宇" => "Tu Xiaoyu",
        "陳賢" => "Chen Xian",
        "郭信駅" => "Guo Xinyi",
        "謝科" => "Xie Ke",
        "陳浩" => "Chen Hao",
        "王沢錦" => "Wang Zejin",
        "楊楷文" => "Yang Kaiwen",
        "劉宇航" => "Liu Yuhang",
        "陳梓健" => "Chen Zijian",
        "陳玉儂" => "Chen Yunong",
        "許一笛" => "Xu Yidi",
        "胡鈺函" => "Hu Yuhan",
        "鄭胥" => "Zheng Xu",
        "伊淩濤" => "Yi Lingtao",
        "韓一洲" => "Han Yizhou",
        "陳正勲" => "Chen Zhengxun",
        "許嘉陽" => "Xu Jiayang",
        "趙晨宇" => "Zhao Chenyu",
        "李沢鋭" => "Li Zerui",
        "陳一純" => "Chen Yichun",
        "葉長欣" => "Ye Changxin",
        "呂立言" => "Lyu Liyan",
        "申真諝" => "Shin Jinseo",
        "姜知勲" => "Kang Jihoon",
        "沈載益" => "Shen Zaiyi",
        "白現宇" => "Baek Hyeonwoo",
        "李元栄" => "Lee Wonyoung",
        "朴常鎭" => "Park Sangjin",
        "崔宰栄" => "Choi Jaeyoung",
        "曹瀟陽" => "Cao Xiaoyang",
        "文儒彬" => "Wen Rubin",
        "文敏鍾" => "Moon Minjong",
        "金炅奐" => "Kim Ilhwan",
        "張強" => "Zhang Qiang",
        "簡靖庭" => "Jian Jingting",
        "賴均輔" => "Lai Junfu",
        "渡辺寛大" => "Watanabe Kandai",
        "崔光戸" => "Choi Kyeongho",
        "黄静遠" => "Huang Jingyuan",
        "陸敏全" => "Lu Minquan",
        "韓宇進" => "Han Woojin",
        "尭瀟童" => "Yao Xiaotong",
        "申載元" => "Shin Jaeweon",
        "成家業" => "Cheng Jiaye",
        "関航太郎" => "Seki Kotaro",
        "盧奕銓" => "Lu Yiquan",
        "鄭宇津" => "Zheng Yujin",
        "仲邑菫" => "Nakamura Sumire",
        "朴鐘勲" => "Park Joonhoon",
        "白賛僖" => "Baek Chanhee",
        "宋知勲" => "Song Jihoon",
        "琴沚玗" => "Gon Jiwoo",
        "金相天" => "Kim Sangcheon",
        "李炫準" => "Lee Hyeonjun",
        "福岡航太朗" => "Fukuoka Kotaro",
        "林漢傑" => "Rin Kanketsu",
        "林傑漢" => "Lin Jiehan",
        "兪俐均" => "Yu Lijun",
        "陳一鳴" => "Chen Yiming",
        "戦鷹" => "Zhan Ying",
        "林聖弈" => "Lin Shengyi",
        "潘陽" => "Pan Yang",
        "趙恵連" => "Cho Hyeyeon",
        "徐靖恩" => "Xu Jingen",
        "呉政娥" => "Oh Jeonga",
        "陳豪鑫" => "Chen Haoxin",
        "王禹程" => "Wang Yucheng",
        "鄭載想" => "Zheng Zaixiang",
        "呂立言" => "Lyu Liyan",
        "肖泽彬" => "Xiao Zebin",
        "何雨涵" => "He Yuhan",
        "陈玉浓" => "Chen Yunong",
        "劉川霆" => "Liu Chuanting",
        "楊子萱" => "Yang Zixuan",
        "白昕卉" => "Bai Xinhui",
        "曺承亜" => "Cho Seungah",
        "劉一芳" => "Liu Yifang",
        "林鈺娗" => "Lin Yuting",
        "金昌勲" => "Kim Changhoon",
        "劉兆哲" => "Liu Zhaozhe",
        "花暢" => "Hua Chang",
        "金凡瑞" => "Kim Beomseo",
        "馬逸超" => "Ma Yichao",
        "王碩" => "Wang Shuo",
        "朴廷恒" => "Park Junghwan",
        "黄春晖" => "Wang Chunhui", # sic
        "刘宇航" => "Liu Yuhang",
        "曹承亚" => "Cho Seungah",
        "羋昱廷" => "Mi Yuting",
        "胡子豪" => "Hu Zihao",
        "赵完珪" => "Cho Wankyu",
        "傅健恒" => "Fu Jianheng",
    )

    merge!(simplenames, manuals)
    simplenames
end

function create_names_db(path_to_json::String; kwargs...)
    text = readline(NAME_DB_JSON);
    json = JSON3.read(text);

    create_names_db(json; kwargs...)
end


text:
#=
# 2D NS eq. Dirichlet bc, CG
=#

### If the Femshop package has already been added, use this line #########
using Femshop # Note: to add the package, first do: ]add "https://github.com/paralab/femshop.git"

### If not, use these four lines (working from the examples directory) ###
# if !@isdefined(Femshop)
#     include("../Femshop.jl");
#     using .Femshop
# end
##########################################################################

init_femshop("NS");

# Optionally generate a log
useLog("NSlog")

# Set up the configuration (order doesn't matter)
domain(2)                   # dimension
solverType(CG)              # Use CG solver (default)
functionSpace(order=1)      # basis polynomial order
nodeType(LOBATTO)           # GLL elemental node arrangement (default)
timeStepper(EULER_IMPLICIT) # time stepper
dt = 0.05;
nsteps = 50;
setSteps(dt, nsteps)        # manually set stepper dt and Nsteps, overriding defaults and interval

# Specify the problem
mesh(QUADMESH, elsperdim=32, bids=4)

u = variable("u")       #Note: many of these variables are needed for
v = variable("v")       #      linearizing the non-linear problem.
uold = variable("uold") #      Eventually, they will be automatically
vold = variable("vold") #      generated and this will become simpler.
p = variable("p")       #
du = variable("du")     #
dv = variable("dv")     #
dp = variable("dp")     #

testSymbol("w")

# Initial conditions
initial(u, "y > 0.99 ? 1 : 0")
initial(uold, "y > 0.99 ? 1 : 0")
initial(du, "0")
initial(v, "0")
initial(vold, "0")
initial(dv, "0")
initial(p, "0")
initial(dp, "0")

# Boundary conditions
boundary(du, 1, DIRICHLET, 0)
boundary(du, 2, DIRICHLET, 0)
boundary(du, 3, DIRICHLET, 0)
boundary(du, 4, DIRICHLET, 0)

boundary(dv, 1, DIRICHLET, 0)
boundary(dv, 2, DIRICHLET, 0)
boundary(dv, 3, DIRICHLET, 0)
boundary(dv, 4, DIRICHLET, 0)

boundary(dp, 1, NO_BC)
boundary(dp, 2, NO_BC)
boundary(dp, 3, NO_BC)
boundary(dp, 4, NO_BC)
# A single Dirichlet reference point must be used to pin values for dp
referencePoint(dp, [0,0], 0)

# Write the weak form
coefficient("mu", 0.01)
coefficient("dtc", dt)
coefficient("h", 1.0 / 32)
coefficient("coe1", 4.0)
coefficient("coe2", 6.0)

parameter("tauM", "1.0 ./ (coe1 ./ dtc ./ dtc+ (u*u+v*v) ./ h ./ h+coe2*mu*mu ./ h ./ h ./ h ./ h) .^ 0.5")
parameter("tauC", "0.1*h*h* ((u*u+v*v) ./ h ./ h+coe2*mu*mu ./ h ./ h ./ h ./ h) .^ 0.5")

weakForm([du, dv, dp], ["w*(du ./ dtc + (u*deriv(du,1)+v*deriv(du,2) + deriv(u,2)*dv)) - deriv(w,1)*dp + mu*dot(grad(w), grad(du)) + tauM*(u*deriv(w,1)+v*deriv(w,2))*(du ./ dtc + (u*deriv(du,1)+v*deriv(du,2)) + deriv(dp,1)) - (w*((u-uold) ./ dtc + (u*deriv(u,1)+v*deriv(u,2))) - deriv(w,1)*p + mu*dot(grad(w), grad(u)) + tauM*(u*deriv(w,1)+v*deriv(w,2))*((u-uold) ./ dtc + (u*deriv(u,1)+v*deriv(u,2)) + deriv(p,1)))",
                        "w*(dv ./ dtc + (u*deriv(dv,1)+v*deriv(dv,2) + deriv(v,1)*du)) - deriv(w,2)*dp + mu*dot(grad(w), grad(dv)) + tauM*(u*deriv(w,1)+v*deriv(w,2))*(dv ./ dtc + (u*deriv(dv,1)+v*deriv(dv,2)) + deriv(dp,2)) - (w*((v-vold) ./ dtc + (u*deriv(v,1)+v*deriv(v,2))) - deriv(w,2)*p + mu*dot(grad(w), grad(v)) + tauM*(u*deriv(w,1)+v*deriv(w,2))*((v-vold) ./ dtc + (u*deriv(v,1)+v*deriv(v,2)) + deriv(p,2)))",
	                     "w*(deriv(du,1)+deriv(dv,2)) + tauC*(deriv(w,1)*( du ./ dtc + (u*deriv(du,1)+v*deriv(du,2)) + deriv(dp,1) ) + deriv(w,2)*( dv ./ dtc + (u*deriv(dv,1)+v*deriv(dv,2)) + deriv(dp,2) )) - (w*(deriv(u,1)+deriv(v,2)) + tauC*(deriv(w,1)*( (u-uold) ./ dtc + (u*deriv(u,1)+v*deriv(u,2)) + deriv(p,1) ) + deriv(w,2)*( (v-vold) ./ dtc + (u*deriv(v,1)+v*deriv(v,2)) + deriv(p,2) )))"])

#exportCode("NScode")
#importCode("NScode")

solve([du, dv, dp], [u, v, p, uold, vold], nonlinear=true);

# solution is stored in the variable's "values"
#using Plots
#pyplot();
#display(plot(Femshop.grid_data.allnodes[1,:], Femshop.grid_data.allnodes[2,:], u.values[:], st = :surface, reuse=false))
#display(plot(Femshop.grid_data.allnodes[1,:], Femshop.grid_data.allnodes[2,:], v.values[:], st = :surface, reuse=false))

# Dump things to the log if desired
log_dump_config();
log_dump_prob();

finalize_femshop() # Finish writing and close any files

ext: jl


text:
module TestGen
using Thrift
using Test

const testdir = dirname(@__FILE__)

function generate()
    @testset "code generation" begin
        @testset "basic" begin
            # generate code
            for (proto_name,type_name) in (("srvcctrl","ServiceControl"), ("proto_tests","ProtoTests"))
                Thrift.generate(proto_name * ".thrift")
                for suffix in ("_constants.jl", "_types.jl", ".jl")
                    @test isfile(joinpath(testdir, "gen-jl", proto_name, proto_name * suffix))
                end
                @test isfile(joinpath(testdir, "gen-jl", proto_name, type_name * ".jl"))
            end

            # add service implementations
            cp(joinpath(testdir, "srvcctrl_impl.jl"), joinpath(testdir, "gen-jl", "srvcctrl", "srvcctrl_impl.jl"); force=true)
            cp(joinpath(testdir, "proto_tests_impl.jl"), joinpath(testdir, "gen-jl", "proto_tests", "proto_tests_impl.jl"); force=true)
        end
        @testset "service method disambiguation" begin
            Thrift.generate("disambiguate_svc_method.thrift")
            touch(joinpath(testdir, "gen-jl", "disambiguate_svc_method", "disambiguate_svc_method_impl.jl"))
            include(joinpath(testdir, "gen-jl", "disambiguate_svc_method", "disambiguate_svc_method.jl"))
            @test disambiguate_svc_method.hello_to_args_Other !== disambiguate_svc_method.hello_to_args_SayHello
        end
    end
end

generate()

end # module TestGen


text:
using TcpInstruments
using Test
using TestSetExtensions
using Unitful

#using Aqua
#Aqua.test_all(TcpInstruments)

const A = u"A"

@testset ExtendedTestSet "TcpInstruments" begin
    @testset "Fake Scope" begin
        f = initialize(TcpInstruments.FakeDSOX4034A)

        data = get_data(f, 1)
        @test data isa TcpInstruments.ScopeData
        @test length(data.time) == length(data.volt)
        @test data.volt[1] isa Unitful.Voltage

        data = get_data(f, [1,2,3,4])
        @test length(data) == 4
    end

    @testset "TcpInstruments.jl" begin
        include("./emulate/test_fake_device.jl")
    end

    @testset "Util Functions" begin
        @testset "split_str_into_host_and_port" begin
            host = "192.168.1.1"
            port = 5055
            host_str = "$host:$port"
            host_out, port_out = TcpInstruments.split_str_into_host_and_port(host_str)
            @test host_out == host
            @test port_out == port

            host_out, port_out = TcpInstruments.split_str_into_host_and_port(host)
            @test host_out == host
            @test port_out == 0
        end
    end
end


text:
#############################################################################
#Copyright (c) 2016-2017 Charles Le Losq
#
#The MIT License (MIT)
#
#Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the #Software without restriction, including without limitation the rights to use, copy, #modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, #and to permit persons to whom the Software is furnished to do so, subject to the #following conditions:
#
#The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
#
#THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, #INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
#
# baseline.jl contains several "baseline" functions. It is directly dependent on functions.jl.
#
#############################################################################

"""

Allows subtracting a baseline under a x y spectrum.

Parameters
----------
x_input : ndarray
	x values.
y_input : ndarray
	y values.
bir : ndarray
	Contain the regions of interest, organised per line.
	For instance, roi = np.array([[100., 200.],[500.,600.]]) will
	define roi between 100 and 200 as well as between 500 and 600.
	Note: This is NOT used by the "als" and "arPLS" algorithms, but still is a requirement when calling the function.
	bir and method probably will become args in a futur iteration of rampy to solve this.
methods : str
	"poly": polynomial fitting, with splinesmooth the degree of the polynomial.
	"unispline": spline with the UnivariateSpline function of Scipy, splinesmooth is
				 the spline smoothing factor (assume equal weight in the present case);
	"gcvspline": spline with the gcvspl.f algorythm, really robust.
				 Spectra must have x, y, ese in it, and splinesmooth is the smoothing factor;
				 For gcvspline, if ese are not provided we assume ese = sqrt(y).
				 Requires the installation of gcvspline with a "pip install gcvspline" call prior to use;
	"exp": exponential background;
	"log": logarythmic background;
	"rubberband": rubberband baseline fitting;
	"als": automatic least square fitting following Eilers and Boelens 2005;
	"arPLS": automatic baseline fit using the algorithm from Baek et al. 2015
			 Baseline correction using asymmetrically reweighted penalized least squares smoothing, Analyst 140: 250-257.

Options
------
polynomial_order : Int
	The degree of the polynomial (0 for a constant), default = 1.
s : Float
	spline smoothing coefficient for the unispline and gcvspline algorithms.
lam : Float
	float, the lambda smoothness parameter for the ALS and ArPLS algorithms. Typical values are between 10**2 to 10**9, default = 10**5.
p : Float
	float, for the ALS algorithm, advised value between 0.001 to 0.1, default = 0.01.
ratio : float
	ratio parameter of the arPLS algorithm. default = 0.01.
niter : Int
	number of iteration of the ALS algorithm, default = 10.
p0_exp : List
	containg the starting parameter for the exp baseline fit with curve_fit. Default = [1.,1.,1.].
p0_log : List
	containg the starting parameter for the log baseline fit with curve_fit. Default = [1.,1.,1.,1.].

Returns
-------
out1 : ndarray
	Contain the corrected signal.
out2 : ndarray
	Contain the baseline.

"""
function baseline(x::Array{Float64},y::Array{Float64},roi::Array{Float64},basetype::AbstractString;polynomial_order=1, s = 1.0, lam = 10^5, p = 0.01, ratio = 0.01, niter = 10, p0_exp = [1.,1.,1.],p0_log =[1.,1.,1.])
	yout, bas = rampy[:baseline](x,y,roi,basetype,polynomial_order=polynomial_order, s=s, lam=lam, niter=niter, p0_exp=p0_exp, p0_log=p0_log)
end


text:
# Base.hash(a::ASTSampleAction) = hash(a.???) # TODO.
Base.hash(a::ASTSeedAction) = hash(a.seed)
function Base.hash(A::Vector{ASTAction})
    if length(A) > 0
        h = hash(A[1])
        for i in 2:length(A)
            h = hash(h, hash(A[i]))
        end
    else
        h = hash([])
    end
    return h
end

Base.hash(s::ASTState) = hash(s.t_index, s.parent, s.action)

function Base.hash(t_index::Int64, parent::Union{Nothing, ASTState}, action::Union{Nothing, ASTAction})
    h = hash(t_index)
    h = hash(h, hash(isnothing(parent) ? nothing : parent.hash))
    h = hash(h, hash(action))
    return h
end

# Base.:(==)(w::ASTSampleAction,v::ASTSampleAction) = w.??? == v.??? # TODO.
Base.:(==)(w::ASTSeedAction,v::ASTSeedAction) = w.seed == v.seed
Base.:(==)(w::ASTState,v::ASTState) = hash(w) == hash(v)
# Base.isequal(w::ASTSampleAction,v::ASTSampleAction) = isequal(w.???, v.???) # TODO.
Base.isequal(w::ASTSeedAction,v::ASTSeedAction) = isequal(w.seed, v.seed)
Base.isequal(w::ASTState,v::ASTState) = hash(w) == hash(v)


text:
doc"""
yield result of `f(x)`
"""
function map_(f::Function; debug=false)
    x = produce()
    while x != :done
        y = _map(x, f)
        if debug
            println(x, " -> ", y)
        end
        x = produce(y)
    end
    return :done
end

doc"""
yield x if f(x) is true
"""
function filter_(f::Function; debug=false)
    x = produce()
    while x != :done
        y = _filter(x, f)
        if debug
            println(x, " -> ", y)
        end
        x = produce(y)
    end
    return :done
end

function reduce{T <: Functor}(f::Function, x0::T; debug=false)
    x = x0
    y = produce()
    while true
        if debug
            print(x, ", ", y)
        end
        x = _reduce(x, y, f)
        if debug
            println(" -> ", x)
        end
        y = produce(x)
    end
end

doc"""
bridge makes a task bridging two functions or tasks `f` and `g` executing sequentially.
"""
function bridge(f::Task, g::Task)
    x = produce()
    while true
        y = consume(f, x)
        z = consume(g, y)
        x = produce(z)
    end
end

function bridge(f::Task, gs::Vector{Task})
    x = produce()
    while true
        y = consume(f, x)
        zs = [consume(g, y) for g in gs]
        x = produce(zs)
    end
end

function bridge(fs::Vector{Task}, g::Task)
    x = produce()
    while true
        ys = [consume(f, x) for f in fs]
        z = consume(g, ys)
        x = produce(z)
    end
end

function bridge(f::Function, g::Task)
    a = init_task(map_, f)
    return bridge(a, g)
end

function bridge(f::Task, g::Function)
    b = init_task(map_, g)
    return bridge(f, b)
end

function bridge(f::Task, gs::Vector{Function})
    bs = [init_task(map_, g) for g in gs]
    return bridge(f, bs)
end

function bridge(f::Function, g::Function)
    a = init_task(map_, f)
    b = init_task(map_, g)
    return bridge(a, b)
end

function bridge(f::Function, gs::Vector{Function})
    a = init_task(map_, f)
    bs = [init_task(map_, g) for g in gs]
    return bridge(a, bs)
end

function bridge(fs::Vector{Function}, g::Function)
    as = [init_task(map_, f) for f in fs]
    b = init_task(map_, g)
    return bridge(as, b)
end

function bridge(f::Task, gs::Vararg)
    g, g_tail = gs[1], gs[2:end]
    t = init_task(bridge, f, g)
    return bridge(t, g_tail...)
end

function bridge(f::Function, gs::Vararg)
    return bridge(init_task(map_, f), gs...)
end


text:
using DataDeps
using DataDepsGenerators
using Test
using ReferenceTests

"""
URLS tend to go out of date fast on github (as we always generate using the latest commit)
so strip them from the recorded reference
"""
function discard_urls(code)
    replace(code, r"http.*?([\s\"\)])" => s"<URL>\1")
end

@testset "fivethirtyeight (folders in 1 repo)" begin
    @testset "538 March" begin
        registration_code = generate(GitHub(), "https://github.com/fivethirtyeight/data/tree/master/march-madness-predictions-2015")

        @testset "Integration Test" begin
            eval(Meta.parse(registration_code))
            @test length(collect(readdir(datadep"March Madness Predictions"))) > 0
        end

        @test_reference "references/538 march.txt" discard_urls(registration_code)
    end

    @test_reference "references/538 steak.txt" discard_urls(replace(
        generate(GitHub(), "https://github.com/fivethirtyeight/data/tree/master/steak-survey"),
        r"Article.*"s=> "")) # Delete most of the text harvested from top level readme, as that will go out of date fast

    @test_reference "references/538 college.txt" discard_urls(generate(GitHub(), "fivethirtyeight/data/tree/master/college-majors"))
end

@testset "BuzzFeedNews (whole repos)" begin
    @testset "Pres Camp" begin
        registration_code = generate(GitHub(), "https://github.com/BuzzFeedNews/presidential-campaign-contributions")

		@testset "Integration Test" begin
			eval(Meta.parse(registration_code)) # evaluate the new code
			@test length(collect(readdir(datadep"Presidential Campaign Contributions"))) > 0
		end

        @test_reference "references/buzzfeed pres-camp.txt" discard_urls(registration_code) # See if new code is same as old
	end

    @test_reference "references/buzzfeed primates.txt" discard_urls(generate(GitHub(), "BuzzFeedNews/2015-07-primates"))
end


text:
import ..Continuations: find_errors

function find_errors(solver::Codim1Solver)
    super = as(solver, ContinuationSolver)
    errors = find_errors(super)
    # TODO: do some Codim1Solver-specific checks
    return errors
end


text:
function LaneReversal(N, link_dict, g, fcoeffs,  carLength=4.572, method="branch_and_bound")

	sort_ld = sorted(link_dict)
	sort_od = sorted(g)
	nLinks = len(G.edges())
	nNodes = len(G.nodes())
	nOD = len(g)
	nPoly = len(fcoeffs)

	m = Model("LaneReversal")
	m.setParam( 'OutputFlag', False )

	x_t = [m.addVar(lb=0, name="x_{"+str(a)+"}") for a in sort_ld]
	x = []
	d = []
	g_id = {}
	w_id = 0
	for (s,t) in sort_od:
		# define flow variables
		x.append([m.addVar(lb=0, name="x_{"+str(a)+","+str(w_id)+"}") for a in sort_ld])
		# set demands
		d.append(np.zeros((nNodes,1)))
		d[w_id][s-1] = g[(s,t)]
		d[w_id][t-1] = -g[(s,t)]
		# add flow conservation and demand comply
		m.update()
		constr = np.dot(N,x[w_id])
		for n in range(nNodes):
			m.addConstr(constr[n] == d[w_id][n][0])
		g_id[w_id] = (s,t)
		w_id +=1

	m.update()
	# second constrains (sum of od flows is overall flow)
	for a in range(nLinks):
		m.addConstr(x_t[a] == sum([x[w][a] for w in g_id.keys()]))
	m.update()

	# define lanes variables
	l = {}
	for a in sort_ld:
		u = link_dict[a][0]
		v = link_dict[a][1]
		l[u,v] = m.addVar(vtype=GRB.INTEGER, lb=0, name="l_{"+str(u)+","+str(v)+"}")
		l[v,u] = m.addVar(vtype=GRB.INTEGER, lb=0, name="l_{"+str(v)+","+str(u)+"}")
		m.addConstr(l[u,v]+l[v,u] <= G[u][v]["lanes"]+G[v][u]["lanes"])

	obj = 0
	for a in sort_ld:
		u = link_dict[a][0]
		v = link_dict[a][1]
		obj += x_t[a] * G[u][v]['t_0'] + (sum([fcoeffs[i]*(x_t[a]^i/(l[u,v]*G[u][v]["length"]*carLength/2)^i)for i in range(nPoly)]))

end

max_stars_repo_name: salomonw/tnet

max_issues_repo_name: salomonw/tnet

max_forks_repo_name: salomonw/tnet


text:
import Base: rand,
             ∈

export LineSegment,
       halfspace_left, halfspace_right,
       constraints_list

"""
    LineSegment{N<:Real} <: AbstractZonotope{N}

Type that represents a line segment in 2D between two points ``p`` and ``q``.

### Fields

- `p` -- first point
- `q` -- second point

### Examples

A line segment along the ``x = y`` diagonal:

```jldoctest linesegment_constructor
julia> s = LineSegment([0., 0], [1., 1.])
LineSegment{Float64}([0.0, 0.0], [1.0, 1.0])
julia> dim(s)
2
```

Use `plot(s)` to plot the extreme points of `s` and the line segment joining
them. Membership test is computed with ∈ (`in`):

```jldoctest linesegment_constructor
julia> [0., 0] ∈ s && [.25, .25] ∈ s && [1., 1] ∈ s && [.5, .25] ∉ s
true
```

We can check the intersection with another line segment, and optionally compute
a witness (which is just the common point in this case):

```jldoctest linesegment_constructor
julia> sn = LineSegment([1., 0], [0., 1.])
LineSegment{Float64}([1.0, 0.0], [0.0, 1.0])
julia> isempty(s ∩ sn)
false
julia> is_intersection_empty(s, sn, true)
(false, [0.5, 0.5])
```
"""
struct LineSegment{N<:Real} <: AbstractZonotope{N}
    p::AbstractVector{N}
    q::AbstractVector{N}

    # default constructor with length constraint
    function LineSegment{N}(p::AbstractVector{N},
                            q::AbstractVector{N}) where {N<:Real}
        @assert length(p) == length(q) == 2 "points for line segments must " *
            "be two-dimensional"
        return new{N}(p, q)
    end
end

isoperationtype(::Type{<:LineSegment}) = false

# convenience constructor without type parameter
LineSegment(p::AbstractVector{N}, q::AbstractVector{N}) where {N<:Real} =
    LineSegment{N}(p, q)

# --- LazySet interface functions ---

"""
    dim(L::LineSegment)::Int

Return the ambient dimension of a line segment.

### Input

- `L` -- line segment

### Output

The ambient dimension of the line segment, which is 2.
"""
function dim(L::LineSegment)::Int
    return 2
end

"""
    σ(d::AbstractVector{N}, L::LineSegment{N}) where {N<:Real}

Return the support vector of a line segment in a given direction.

### Input

- `d` -- direction
- `L` -- line segment

### Output

The support vector in the given direction.

### Algorithm

If the angle between the vector ``q - p`` and ``d`` is bigger than 90° and less
than 270° (measured in counter-clockwise order), the result is ``p``, otherwise
it is ``q``.
If the angle is exactly 90° or 270°, or if the direction has norm zero, this
implementation returns ``q``.
"""
function σ(d::AbstractVector{N}, L::LineSegment{N}) where {N<:Real}
    return sign(dot(L.q - L.p, d)) >= 0 ? L.q : L.p
end

"""
    an_element(L::LineSegment{N}) where {N<:Real}

Return some element of a line segment.

### Input

- `L` -- line segment

### Output

The first vertex of the line segment.
"""
function an_element(L::LineSegment{N}) where {N<:Real}
    return L.p
end

"""
    ∈(x::AbstractVector{N}, L::LineSegment{N})::Bool where {N<:Real}

Check whether a given point is contained in a line segment.

### Input

- `x` -- point/vector
- `L` -- line segment

### Output

`true` iff ``x ∈ L``.

### Algorithm

Let ``L = (p, q)`` be the line segment with extremes ``p`` and ``q``, and let
``x`` be the given point.

1. A necessary condition for ``x ∈ (p, q)`` is that the three points are aligned,
   thus their cross product should be zero.
2. It remains to check that ``x`` belongs to the box approximation of ``L``.
   This amounts to comparing each coordinate with those of the extremes ``p``
   and ``q``.

### Notes

The algorithm is inspired from [here](https://stackoverflow.com/a/328110).
"""
function ∈(x::AbstractVector{N}, L::LineSegment{N})::Bool where {N<:Real}
    @assert length(x) == dim(L)
    # check if the point is on the line through the line segment
    if (x[2] - L.p[2]) * (L.q[1] - L.p[1]) -
            (x[1] - L.p[1]) * (L.q[2] - L.p[2]) != 0
        return false
    end
    # check if the point is inside the box approximation of the line segment
    return min(L.p[1], L.q[1]) <= x[1] <= max(L.p[1], L.q[1]) &&
           min(L.p[2], L.q[2]) <= x[2] <= max(L.p[2], L.q[2])
end

# --- AbstractCentrallySymmetric interface functions ---

"""
    center(L::LineSegment{N})::Vector{N} where {N<:Real}

Return the center of a line segment.

### Input

- `L` -- line segment

### Output

The center of the line segment.
"""
function center(L::LineSegment{N})::Vector{N} where {N<:Real}
    return L.p + (L.q - L.p) / 2
end

# --- AbstractZonotope interface functions ---

"""
   genmat(L::LineSegment)

Return the generator matrix of a line segment.

### Input

- `L` -- line segment

### Output

A matrix with a single column representing the generator of `L`.
"""
function genmat(L::LineSegment)
    return genmat_fallback(L)
end

"""
    generators(L::LineSegment{N}) where {N<:Real}

Return an iterator over the (single) generator of a line segment.

### Input

- `L` -- line segment

### Output

A one-element iterator over the generator of `L`.
"""
function generators(L::LineSegment{N}) where {N<:Real}
    if L.p == L.q
        # degenerate line segment has no generators
        return EmptyGeneratorIterator{N}()
    end
    return [L.p - center(L)]
end

# --- AbstractPolytope interface functions ---

"""
    vertices_list(L::LineSegment{N}
                 )::Vector{<:AbstractVector{N}} where {N<:Real}

Return the list of vertices of a line segment.

### Input

- `L` -- line segment

### Output

The list of end points of the line segment.
"""
function vertices_list(L::LineSegment{N}
                      )::Vector{<:AbstractVector{N}} where {N<:Real}
    return [L.p, L.q]
end

# --- LazySet interface functions ---

"""
    rand(::Type{LineSegment}; [N]::Type{<:Real}=Float64, [dim]::Int=2,
         [rng]::AbstractRNG=GLOBAL_RNG, [seed]::Union{Int, Nothing}=nothing
        )::LineSegment{N}

Create a random line segment.

### Input

- `LineSegment` -- type for dispatch
- `N`           -- (optional, default: `Float64`) numeric type
- `dim`         -- (optional, default: 2) dimension
- `rng`         -- (optional, default: `GLOBAL_RNG`) random number generator
- `seed`        -- (optional, default: `nothing`) seed for reseeding

### Output

A random line segment.

### Algorithm

All numbers are normally distributed with mean 0 and standard deviation 1.
"""
function rand(::Type{LineSegment};
              N::Type{<:Real}=Float64,
              dim::Int=2,
              rng::AbstractRNG=GLOBAL_RNG,
              seed::Union{Int, Nothing}=nothing
             )::LineSegment{N}
    @assert dim == 2 "cannot create a random LineSegment of dimension $dim"
    rng = reseed(rng, seed)
    p = randn(rng, N, dim)
    q = randn(rng, N, dim)
    return LineSegment(p, q)
end

# --- LineSegment functions ---

"""
    halfspace_left(L::LineSegment)

Return a half-space describing the 'left' of a two-dimensional line segment
through two points.

### Input

 - `L` -- line segment

### Output

The half-space whose boundary goes through the two points `p` and `q` and which
describes the left-hand side of the directed line segment `pq`.
"""
halfspace_left(L::LineSegment) = halfspace_left(L.p, L.q)

"""
    halfspace_right(L::LineSegment)

Return a half-space describing the 'right' of a two-dimensional line segment
through two points.

### Input

 - `L` -- line segment

### Output

The half-space whose boundary goes through the two points `p` and `q` and which
describes the right-hand side of the directed line segment `pq`.
"""
halfspace_right(L::LineSegment) = halfspace_right(L.p, L.q)

"""
    constraints_list(L::LineSegment{N}) where {N<:Real}

Return the list of constraints defining a line segment in 2D.

### Input

- `L` -- line segment

### Output

A vector of constraints that define the line segment.

### Algorithm

``L`` is defined by 4 constraints. In this algorithm, the first two constraints
are returned by ``halfspace_right`` and ``halfspace_left``, and the other two
are obtained by considering the vector normal to the line segment that passes
through each opposite vertex.

### Notes

This function returns a vector of halfspaces. It does not return equality
constraints.
"""
function constraints_list(L::LineSegment{N}) where {N<:Real}
    p, q = L.p, L.q
    d = [(p[2]-q[2]), (q[1]-p[1])]
    return [halfspace_left(L), halfspace_right(L),
            halfspace_right(p, p + d), halfspace_left(q, q + d)]
end

"""
    translate(L::LineSegment{N}, v::AbstractVector{N}) where {N<:Real}

Translate (i.e., shift) a line segment by a given vector.

### Input

- `L` -- line segment
- `v` -- translation vector

### Output

A translated line segment.

### Algorithm

We add the vector to both defining points of the line segment.
"""
function translate(L::LineSegment{N}, v::AbstractVector{N}) where {N<:Real}
    @assert length(v) == dim(L) "cannot translate a $(dim(L))-dimensional " *
                                "set by a $(length(v))-dimensional vector"
    return LineSegment(L.p + v, L.q + v)
end

"""
    plot_recipe(L::LineSegment{N}, [ε]::N=zero(N)) where {N<:Real}

Convert a line segment to a pair `(x, y)` of points for plotting.

### Input

- `L` -- line segment
- `ε` -- (optional, default: `0`) ignored, used for dispatch

### Output

A pair `(x, y)` of two points that can be plotted.
"""
function plot_recipe(L::LineSegment{N}, ε::N=zero(N)) where {N<:Real}
    return [L.p[1], L.q[1]], [L.p[2], L.q[2]]
end


text:
struct ByteSet <: AbstractSet{UInt8}
    data::NTuple{4, UInt64}
    ByteSet(x::NTuple{4, UInt64}) = new(x)
end

ByteSet() = ByteSet((UInt64(0), UInt64(0), UInt64(0), UInt64(0)))
Base.length(s::ByteSet) = mapreduce(count_ones, +, s.data)
Base.isempty(s::ByteSet) = s === ByteSet()

function ByteSet(it)
    a = b = c = d = UInt64(0)
    for i in it
        vi = convert(UInt8, i)
        if vi < 0x40
            a |= UInt(1) << ((vi - 0x00) & 0x3f)
        elseif vi < 0x80
            b |= UInt(1) << ((vi - 0x40) & 0x3f)
        elseif vi < 0xc0
            c |= UInt(1) << ((vi - 0x80) & 0x3f)
        else
            d |= UInt(1) << ((vi - 0xc0) & 0x3f)
        end
    end
    ByteSet((a, b, c, d))
end

function Base.minimum(s::ByteSet)
    isempty(s) && Base._empty_reduce_error()
    (iterate(s)::Tuple{UInt8, Any})[1]
end

function Base.maximum(s::ByteSet)
    isempty(s) && Base._empty_reduce_error()
    offset = 0x03 * UInt8(64)
    for i in 0:3
        @inbounds bits = s.data[4 - i]
        lz = leading_zeros(bits)
        if lz != 64
            return offset + UInt8(64) - UInt8(lz) - 0x01
        end
        offset -= UInt8(64)
    end
end

function Base.in(byte::UInt8, s::ByteSet)
    i, o = divrem(byte, UInt8(64))
    @inbounds !(iszero(s.data[i & 0x03 + 0x01] >>> (o & 0x3f) & UInt(1)))
end

function Base.iterate(s::ByteSet, state=UInt(0))
    ioffset, offset = divrem(state, UInt(64))
    ioffset == 4 && return nothing
    bits = @inbounds s.data[ioffset & 0x03 + 0x01] >>> offset
    tz = trailing_zeros(bits)
    while tz == 64
        ioffset += 1
        offset = UInt(0)
        ioffset == 4 && return nothing
        bits = s.data[ioffset & 0x03 + 0x01]
        tz = trailing_zeros(bits)
    end
    result = (64 * ioffset + offset + tz) % UInt8
    result, UInt(result) + UInt(1)
end

function Base.:~(s::ByteSet)
    a, b, c, d = s.data
    ByteSet((~a, ~b, ~c, ~d))
end

is_contiguous(s::ByteSet) = isempty(s) || (maximum(s)::UInt8 - minimum(s)::UInt8 + 1 == length(s))


text:
# test functions for the 5taxon networks
# functions used in tests_5taxon.jl
# Claudia September 2014
# changed to new reparametrizations: bad diamondI,II and bad triangle I,II
# Claudia November 2014

# Case C bad triangle II
function testCaseC(net::HybridNetwork)
    n=searchHybridNode(net);
    node = n[1];
    node.k != 3 ? error("k diff than 3") : nothing
    node.isVeryBadTriangle ? nothing : error("does not know it is very bad triangle")
    node.isExtBadTriangle ? error("thinks it is extremely bad triangle") : nothing
    net.hasVeryBadTriangle ? nothing : error("net does not know it has very bad triangle")
    net.numBad == 0 ? nothing : error("net.numBad should be 0")
    net.numHybrids != 1 ? error("should have 1 hybrid, but net.numHybrids is $(net.numHybrids): $([n.number for n in net.hybrid])") : nothing
end

# Case F bad diamond
function testCaseF(net::HybridNetwork)
    net.visited = [e.istIdentifiable for e in net.edge];
    n = searchHybridNode(net);
    node = n[1];
    node.k != 4 ? error("k diff than 4") : nothing
    edge5 = getIndexEdge(5,net);
    edge9 = getIndexEdge(9,net);
    edge11 = getIndexEdge(11,net);
    edge12 = getIndexEdge(12,net);
    edge15 = getIndexEdge(15,net);
    edge14 = getIndexEdge(14,net);
    node1 = getIndexNode(1,net);
    node5 = getIndexNode(5,net);
    node11 = getIndexNode(11,net);
    node12 = getIndexNode(12,net);
    (net.edge[edge5].inCycle != node.number || net.edge[edge11].inCycle != node.number || net.edge[edge12].inCycle != node.number || net.edge[edge15].inCycle != node.number ) ? error("edges not correctly in cycle") : nothing
    (net.node[node1].inCycle  != node.number || net.node[node5].inCycle  != node.number || net.node[node11].inCycle  != node.number || net.node[node12].inCycle  != node.number) ? error("nodes 1,5,11,12 not correctly in cycle") : nothing
    !node.isBadDiamondI ? error("does not know it is bad diamondI") : nothing
    node.isBadDiamondII ? error("thinks it is bad diamond II") : nothing
    net.edge[edge14].containRoot ? error("edge can contain root") : nothing
    (!net.edge[edge11].hybrid || !net.edge[edge11].isMajor) ? error("edge 11 is not hybrid or major") : nothing
    net.node[node12].gammaz != net.edge[edge15].gamma*net.edge[edge12].z ? error("node 12 gammaz not correctly calculated") : nothing
    net.node[1].gammaz != net.edge[edge11].gamma*net.edge[edge5].z ? error("node 11 gammaz not correctly calculated") : nothing
    !net.edge[edge9].istIdentifiable ? error("edge9 not identifiable") : nothing
    net.visited[edge9] = false;
    !all([!id for id in net.visited]) ? error("edges not identifiable as identifiable") : nothing
    net.numHybrids != 1 ? error("should have 1 hybrid, but net.numHybrids is $(net.numHybrids): $([n.number for n in net.hybrid])") : nothing
    net.numBad == 1 ? nothing : error("net.numBad should be 1")
end

# Case G
function testCaseG(net::HybridNetwork)
    net.visited = [e.istIdentifiable for e in net.edge];
    n = searchHybridNode(net);
    node = n[1];
    node.k != 4 ? error("k diff than 4") : nothing
    edge9 = getIndexEdge(9,net);
    edge12 = getIndexEdge(12,net);
    edge15 = getIndexEdge(15,net);
    edge5 = getIndexEdge(5,net);
    edge13 = getIndexEdge(13,net);
    edge14 = getIndexEdge(14,net);
    node9 = getIndexNode(9,net);
    node5 = getIndexNode(5,net);
    node11 = getIndexNode(11,net);
    node12 = getIndexNode(12,net);
    (net.edge[edge9].inCycle != node.number || net.edge[edge12].inCycle != node.number || net.edge[edge13].inCycle != node.number || net.edge[edge15].inCycle != node.number ) ? error("edges not correctly in cycle") : nothing
    (net.node[node9].inCycle  != node.number || net.node[node5].inCycle  != node.number || net.node[node11].inCycle  != node.number || net.node[node12].inCycle  != node.number) ? error("nodes not correctly in cycle") : nothing
    (node.isBadDiamondI || node.isBadDiamondII) ? error("thinks it is bad diamond") : nothing
    net.edge[edge14].containRoot ? error("14 can contain root") : nothing
    (!net.edge[edge12].hybrid || !net.edge[edge12].isMajor) ? error("edge 12 is not hybrid or major") : nothing
    net.node[node12].gammaz != -1 ? error("node 12 gammaz should be -1") : nothing
    (!net.edge[edge9].istIdentifiable || !net.edge[edge5].istIdentifiable || !net.edge[edge13].istIdentifiable) ? error("edge9,5,13not identifiable") : nothing
    net.visited[edge9] = false;
    net.visited[edge5] = false;
    net.visited[edge13] = false;
    !all([!id for id in net.visited]) ? error("edges not identifiable as identifiable") : nothing
    net.numHybrids != 1 ? error("should have 1 hybrid, but net.numHybrids is $(net.numHybrids): $([n.number for n in net.hybrid])") : nothing
end

# Case H
function testCaseH(net::HybridNetwork)
    net.visited = [e.istIdentifiable for e in net.edge];
    n = searchHybridNode(net);
    node = n[1];
    node.k != 4 ? error("k diff than 4") : nothing
    edge9 = getIndexEdge(9,net);
    edge14 = getIndexEdge(14,net);
    edge15 = getIndexEdge(15,net);
    edge13 = getIndexEdge(13,net);
    edge5 = getIndexEdge(5,net);
    edge8 = getIndexEdge(8,net);
    node9 = getIndexNode(9,net);
    node5 = getIndexNode(5,net);
    node11 = getIndexNode(11,net);
    node12 = getIndexNode(12,net);
    (net.edge[edge9].inCycle != node.number || net.edge[edge5].inCycle != node.number || net.edge[edge14].inCycle != node.number || net.edge[edge15].inCycle != node.number ) ? error("edges not correctly in cycle") : nothing
    (net.node[node9].inCycle  != node.number || net.node[node5].inCycle  != node.number || net.node[node11].inCycle  != node.number || net.node[node12].inCycle  != node.number) ? error("nodes 9,5,11,12 not correctly in cycle") : nothing
    (node.isBadDiamondI || node.isBadDiamondII )? error("thinks it is bad diamond") : nothing
    net.edge[edge8].containRoot ? error("8 can contain root") : nothing
    (!net.edge[edge14].hybrid || !net.edge[edge14].isMajor) ? error("edge 14 is not hybrid or major") : nothing
    net.node[node12].gammaz != -1 ? error("node 12 gammaz should be -1") : nothing
    (!net.edge[edge9].istIdentifiable || !net.edge[edge5].istIdentifiable || !net.edge[edge13].istIdentifiable) ? error("edge9,5,13not identifiable") : nothing
    net.visited[edge9] = false;
    net.visited[edge5] = false;
    net.visited[edge13] = false;
    !all([!id for id in net.visited]) ? error("edges not identifiable as identifiable") : nothing
    net.numHybrids != 1 ? error("should have 1 hybrid, but net.numHybrids is $(net.numHybrids): $([n.number for n in net.hybrid])") : nothing
end

# Case J
function testCaseJ(net::HybridNetwork)
    net.visited = [e.istIdentifiable for e in net.edge];
    n = searchHybridNode(net);
    node = n[1];
    node.k != 5 ? error("k diff than 5") : nothing
    edge9 = getIndexEdge(9,net);
    edge10 = getIndexEdge(10,net);
    edge5 = getIndexEdge(5,net);
    edge15 = getIndexEdge(15,net);
    edge6 = getIndexEdge(6,net);
    edge14 = getIndexEdge(14,net);
    node9 = getIndexNode(9,net);
    node1 = getIndexNode(1,net);
    node5 = getIndexNode(5,net);
    node11 = getIndexNode(11,net);
    node12 = getIndexNode(12,net);
    (net.edge[edge9].inCycle != node.number || net.edge[edge10].inCycle != node.number || net.edge[edge6].inCycle != node.number || net.edge[edge15].inCycle != node.number || net.edge[edge5].inCycle != node.number) ? error("edges not correctly in cycle") : nothing
    (net.node[node9].inCycle  != node.number || net.node[node5].inCycle  != node.number || net.node[node11].inCycle  != node.number || net.node[node12].inCycle  != node.number || net.node[node1].inCycle) != node.number ? error("nodes 9,5,11,12,1 not correctly in cycle") : nothing
    net.edge[edge14].containRoot ? error("14 can contain root") : nothing
    (!net.edge[edge6].hybrid || !net.edge[edge6].isMajor) ? error("edge 6 is not hybrid or major") : nothing
    net.node[node12].gammaz != -1 ? error("node 12 gammaz should be -1") : nothing
    (!net.edge[edge9].istIdentifiable || !net.edge[edge5].istIdentifiable || !net.edge[edge10].istIdentifiable) ? error("edge9,5,10not identifiable") : nothing
    net.visited[edge9] = false;
    net.visited[edge5] = false;
    net.visited[edge10] = false;
    !all([!id for id in net.visited]) ? error("edges not identifiable as identifiable") : nothing
    net.numHybrids != 1 ? error("should have 1 hybrid, but net.numHybrids is $(net.numHybrids): $([n.number for n in net.hybrid])") : nothing
end

# Case D bad triangle I
function testCaseD(net::HybridNetwork)
    net.visited = [e.istIdentifiable for e in net.edge];
    n = searchHybridNode(net);
    node = n[1];
    node.k != 3 ? error("k diff than 3") : nothing
    node.isVeryBadTriangle ? nothing : error("does not know it is very bad triangle")
    node.isExtBadTriangle ? error("thinks it is extremely bad triangle") : nothing
    net.hasVeryBadTriangle ? nothing : error("net does not know it has very bad triangle")
    net.numBad == 0 ? nothing : error("net.numBad should be 0")
    net.numHybrids != 1 ? error("should have 1 hybrid, but net.numHybrids is $(net.numHybrids): $([n.number for n in net.hybrid])") : nothing
end

# Case E bad triangle I
function testCaseE(net::HybridNetwork)
    net.visited = [e.istIdentifiable for e in net.edge];
    n = searchHybridNode(net);
    node = n[1];
    node.k != 3 ? error("k diff than 3") : nothing
    node.isVeryBadTriangle ? nothing : error("does not know it is very bad triangle")
    node.isExtBadTriangle ? error("thinks it is extremely bad triangle") : nothing
    net.hasVeryBadTriangle ? nothing : error("net does not know it has very bad triangle")
    net.numBad == 0 ? nothing : error("net.numBad should be 0")
    net.numHybrids != 1 ? error("should have 1 hybrid, but net.numHybrids is $(net.numHybrids): $([n.number for n in net.hybrid])") : nothing
end

# Case I bad diamond II
function testCaseI(net::HybridNetwork)
    n = searchHybridNode(net);
    node = n[1];
    node.isBadDiamondII ? nothing : error("does not know it is bad diamond II")
    !node.isBadDiamondI ? nothing : error("thinks it is bad diamond I")
    net.numHybrids != 1 ? error("should have 1 hybrid, but net.numHybrids is $(net.numHybrids): $([n.number for n in net.hybrid])") : nothing
    net.visited = [e.istIdentifiable for e in net.edge];
    edge14 = getIndexEdge(14,net);
    edge10 = getIndexEdge(10,net);
    edge5 = getIndexEdge(5,net);
    edge15 = getIndexEdge(15,net);
    edge9 = getIndexEdge(9,net);
    edge11 = getIndexEdge(11,net);
    edge8 = getIndexEdge(8,net);
    node1 = getIndexNode(1,net);
    node12 = getIndexNode(12,net);
    node5 = getIndexNode(5,net);
    node11 = getIndexNode(11,net);
    (net.edge[edge5].inCycle != node.number || net.edge[edge9].inCycle != node.number || net.edge[edge11].inCycle != node.number || net.edge[edge15].inCycle != node.number ) ? error("edges not correctly in cycle") : nothing
    (net.node[node1].inCycle  != node.number || net.node[node12].inCycle  != node.number || net.node[node5].inCycle  != node.number || net.node[node11].inCycle  != node.number) ? error("nodes 1,5,11,12 not correctly in cycle") : nothing
    (net.edge[edge14].containRoot || net.edge[edge10].containRoot || net.edge[edge8].containRoot) ? error("edges can contain root and shouldn't") : nothing
    (!net.edge[edge9].hybrid || !net.edge[edge9].isMajor) ? error("edge 4 is not hybrid or major") : nothing
    net.edge[edge14].length != 0 ? error("edges should have length 0") : nothing
    net.edge[edge14].istIdentifiable ? error("edge14 identifiable and should not") : nothing
    net.visited[edge9] = false;
    net.visited[edge5] = false;
    net.visited[edge11] = false;
    net.visited[edge15] = false;
    !all([!id for id in net.visited]) ? error("edges not identifiable as identifiable") : nothing
    net.numBad == 0 ? nothing : error("net.numBad should be 0")
end

# tree example
# warning: if added hybridization (bad diamond/triangle) and then delete,
#          original edge lengths cannot be recovered
function testTree(net::HybridNetwork)
    !all([!e.hybrid for e in net.edge]) ? error("some edge is still hybrid") : nothing
    !all([!e.hybrid for e in net.node]) ? error("some node is still hybrid") : nothing
    !all([!e.hasHybEdge for e in net.node]) ? error("some node has hybrid edge") : nothing
    !all([e.isMajor for e in net.edge]) ? error("some edge is not major") : nothing
    !all([e.containRoot for e in net.edge]) ? error("some edge cannot contain root") : nothing
    edge9 = getIndexEdge(9,net);
    edge5 = getIndexEdge(5,net);
    (!net.edge[edge9].istIdentifiable || !net.edge[edge5].istIdentifiable) ? error("edge9,5 not identifiable") : nothing
    net.visited = [e.istIdentifiable for e in net.edge];
    net.visited[edge9] = false;
    net.visited[edge5] = false;
    !all([!id for id in net.visited]) ? error("edges not identifiable as identifiable") : nothing
    isempty(net.hybrid) ? nothing : error("something in net.hybrid")
    net.numHybrids != 0 ? error("should have 0 hybrid, but net.numHybrids is $(net.numHybrids)") : nothing
    #edge11 = getIndexEdge(11,net);
    #edge12 = getIndexEdge(12,net);
    #net.edge[edge11].length != 1.5 ? error("edge length for 11 is wrong") : nothing
    #net.edge[edge12].length != 0.2 ? error("edge length for 12 is wrong") : nothing
end


text:
using FastGaussQuadrature

function getCoefficientsPGLRK(s::Int)
    @assert s ≥ 2

    # order
    o = 2s

    # obtain Gauss-Legendre nodes and weights
    gl = gausslegendre(s)

    # scale from [-1,+1] to [0,1]
    c = (gl[1] .+ 1) ./ 2
    b = gl[2] ./ 2

    T = eltype(c)

    local ξᵢ::T

    leg_basis = LegendreBasis(T,s)

    a = zeros(T, s, s)
    t = zeros(T, s, s)
    P = zeros(T, s, s)
    X = zeros(T, s, s)
    W = zeros(T, s, s)

    for j in 1:s
        for i in 1:s
            P[i,j] = evaluate(leg_basis, j, c[i])
        end
    end

    Q = inv(P)

    X[1,1] = 1/2

    for i in 1:s-1
        ξᵢ = one(T)/sqrt(convert(T, 4i^2-1))/2
        X[i+1,i] = +ξᵢ
        X[i,i+1] = -ξᵢ
    end

    W[s, s-1] = +1
    W[s-1, s] = -1

    simd_mult!(t, X, Q)
    simd_mult!(a, P, t)

    CoefficientsPGLRK(Symbol("pglrk", s), o, a, b, c, P, X, W)
end


text:
using SDtoolbox, Test

@testset "SDtoolbox.jl" begin
    include("CJspeed.jl")
    @test rel_err < 1e-4
    @test ref_err < 1e-2

    include("cv.jl")
    @test err > 0.99

    include("znd.jl")
    @test err_T > 0.99 && err_therm > 0.99

    include("cell_size.jl")
    @test 9e-3 <= λ <= 10e-3
end


text:
__precompile__()
# WOS Fractional Laplacian
# The algorithms are implemented here.
########################
module KOS_FIELD_SOLVE
#######################
using Distributions
using Interpolations
#
using KOS
using KOSHELPER
################################
export init, get_solution, wos_ev
################################
function __init__()
end
##########################
function ann_ml_params(p, dict, G, no_samps, report_flag)
  # Analyse working of ML
  # compute variance rates, cputimes, total runtimes
  # and fastest start level
  no_levels=dict["no_levels"];  l0=dict["l0"]
  #
  if report_flag
    println("====\nAnalyse multilevel params with lmax=", no_levels)
  end
  #############################################################
  # get variances and costs
  ##############################################################
  cvariance=zeros(no_levels);  ccost=ones(no_levels)
  variance=zeros(no_levels);  cost=ones(no_levels)
  #
  # coarse level
  for ell=l0:no_levels
    coarse_level=intersect(G.int,1:G.get_level_ind(ell)[end])
    no_samps=Int64(ceil(dict["no_samps"]))
    # initialise
    tic()
    v,no_samps=uncoupled_factor_par2(no_samps, G, p, dict,coarse_level)
    cvariance[ell]=G.sq_int(sqrt.(v),l0)
    ccost[ell]    =toq()/(no_samps)
  end
  # coupled levels
  for ell=(l0+1):(no_levels)
    finest_level_int=intersect(G.int,G.get_level_ind(ell))
    correction_level=intersect(G.int,1:G.get_level_ind(ell)[end])
    no_samps=Int64(ceil(dict["no_samps"]))
    tic()
    v,no_samps=coupled_factor_par2(no_samps, G, p, dict,
       correction_level, finest_level_int)
    variance[ell]=G.sq_int(sqrt.(v),ell)
    cost[ell]=toq()/no_samps
  end # loop over levels
  #
  if report_flag
    println("===Completed ml params.")
  end
  ##################################################################
  ## get ratios for variance
  beta_ratio=log.(variance[l0+1:end-1]./variance[l0+2:end])/log(2)
  if report_flag
    println("beta ML param (base 2) ",beta_ratio)
    coord_print(1:length(variance),variance)
    coord_print(2*2.^(-1.0*collect(l0+1:no_levels) ), variance[l0+1:no_levels])
    println("time per coarse uncoupled sample=",ccost)
    println("                  coupled sample=",cost)
    println("variance for uncoupled =",cvariance)
    println("              coupled  =",variance)
  end
  #################################################################### compute cpu times
  level_samps=zeros(Int64,no_levels);  total_time=zeros(no_levels)
  #
  cputimes=zeros(no_levels-l0+1,2); tols=zeros(no_levels-l0+1)
  for my_no_levels=l0:no_levels # outer loop (max levels)
    my_k=1+my_no_levels-l0
    hfine=2.0^(-my_no_levels); my_tol=hfine^2;
    println("max_levels=", my_no_levels,", hfine=", hfine,"; tolerance=",my_tol);
    for ell=my_no_levels:-1:l0 # inner loop coarse level
      r=ell+1:my_no_levels
      level_samps[1:ell]*=0
      if ell<my_no_levels
        myconst=(sqrt(cvariance[ell]*ccost[ell])
                +sum(sqrt.(variance[r].*cost[r]))
                )*(2/my_tol^2)
        level_samps[r]  =Int64.(ceil.(sqrt.(variance[r]./cost[r]) *myconst))
      else # vanilla
        myconst=sqrt(cvariance[ell]*ccost[ell])*(2/my_tol^2)
      end
      level_samps[ell]=Int64(ceil(sqrt(cvariance[ell]/ccost[ell])*myconst) )
      #
      mycost=deepcopy(cost); mycost[ell]=ccost[ell]
      myvars=deepcopy(variance); myvars[ell]=cvariance[ell]
      total_time[ell]=sum(mycost.*level_samps)
      #
      if report_flag
        println("ell=",ell,"   level_samps=",level_samps[ell:my_no_levels])
        println("             variances=",myvars[ell:my_no_levels])
        println("          time per level=",(mycost.*level_samps)[ell:my_no_levels],
              " total time=",total_time[ell])
      end # end if
    end # end inner loop
    # find quickest option
    mint,min_ell=findmin(total_time[l0:my_no_levels]);
    ell=min_ell+l0-1;   r=(ell+1):my_no_levels
    if ell<my_no_levels
        myconst=(sqrt(cvariance[ell]*ccost[ell])
                +sum(sqrt.(variance[r].*cost[r]))
                )*(2/my_tol^2)
      else
        myconst=sqrt(cvariance[ell]*ccost[ell])*(2/my_tol^2)
      end
    level_samps[1:ell]*=0
    level_samps[ell]=Int64(ceil(sqrt(cvariance[ell]/ccost[ell])*myconst) )
    level_samps[r]=Int64.(ceil.(sqrt.(variance[r]./cost[r]) *myconst))
    if report_flag
      println(total_time )
      println("choose ell0=",ell," level_samps",
              level_samps[ell:my_no_levels])
      cputimes[my_k,1]=total_time[ell];
      tols[my_k]=(2.0^(-my_no_levels))^2;
      cputimes[my_k,2]=total_time[my_no_levels]
      println("% ML")
      coord_print(tols,cputimes[:,1])
      println("% Vanilla")
      coord_print(tols,cputimes[:,2])
    end
  end
  #
end
##############################################################
function coord_print(x,y)
  # convenient output of co-ordinates for including in LaTeX
  n=length(x)
  println("coordinates {")
  for i=1:n
    print("(",x[i],",",y[i],") ")
  end
  println("\n};")
end
#####################################################
function ann_ml_coupling(p, dict, G,  report_flag)
  # Multilevel field solve
  no_levels=dict["no_levels"];  l0=2
  #
  if report_flag
    println("====\nAnalyse multilevel params with lmax=", no_levels,
            " with tolerance=",dict["tol"]," and ",dict["no_samps"]," samples")
    println("No workers",nworkers())
    no_samps=dict["no_samps"]
  else
    no_samps=3
  end
  variance=zeros(no_levels);  cost=ones(no_levels)
  #######################################
  # coupled levels
  for ell=(l0+1):(no_levels-1)
    finest_level_int=intersect(G.int,G.get_level_ind(ell))
    correction_level=intersect(G.int,1:G.get_level_ind(ell)[end])
    tic()
    v,no_samps=coupled_factor_par2(no_samps, G, p, dict,
       correction_level, finest_level_int)
    variance[ell]=G.sq_int(sqrt.(v),ell)
    cost[ell]=toq()/no_samps
    println(ell," variance ", variance[ell], ", cost ", cost[ell])
  end # loop over levels
  #
  beta_ratio=log.(variance[l0+1:end-1]./variance[l0+2:end])/log(2)
  hv=2*2.^(-1.0*collect(l0+1:(no_levels-1)) ); vv=variance[l0+1:(no_levels-1)];
  #
  if report_flag
    println("beta ML param (base 2) ",beta_ratio)
    #coord_print(1:length(variance),variance)
    coord_print(hv, vv,p.alpha,true)
    println("time per coupled sample=",cost)
    println(" coupled variance=",variance)
  end
  ell=0
  #
  level_samps=zeros(Int64,no_levels);  total_time=zeros(no_levels)
  return ell, level_samps
end
#######################################
function coord_print(x,y,alpha,fit_flag)
  # convenient output of co-ordinates for including in LaTeX
  # include linear fit
  n=length(x)
  #
  println("coordinates {")
  for i=1:n
    print("(",x[i],",",y[i],") ")
  end
  println("\n}; % alpha=", alpha)
  #
  if (fit_flag)
    lft=linreg(log.(x),log.(y))
    yy=exp(lft[1])*x.^lft[2];
    println(lft)
    println("coordinates {")
    for i=1:n
      print("(",x[i],",",yy[i],") ")
    end
    println("\n};% slope ",lft[2], " alp/(alp+1)=", alpha/(1+alpha))
  end
end
#####################################################
function get_ml_params(p, dict, G, l0, no_samps, report_flag)
  # Multilevel field solve
  # find optimum parametrs (l0,L)
  no_levels=dict["no_levels"];
  #
  if report_flag
    println("====\nGet multilevel params with l0=",l0,
            " and lmax=", no_levels,
            " with tolerance=",dict["tol"])
  end
  # scale constant for stopping criteria (K_C)
  cvariance=zeros(no_levels);  ccost=ones(no_levels)
  variance=zeros(no_levels);  cost=ones(no_levels)
  ########################################
  # coarse level
  ell=l0
  coarse_level=intersect(G.int,1:G.get_level_ind(ell)[end])
  # initialise
  tic()
  v,no_samps=uncoupled_factor_par2(no_samps, G, p, dict,
       coarse_level)
  cvariance=G.sq_int(sqrt.(v),l0)
  ccost    =toq()/(no_samps)
  #
  #######################################
  # coupled levels
  for ell=(l0+1):(no_levels)
    finest_level_int=intersect(G.int,G.get_level_ind(ell))
    correction_level=intersect(G.int,1:G.get_level_ind(ell)[end])
    tic()
    v,no_samps=coupled_factor_par2(no_samps, G, p, dict,
       correction_level, finest_level_int)
    variance[ell]=G.sq_int(sqrt.(v),ell)
    cost[ell]=toq()/no_samps
  end # loop over correction levels
  #
  if report_flag
    println("===Completed ml params.")
  end
  #
  beta_ratio=log.(variance[l0+1:end-1]./variance[l0+2:end])/log(2)
  if report_flag
    println("beta ML param (base 2) ",beta_ratio)
    println("time per coarse uncoupled sample=",ccost)
    println("time per coupled sample=",cost)
    println(" coarse uncoupled variance=",cvariance)
    println(" coupled variance=",variance)
  end
  #
  level_samps=zeros(Int64,no_levels);  total_time=zeros(no_levels)
  #
  ell=l0
  r=ell+1:no_levels
  level_samps[1:ell]*=0
  if ell<no_levels
      myconst=(sqrt(cvariance*ccost)
                +sum(sqrt.(variance[r].*cost[r]))
                )*(2/dict["tol"]^2)
      level_samps[r]  =Int64.(ceil.(sqrt.(variance[r]./cost[r]) *myconst))
      level_samps[r]=min.(level_samps[r],dict["max_samples"])
  else
      myconst=sqrt(cvariance*ccost)*(2/dict["tol"]^2)
  end
  level_samps[ell]=Int64(ceil(sqrt(cvariance/ccost)*myconst) )
  level_samps[ell]=min.(level_samps[ell],dict["max_samples"])
    #
  mycost=deepcopy(cost); mycost[ell]=ccost
  total_time[ell]=sum(mycost.*level_samps)
  #
  if report_flag
     println("ell=",ell,"   level_samps=",level_samps)
     println("          time per level=",mycost.*level_samps,
              " total time=",total_time[ell])
  end
  #
  return level_samps
end
#####################################################
######################################################
function wos_field_ml(p, dict, G)
  # Multilevel field solve
  no_levels=dict["no_levels"]; l0=dict["l0"]
  # get key WoS functions
  level_samps=get_ml_params(p,dict,G,l0,600,false)
  #
  if haskey(dict,"report")
    println("====\nMultilevel field solve with l0=",l0,
            " and lmax=", no_levels,
            " with tolerance=",dict["tol"]  )
  end
  ########################################
  # coarse level
  coarse_level=intersect(G.int,1:G.get_level_ind(l0)[end])
  if haskey(dict,"report")
      println("level=",l0,
              ", points on level=",length(coarse_level),
              ", with samples=", level_samps[l0],
              "." )
  end
  tic()
  out,no_samps=uncoupled_factor_par(level_samps[l0], G, p, dict,
        coarse_level)
  elapsed=toq(); println("elapsed_time=",elapsed)
  #######################################
  # coupled levels
  for ell=(l0+1):no_levels
    # get indices
    finest_level_int=intersect(G.int,G.get_level_ind(ell))
    correction_level=intersect(G.int,1:G.get_level_ind(ell)[end])
    #
    if haskey(dict,"report")
      println("level=",ell,
              ", points on level=",length(finest_level_int),
              ", with samples=", level_samps[ell],".")
    end
    #
    tic()
    delta_sum,no_samps=coupled_factor_par(level_samps[ell],
                     G, p, dict,
                    correction_level, finest_level_int)
    elapsed=toq(); println("elapsed_time=",elapsed)
    #
    for i=finest_level_int # update out with correction
         pa=G.get_parents(i)
         out[i]=mean(out[pa])+delta_sum[i]
    end
    #
    if haskey(dict,"report")
      println("level=",ell,
              ", H=",G.scale(ell)^2, # length scale for this lev,
              ", points on level=",length(finest_level_int))
    end
  end # loop over correction levels
  #
  #
  if haskey(dict,"report")
    println("===Completed multilevel solve.")
  end
  return out
end
###########################################
###########################################
function coupled_factor_par(no_samps, G, p, dict,
       correction_level, finest_level_int)
  # compute contribution form a ML correction levels
  pts=G.get_pt_all(correction_level)
  N=G.n_int+G.n_ext; pa=zeros(Int64,N,2)
  #
  for i=finest_level_int # get parents
      pa[i,:] = G.get_parents(i)
  end
  #
  np=nworkers() # no processors
  delta_sum=SharedArray{Float64}((N,np) )
  my_rhs_ball=Array{Any}(np)
  my_step=Array{Any}(np)
  #
  for i=1:np
    my_step[i],my_rhs_ball[i],dist_fn,ext_fn=setup_run_wos(p,dict)
  end
  dist_fn=p.d; ext_fn=p.fn_ext
  #
  no_samps=ceil(no_samps/np)
  @sync begin
    for k in 1:np
      j=workers()[k]
      @async  remotecall_wait(chunkchunk!,j,
                    delta_sum, no_samps,  my_step[k],
                    my_rhs_ball[k], dist_fn, ext_fn,
                    finest_level_int,
                    correction_level, pa, pts)
    end # proc loop
  end # @sync
  if haskey(dict,"report")
    print(".")
  end
  return mean(delta_sum,2)/no_samps, no_samps*np
end
###########################################
function coupled_factor_par2(no_samps, G, p, dict,
       correction_level, finest_level_int)
  # variant of above with variance output
  pts=G.get_pt_all(correction_level)
  N=G.n_int+G.n_ext; pa=zeros(Int64,N,2)
  #
  for i=finest_level_int # get parents
      pa[i,:] = G.get_parents(i)
  end
  #
  np=nworkers()# no processors for parallel
  delta_sum=SharedArray{Float64}((N,np) )
  delta_sq_sum=SharedArray{Float64}((N,np) )
  #tmp=SharedArray{Float64}((N,np) )
  my_rhs_ball=Array{Any}(np)
  my_step=Array{Any}(np)
  #
  for i=1:np
    my_step[i],my_rhs_ball[i],dist_fn,ext_fn=setup_run_wos(p,dict)
  end
  dist_fn=p.d; ext_fn=p.fn_ext
  #
  #  if haskey(dict,"report")
  #   print("Start parallel comp on ", workers())
  # end
  #
  no_samps=ceil(no_samps/np)
  println("no_samps=", no_samps,", length correction_level=",length(correction_level))
  @sync begin
    for k in 1:nworkers()
      j=workers()[k]
      @async  remotecall_wait(chunkchunk2!,j,
                    delta_sum, delta_sq_sum, #tmp,
                    no_samps,
                    my_step[k],
                    my_rhs_ball[k], dist_fn, ext_fn,
                    finest_level_int,
                    correction_level, pa, pts)
    end # proc loop
  end # @sync
  if haskey(dict,"report")
    print(".")
  end
  m=mean(delta_sum,2)/no_samps
  m2=mean(delta_sq_sum,2)/no_samps
  return m2-m.^2, no_samps*np
end
###########################################
###########################################
function uncoupled_factor_par(no_samps, G, p, dict,
        coarse_level)
  # compute uncoupled contribution to ML
  pts=G.get_pt_all(coarse_level)
  N=G.n_int+G.n_ext;
  #
  np=nworkers()# no processors for parallel
  delta_sum=SharedArray{Float64}((N,np) )
  my_rhs_ball=Array{Any}(np)
  my_step=Array{Any}(np)
  #
  for i=1:np
    my_step[i],my_rhs_ball[i],dist_fn,ext_fn=setup_run_wos(p,dict)
  end
  dist_fn=p.d; ext_fn=p.fn_ext
  #
  no_samps=ceil(no_samps/np)
  @sync begin
    for k in 1:nworkers()
      j=workers()[k]
      @async  remotecall_wait(chunkchunkchunk!, j,
                    delta_sum,  no_samps,
                    my_step[k], my_rhs_ball[k], dist_fn, ext_fn,
                    coarse_level, pts)
    end # proc loop
  end # @sync
  if haskey(dict,"report")
    print(".")
  end
  return mean(delta_sum,2)/no_samps,no_samps*np
end
###########################################
###########################################
function uncoupled_factor_par2(no_samps, G, p, dict,
        coarse_level)
  # variant of above with variance output
  pts=G.get_pt_all(coarse_level)
  N=G.n_int+G.n_ext;
  #
  np=nworkers()# no processors for parallel
  delta_sum=SharedArray{Float64}((N,np) )
  delta_sq_sum=SharedArray{Float64}((N,np) )
  my_rhs_ball=Array{Any}(np)
  my_step=Array{Any}(np)
  #
  for i=1:np
    my_step[i],my_rhs_ball[i],dist_fn,ext_fn=setup_run_wos(p,dict)
  end
  dist_fn=p.d; ext_fn=p.fn_ext
  #
  no_samps=ceil(no_samps/np)
  @sync begin
    for k in 1:nworkers()
      j=workers()[k]
      @async  remotecall_wait(chunkchunkchunk2!,j,
                    delta_sum, delta_sq_sum,  no_samps,
                    my_step[k], my_rhs_ball[k], dist_fn, ext_fn,
                    coarse_level, pts)
    end # proc loop
  end # @sync
  m=mean(delta_sum,2)/no_samps
  m2=mean(delta_sq_sum,2)/no_samps
  if haskey(dict,"report")
    print(".")
  end
  return abs.(m2-m.^2), no_samps*np
end
###########################
###########################
function get_solution(p, dict, G)
  # driver routine for access to WOS field solve
  tic()
  if  dict["solve_method"]=="simple"
    l0=dict["l0"]
    dict["l0"]=dict["no_levels"]
    out=wos_field_ml(p, dict, G)
    dict["l0"]=l0
  elseif dict["solve_method"]=="ml_analyse"
    # find bests ML paramters
    #ann_ml_params(p, dict, G,3,true)
    ann_ml_params(p, dict, G,dict["no_samps"],true)
    out=zeros(G.n_int)
  elseif dict["solve_method"]=="couple_analyse"
    ann_ml_coupling(p, dict, G,true)
    #
    out=zeros(G.n_int)
  else
    out=wos_field_ml(p, dict, G)
  end
  elapsed_time=toq()

  return out[:],elapsed_time
end
######################
######################
function  wos_ev(p::KOS.Params, dict, G, arnoldi, lam_old)
  # drive routine for access to Arnoldi-method based for finding eigenvalues
  #
  if haskey(dict,"report")
    println("====\nStarting Arnoldi iteration",
    " lam_tol=",dict["lam_tol"], " and ",dict["no_arnoldi_its"]," iterations.")
  end
  tic()
  #
  my_step,my_rhs_ball,dist_fn,ext_fn=setup_run_wos(p,dict)
  #
  its=arnoldi.k; eig_res=1.; delta_lam=1.; sep=1.
  while (its<dict["no_arnoldi_its"]
    && eig_res>0
    #    && delta_lam>dict["lam_tol"]
    )

    its+=1
    # update rhs with interpolant
    smooth_interpolate!(G,dict,p,arnoldi.rhs,1.)
    my_rhs_ball=update_rhs(p,dict)
    # set tolerance
    if dict["relax"]==true && its>1
      relax_factor=max(1,sep/eig_res)
    else
      relax_factor=1.
    end
    dict["tol"]=relax_factor*dict["lam_tol"]/(dict["no_arnoldi_its"]-1)/dict["B"]
    # WoS solve (inner iteration)
    out,=get_solution(p, dict, G)
    # Arnoldi
    lambda,eig_res,sep=Add_Arnoldi!(arnoldi,out[G.int]) # sep=0.2 artificial
    lam_inv=1./lambda
    delta_lam=abs(lam_inv-lam_old); lam_old=lam_inv
    #
    if eig_res>0 && haskey(dict,"report")
      println("\nArnoldi it=",its,
              ", lam_inv=",lam_inv,
              ", eig_res=",eig_res,
              ", delta_lam=",delta_lam,", sep=",sep,".")
    end
  end
  lam,evec,all_lam,sep=Final_Arnoldi(arnoldi)
  println("All eigenvalues=",1./all_lam)
  el_time=toq()
  println("Elapsed time for wos_ev() is ", el_time)
  return (1./real(lam)),evec,delta_lam,arnoldi,el_time,eig_res
end
#
end# end of module kos_field_solve


text:
using YAXArrays
using Test

include("tools.jl")
include("Cubes/axes.jl")
include("Cubes/cubes.jl")
include("Cubes/transformedcubes.jl")

include("Datasets/datasets.jl")

include("DAT/PickAxisArray.jl")
include("DAT/MovingWindow.jl")
include("DAT/tablestats.jl")
include("DAT/mapcube.jl")
include("DAT/DAT.jl")
include("DAT/loopchunks.jl")

max_stars_repo_name: meggart/YAXArrays

max_forks_repo_name: meggart/YAXArrays


text:
#=
Created on 02/02/2021 14:28:02
Last update: -

@author: Michiel Stock
michielfmstock@gmail.com

A simple Dungeons and Dragons interface to illustrate the (composite) type system.
=#

# First we define all the different types we want to use.

abstract type Creature end

abstract type Humanoid <: Creature end

#=
Notice that we want to use a mutable struct for the characters below since
the values of their field will be changed throughout the game.
=#

mutable struct Human <: Humanoid
    name::String
    hp::Int
    level::Int
    gold::Int
end

Human(name) = Human(name, 10, 1, rand(20:50))

mutable struct Elf <: Humanoid
    name::String
    hp::Int
    level::Int
    gold::Int
    function Elf(name)
        return new(name, 12, 1, rand(10:40))
    end
end

mutable struct Gnome <: Humanoid
    name::String
    hp::Int
    level::Int
    spells::Array{String}
    gold::Int
end

abstract type Monster <: Creature end

#=
Below you can find an example of when you would use the Union{} function.
This makes is so that all types that fall under the Float64 or Int type will be accepted.
=#

mutable struct Dragon{T<:Union{Float64,Int}} <: Creature
    hp::T
    strength::Int
    magic::Int
end

mutable struct GelatinousCube{T<:Union{Float64,Int}} <: Creature
    hp::T
    strength::Int
    object::String
end

# Lastly we define the functions we need.

languages(::Humanoid) = "Common Tongue"
languages(::Elf) = "Common Tongue, Elvish"

encounter(::Creature, ::Creature) = nothing
encounter(::Humanoid, ::Monster) = "run!"
encounter(::Monster, ::Humanoid) = "attack!"
encounter(::Humanoid, ::Humanoid) = "talk"
encounter(::Gnome, ::Humanoid) = "steals gold..."

Base.size(::Creature) = "medium"
Base.size(::Gnome) = "small"
Base.size(::Dragon) = "large"

damage(opponent, dragon::Dragon{Int}) = max(0, opponent.level - dragon.strength)
damage(opponent, dragon::Dragon{Float64}) = max(0.0, (opponent.level - dragon.strength) / dragon.strength)


text:
"""
Test for AVR model : EXST1 available in PSS/e
This case study defines a four bus system with an infinite bus in 1,
a GENSAL in bus 2 and a constant impedance load in bus 3
The GENSAL machine has the EXST1 and a HYGOV.
The disturbance is the outage of one line between buses 1 and 4
"""
##################################################
############### LOAD DATA ########################
##################################################

include(joinpath(TEST_FILES_DIR, "data_tests/test39.jl"))

##################################################
############### SOLVE PROBLEM ####################
##################################################

csv_file = joinpath(TEST_FILES_DIR, "benchmarks/psse//EXST1/results_PSSe.csv")

@testset "Test 39 EXST1 ResidualModel" begin
    path = joinpath(pwd(), "test-psse-exst1")
    !isdir(path) && mkdir(path)
    try
        # Define Simulation Problem
        sim = Simulation(
            ResidualModel,
            sys, #system
            path,
            (0.0, 20.0), #time span
            BranchTrip(1.0, Line, "BUS1-BUS4-i_1"), #Type of Fault
        )

        # Test Initial Condition
        diff_val = [0.0]
        res = get_init_values_for_comparison(sim)
        for (k, v) in test39_x0_init
            diff_val[1] += LinearAlgebra.norm(res[k] - v)
        end

        @test diff_val[1] < 1e-3

        # Obtain small signal results for initial conditions
        small_sig = small_signal_analysis(sim)
        eigs = small_sig.eigenvalues
        @test small_sig.stable

        # Test Eigenvalues
        @test LinearAlgebra.norm(eigs - test39_eigvals) < 1e-3

        # Solve problem
        @test execute!(sim, IDA(), dtmax = 0.005, saveat = 0.005) ==
              PSID.SIMULATION_FINALIZED
        results = read_results(sim)

        # Obtain results

        t_psid, v2_psid = get_voltage_magnitude_series(results, 2)
        _, v3_psid = get_voltage_magnitude_series(results, 3)
        _, ω_psid = get_state_series(results, ("generator-2-1", :ω))

        # Obtain PSSE results
        M = get_csv_data(csv_file)
        t_psse = M[:, 1]
        v1_psse = M[:, 2]
        v2_psse = M[:, 3]
        v3_psse = M[:, 4]
        v4_psse = M[:, 5]
        ω_psse = M[:, 6] .+ 1.0
        efd_psse = M[:, 7]

        # Test Transient Simulation Results

        @test LinearAlgebra.norm(t_psid - round.(t_psse, digits = 3)) == 0.0
        @test LinearAlgebra.norm(v2_psid - v2_psse, Inf) <= 1e-3
        @test LinearAlgebra.norm(v3_psid - v3_psse, Inf) <= 1e-3
        @test LinearAlgebra.norm(ω_psid - ω_psse, Inf) <= 1e-3

    finally
        @info("removing test files")
        rm(path, force = true, recursive = true)
    end
end


text:
module Wrapper
using Libdl
using LinearAlgebra
using OffsetArrays
using ..Mie
using ..SpecialFunctions

# Constants <=> numconstants

function init!(mstm, notd::Int64)
    ccall(Libdl.dlsym(mstm, :__numconstants_MOD_init), Cvoid, (Ref{Int32},), convert(Int32, notd))

    return
end

function bcof(mstm, notd::Int64)
    nbc = 6notd + 6
    bcof_addr_ptr = cglobal(Libdl.dlsym(mstm, :__numconstants_MOD_bcof), UInt64)
    bcof_addr = convert(Ptr{Float64}, unsafe_load(bcof_addr_ptr))
    bcof = unsafe_wrap(Array{Float64,2}, bcof_addr, (nbc + 1, nbc + 1))

    return OffsetArray(bcof, 0:nbc, 0:nbc)
end

function vwh_coef(mstm, notd::Int64)
    vwh_coef_addr_ptr = cglobal(Libdl.dlsym(mstm, :__numconstants_MOD_vwh_coef), UInt64)
    vwh_coef_addr = convert(Ptr{Float64}, unsafe_load(vwh_coef_addr_ptr))
    vwh_coef = unsafe_wrap(Array{Float64}, vwh_coef_addr, ((2notd + 1) * notd * 3 * 3,))

    return OffsetArray(reshape(vwh_coef, 2notd + 1, notd, 3, 3), (-notd):notd, 1:notd, -1:1, -1:1)
end

# SpecialFunctions <=> specialfuncs

function ricbessel(mstm, n::Int64, ds::Float64, ϵ::Float64)
    Ψ = zeros(n + 1)
    ccall(
        Libdl.dlsym(mstm, :__specialfuncs_MOD_ricbessel),
        Cvoid,
        (Ref{Int32}, Ref{Float64}, Ref{Float64}, Ref{Int64}, Ptr{Float64}),
        convert(Int32, n),
        ds,
        ϵ,
        0,
        Ψ,
    )
    return OffsetArray(Ψ, 0:n)
end

function richankel(mstm, n::Int64, ds::Float64)
    Ξ = zeros(ComplexF64, n + 1)
    ccall(
        Libdl.dlsym(mstm, :__specialfuncs_MOD_richankel),
        Cvoid,
        (Ref{Int32}, Ref{Float64}, Ptr{ComplexF64}),
        convert(Int32, n),
        ds,
        Ξ,
    )
    return OffsetArray(Ξ, 0:n)
end

function cricbessel(mstm, n::Int64, ds::ComplexF64)
    Ψ = zeros(ComplexF64, n + 1)

    ccall(
        Libdl.dlsym(mstm, :__specialfuncs_MOD_cricbessel),
        Cvoid,
        (Ref{Int32}, Ref{ComplexF64}, Ptr{ComplexF64}),
        convert(Int32, n),
        ds,
        Ψ,
    )
    return OffsetArray(Ψ, 0:n)
end

function crichankel(mstm, n::Int64, ds::ComplexF64)
    Ξ = zeros(ComplexF64, n + 1)
    ccall(
        Libdl.dlsym(mstm, :__specialfuncs_MOD_crichankel),
        Cvoid,
        (Ref{Int32}, Ref{ComplexF64}, Ptr{ComplexF64}),
        convert(Int32, n),
        ds,
        Ξ,
    )
    return OffsetArray(Ξ, 0:n)
end

function cspherebessel(mstm, n::Int64, z::ComplexF64)
    csj = zeros(ComplexF64, n + 1)
    csy = zeros(ComplexF64, n + 1)
    ccall(
        Libdl.dlsym(mstm, :__specialfuncs_MOD_cspherebessel),
        Cvoid,
        (Ref{Int32}, Ref{ComplexF64}, Ptr{ComplexF64}, Ptr{ComplexF64}),
        convert(Int32, n),
        z,
        csj,
        csy,
    )
    return OffsetArray(csj, 0:n), OffsetArray(csy, 0:n)
end

function vcfunc(mstm, m::Int64, n::Int64, k::Int64, l::Int64)
    init!(mstm, n + l)

    vcn = zeros(n + l + 1)
    ccall(
        Libdl.dlsym(mstm, :__specialfuncs_MOD_vcfunc),
        Cvoid,
        (Ref{Int32}, Ref{Int32}, Ref{Int32}, Ref{Int32}, Ptr{Float64}),
        convert(Int32, m),
        convert(Int32, n),
        convert(Int32, k),
        convert(Int32, l),
        vcn,
    )
    return OffsetArray(vcn, 0:(n + l))
end

function normalizedlegendre(mstm, cbe::Float64, mmax::Int64, nmax::Int64)
    init!(mstm, mmax + nmax)

    dc = zeros(2mmax + 1, nmax + 1)
    ccall(
        Libdl.dlsym(mstm, :__specialfuncs_MOD_normalizedlegendre),
        Cvoid,
        (Ref{Float64}, Ref{Int64}, Ref{Int64}, Ptr{Float64}),
        cbe,
        mmax,
        nmax,
        dc,
    )
    return OffsetArray(dc, (-mmax):mmax, 0:nmax)
end

function rotcoef(mstm, cbe::Float64, kmax::Int64, nmax::Int64)
    init!(mstm, kmax + nmax)

    dc = zeros(2kmax + 1, nmax * (nmax + 2) + 1)
    ccall(
        Libdl.dlsym(mstm, :__specialfuncs_MOD_rotcoef),
        Cvoid,
        (Ref{Float64}, Ref{Int64}, Ref{Int64}, Ptr{Float64}),
        cbe,
        kmax,
        nmax,
        dc,
    )
    return OffsetArray(dc, (-kmax):kmax, 0:(nmax * (nmax + 2)))
end

function taufunc(mstm, cb::Float64, nmax::Int64)
    init!(mstm, nmax + 1)

    τ = zeros(nmax + 2, nmax, 2)
    ccall(Libdl.dlsym(mstm, :__specialfuncs_MOD_taufunc), Cvoid, (Ref{Float64}, Ref{Int64}, Ptr{Float64}), cb, nmax, τ)
    return OffsetArray(τ, 0:(nmax + 1), 1:nmax, 1:2)
end

function pifunc(mstm, cb::Float64, ephi::ComplexF64, nmax::Int64, ndim::Int64)
    init!(mstm, nmax + 1)

    π_vec = zeros(ComplexF64, ndim + 2, ndim, 2)
    ccall(
        Libdl.dlsym(mstm, :__specialfuncs_MOD_pifunc),
        Cvoid,
        (Ref{Float64}, Ref{ComplexF64}, Ref{Int32}, Ref{Int32}, Ptr{ComplexF64}),
        cb,
        ephi,
        convert(Int32, nmax),
        convert(Int32, ndim),
        π_vec,
    )
    return OffsetArray(π_vec, 0:(ndim + 1), 1:ndim, 1:2)
end

function planewavecoef(mstm, α::Float64, β::Float64, nodr::Int64)
    pmnp0 = zeros(ComplexF64, nodr + 2, nodr, 2, 2)
    ccall(
        Libdl.dlsym(mstm, :__specialfuncs_MOD_planewavecoef),
        Cvoid,
        (Ref{Float64}, Ref{Float64}, Ref{Int32}, Ptr{ComplexF64}),
        α,
        β,
        convert(Int32, nodr),
        pmnp0,
    )
    return OffsetArray(pmnp0, 0:(nodr + 1), 1:nodr, 1:2, 1:2)
end

function gaussianbeamcoef(mstm, α::Float64, β::Float64, cbeam::Float64, nodr::Int64)
    pmnp0 = zeros(ComplexF64, nodr + 2, nodr, 2, 2)
    ccall(
        Libdl.dlsym(mstm, :__specialfuncs_MOD_gaussianbeamcoef),
        Cvoid,
        (Ref{Float64}, Ref{Float64}, Ref{Float64}, Ref{Int32}, Ptr{ComplexF64}),
        α,
        β,
        cbeam,
        convert(Int32, nodr),
        pmnp0,
    )
    return OffsetArray(pmnp0, 0:(nodr + 1), 1:nodr, 1:2, 1:2)
end

function sphereplanewavecoef(
    mstm,
    nodr::Array{Int64,1},
    α::Float64,
    β::Float64,
    rpos::Array{Float64,2},
    hostsphere::Array{Int64,1},
    numberfieldexp::Array{Int64,1},
    rimedium::Array{ComplexF64,1},
)
    nsphere = length(nodr)
    nodrmax = maximum(nodr)

    # Valid the input
    @assert length(hostsphere) == nsphere
    @assert length(numberfieldexp) == nsphere
    @assert size(rpos) == (3, nsphere)
    @assert length(rimedium) == 2

    neqns = 2 * sum(numberfieldexp[i] * nodr[i] * (nodr[i] + 2) for i in 1:nsphere)
    pmnp = zeros(ComplexF64, 2neqns)

    ccall(
        Libdl.dlsym(mstm, :__specialfuncs_MOD_sphereplanewavecoef),
        Cvoid,
        (
            Ref{Int32},
            Ref{Int32},
            Ptr{Int32},
            Ref{Int32},
            Ref{Float64},
            Ref{Float64},
            Ptr{Float64},
            Ptr{Int32},
            Ptr{Int32},
            Ptr{ComplexF64},
            Ptr{ComplexF64},
        ),
        convert(Int32, nsphere),
        convert(Int32, neqns),
        convert(Array{Int32,1}, nodr),
        convert(Int32, nodrmax),
        α,
        β,
        rpos,
        convert(Array{Int32,1}, hostsphere),
        convert(Array{Int32,1}, numberfieldexp),
        rimedium,
        pmnp,
    )

    return pmnp
end

function axialtrancoefrecurrence(mstm, itype::Int64, r::Float64, ri::Array{ComplexF64,1}, nmax::Int64, lmax::Int64)
    @assert itype == 1 || itype == 3
    @assert length(ri) == 2

    init!(mstm, nmax)
    ndim = SpecialFunctions.atcdim(nmax, lmax)
    ac = zeros(ComplexF64, ndim)

    ccall(
        Libdl.dlsym(mstm, :__specialfuncs_MOD_axialtrancoefrecurrence),
        Cvoid,
        (Ref{Int32}, Ref{Float64}, Ptr{ComplexF64}, Ref{Int32}, Ref{Int32}, Ref{Int32}, Ptr{ComplexF64}),
        convert(Int32, itype),
        r,
        ri,
        convert(Int32, nmax),
        convert(Int32, lmax),
        convert(Int32, ndim),
        ac,
    )

    return ac
end

function tranordertest(mstm, r::Float64, ri::ComplexF64, lmax::Int64, ϵ::Float64)
    nmax = Ref{Int32}(0)

    ccall(
        Libdl.dlsym(mstm, :__specialfuncs_MOD_tranordertest),
        Cvoid,
        (Ref{Float64}, Ref{ComplexF64}, Ref{Int32}, Ref{Float64}, Ref{Int32}),
        r,
        ri,
        convert(Int32, lmax),
        ϵ,
        nmax,
    )

    return convert(Int64, nmax.x)
end

function gentrancoef(
    mstm,
    itype::Int64,
    xptran::Array{Float64,1},
    ri::Array{ComplexF64,1},
    nrow0::Int64,
    nrow1::Int64,
    ncol0::Int64,
    ncol1::Int64,
    iaddrow0::Int64,
    iaddcol0::Int64,
)
    @assert itype == 1 || itype == 3
    @assert length(xptran) == 3
    @assert length(ri) == 2
    @assert 1 <= nrow0 <= nrow1
    @assert 1 <= ncol0 <= ncol1

    ac = zeros(
        ComplexF64,
        2,
        nrow1 * (nrow1 + 2) - (nrow0 - 1) * (nrow0 - 1) + iaddrow0,
        ncol1 * (ncol1 + 2) - (ncol0 - 1) * (ncol0 + 1) + iaddcol0,
    )

    ccall(
        Libdl.dlsym(mstm, :__specialfuncs_MOD_gentrancoef),
        Cvoid,
        (
            Ref{Int32},
            Ptr{Float64},
            Ptr{ComplexF64},
            Ref{Int32},
            Ref{Int32},
            Ref{Int32},
            Ref{Int32},
            Ref{Int32},
            Ref{Int32},
            Ptr{ComplexF64},
        ),
        convert(Int32, itype),
        xptran,
        ri,
        convert(Int32, nrow0),
        convert(Int32, nrow1),
        convert(Int32, ncol0),
        convert(Int32, ncol1),
        convert(Int32, iaddrow0),
        convert(Int32, iaddcol0),
        ac,
    )

    return ac
end

function eulerrotation(mstm, xp::Array{Float64,1}, eulerangf::Array{Float64,1}, dir::Int64)
    xprot = zeros(3)

    ccall(
        Libdl.dlsym(mstm, :__specialfuncs_MOD_eulerrotation),
        Cvoid,
        (Ptr{Float64}, Ptr{Float64}, Ref{Int32}, Ptr{Float64}),
        xp,
        eulerangf,
        convert(Int32, dir),
        xprot,
    )

    return xprot
end

function planewavetruncationorder(mstm, r::Float64, rimedium::Array{ComplexF64,1}, ϵ::Float64)
    nodr = Ref{Int32}(0)

    ccall(
        Libdl.dlsym(mstm, :__specialfuncs_MOD_planewavetruncationorder),
        Cvoid,
        (Ref{Float64}, Ptr{ComplexF64}, Ref{Float64}, Ref{Int32}),
        r,
        rimedium,
        ϵ,
        nodr,
    )

    return convert(Int64, nodr.x)
end

function vwhcalc(mstm, rpos::Array{Float64,1}, ri::Array{ComplexF64,1}, nodr::Int64, itype::Int64)
    @assert itype == 1 || itype == 3
    @assert length(rpos) == 3
    @assert length(ri) == 2

    vwh = zeros(ComplexF64, 3, 2, nodr * (nodr + 2))

    ccall(
        Libdl.dlsym(mstm, :__specialfuncs_MOD_vwhcalc),
        Cvoid,
        (Ptr{Float64}, Ptr{ComplexF64}, Ref{Int32}, Ref{Int32}, Ptr{ComplexF64}),
        rpos,
        ri,
        convert(Int32, nodr),
        convert(Int32, itype),
        vwh,
    )

    return vwh
end

function vwhaxialcalc(mstm, rpos::Array{Float64,1}, ri::Array{ComplexF64,1}, nodr::Int64, itype::Int64)
    @assert itype == 1 || itype == 3
    @assert length(rpos) == 3
    @assert length(ri) == 2

    vwh = zeros(ComplexF64, 3, 2, 2, nodr)

    ccall(
        Libdl.dlsym(mstm, :__specialfuncs_MOD_vwhaxialcalc),
        Cvoid,
        (Ptr{Float64}, Ptr{ComplexF64}, Ref{Int32}, Ref{Int32}, Ptr{ComplexF64}),
        rpos,
        ri,
        convert(Int32, nodr),
        convert(Int32, itype),
        vwh,
    )

    return vwh
end

function twobytwoinverse(mstm, mat::Array{ComplexF64,2})
    @assert size(mat) == (2, 2)
    @assert det(mat) != 0

    imat = zeros(ComplexF64, 2, 2)

    ccall(Libdl.dlsym(mstm, :__specialfuncs_MOD_twobytwoinverse), Cvoid, (Ptr{ComplexF64}, Ptr{ComplexF64}), mat, imat)

    return imat
end

# Mie <=> miedata

function miecoefcalc(
    mstm,
    xsp::Array{Float64,1},
    ri::OffsetArray{ComplexF64,2},
    hostsphere::Array{Int64,1},
    numberfieldexp::Array{Int64,1},
    qeps::Float64,
)
    nsphere = length(xsp)
    @assert size(ri) == (2, nsphere + 1)
    @assert length(hostsphere) == nsphere
    @assert length(numberfieldexp) == nsphere

    hostsphere = convert(Array{Int32,1}, hostsphere)
    numberfieldexp = convert(Array{Int32,1}, numberfieldexp)

    tmatrix_file = Ptr{String}()
    ccall(
        Libdl.dlsym(mstm, :__miecoefdata_MOD_miecoefcalc),
        Cvoid,
        (Ref{Int32}, Ptr{Float64}, Ptr{ComplexF64}, Ptr{Int32}, Ptr{Int32}, Ref{Float64}, Ptr{String}),
        convert(Int32, nsphere),
        xsp,
        ri,
        hostsphere,
        numberfieldexp,
        qeps,
        tmatrix_file,
    )

    return getmiedataall(mstm, nsphere)
end

function mieoa(
    mstm,
    x::Float64,
    ri::Array{ComplexF64,1},
    nodr0::Int64,
    qeps::Float64;
    ri_medium::Union{Array{ComplexF64,1},Nothing} = nothing,
)
    @assert length(ri) == 2
    @assert isnothing(ri_medium) || length(ri_medium) == 2

    qext = Ref{Float64}(0.0)
    qsca = Ref{Float64}(0.0)
    qabs = Ref{Float64}(0.0)
    nodr = Ref{Int32}(nodr0)

    # FIXME: present check does not work normally, pass dummy arrays as a workaround
    anp_mie = zeros(ComplexF64, 2, 2, 100)
    dnp_mie = zeros(ComplexF64, 2, 2, 100)
    unp_mie = zeros(ComplexF64, 2, 2, 100)
    vnp_mie = zeros(ComplexF64, 2, 2, 100)
    cnp_mie = zeros(ComplexF64, 2, 2, 100)
    if isnothing(ri_medium)
        ri_medium = [1.0 + 0.0im, 1.0 + 0.0im]
    end
    anp_inv_mie = zeros(ComplexF64, 2, 2, 100)

    ccall(
        Libdl.dlsym(mstm, :__miecoefdata_MOD_mieoa),
        Cvoid,
        (
            Ref{Float64},
            Ptr{ComplexF64},
            Ref{Int32},
            Ref{Float64},
            Ref{Float64},
            Ref{Float64},
            Ref{Float64},
            Ptr{ComplexF64},
            Ptr{ComplexF64},
            Ptr{ComplexF64},
            Ptr{ComplexF64},
            Ptr{ComplexF64},
            Ptr{ComplexF64},
            Ptr{ComplexF64},
        ),
        x,
        ri,
        nodr,
        qeps,
        qext,
        qsca,
        qabs,
        anp_mie,
        dnp_mie,
        unp_mie,
        vnp_mie,
        cnp_mie,
        ri_medium,
        anp_inv_mie,
    )

    return qext[], qsca[], qabs[], convert(Int64, nodr[])
end

function getmiedataall(mstm, nsphere::Int64)::MieData
    # Do not know termstot at first, so we set an's and cn's sizes to a large number.

    MAXIMUM_TERM = 100_000
    number_eqns = Ref{Int32}(0)
    max_order = Ref{Int32}(0)
    order = zeros(Int32, nsphere)
    offset = zeros(Int32, nsphere + 1)
    block = zeros(Int32, nsphere)
    block_offset = zeros(Int32, nsphere + 1)
    number_field_expansions = zeros(Int32, nsphere)
    is_optically_active = zeros(Int32, nsphere)
    qext = zeros(nsphere)
    qabs = zeros(nsphere)
    an = zeros(ComplexF64, MAXIMUM_TERM)
    cn = zeros(ComplexF64, MAXIMUM_TERM)

    ccall(
        Libdl.dlsym(mstm, :__miecoefdata_MOD_getmiedataallreq),
        Cvoid,
        (
            Ref{Int32},
            Ptr{Int32},
            Ptr{Int32},
            Ptr{Int32},
            Ptr{Int32},
            Ptr{Float64},
            Ptr{Float64},
            Ptr{ComplexF64},
            Ptr{ComplexF64},
            Ref{Int32},
            Ref{Int32},
            Ptr{Int32},
            Ptr{Int32},
        ),
        convert(Int32, nsphere),
        order,
        block,
        offset,
        block_offset,
        qext,
        qabs,
        an,
        cn,
        number_eqns,
        max_order,
        is_optically_active,
        number_field_expansions,
    )

    mie = MieData()
    mie.number_spheres = nsphere
    mie.block = convert(Array{Int64,1}, block)
    mie.offset = convert(Array{Int64,1}, offset)
    mie.block_offset = convert(Array{Int64,1}, block_offset)
    mie.is_optically_active = convert(Array{Bool, 1}, is_optically_active)
    mie.order = convert(Array{Int64,1}, order)
    mie.max_order = convert(Int64, max_order.x)
    mie.number_field_expansions = convert(Array{Int64,1}, number_field_expansions)
    mie.number_eqns = convert(Int32, number_eqns.x)
    mie.qext = qext
    mie.qabs = qabs

    ntermstot = offset[nsphere + 1]
    mie.an = an[1:ntermstot]
    mie.cn = cn[1:ntermstot]
    return mie
end

function onemiecoeffmult(mstm, i::Int64, nodr::Int64, cx::OffsetArray{ComplexF64,3}, mie_coefficent::Char = 'a')
    cy = OffsetArray(zeros(ComplexF64, nodr + 2, nodr, 2), 0:(nodr + 1), 1:nodr, 1:2)

    ccall(
        Libdl.dlsym(mstm, :__miecoefdata_MOD_onemiecoeffmult),
        Cvoid,
        (Ref{Int32}, Ref{Int32}, Ptr{ComplexF64}, Ptr{ComplexF64}, Ref{Cchar}),
        convert(Int32, i),
        convert(Int32, nodr),
        cx,
        cy,
        mie_coefficent,
    )

    return cy
end

function multmiecoeffmult(
    mstm,
    neqns::Int64,
    nrhs::Int64,
    idir::Int64,
    ain::Array{ComplexF64,1},
    rhs_list::Union{Array{Bool,1},Nothing} = nothing,
)
    @assert length(ain) == neqns * nrhs
    @assert isnothing(rhs_list) || length(rhs_list) == nrhs

    if isnothing(rhs_list)
        rhs_list = fill(true, nrhs)
    end

    aout = zeros(ComplexF64, neqns * nrhs)

    ccall(
        Libdl.dlsym(mstm, :__miecoefdata_MOD_multmiecoeffmult),
        Cvoid,
        (Ref{Int32}, Ref{Int32}, Ref{Int32}, Ptr{ComplexF64}, Ptr{ComplexF64}, Ptr{Int32}),
        convert(Int32, neqns),
        convert(Int32, nrhs),
        convert(Int32, idir),
        ain,
        aout,
        convert(Array{Int32,1}, rhs_list), # Fortran bool takes 4 bytes by default !!!
    )

    return aout
end

# # ? <=> spheredata

# function getspheredata(mstm)
#     # Do not know nspheres at first, so we set the arrays' sizes to be a large number.

#     MAXIMUM_NSPHERE = 100_000
#     number_spheres = Ref{Int32}(0)
#     sphere_size_parameters = zeros(MAXIMUM_NSPHERE)
#     sphere_positions = zeros(3, MAXIMUM_NSPHERE)
#     sphere_refractive_indices = OffsetArray(zeros(ComplexF64, 2, MAXIMUM_NSPHERE + 1), 1:2, 0:MAXIMUM_NSPHERE)
#     volume_size_parameter = Ref{Float64}(0.0)
#     host_spheres = zeros(Int32, MAXIMUM_NSPHERE)
#     tmatrix_on_file = fill(false, MAXIMUM_NSPHERE)
#     tmatrix_file = Ptr{String}()

#     ccall(
#         Libdl.dlsym(mstm, :__spheredata_MOD_getspheredata),
#         Cvoid,
#         (Ref{Int32}, Ptr{Float64}, Ptr{Float64}, Ptr{ComplexF64}, Ref{Float64}, Ptr{Int32}, Ptr{Bool}, Ptr{String}),
#         number_spheres,
#         Ptr{Float64}(),
#         Ptr{Float64}(),
#         Ptr{ComplexF64}(),
#         volume_size_parameter,
#         Ptr{Int32}(),
#         Ptr{Bool}(),
#         Ptr{String}(),
#         # sphere_size_parameters,
#         # sphere_positions,
#         # sphere_refractive_indices,
#         # volume_size_parameter,
#         # host_spheres,
#         # tmatrix_on_file,
#         # tmatrix_file,
#     )

#     number_spheres = convert(Int64, number_spheres.x)
#     return number_spheres
# end

end # module Wrapper


text:
############################################################################
#
#  Actuator.jl
#
#  Purpose:    Nonlinear Actuator model
#
#  Author:     Edward Daughtery
#
#  Date:       05 August 2016
#
#
############################################################################

using Aerospace

############################################################################
## Actuator Initialization
function Initialize_Act(dt::Float64,
                        wn::Float64,
                        zeta::Float64,
                        pos_lim_pos::Float64,
                        pos_lim_neg::Float64,
                        rate_lim_pos::Float64,
                        rate_lim_neg::Float64
                        )

    ## Time, Position, and Rate Initialization
    time = 0.0;
    position = 0.0;
    rate = 0.0;

    ## Creating Actuator Structure
    actout = Actuator(time, dt, wn, zeta, position, rate,
                      pos_lim_pos,
                      pos_lim_neg,
                      rate_lim_pos,
                      rate_lim_neg
                      );

    return actout;
end

############################################################################
## Actuator Update
function Update_Act(act::Actuator, position_cmd::Float64)

    act_out = act;

    ## Variable Initialization
    x1 = act.position;
    position = x1;
    x2 = act.rate;
    rate = x2;
    wn = act.wn;
    zeta = act.zeta;
    pos_lim_pos = act.pos_lim_pos;
    pos_lim_neg = act.pos_lim_neg;
    rate_lim_pos = act.pos_lim_pos;
    rate_lim_neg = act.rate_lim_neg;
    dt = act.dt;
    time = act.time;

    ## Limiting Command
    if (position_cmd > pos_lim_pos )
        position_cmd = pos_lim_pos;
    end

    if (position_cmd < pos_lim_neg )
        position_cmd = pos_lim_neg;
    end

    wn2 = wn*wn;
    d2 = 2.0*wn*zeta;

    ## Integrating to get X2
    dx2 = wn2*(position_cmd - position) - d2*rate;
    x2p = x2 + dt*dx2;

    ## Limiting Rate
    if (x2p > rate_lim_pos  )
        x2p = rate_lim_pos;
    end

    if (x2p < rate_lim_neg )
        x2p = rate_lim_neg;
    end

    ## Integrating to get X1
    dx1 = x2p;
    x1 = x1 + dt*dx1;

    ## Limiting Command
    if (x1 > pos_lim_pos )
        x1 = pos_lim_pos;
        x2 = 0.0;
    elseif (x1 < pos_lim_neg )
        x1 = pos_lim_neg;
        x2 = 0.0;
    else
        x2 = x2p;
    end

    ## Updating Time and Member data
    time = time + dt;
    position = x1;
    rate = x2;
    act_out.time = time;
    act_out.position = position;
    act_out.rate = rate;

    return act_out;
end