text:
# This file is part of EndemicModel project, which is licensed under BDS-3-Clause.
# See file LICENSE.md for more information.

# CSSE COVID-19 Repository

const CSSE_COLUMNS = (;
    Symbol("Country/Region") => :country,
    Symbol("Province/State") => :state,
    Symbol("Lat") => :latitude,
    Symbol("Long") => :longitude,
)

function _csse_url(tablename::String)
    "https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/" *
    "csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_" *
    tablename *
    "_global.csv"
end

const CSSE_TABLES = (
    confirmed = _csse_url("confirmed"),
    deaths = _csse_url("deaths"),
    recovered = _csse_url("recovered"),
)

"""
Novel Coronavirus (COVID-19) Cases, provided by JHU CSSE
https://github.com/CSSEGISandData/COVID-19
https://systems.jhu.edu/research/public-health/ncov/
"""
struct CsseSource <: AbstractDataSource
    function CsseSource()
        new()
    end
end
CsseSource(::Nothing; kwargs...) = CsseSource()

DATA_SOURCES[:csse] = CsseSource

function _import_data(source::CsseSource; kwargs...)
    result = DataDict()
    for (dfsym, url) ∈ pairs(CSSE_TABLES)
        file = import_data(DownloadSource(url))
        df = csv_read(file; copycols = true)

        newcolnames = [
            if all(c -> c ∈ "0123456789/", string(colname))
                date_value = Date(string(colname), "m/d/y") + Year(2000)
                Symbol(string(date_value))
            else
                CSSE_COLUMNS[Symbol(colname)]
            end for colname ∈ names(df)
        ]
        rename!(df, newcolnames)
        df[!, :country] = correct_countries!(df.country)
        result[dfsym] = df
    end
    result
end

function import_data(source::CsseSource; kwargs...)
    @debug "Importing CSSE data."
    raw_data = _import_data(source)
    key_columns = [:country, :state, :latitude, :longitude]
    dataframes = []

    for (main_col_name, df) ∈ raw_data
        dates = []
        date_syms = []
        for colsym ∈ names(df)
            colstr = string(colsym)
            if all(c -> c ∈ "0123456789-", colstr)
                push!(dates, Date(colstr))
                push!(date_syms, colsym)
            end
        end
        datadf = combine(groupby(df, key_columns)) do subdf
            col_values = [sum(subdf[!, date_sym]) for date_sym ∈ date_syms]
            DataFrame(; date = dates, main_col_name => col_values)
        end
        push!(dataframes, datadf)
    end
    parsedata = innerjoin(dataframes...; on = key_columns ∪ [:date]) |> DataFrame
    sort!(parsedata, [:country, :state, :date])
    DataDict(:world => parsedata)
end

"""
Number of tests per country from Our World in Data.
https://ourworldindata.org/covid-testing
"""
struct TotalTestsOWID <: AbstractDataSource
    function TotalTestsOWID()
        new()
    end
end
TotalTestsOWID(::Nothing; kwargs...) = TotalTestsOWID()

DATA_SOURCES[:total_tests] = TotalTestsOWID

const TESTING_OWID_URL = "https://github.com/owid/covid-19-data/raw/master/public/data/testing/covid-testing-all-observations.csv"

function import_data(source::TotalTestsOWID; url = TESTING_OWID_URL, kwargs...)
    @debug "Importing testing data from Our World In Data."
    file = import_data(DownloadSource(url))
    col_select = [:Entity, :Date, Symbol("Cumulative total")]
    data::DataFrame = csv_read(file; copycols = true, select = col_select)
    rename!(data, [:country, :date, :total_tests])
    country_split = map(name -> split(name, " - "), data.country)
    data[!, :country] = map(row -> row[1], country_split) |> correct_countries!
    test_kinds = map(row -> row[2], country_split) |> correct_test_kind!
    insertcols!(data, 3, :test_kind => test_kinds)
    sort!(data, [:country, :date])

    rows = collect(eachrow(data))
    nextrows = rows[2:end]
    lastrow = rows[end]
    nlastrow = DataFrame(;
        country = [lastrow.country],
        date = [today()],
        test_kind = [lastrow.test_kind],
        total_tests = [lastrow.total_tests],
    )
    push!(nextrows, nlastrow[1, :])
    data = copy(data)

    for (row, nextrow) ∈ zip(rows, nextrows)
        if row.country == nextrow.country
            dtend = nextrow.date - Day(1)
        else
            dtend = today() - Day(1)
        end

        for dt ∈ (row.date+Day(1)):Day(1):dtend
            newrow = (;
                country = row.country,
                date = dt,
                test_kind = row.test_kind,
                total_tests = row.total_tests,
            )
            push!(data, newrow)
        end
    end
    sort!(data, [:country, :test_kind, :date])
    data
end

# Brasil.io COVID-19 database

const BRASIL_IO_URL = "https://data.brasil.io/dataset/covid19/caso.csv.gz"

"""
Novel Coronavirus (COVID-19) Cases in Brazil, provided by Brasil.io
https://brasil.io/dataset/covid19/caso_full/
"""
struct BrasilIo <: AbstractDataSource
    function BrasilIo()
        new()
    end
end
BrasilIo(::Nothing; kwargs...) = BrasilIo()

DATA_SOURCES[:brasil_io] = BrasilIo

const BRAZIL_POP_FACTOR = 1.01

function import_data(source::BrasilIo; url = BRASIL_IO_URL, kwargs...)
    @debug "Importing data from Brasil.io."
    col_drop = [:confirmed_per_100k_inhabitants, :death_rate, :is_last]
    gzfile = import_data(DownloadSource(url))
    file = gzopen(string(gzfile))
    fstring = read(file)
    close(file)
    raw_data::DataFrame = csv_read(
        IOBuffer(fstring);
        copycols = true,
        truestrings = ["True"],
        falsestrings = ["False"],
        drop = col_drop,
    )
    raw_data[(raw_data.deaths.===missing), :deaths] .= 0
    raw_data[(raw_data.city.===missing), :city] .= ""
    sort!(raw_data, [:state, :city, :date])

    #rename!(raw_data, "estimated_population_2019" => "estimated_population")
    raw_data.estimated_population =
        round.(OptInt, raw_data.estimated_population .* BRAZIL_POP_FACTOR)

    estimates = get_input(
        "brazil_estimates_per_state";
        types = Dict(:infected_per_confirmed => Float),
    )
    est_sel = select(estimates, :state, :infected_per_confirmed)
    raw_data = leftjoin(raw_data, est_sel; on = :state)
    indexes = .!ismissing.(raw_data.infected_per_confirmed)
    ipc, conf = raw_data.infected_per_confirmed[indexes], raw_data.confirmed[indexes]
    avg_ipc = est_sel[est_sel.state.=="total", :infected_per_confirmed][1]
    raw_data.infected_per_confirmed[.!indexes] .= avg_ipc

    ipc = raw_data.infected_per_confirmed
    infected = round.(Int, ipc .* raw_data.confirmed)
    insertcols!(raw_data, :infected => infected)

    states = @where(raw_data, :place_type .== "state")
    cities = @where(raw_data, :place_type .== "city")

    select!(states, Not([:city, :place_type]))
    select!(cities, Not([:place_type]))

    total = select(states, :date, :confirmed, :deaths, :infected)
    total = combine(groupby(total, :date)) do subdf
        (;
            confirmed = sum(subdf.confirmed),
            deaths = sum(subdf.deaths),
            infected = sum(subdf.infected),
        )
    end

    DataDict(:states => states, :cities => cities, :total => total, :estimates => estimates)
end


text:
"""
# Module LinuxKernelPerfAPI

This module is an autogenerated set of constants and enums defined
mostly by C preprocessor macros in the Linux kernel header files.

These definitions are kernel version and architecture dependent.

Autogeneration from Linux 4.4.49 header files took place on Apr 23 2017 at 15:33:40.
"""
module LinuxKernelPerfAPI

const LINUX_VERSION_CODE = 263217
const LINUX_KERNEL_VERSION = v"4.4.49"

# syscall ids
const __NR_perf_event_open = 298
const __NR_getcpu = 309
const __NR_getpid = 39
const __NR_gettid = 186

@enum( perf_type_id::UInt32
     , PERF_TYPE_HARDWARE = 0
     , PERF_TYPE_SOFTWARE = 1
     , PERF_TYPE_TRACEPOINT = 2
     , PERF_TYPE_HW_CACHE = 3
     , PERF_TYPE_RAW = 4
     , PERF_TYPE_BREAKPOINT = 5
     , PERF_TYPE_MAX = 6
     )

@enum( perf_hw_id::UInt32
     , PERF_COUNT_HW_CPU_CYCLES = 0
     , PERF_COUNT_HW_INSTRUCTIONS = 1
     , PERF_COUNT_HW_CACHE_REFERENCES = 2
     , PERF_COUNT_HW_CACHE_MISSES = 3
     , PERF_COUNT_HW_BRANCH_INSTRUCTIONS = 4
     , PERF_COUNT_HW_BRANCH_MISSES = 5
     , PERF_COUNT_HW_BUS_CYCLES = 6
     , PERF_COUNT_HW_STALLED_CYCLES_FRONTEND = 7
     , PERF_COUNT_HW_STALLED_CYCLES_BACKEND = 8
     , PERF_COUNT_HW_REF_CPU_CYCLES = 9
     , PERF_COUNT_HW_MAX = 10
     )

@enum( perf_hw_cache_id::UInt32
     , PERF_COUNT_HW_CACHE_L1D = 0
     , PERF_COUNT_HW_CACHE_L1I = 1
     , PERF_COUNT_HW_CACHE_LL = 2
     , PERF_COUNT_HW_CACHE_DTLB = 3
     , PERF_COUNT_HW_CACHE_ITLB = 4
     , PERF_COUNT_HW_CACHE_BPU = 5
     , PERF_COUNT_HW_CACHE_NODE = 6
     , PERF_COUNT_HW_CACHE_MAX = 7
     )

@enum( perf_hw_cache_op_id::UInt32
     , PERF_COUNT_HW_CACHE_OP_READ = 0
     , PERF_COUNT_HW_CACHE_OP_WRITE = 1
     , PERF_COUNT_HW_CACHE_OP_PREFETCH = 2
     , PERF_COUNT_HW_CACHE_OP_MAX = 3
     )

@enum( perf_hw_cache_op_result_id::UInt32
     , PERF_COUNT_HW_CACHE_RESULT_ACCESS = 0
     , PERF_COUNT_HW_CACHE_RESULT_MISS = 1
     , PERF_COUNT_HW_CACHE_RESULT_MAX = 2
     )

@enum( perf_sw_ids::UInt32
     , PERF_COUNT_SW_CPU_CLOCK = 0
     , PERF_COUNT_SW_TASK_CLOCK = 1
     , PERF_COUNT_SW_PAGE_FAULTS = 2
     , PERF_COUNT_SW_CONTEXT_SWITCHES = 3
     , PERF_COUNT_SW_CPU_MIGRATIONS = 4
     , PERF_COUNT_SW_PAGE_FAULTS_MIN = 5
     , PERF_COUNT_SW_PAGE_FAULTS_MAJ = 6
     , PERF_COUNT_SW_ALIGNMENT_FAULTS = 7
     , PERF_COUNT_SW_EMULATION_FAULTS = 8
     , PERF_COUNT_SW_DUMMY = 9
     , PERF_COUNT_SW_BPF_OUTPUT = 10
     , PERF_COUNT_SW_MAX = 11
     )

@enum( perf_event_sample_format::UInt32
     , PERF_SAMPLE_IP = 1
     , PERF_SAMPLE_TID = 2
     , PERF_SAMPLE_TIME = 4
     , PERF_SAMPLE_ADDR = 8
     , PERF_SAMPLE_READ = 16
     , PERF_SAMPLE_CALLCHAIN = 32
     , PERF_SAMPLE_ID = 64
     , PERF_SAMPLE_CPU = 128
     , PERF_SAMPLE_PERIOD = 256
     , PERF_SAMPLE_STREAM_ID = 512
     , PERF_SAMPLE_RAW = 1024
     , PERF_SAMPLE_BRANCH_STACK = 2048
     , PERF_SAMPLE_REGS_USER = 4096
     , PERF_SAMPLE_STACK_USER = 8192
     , PERF_SAMPLE_WEIGHT = 16384
     , PERF_SAMPLE_DATA_SRC = 32768
     , PERF_SAMPLE_IDENTIFIER = 65536
     , PERF_SAMPLE_TRANSACTION = 131072
     , PERF_SAMPLE_REGS_INTR = 262144
     , PERF_SAMPLE_MAX = 524288
     )

@enum( perf_branch_sample_type_shift::UInt32
     , PERF_SAMPLE_BRANCH_USER_SHIFT = 0
     , PERF_SAMPLE_BRANCH_KERNEL_SHIFT = 1
     , PERF_SAMPLE_BRANCH_HV_SHIFT = 2
     , PERF_SAMPLE_BRANCH_ANY_SHIFT = 3
     , PERF_SAMPLE_BRANCH_ANY_CALL_SHIFT = 4
     , PERF_SAMPLE_BRANCH_ANY_RETURN_SHIFT = 5
     , PERF_SAMPLE_BRANCH_IND_CALL_SHIFT = 6
     , PERF_SAMPLE_BRANCH_ABORT_TX_SHIFT = 7
     , PERF_SAMPLE_BRANCH_IN_TX_SHIFT = 8
     , PERF_SAMPLE_BRANCH_NO_TX_SHIFT = 9
     , PERF_SAMPLE_BRANCH_COND_SHIFT = 10
     , PERF_SAMPLE_BRANCH_CALL_STACK_SHIFT = 11
     , PERF_SAMPLE_BRANCH_IND_JUMP_SHIFT = 12
     , PERF_SAMPLE_BRANCH_CALL_SHIFT = 13
     , PERF_SAMPLE_BRANCH_MAX_SHIFT = 14
     )

@enum( perf_branch_sample_type::UInt32
     , PERF_SAMPLE_BRANCH_USER = 1
     , PERF_SAMPLE_BRANCH_KERNEL = 2
     , PERF_SAMPLE_BRANCH_HV = 4
     , PERF_SAMPLE_BRANCH_ANY = 8
     , PERF_SAMPLE_BRANCH_ANY_CALL = 16
     , PERF_SAMPLE_BRANCH_ANY_RETURN = 32
     , PERF_SAMPLE_BRANCH_IND_CALL = 64
     , PERF_SAMPLE_BRANCH_ABORT_TX = 128
     , PERF_SAMPLE_BRANCH_IN_TX = 256
     , PERF_SAMPLE_BRANCH_NO_TX = 512
     , PERF_SAMPLE_BRANCH_COND = 1024
     , PERF_SAMPLE_BRANCH_CALL_STACK = 2048
     , PERF_SAMPLE_BRANCH_IND_JUMP = 4096
     , PERF_SAMPLE_BRANCH_CALL = 8192
     , PERF_SAMPLE_BRANCH_MAX = 16384
     )

@enum( perf_sample_regs_abi::UInt32
     , PERF_SAMPLE_REGS_ABI_NONE = 0
     , PERF_SAMPLE_REGS_ABI_32 = 1
     , PERF_SAMPLE_REGS_ABI_64 = 2
     )

# this is in the original header files _not_ an enum !
const PERF_TXN_ABORT_MASK_size = 8
@enum( perf_txn::UInt64
     , PERF_TXN_ELISION = 1
     , PERF_TXN_TRANSACTION = 2
     , PERF_TXN_SYNC = 4
     , PERF_TXN_ASYNC = 8
     , PERF_TXN_RETRY = 16
     , PERF_TXN_CONFLICT = 32
     , PERF_TXN_CAPACITY_WRITE = 64
     , PERF_TXN_CAPACITY_READ = 128
     , PERF_TXN_MAX = 256
     , PERF_TXN_ABORT_MASK = 18446744069414584320
     , PERF_TXN_ABORT_SHIFT = 32
     )

@enum( perf_event_read_format::UInt32
     , PERF_FORMAT_TOTAL_TIME_ENABLED = 1
     , PERF_FORMAT_TOTAL_TIME_RUNNING = 2
     , PERF_FORMAT_ID = 4
     , PERF_FORMAT_GROUP = 8
     , PERF_FORMAT_MAX = 16
     )

# sizes of `perf_event_attr` for various kernel API versions
const PERF_ATTR_SIZE_VER0 = 64
const PERF_ATTR_SIZE_VER1 = 72
const PERF_ATTR_SIZE_VER2 = 80
const PERF_ATTR_SIZE_VER3 = 96
const PERF_ATTR_SIZE_VER4 = 104
const PERF_ATTR_SIZE_VER5 = 112

# size of current implementation of `perf_event_attr`
const perf_event_attr_size = 112

# `ioctl` ids applicable to perf event file descriptors
const PERF_EVENT_IOC_ENABLE = 9216
const PERF_EVENT_IOC_DISABLE = 9217
const PERF_EVENT_IOC_REFRESH = 9218
const PERF_EVENT_IOC_RESET = 9219
const PERF_EVENT_IOC_PERIOD = 1074275332
const PERF_EVENT_IOC_SET_OUTPUT = 9221
const PERF_EVENT_IOC_SET_FILTER = 1074275334
const PERF_EVENT_IOC_ID = 2148017159
const PERF_EVENT_IOC_SET_BPF = 1074013192

@enum( perf_event_ioc_flags::UInt32
     , PERF_IOC_FLAG_GROUP = 1
     )

const PERF_RECORD_MISC_CPUMODE_MASK = 7
const PERF_RECORD_MISC_CPUMODE_UNKNOWN = 0
const PERF_RECORD_MISC_KERNEL = 1
const PERF_RECORD_MISC_USER = 2
const PERF_RECORD_MISC_HYPERVISOR = 3
const PERF_RECORD_MISC_GUEST_KERNEL = 4
const PERF_RECORD_MISC_GUEST_USER = 5
const PERF_RECORD_MISC_PROC_MAP_PARSE_TIMEOUT = 4096
const PERF_RECORD_MISC_MMAP_DATA = 8192
const PERF_RECORD_MISC_COMM_EXEC = 8192
const PERF_RECORD_MISC_SWITCH_OUT = 8192
const PERF_RECORD_MISC_EXACT_IP = 16384
const PERF_RECORD_MISC_EXT_RESERVED = 32768
@enum( perf_event_type::UInt32
     , PERF_RECORD_MMAP = 1
     , PERF_RECORD_LOST = 2
     , PERF_RECORD_COMM = 3
     , PERF_RECORD_EXIT = 4
     , PERF_RECORD_THROTTLE = 5
     , PERF_RECORD_UNTHROTTLE = 6
     , PERF_RECORD_FORK = 7
     , PERF_RECORD_READ = 8
     , PERF_RECORD_SAMPLE = 9
     , PERF_RECORD_MMAP2 = 10
     , PERF_RECORD_AUX = 11
     , PERF_RECORD_ITRACE_START = 12
     , PERF_RECORD_LOST_SAMPLES = 13
     , PERF_RECORD_SWITCH = 14
     , PERF_RECORD_SWITCH_CPU_WIDE = 15
     , PERF_RECORD_MAX = 16
     )

const PERF_MAX_STACK_DEPTH = 127

@enum( perf_callchain_context::UInt64
     , PERF_CONTEXT_HV = 18446744073709551584
     , PERF_CONTEXT_KERNEL = 18446744073709551488
     , PERF_CONTEXT_USER = 18446744073709551104
     , PERF_CONTEXT_GUEST = 18446744073709549568
     , PERF_CONTEXT_GUEST_KERNEL = 18446744073709549440
     , PERF_CONTEXT_GUEST_USER = 18446744073709549056
     , PERF_CONTEXT_MAX = 18446744073709547521
     )

const PERF_AUX_FLAG_TRUNCATED = 1
const PERF_AUX_FLAG_OVERWRITE = 2

const PERF_FLAG_FD_NO_GROUP = 1
const PERF_FLAG_FD_OUTPUT = 2
const PERF_FLAG_PID_CGROUP = 4
const PERF_FLAG_FD_CLOEXEC = 8

const PERF_MEM_OP_NA = 1
const PERF_MEM_OP_LOAD = 2
const PERF_MEM_OP_STORE = 4
const PERF_MEM_OP_PFETCH = 8
const PERF_MEM_OP_EXEC = 16
const PERF_MEM_OP_SHIFT = 0

const PERF_MEM_LVL_NA = 1
const PERF_MEM_LVL_HIT = 2
const PERF_MEM_LVL_MISS = 4
const PERF_MEM_LVL_L1 = 8
const PERF_MEM_LVL_LFB = 16
const PERF_MEM_LVL_L2 = 32
const PERF_MEM_LVL_L3 = 64
const PERF_MEM_LVL_LOC_RAM = 128
const PERF_MEM_LVL_REM_RAM1 = 256
const PERF_MEM_LVL_REM_RAM2 = 512
const PERF_MEM_LVL_REM_CCE1 = 1024
const PERF_MEM_LVL_REM_CCE2 = 2048
const PERF_MEM_LVL_IO = 4096
const PERF_MEM_LVL_UNC = 8192
const PERF_MEM_LVL_SHIFT = 5

const PERF_MEM_SNOOP_NA = 1
const PERF_MEM_SNOOP_NONE = 2
const PERF_MEM_SNOOP_HIT = 4
const PERF_MEM_SNOOP_MISS = 8
const PERF_MEM_SNOOP_HITM = 16
const PERF_MEM_SNOOP_SHIFT = 19

const PERF_MEM_LOCK_NA = 1
const PERF_MEM_LOCK_LOCKED = 2
const PERF_MEM_LOCK_SHIFT = 24

const PERF_MEM_TLB_NA = 1
const PERF_MEM_TLB_HIT = 2
const PERF_MEM_TLB_MISS = 4
const PERF_MEM_TLB_L1 = 8
const PERF_MEM_TLB_L2 = 16
const PERF_MEM_TLB_WK = 32
const PERF_MEM_TLB_OS = 64
const PERF_MEM_TLB_SHIFT = 26

# equivalent type for `pid_t`
const PID_T = Int32

end  # autogenerated module LinuxKernelPerfAPI

ext: jl


text:
using StaticArrays
using ConstrainedDynamics

# Parameters
ex = [1.;0.;0.]
ey = [0.;1.;0.]
ez = [0.;0.;1.]

l1 = 2.0
x, y = .1, .1
r = l1 / 4
h = r / 10
box = Box(x, l1, y, l1, color = RGBA(1., 1., 0.))
cyl = Cylinder(r, h, 2 * r, color = RGBA(1., 0., 0.))

p0 = [0;l1 / 2;0]
p1 = -p0
p2 = [0;r;0]
p3 = -p2

# Initial orientation
q1 = Quaternion(RotY(pi / 2))

# Links
origin = Origin{Float64}()
link1 = Body(box)
link2 = Body(cyl)

# Constraints
joint0to12 = EqualityConstraint(CylindricalFree(origin, link1, zeros(3), p0, ey), Revolute(origin, link2, p3 + p1, zeros(3), ez))
joint1to2 = EqualityConstraint(Cylindrical(link1, link2, p1, p3, ez))

links = [link1; link2]
constraints = [joint0to12;joint1to2]
shapes = [box,cyl]

mech = Mechanism(origin, links, constraints, shapes = shapes)
setPosition!(mech,link1,x = -p0)
setPosition!(mech,link2,x = p3 + p1)

a1 = nothing # SVector{0,Float64}() # SVector{1,Float64}(0)
a2 = nothing # SVector{0,Float64}() # SVector{3,Float64}(0,0,0)
a3 = nothing # SVector{0,Float64}() # SVector{0,Float64}()
a4 = 0.3 # SVector{1,Float64}(0.3)

a1 = SVector{0,Float64}() # SVector{1,Float64}(0)
a2 = SVector{0,Float64}() # SVector{3,Float64}(0,0,0)
a3 = SVector{0,Float64}() # SVector{0,Float64}()
a4 = 0.3 # SVector{1,Float64}(0.3)

a = [[a1];[a2];[a3];[a4]]
# a = [a1;a2;a3;a4]

setForce!(a,joint0to12,mech)
# setForce!(mech,link2,τ=[0;0;0.3])

simulate!(mech,save = true)
visualize!(mech)


text:
function num_divisors(n::Int64)
	sqrtn::Int64=floor(Int64,sqrt(n))
	result=0

	for i=1:(sqrtn+1)
		 if n%i==0
			result+=2
		end
	end

	if sqrtn^2==n
		result-=1
	end
    return result
end
"""
#print( num_divisors(100))
"""


text:
for sense in (MOI.MIN_SENSE, MOI.MAX_SENSE), offset in (0, 6 // 5)
    o = iszero(offset) ? "" : "_offset"
    f_unbd_name = Symbol("test_unbounded_$(Symbol(sense))$(o)")
    f_infeas_name = Symbol("test_infeasible_$(Symbol(sense))$(o)")
    f_infeas_affine_name = Symbol("test_infeasible_affine_$(Symbol(sense))$(o)")
    @eval begin
        """
        Test the infeasibility certificates of an unbounded model.
        """
        function $(f_unbd_name)(
            model::MOI.ModelLike,
            config::Config{T},
        ) where {T}
            @requires _supports(config, MOI.optimize!)
            x = MOI.add_variable(model)
            MOI.set(model, MOI.ObjectiveSense(), $sense)
            f = T(11 // 5) * x + T($offset)
            MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f)
            c = if $sense == MOI.MIN_SENSE
                MOI.add_constraint(
                    model,
                    T(13 // 10) * x,
                    MOI.LessThan(T(11 // 10)),
                )
            else
                MOI.add_constraint(
                    model,
                    T(13 // 10) * x,
                    MOI.GreaterThan(T(11 // 10)),
                )
            end
            MOI.optimize!(model)
            @requires(
                MOI.get(model, MOI.TerminationStatus()) == MOI.DUAL_INFEASIBLE,
            )
            @requires(
                MOI.get(model, MOI.PrimalStatus()) ==
                MOI.INFEASIBILITY_CERTIFICATE,
            )
            d = MOI.get(model, MOI.VariablePrimal(), x)
            obj = MOI.get(model, MOI.ObjectiveValue())
            if $sense == MOI.MIN_SENSE
                @test d < -config.atol
                @test obj < -config.atol
            else
                @test d > config.atol
                @test obj > config.atol
            end
            @test isapprox(T(11 // 5) * d, obj, config)
            if _supports(config, MOI.ConstraintPrimal)
                Ad = MOI.get(model, MOI.ConstraintPrimal(), c)
                @test isapprox(T(13 // 10) * d, Ad, config)
            end
            return
        end
        version_added(::typeof($(f_unbd_name))) = v"0.10.6"
        function setup_test(
            ::typeof($(f_unbd_name)),
            mock::MOIU.MockOptimizer,
            config::Config{T},
        ) where {T}
            MOIU.set_mock_optimize!(
                mock,
                (mock::MOIU.MockOptimizer) -> begin
                    MOIU.mock_optimize!(
                        mock,
                        MOI.DUAL_INFEASIBLE,
                        (
                            MOI.INFEASIBILITY_CERTIFICATE,
                            [$sense == MOI.MIN_SENSE ? T(-1) : T(1)],
                        ),
                        MOI.NO_SOLUTION,
                    )
                end,
            )
            return
        end

        """
        Test the infeasibility certificates of an infeasible model.
        """
        function $(f_infeas_name)(
            model::MOI.ModelLike,
            config::Config{T},
        ) where {T}
            @requires _supports(config, MOI.optimize!)
            x = MOI.add_variable(model)
            MOI.set(model, MOI.ObjectiveSense(), $sense)
            f = T(11 // 5) * x + T($offset)
            MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f)
            cl = MOI.add_constraint(model, x, MOI.GreaterThan(T(5 // 2)))
            cu = MOI.add_constraint(model, T(1) * x, MOI.LessThan(T(7 // 5)))
            MOI.optimize!(model)
            @requires(
                MOI.get(model, MOI.TerminationStatus()) ==
                config.infeasible_status,
            )
            @requires(
                MOI.get(model, MOI.DualStatus()) ==
                MOI.INFEASIBILITY_CERTIFICATE,
            )
            dl = MOI.get(model, MOI.ConstraintDual(), cl)
            du = MOI.get(model, MOI.ConstraintDual(), cu)
            @test dl > config.atol
            @test du < -config.atol
            @test isapprox(dl + du, T(0), config)
            if _supports(config, MOI.DualObjectiveValue)
                obj = MOI.get(model, MOI.DualObjectiveValue())
                target = T(7 // 5) * dl + T(5 // 2) * du
                if $sense == MOI.MIN_SENSE
                    @test obj > config.atol
                    @test isapprox(-target, obj, config)
                else
                    @test obj < -config.atol
                    @test isapprox(target, obj, config)
                end
            end
            return
        end
        version_added(::typeof($(f_infeas_name))) = v"0.10.6"
        function setup_test(
            ::typeof($(f_infeas_name)),
            mock::MOIU.MockOptimizer,
            config::Config{T},
        ) where {T}
            MOIU.set_mock_optimize!(
                mock,
                (mock::MOIU.MockOptimizer) -> begin
                    MOIU.mock_optimize!(
                        mock,
                        MOI.INFEASIBLE,
                        MOI.NO_SOLUTION,
                        MOI.INFEASIBILITY_CERTIFICATE,
                        (MOI.VariableIndex, MOI.GreaterThan{T}) => T[1],
                        (MOI.ScalarAffineFunction{T}, MOI.LessThan{T}) =>
                            T[-1],
                    )
                end,
            )
            return
        end

        """
        Test the infeasibility certificates of an infeasible model with affine
        constraints.
        """
        function $(f_infeas_affine_name)(
            model::MOI.ModelLike,
            config::Config{T},
        ) where {T}
            @requires _supports(config, MOI.optimize!)
            x = MOI.add_variable(model)
            MOI.set(model, MOI.ObjectiveSense(), $sense)
            f = T(11 // 5) * x + T($offset)
            MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f)
            cl = MOI.add_constraint(model, T(1) * x, MOI.GreaterThan(T(5 // 2)))
            cu = MOI.add_constraint(model, x, MOI.LessThan(T(7 // 5)))
            MOI.optimize!(model)
            @requires(
                MOI.get(model, MOI.TerminationStatus()) ==
                config.infeasible_status,
            )
            @requires(
                MOI.get(model, MOI.DualStatus()) ==
                MOI.INFEASIBILITY_CERTIFICATE,
            )
            dl = MOI.get(model, MOI.ConstraintDual(), cl)
            du = MOI.get(model, MOI.ConstraintDual(), cu)
            @test dl > config.atol
            @test du < -config.atol
            @test isapprox(dl + du, T(0), config)
            if _supports(config, MOI.DualObjectiveValue)
                obj = MOI.get(model, MOI.DualObjectiveValue())
                target = T(7 // 5) * dl + T(5 // 2) * du
                if $sense == MOI.MIN_SENSE
                    @test obj > config.atol
                    @test isapprox(-target, obj, config)
                else
                    @test obj < -config.atol
                    @test isapprox(target, obj, config)
                end
            end
            return
        end
        version_added(::typeof($(f_infeas_affine_name))) = v"0.10.6"
        function setup_test(
            ::typeof($(f_infeas_affine_name)),
            mock::MOIU.MockOptimizer,
            config::Config{T},
        ) where {T}
            MOIU.set_mock_optimize!(
                mock,
                (mock::MOIU.MockOptimizer) -> begin
                    MOIU.mock_optimize!(
                        mock,
                        MOI.INFEASIBLE,
                        MOI.NO_SOLUTION,
                        MOI.INFEASIBILITY_CERTIFICATE,
                        (MOI.ScalarAffineFunction{T}, MOI.GreaterThan{T}) =>
                            T[1],
                        (MOI.VariableIndex, MOI.LessThan{T}) => T[-1],
                    )
                end,
            )
            return
        end
    end
end


text:
using PyPlot

A = [1.0 2.0
     2.0 1.0]

evals, evecs = eig(A)
evecs =(evecs[:,1], evecs[:,2])

fig, ax = subplots()
for spine in ["left", "bottom"]
    ax[:spines][spine][:set_position]("zero")
end

for spine in ["right", "top"]
    ax[:spines][spine][:set_color]("none")
end
ax[:grid](alpha=0.4)

xmin, xmax = -3, 3
ymin, ymax = -3, 3
ax[:set_xlim](xmin, xmax)
ax[:set_ylim](ymin, ymax)

for v in evecs
    # Plot each eigenvector
    ax[:annotate](" ", xy=v, xytext=[0, 0],
                  arrowprops={"facecolor"=>"blue",
                              "shrink"=>0,
                              "alpha"=>0.6,
                              "width"=>0.5})

    # Plot the image of each eigenvector
    v = A * v
    ax[:annotate](" ", xy=v, xytext=[0, 0],
                  arrowprops={"facecolor"=>"red",
                              "shrink"=>0,
                              "alpha"=>0.6,
                              "width"=>0.5})
end

x = linspace(xmin, xmax, 3)
for v in evecs
    a = v[2] / v[1]
    ax[:plot](x, a .* x, "b-", lw=0.4)
end


text:
module MultiAgentPathFinding

## Dependencies

using DataFrames
using DataGraphs
using DataStructures
using FillArrays
using Graphs
using LinearAlgebra
using OffsetArrays
using ProgressMeter
using PythonCall
using Random
using Requires
using Statistics
using StatsBase
using SparseArrays
using UnicodePlots

## Includes

include("mapf.jl")

include("utils/eval_sol.jl")
include("utils/conflicts.jl")
include("utils/reservation.jl")
include("utils/vectorize.jl")

include("paths/dijkstra.jl")
include("paths/temporal_astar.jl")
include("paths/cooperative_astar.jl")
include("paths/independent_shortest_paths.jl")

include("exact_methods/conflict_based_search.jl")

include("local_search/neighborhood.jl")
include("local_search/large_neighborhood_search.jl")
include("local_search/feasibility_search.jl")
include("local_search/permutation_search.jl")

include("learning/features_agents.jl")
include("learning/features_edges.jl")

include("datasets/flatland/constants.jl")
include("datasets/flatland/agent.jl")
include("datasets/flatland/graph.jl")
include("datasets/flatland/utils.jl")
include("datasets/flatland/mapf.jl")

include("datasets/benchmarks/read.jl")
include("datasets/benchmarks/graph.jl")
include("datasets/benchmarks/mapf.jl")

## Exports

export Path, Solution
export MAPF, nb_agents

export flowtime, max_time
export VectorPriorityQueue
export find_conflict, conflict_exists, count_conflicts
export is_feasible
export path_to_vec, solution_to_mat, solution_to_mat2

export my_dijkstra_shortest_paths
export temporal_astar
export independent_astar, independent_dijkstra, independent_topological_sort
export compute_reservation
export cooperative_astar!, cooperative_astar

export conflict_based_search

export local_search_permutations
export feasibility_search!, feasibility_search
export large_neighborhood_search, large_neighborhood_search!

export agents_embedding
export edges_embedding

export flatland_mapf

export read_map, read_scenario
export GridGraph, shortest_path_grid
export benchmark_mapf

## Conditional dependencies

function __init__()
    @require GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a" begin
        using .GLMakie
        include("flatland/plot.jl")
        export plot_flatland_graph, flatland_agent_coords
    end
    @require Images = "916415d5-f1e6-5110-898d-aaa5f9f070e0" begin
        using .Images
        include("benchmarks/plot.jl")
        export display_map
    end
    @require JuMP = "4076af6c-e467-56ae-b986-b466b2749572" begin
        @require SCIP = "82193955-e24f-5292-bf16-6f2c5261a85f" begin
            using .JuMP, .SCIP
            include("exact_methods/linear_program.jl")
            export solve_lp
        end
    end
end

end # module


text:
module KalmanWrapper

using JavaCall

const AbstractKalmanFilter = @jimport org.hipparchus.filtering.kalman.AbstractKalmanFilter
const Class = @jimport java.lang.Class
const JString = @jimport java.lang.String
const KalmanFilter = @jimport org.hipparchus.filtering.kalman.KalmanFilter
const Measurement = @jimport org.hipparchus.filtering.kalman.Measurement
const Object = @jimport java.lang.Object
const ProcessEstimate = @jimport org.hipparchus.filtering.kalman.ProcessEstimate
const RealMatrix = @jimport org.hipparchus.linear.RealMatrix
const RealVector = @jimport org.hipparchus.linear.RealVector

export AbstractKalmanFilter
export KalmanFilter
export Measurement
export ProcessEstimate
export equals
export estimation_step
export get_class
export get_corrected
export get_covariance
export get_innovation_covariance
export get_kalman_gain
export get_measurement_jacobian
export get_predicted
export get_state
export get_state_transition_matrix
export get_time
export get_value
export hash_code
export notify
export notify_all
export to_string
export wait

include("ExtendedWrapper/ExtendedWrapper.jl")
include("LinearWrapper/LinearWrapper.jl")
include("abstract_kalman_filter.jl")
include("kalman_filter.jl")
include("measurement.jl")
include("process_estimate.jl")

end


text:
using LinearAlgebra
using SparseArrays
using Plots

"This code solves the 2D Laplace's equation with nonzero Dirichlet and Neumann
boundary conditions using a finite difference method."

m = 50 # number of points
x = LinRange(-1,1,m+2)
y = LinRange(-1,1,m+2)
xint = x[2:end-1]
yint = y[2:end-1]
h = x[2]-x[1]

T = 2.0 # final time
dt = .01 # timestep

id(i,j,m) = i + (j-1)*m # return global index from i,j
function FD_matrix_2D(h,m)
    A = spzeros(m*m,m*m)
    for i = 1:m, j = 1:m # loop thru x,y indices
        A[id(i,j,m),id(i,j,m)] = 4/h^2
        if i > 1 # avoids leftmost line of nodes
            A[id(i,j,m),id(i-1,j,m)] = -1/h^2 # x-derivative
        end
        if i < m # avoids rightmost line of nodes
            A[id(i,j,m),id(i+1,j,m)] = -1/h^2
        end
        if j > 1
            A[id(i,j,m),id(i,j-1,m)] = -1/h^2 # y-derivative
        end
        if j < m
            A[id(i,j,m),id(i,j+1,m)] = -1/h^2
        end
    end
    return A
end

A = FD_matrix_2D(h,m)

# initial condition
u0(x,y) = exp(-100*(x^2+y^2))
u0(x,y) = 0.0

# forcing
f(x,y,t) = exp(-25*((x-.25)^2 + (y-.75)^2))*exp(2*t)*(t < 1.5)
# f(x,y,t) = 0.0

F(t) = vec(f.(xint,yint',t)) # zero Dirichlet BCs
u = vec(u0.(xint,yint'))
interval = 1
Nsteps = ceil(Int,T / dt)
dt = T/Nsteps

Bsparse = (I + .5*dt*A)
#B = cholesky(I + .5*dt*A)
B = lu(I + .5*dt*A)

plot()
@gif for i = 1:Nsteps
    t = i*dt
    Favg = .5*(F(t) + F(t+dt))
    u .= B\(u + .5*dt*(Favg-A*u))
    if i%interval==0
        println("on timestep $i / $Nsteps")
        contourf(xint,yint,u,clims=(0,.5)) #vec(reshape(u,m,m)'))
    end
end every interval


text:
# A Timeline is a container for many tweens, which can be played either sequentially,
# or in parallel. A timeline can contain other timelines as well.
mutable struct TimeLine <: AbstractTween
    tween::Tween

    children::Array{AbstractTween,1}

    current::TimeLine
    parent::TimeLine

    mode::Int64
    built::Bool

    function TimeLine(type::Int64)
        o = new()
        o.tween = Tween()
        o.children = Array{AbstractTween,1}()
        o.mode = type
        o.built = false
        o
    end
end

# Called if the TimeLine is managing itself.
function build!(timeline::TimeLine)
    if timeline.built
        return
    end

    timeline.tween.base.duration = 0

    for child in timeline.children
        if child.tween.base.repeat_cnt < 0
            throw("You can't push an object with infinite repetitions in a TimeLine")
        end

        build!(child)

        if timeline.mode == TIMELINE_SEQUENCE
            t_delay = timeline.tween.base.duration
            timeline.tween.base.duration += full_duration(child.base)
            child.base.delay += t_delay
        elseif timeline.mode == TIMELINE_PARALLEL
            timeline.tween.base.duration = max(timeline.tween.base.duration, full_duration(child.tween.base))
        end
    end

    timeline.built = true;
end

function start!(timeline::TimeLine, man::TweenManager)
    build!(timeline)

    timeline.tween.base.current_time = 0
    timeline.tween.base.started = true

    for child in timeline.children
        start!(child, man)
    end
end

function free!(timeline::TimeLine)
    for child in timeline.children
        free!(child)
    end

    empty!(timeline.children)

    free!(timeline)
end

function contains_target_field(timeline::TimeLine, field_id::Int64)
    for child in timeline.children
        if child.tween.base.field_id == field_id
            return true
        end
    end

    false
end

function update_override!(timeline::TimeLine, step::Int64, lastStep::Int64, isIterationStep::Bool, delta::Float64)
    if !isIterationStep && step > lastStep
        dt = if is_reverse(lastStep)
            -delta - 1.0
        else
            delta + 1.0
        end

        for child in timeline.children
            update!(child, dt)
        end

        return
    end

    if !isIterationStep && step < lastStep
        dt = if is_reverse(lastStep)
            delta + 1.0
        else
            -delta - 1.0
        end

        # Interate in reverse
        for child in reverse!(timeline.children)
            update!(child, dt)
        end

        return
    end

    @assert(isIterationStep, "isIterationStep was false")

    if step > lastStep
        if is_reverse(step)
            force_end_values!(timeline)
            for child in reverse!(timeline.children)
                update!(child, delta)
            end
        else
            force_start_values!(timeline)
            for child in timeline.children
                update!(child, delta)
            end
        end
    elseif step < lastStep
        if is_reverse(step)
            force_start_values!(timeline)
            for child in reverse!(timeline.children)
                update!(child, delta)
            end
        else
            force_end_values!(timeline)
            for child in reverse!(timeline.children)
                update!(child, delta)
            end
        end
    else
        dt = if is_reverse(lastStep)
            -delta
        else
            delta
        end

        if delta >= 0.0
            for child in timeline.children
                update!(child, dt)
            end
        else
            for child in reverse!(timeline.children)
                update!(child, dt)
            end
        end
    end
end

function force_start_values!(timeline::TimeLine)
    for child in timeline.children
        force_to_start!(child.base, timeline)
    end
end

function force_end_values!(timeline::TimeLine)
    for child in timeline.children
        force_to_end!(child.base, timeline)
    end
end

max_stars_repo_name: wdevore/Ranger-Julia-SDL

max_issues_repo_name: wdevore/Ranger-Julia-SDL

max_forks_repo_name: wdevore/Ranger-Julia-SDL


text:
using Plots
using Random: randn;

@testset "plot wavelet transformed" begin
    X = randn(1024);
    wt = ContinuousWaveletTransform(CauchyWavelet(10), LinRange(2,200,64))
    wX = transform(wt, X)
    @test_nowarn contourf(wX)
end

@testset "plot coherence" begin
    N, M = 1024, 10;
    X = zeros(Float64, N, M)
    Y = zeros(Float64, N, M)
    for i in 1:M
        X[:,i] = map(t -> sin(2π*t/50), 1:N) + randn(N)
        Y[:,i] = map(t -> sin(2π*t/50), 1:N) + randn(N)
    end
    wt = ContinuousWaveletTransform(CauchyWavelet(10), LinRange(2,200,64))
    C = coherence(wt, X,Y)
    @test_nowarn contourf(C);
end


text:
# Exploration of Boos Stanski Chapter 1 - consulting example

# You are trying to model maximum flooding level Q
# while A is floodplane
# and k and η are parameters to be estimated

# literature suggests over 100 years the maximum level is:
Q(k, A, η) = k * A^(η - 1)

# Take the log to linearlize
logQ(k, A, η) = log(k) + (η - 1)*log(A)

# Verify
round(logQ(5, 100, .5), digits=5) == round(log(Q(5, 100, .5)), digits=5)

# The problem is that we don't actually see any observation which spans 100 years

# Example is of 140 gauging stations - 140 areas
# Stations have between 6 and 83 years of observations.

# The expected value E for subset of observations X₁ which is smaller than subset
# X₂ which is smaller than the full set X₁₀₀ is going to have the following property:
# E(max(X₁)) < E(max(X₂)) < E(max(X₁₀₀))

# What we would like to know is max(X₁₀₀)

# But the maximum of any of our observed values will be biased downwards since we
# never observe 100 years.

# Authors suggest using an estimate of the the joint likelihood that any observation
# will be below the median using the following identity allowing unobserved and observed
# X_i through X_100 to be ordered
# P(X_100 ≤ t) = P(X_1 ≤, ... X_100 ≤ t) = Π P(X_i < t) = [F(t)]^100
# by iid P(X_100 ≤ t) => [F(t)]^100

# P(X_100 ≤ t) = [F(t)]^100 = 1/2
# → P(X_100 ≤ t) = F(t₀) = (1/2)^.01 → t₀ = .993

# Using Extreme value function
# P(X_i <= t) = exp(-exp(-(t-μ)/σ))

Q993(μ, σ) = σ*(-log(-log(.993))) + μ

# The Authors suggest that we can estimate μ and σ then plug them into the function
# to get the solution.

# lets try it

using Random, Distributions

waterlevel = rand(Poisson(50), 100)
observe = rand(Binomial(1, .25), 100)

μhat = mean(waterlevel[observe .== 1])
σhat = var(waterlevel[observe .== 1])^.5

Q993(μhat, σhat)
maximum(waterlevel)


text:
module TypesTest

using MAFE, Test

function make_test_history(nMig)
    if nMig == 0
        return MigrationHistory(Vector{Int}());
    else
        migV = Vector{Migration}(undef, nMig);
        sy = 1970;
        ey = 1973;
        for i1 = 1 : nMig
            migV[i1] = Migration(sy, ey, 90394 + i1);
            sy = ey + 4;
            ey = sy + 3;
        end
        return MigrationHistory(migV)
    end
end

@testset "Types" begin
    mh = make_test_history(0);
    @test n_spells(mh) == 0
    syV, eyV = spell_years(mh);
    @test isempty(syV)
    @test isempty(eyV)
    @test spell_from_year_range(mh, 1970, 1980) == 0

    nMig = 4;
    mh = make_test_history(nMig);
    @test n_spells(mh) == 4
    syV, eyV = spell_years(mh);
    @test length(syV) == length(eyV) == nMig
    @test start_year(mh, 2) == syV[2]
    @test end_year(mh, 1) == eyV[1]

    @test spell_from_year_range(mh, 1940, 1950) == 0
    @test spell_from_year_range(mh, eyV[nMig], eyV[nMig]+2) == nMig + 1
    @test spell_from_year_range(mh, syV[2], eyV[2]) == 2
    @test spell_from_year_range(mh, syV[2] - 1, eyV[2]) < 0
    @test spell_from_year_range(mh, syV[2], eyV[2] + 1) < 0
end

end


text:
using LinearDynamicsModels
using DifferentialDynamicsModels
using LinearAlgebra
using StaticArrays
using Test

@maintain_type struct RoboState{T} <: FieldVector{2,T}
    x::T
    y::T
end

struct TestLinearDynamics{Dx,Du,TA<:StaticMatrix{Dx,Dx},TB<:StaticMatrix{Dx,Du},Tc<:StaticVector{Dx}} <: DifferentialDynamics
    A::TA
    B::TB
    c::Tc
end
(f::TestLinearDynamics{Dx,Du})(x::StaticVector{Dx}, u::StaticVector{Du}) where {Dx,Du} = f.A*x + f.B*u + f.c

Base.isapprox(x, y; atol=0) = all(isapprox(getfield(x, i), getfield(y, i), atol=atol) for i in 1:fieldcount(typeof(x)))

DI = DoubleIntegratorDynamics(2)
J  = TimePlusQuadraticControl(SMatrix{2,2,Rational{Int}}(I))
@test_throws ErrorException SteeringBVP(DI, J, compile=Val(true))
using SymPy
bvp_compiled = SteeringBVP(DI, J, compile=Val(true))
bvp_adhoc    = SteeringBVP(DI, J, compile=Val(false))

for T in (Float32, Float64)
    # LinearDynamics vs. Ad Hoc Single Integrator (from DifferentialDynamicsModels)
    SI_linear = NIntegratorDynamics(1, 2)
    SI_adhoc  = SingleIntegratorDynamics{2}()
    @test state_dim(SI_linear) == 2
    @test control_dim(SI_linear) == 2

    x0 = rand(RoboState{T})
    sc = StepControl(.5, rand(SVector{2,T}))
    rc = RampControl(.5, rand(SVector{2,T}), rand(SVector{2,T}))
    @test propagate(SI_linear, x0, sc) ≈ propagate(SI_adhoc, x0, sc)
    @test propagate(SI_linear, x0, sc) isa RoboState
    @test propagate(SI_linear, x0, rc) ≈ propagate(SI_adhoc, x0, rc)
    @test propagate(SI_linear, x0, rc) isa RoboState

    for dyn in (DoubleIntegratorDynamics(2),
                LinearDynamics(rand(SMatrix{4,4}), rand(SMatrix{4,2}), rand(SVector{4})))
        dyn_test = TestLinearDynamics(dyn.A, dyn.B, dyn.c)
        @test state_dim(dyn) == 4
        @test control_dim(dyn) == 2

        x0 = rand(SVector{4,T})
        u0 = sc.u
        # Exact vs. ODE Propagation
        @test propagate(dyn, x0, sc) ≈ propagate(dyn_test, x0, sc) atol=1e-3
        @test propagate(dyn, x0, rc) ≈ propagate(dyn_test, x0, rc) atol=1e-3

        # Exact vs. ODE Linearization (+ Reduction)
        X, U = SVector(1,3), SVector(1)
        for u_sc_or_rc in (u0, sc, rc)
            @test linearize(dyn, x0, u_sc_or_rc) ≈ linearize(dyn_test, x0, u_sc_or_rc) atol=1e-3
            @test linearize(dyn, x0, u_sc_or_rc, keep_control_dims=U) ≈
                  linearize(dyn_test, x0, u_sc_or_rc, keep_control_dims=U) atol=1e-3
            @test linearize(dyn, x0, u_sc_or_rc, keep_state_dims=X) ≈
                  linearize(dyn_test, x0, u_sc_or_rc, keep_state_dims=X) atol=1e-3
            @test linearize(dyn, x0, u_sc_or_rc, keep_state_dims=X, keep_control_dims=U) ≈
                  linearize(dyn_test, x0, u_sc_or_rc, keep_state_dims=X, keep_control_dims=U) atol=1e-3
        end
        @test propagate(dyn, x0, sc) ≈ linearize(dyn, x0, sc)(x0, sc)
        @test propagate(dyn, x0, rc) ≈ linearize(dyn, x0, rc)(x0, rc)
    end

    # Compiled vs Ad Hoc Steering
    x0 = rand(SVector{4,T})
    xf = rand(SVector{4,T})

    sol_compiled = bvp_compiled(x0, xf, 10.0)
    sol_adhoc    = bvp_adhoc(x0, xf, 10.0)
    @test sol_compiled.cost ≈ sol_adhoc.cost atol=1e-2
    @test duration(sol_compiled.controls) ≈ duration(sol_adhoc.controls) atol=1e-2
    @test waypoints(DI, x0, sol_compiled.controls, 10) ≈ waypoints(DI, x0, sol_adhoc.controls, 10) atol=1e-2
    @test instantaneous_control(sol_compiled.controls, duration(sol_compiled.controls)/2) ≈
          instantaneous_control(sol_adhoc.controls, duration(sol_adhoc.controls)/2) atol=1e-2
end


text:
type Bias <: Op; b; Bias(b::KUparam)=new(b); end

Bias(d1, d...; o...)=Bias(KUparam(d1, d...; o...))
Bias(; o...)=Bias(KUparam(0; o...))

params(l::Bias)=Any[l.b]
ninputs(::Bias)=1
ysize(::Bias,x)=size(x)
overwrites(::Bias)=true
back_reads_x(::Bias)=false
back_reads_y(::Bias)=false

function forw(l::Bias, x; y=x, o...)
    (b,x,y)=initforw(l,x,y;o...)
    biasforw(b,x,y)
    return y
end

function initforw(l::Bias, x, y; train=true, o...)
    issimilar(x,y) || error("Output mismatch")
    nb = size(x, ndims(x)==1 ? 1 : ndims(x)-1)
    if isempty(l.b)
        nz(l.b,:init,nothing) || (l.b.init = fill!; l.b.initp = 0)
        init(l.b, eltype(x), (nb,))
    end
    length(l.b) == nb || error("length mismatch")
    eltype(l.b) == eltype(x) || error("eltype mismatch")
    return ((!train && nz(l.b, :average, false)) ?
            (l.b.avg, x, y) :
            (l.b.arr, x, y))
end

function back(l::Bias, dy; dx=dy, incr=false, returndx=true, o...)
    # display((:back00,incr,vecnorm0(dy),vecnorm0(params(l))))
    initback(l, dy, incr)
    # display((:back01,incr,vecnorm0(dy),vecnorm0(params(l))))
    if incr
        biasback(l.b.inc, dy)
        axpy!(1, l.b.inc, l.b.diff)
    else
        biasback(l.b.diff, dy)
    end
    if returndx
        issimilar(dx,dy) || error("Gradient mismatch")
        (dx===dy ? dx : copy!(dx,dy))
    end
    # display((:back02,incr,vecnorm0(dy),vecnorm0(params(l))))
end

function initback(l::Bias, dy, incr)
    nb = size(dy, ndims(dy)==1 ? 1 : ndims(dy)-1)
    length(l.b) == nb || error("length mismatch")
    eltype(l.b) == eltype(dy) || error("eltype mismatch")
    similar!(l.b, :diff, l.b.arr; fill=0)
    incr && similar!(l.b, :inc, l.b.arr)
end

# We are implementing the CUDNN_ADD_SAME_C mode of cudnn:
# In this mode if x has dimensions (X1,X2,...,C,N) then
# bias has length=C.

biasforw(b::Vector, x::Vector, y::Vector=x)=(for i=1:length(y); y[i] = x[i] + b[i]; end; y)
biasback(db::Vector, dy::Vector)=(for i=1:length(dy); db[i]=dy[i]; end)

biasforw(b::Vector, x::Array, y::Array=x)=(c=ndims(x)-1; for i=1:length(y); y[i] = x[i] + b[ind2sub(size(x),i)[c]]; end; y)
biasback(db::Vector, dy::Array)=(c=ndims(dy)-1; fill!(db, zero(eltype(db))); for i=1:length(dy); db[ind2sub(size(dy),i)[c]] += dy[i]; end)

biasforw(b, x::KUdense, y::KUdense=x)=(biasforw(b,x.arr,y.arr); y)
biasback(db, dy::KUdense)=biasback(db, dy.arr)

GPU && (biasforw(b::CudaArray, x::CudaArray, y::CudaArray=x)=(y===x||copy!(y,x);cudnnAddTensor(b, y; mode=CUDNN_ADD_SAME_C)))
GPU && (biasback(db::CudaArray, dy::CudaArray)=cudnnConvolutionBackwardBias(dy, db))


text:
using ControlSystems
using DeterministicPolicyGradient
try
close("all")
catch
end
# const G = tf([1,-1],[1,0.01,1])
const G = tf(0.1,[1,0.01,0.1])

const T = 20
const t = 1:T
const m = 1
const n = 2
const p = 5
x0 = [0,1.]

# x,u = lsim(G, (i,s)->[-0.2s[2]],t,[0,1])[3:4]# Kolla att statesen är de förväntade
# lsimplot(G, (i,s)->[-0.2s[2]],t,[0,1])
# J(x,u)
# lsimplot(G, (i,s)->[0], t,x0), plot!(title="Before")
# gui()

# Cost function             ==========================================
const Q1 = diagm([5,1])
const Q2 = 1eye(m)
# r(s,a) = -sum(abs(Q1*s)) - sum(abs(Q2*a))
r(s,a,i) = (-(s)'Q1*(s) - a'Q2*a)[1]
L = lqr(ss(G),Q1,Q2)
μlqr(x) = -L*x

# Initialize solver options ==========================================
σβ      = 2
αΘ      = 0.005
αw      = 0.001
αv      = 0.001
αu      = 0.01
γ       = 0.999
τ       = 0.01
iters   = 30_000
critic_update = :rls # Choose between :rls and :gradient
λrls    = 0.99999
stepreduce_interval = 5000
stepreduce_factor = 0.995
hold_actor = 1000
opts = DPGopts(σβ,αΘ,αw,αv,αu,γ,τ,iters,m,critic_update,λrls,stepreduce_interval,stepreduce_factor,hold_actor)

# Initialize functions      ==========================================
cp = linspace(-5,5,p)
cv = linspace(-5,5,p)

grid1 = meshgrid(cp,cv)
const c1 = [grid1[1][:] grid1[2][:]]
P = size(c1,1)

function ϕ(s)
    a = exp(-1/(2*2)*sum((s'.-c1).^2,2))#exp(-1/2*(s.-c).^2)
    a ./= sum(a)
    a
    # [a; s[2]]
end

μ(s,Θ,t)        = Θ'ϕ(s)
∇μ(s)           = ϕ(s)
β(s,Θ,noise,i)  = Θ'ϕ(s) + noise[i]
ϕ(s,a,Θ)        = ∇μ(s)*(a-μ(s,Θ,t))
V(s,v)          = v'ϕ(s) # It's a good idea to bias V to some mean of the final landscape
Q(s,a,v,w,Θ,t)  = (ϕ(s,a,Θ)'w + V(s,v))[1]

function gradients(s1,s,a1,a,Θ,w,v,t)
    ∇μ = ϕ(s)
    ∇aQ = ∇μ'w
    ∇wQ = ∇μ*(a-Θ'∇μ)
    ∇vQ = ∇μ
    ∇aQ, ∇wQ, ∇vQ, ∇μ
end
simulate(Θ,x0, noise) = lsim(G, (t,s)->β(s,Θ,noise,t), t, x0)[3:4]
simulate(Θ,x0) = lsim(G, (t,s)->μ(s,Θ,t), t, x0)[3:4]
exploration(σβ) = filt(ones(5),[5],σβ*randn(T))
funs            = DPGfuns(μ,Q, gradients, simulate, exploration, r)

Θ               = zeros(P,m) # Weights
w               = 0.001randn(P)
v               = 0.001randn(P)
initial_state   = DPGstate(Θ,w,v)

# Train the actor and critic ==========================================
cost, Θ, w, v = dpg(opts, funs, initial_state, x0)

plot(cost[1:100:end])
gui()

lsimplot(G, (i,s)->[μ(s,Θ)], linspace(0,T,200),x0), plot!(title="After")
gui()
lsimplot(G, (i,s)->[μlqr(s)], linspace(0,T,200),x0), plot!(title="LQR")
gui()

Qpv_map = Float64[Q([sv,sp],0,v,w,Θ)[1] for sp in linspace(-5,5,100), sv in linspace(-5,5,100)]
PyPlot.matshow(Qpv_map); PyPlot.colorbar()
PyPlot.title("Q map");PyPlot.xlabel("Velocity");PyPlot.ylabel("Position")
PyPlot.yticks(0:10:100,-5:5);PyPlot.xticks(0:10:100,-5:5)

Qpa_map = Float64[Q([0,sp],a,v,w,Θ)[1] for sp in linspace(-5,5,100), a in linspace(-10,10,100)]
PyPlot.matshow(Qpa_map); PyPlot.colorbar()
PyPlot.title("Q map");PyPlot.xlabel("Control");PyPlot.ylabel("Position")
PyPlot.yticks(0:10:100,-5:5);PyPlot.xticks(0:10:100,-10:2:10)

control_map = Float64[μ([sv,sp],Θ)[1] for sp in linspace(-5,5,100), sv in linspace(-5,5,100)]
PyPlot.matshow(control_map); PyPlot.colorbar()
PyPlot.title("Control map");PyPlot.xlabel("Velocity");PyPlot.ylabel("Position")
PyPlot.yticks(0:10:100,-5:5);PyPlot.xticks(0:10:100,-5:5)

control_map = Float64[μlqr([sv,sp])[1] for sp in linspace(-5,5,100), sv in linspace(-5,5,100)]
PyPlot.matshow(control_map); PyPlot.colorbar()
PyPlot.title("Control map LQR");PyPlot.xlabel("Velocity");PyPlot.ylabel("Position")
PyPlot.yticks(0:10:100,-5:5);PyPlot.xticks(0:10:100,-5:5)

show()


text:
"""`bisquareWM(x, k, center, tol)`

Return the bisquare weighted mean, a location estimator.
The bisquare weighted mean of the data `x` with a k-value of `k`.
A sensible choice of `k` is 3 to 5 times the rms width or 1.3 to 2 times the
full width at half max of a peak.  For strictly Gaussian data, the choices of
k= 3.14, 3.88, and 4.68 times sigma will be 80%, 90%, and 95% efficient.
Data values a distance of more than `k` from the weighted mean are given no
weight.

`center` is used as an initial guess at the weighted mean; if omitted,
then the data median will be used.

The answer is found iteratively, revised until it changes by less than `tol`.  If
`tol` is omitted, then `tol` will use 1e-5 times the median absolute
deviation of `x` about its median."""
function bisquareWM(x::AbstractArray{T}, k::Real, center::Real, tol::Real) where {T <: Real}
    for _iteration in 1:100
        weights = (1.0-((x-center)/k).^2).^2
        weights[abs.(x-center).>k] = 0.0
        newcenter = dot(weights, x)/sum(weights)
        if abs(newcenter - center)<tol
            return newcenter
        end
        center = newcenter
    end
    error("bisquareWM used too many iterations.\nConsider using higher `tol` or better `center`, or change to trimean(x).")
end

bisquareWM(x::AbstractArray{T}, k::Real, center::Real) where {T<:Real} =
        bisquareWM(x, k, center, 1e-5*StatsBase.mad(float(x)))
bisquareWM(x::AbstractArray{T}, k::Real) where {T<:Real} = bisquareWM(x, k, median(x))

"""huberWM(x, k, center, tol)`

Return Huber's weighted mean, a location estimator.
Huber's weighted mean of the data `x` with a k-value of `k`.
A sensible choice of `k` is 1 to 1.5 times the rms width or 0.4 to 0.6 times the
full width at half max of a peak.  For strictly Gaussian data, the choices of
k=1.0 and 1.4 sigma give ...
Data values a distance of more than `k` from the weighted mean are given no
weight.

`center` is used as an initial guess at the weighted mean.
If `center` is None, then the data median will be used.

The answer is found iteratively, revised until it changes by less than `tol`.  If
`tol` is None (the default), then `tol` will use 1e-5 times the median absolute
deviation of `x` about its median."""
function huberWM(x::AbstractArray{T}, k::Real, center::Real, tol::Real) where {T <: Real}
    for _iteration = 1:100
        weights = float(k)./abs.(x-center)
        weights[weights .> 1.0] = 1.0
        newcenter = dot(weights, x)/sum(weights)
        if abs(newcenter - center) < tol
            return newcenter
        end
        center = newcenter
    end
    error("huberWM used too many iterations.\nConsider using higher `tol` or better `center`, or change to trimean(x).")
end

huberWM(x::AbstractArray{T}, k::Real, center::Real) where {T<:Real} =
        huberWM(x, k, center, 1e-5*StatsBase.mad(float(x)))
huberWM(x::AbstractArray{T}, k::Real) where {T<:Real} = huberWM(x, k, median(x))

"""`trimean(x)`

Return Tukey's trimean, a location estimator, for a data set `x`.

If (q1,q2,q3) are the quartiles (i.e., the 25%ile, median, and 75 %ile),
the trimean is (q1+q3)/4 + q2/2. Uses `idealf(x)` for finding quartiles.
"""
function trimean(x::AbstractArray{T}) where {T <: Real}
    mean([mean(idealf(x)), median(x)])
end


text:
immutable DropoutWork
    states
    statessize
    reservespace
    reservesize
end

function dropout(x, droprate)
    h = handle(x)
    y = similar(x)
    p = Ptr{Void}[0]
    cudnnCreateDropoutDescriptor(p)
    dropoutdesc = p[1]

    statessize_p = Cint[0]
    cudnnDropoutGetStatesSize(h, statessize_p)
    statessize = statessize_p[1]
    states = CuArray{Int8}(Int(statessize))

    xdesc = tensor_desc(x)
    ydesc = tensor_desc(y)
    reservesize_p = Cint[0]
    cudnnDropoutGetReserveSpaceSize(xdesc, reservesize_p)
    reservesize = reservesize_p[1]
    reservespace = CuArray{Int8}(Int(reservesize))

    states_p = Ptr{Void}[0]
    cudnnSetDropoutDescriptor(dropoutdesc, h, Cfloat(droprate), states, statessize, 0)
    cudnnDropoutForward(h, dropoutdesc, xdesc, x, ydesc, y, reservespace, reservesize)

    cudnnDestroyDropoutDescriptor(dropoutdesc)
    cudnnDestroyTensorDescriptor(xdesc)
    cudnnDestroyTensorDescriptor(ydesc)
    y, DropoutWork(states, statessize, reservespace, reservesize)
end

function ∇dropout!(dy, dx, droprate, work::DropoutWork)
    h = handle(dy)
    p = Ptr{Void}[0]
    cudnnCreateDropoutDescriptor(p)
    dropoutdesc = p[1]

    dydesc = tensor_desc(dy)
    dxdesc = tensor_desc(dx)
    cudnnSetDropoutDescriptor(dropoutdesc, h, Cfloat(droprate), work.states, work.statessize, 0)
    cudnnDropoutBackward(h, dropoutdesc, dydesc, dy, dxdesc, dx, work.reservespace, work.reservesize)

    cudnnDestroyDropoutDescriptor(dropoutdesc)
    cudnnDestroyTensorDescriptor(dydesc)
    cudnnDestroyTensorDescriptor(dxdesc)
    dx
end


text:
import gmsh

gmsh.initialize()

gmsh.model.add("t16")

gmsh.model.occ.addBox(0,0,0, 1,1,1, 1)
gmsh.model.occ.addBox(0,0,0, 0.5,0.5,0.5, 2)

gmsh.model.occ.cut([(3,1)], [(3,2)], 3)

x = 0; y = 0.75; z = 0; r = 0.09

holes = []
for t in 1:5
    global x, z
    x += 0.166
    z += 0.166
    gmsh.model.occ.addSphere(x,y,z,r, 3 + t)
    t = (3, 3 + t)
    push!(holes, t)
end

ov = gmsh.model.occ.fragment([(3,3)], holes)
gmsh.model.occ.synchronize()

lcar1 = .1
lcar2 = .0005
lcar3 = .055

ov = gmsh.model.getEntities(0);
gmsh.model.mesh.setSize(ov, lcar1);

ov = gmsh.model.getBoundary(holes, false, false, true);
gmsh.model.mesh.setSize(ov, lcar3);

eps = 1e-3
ov = gmsh.model.getEntitiesInBoundingBox(0.5-eps, 0.5-eps, 0.5-eps,
                                         0.5+eps, 0.5+eps, 0.5+eps, 0)
gmsh.model.mesh.setSize(ov, lcar2)

gmsh.model.mesh.generate(3)

gmsh.write("t16.msh")

gmsh.finalize()

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl


text:
module MaxwellBenchmarks

using BenchmarkTools

const SUITE = BenchmarkGroup()

include(joinpath(dirname(@__FILE__),"..","src","MaxwellSolver","PoissonFourier1D.jl"))
import .PoissonFourier1D

Δx = Dict()
ρ = Dict()
e = Dict()
solve! = Dict()

for N in [100, 1000, 10000]
    Δx[N] = 0.01
    ρ[N] = [sin(2π*j*Δx[N]) for j = 0:N]
    e[N] = similar(ρ[N])
    solve![N] = PoissonFourier1D.generate_solver(ρ[N], Δx[N])
    SUITE["poisson $N"] = @benchmarkable solve![$N]($(e[N]), $(ρ[N]))
end

end # module


text:
@testset "sinkhorn" begin
    a = Float64.(1:5)/sum(1:5)
    b = Float64.(1:10)/sum(1:10)
    m = zeros(5, 10)
    for i = 1:5
        for j = 1:10
            m[i,j] = 1/(i+j)
        end
    end
    reg = 1.0
    iter = 100
    tol = 1e-10
    loss = sinkhorn(a, b, m)
    @test run(sess, loss) ≈ 0.10269121
    # @test_nowarn loss = sinkhorn(a, b, m, method = "lp")
end

@testset "dist" begin
    a = rand(10,3)
    b = rand(20,3)
    for order in [1,2,3,4,5]
        M = zeros(10, 20)
        for i = 1:10
            for j = 1:20
                M[i,j] = norm(a[i,:] - b[j,:], order)
            end
        end
        m = ot_dist(a, b, order)
        @test m≈M
        m = ot_dist(a, constant(b), order)
        @test run(sess, m)≈M
    end
end

# @testset "fastdtw" begin
#     Sample = Float64[1,2,3,5,5,5,6]
#     Test = Float64[1,1,2,2,3,5]
#     u, p = dtw(Sample, Test, true)
#     @test run(sess, u)≈1.0
# end


text:
"""
    struct ScaledDiagonallyDominantBridge{T} <: MOI.Bridges.Variable.AbstractBridge
        side_dimension::Int
        variables::Vector{Vector{MOI.VariableIndex}}
        constraints::Vector{MOI.ConstraintIndex{
            MOI.VectorOfVariables, SOS.PositiveSemidefinite2x2ConeTriangle}}
    end

A matrix is SDD iff it is the sum of psd matrices Mij that are zero except
for entries ii, ij and jj [Lemma 9, AM17]. This bridge substitute the
constrained variables in [`SumOfSquares.ScaledDiagonallyDominantConeTriangle`](@ref)
into a sum of constrained variables in [`SumOfSquares.PositiveSemidefinite2x2ConeTriangle`](@ref).

[AM17] Ahmadi, A. A. & Majumdar, A.
*DSOS and SDSOS Optimization: More Tractable Alternatives to Sum of Squares and Semidefinite Optimization*
ArXiv e-prints, **2017**.
"""
struct ScaledDiagonallyDominantBridge{T} <: MOIB.Variable.AbstractBridge
    side_dimension::Int
    variables::Vector{Vector{MOI.VariableIndex}}
    constraints::Vector{MOI.ConstraintIndex{
        MOI.VectorOfVariables, SOS.PositiveSemidefinite2x2ConeTriangle}}
end

function MOIB.Variable.bridge_constrained_variable(
    ::Type{ScaledDiagonallyDominantBridge{T}}, model::MOI.ModelLike,
    s::SOS.ScaledDiagonallyDominantConeTriangle) where T
    n = s.side_dimension
    if n <= 1
        error("The bridges does not work with 1, `matrix_cone` should have returned `Nonnegatives` instead.")
    end
    N = MOI.dimension(MOI.PositiveSemidefiniteConeTriangle(n - 1))
    variables = Vector{Vector{MOI.VariableIndex}}(undef, N)
    constraints = Vector{MOI.ConstraintIndex{MOI.VectorOfVariables, SOS.PositiveSemidefinite2x2ConeTriangle}}(undef, N)
    k = 0
    for j in 1:n
        for i in 1:(j-1)
            k += 1
            # PSD constraints on 2x2 matrices are SOC representable
            variables[k], constraints[k] = MOI.add_constrained_variables(model, SOS.PositiveSemidefinite2x2ConeTriangle())
        end
    end
    return ScaledDiagonallyDominantBridge{T}(n, variables, constraints)
end

function MOIB.Variable.supports_constrained_variable(
    ::Type{<:ScaledDiagonallyDominantBridge},
    ::Type{SOS.ScaledDiagonallyDominantConeTriangle})
    return true
end
function MOIB.added_constrained_variable_types(
    ::Type{<:ScaledDiagonallyDominantBridge})
    return [(SOS.PositiveSemidefinite2x2ConeTriangle,)]
end
function MOIB.added_constraint_types(
    ::Type{<:ScaledDiagonallyDominantBridge})
    return Tuple{DataType, DataType}[]
end

# Attributes, Bridge acting as a model
function MOI.get(bridge::ScaledDiagonallyDominantBridge,
                 ::MOI.NumberOfVariables)
    return 3length(bridge.variables)
end
function MOI.get(bridge::ScaledDiagonallyDominantBridge,
                 ::MOI.ListOfVariableIndices)
    return Iterators.flatten(bridge.variables)
end
function MOI.get(bridge::ScaledDiagonallyDominantBridge,
                 ::MOI.NumberOfConstraints{MOI.VectorOfVariables, SOS.PositiveSemidefinite2x2ConeTriangle})
    return length(bridge.constraints)
end
function MOI.get(bridge::ScaledDiagonallyDominantBridge,
                 ::MOI.ListOfConstraintIndices{MOI.VectorOfVariables, SOS.PositiveSemidefinite2x2ConeTriangle})
    return bridge.constraints
end

# Indices
function MOI.delete(model::MOI.ModelLike,
                    bridge::ScaledDiagonallyDominantBridge)
    for variables in bridge.variables
        MOI.delete(model, variables)
    end
end

# Attributes, Bridge acting as a constraint

function MOI.get(::MOI.ModelLike, ::MOI.ConstraintSet,
                 bridge::ScaledDiagonallyDominantBridge)
    return SOS.ScaledDiagonallyDominantConeTriangle(bridge.side_dimension)
end

# TODO ConstraintPrimal, ConstraintDual

trimap(i, j) = div(j * (j - 1), 2) + i

function MOI.get(model::MOI.ModelLike, attr::MOI.VariablePrimal,
                 bridge::ScaledDiagonallyDominantBridge{T}, i::MOIB.IndexInVector) where T
    i, j = matrix_indices(i.value)
    if i == j
        value = zero(T)
        for k in 1:(i - 1)
            idx = offdiag_vector_index(k, i)
            value += MOI.get(model, attr, bridge.variables[idx][3])
        end
        for k in (i + 1):bridge.side_dimension
            idx = offdiag_vector_index(i, k)
            value += MOI.get(model, attr, bridge.variables[idx][1])
        end
        return value
    else
        idx = offdiag_vector_index(i, j)
        return MOI.get(model, attr, bridge.variables[idx][2])
    end
end

function MOIB.bridged_function(bridge::ScaledDiagonallyDominantBridge{T},
                               i::MOIB.IndexInVector) where T
    i, j = matrix_indices(i.value)
    if i == j
        func = zero(MOI.ScalarAffineFunction{T})
        for k in 1:(i - 1)
            idx = offdiag_vector_index(k, i)
            MOIU.operate!(+, T, func, bridge.variables[idx][3])
        end
        for k in (i + 1):bridge.side_dimension
            idx = offdiag_vector_index(i, k)
            MOIU.operate!(+, T, func, bridge.variables[idx][1])
        end
        return func
    else
        idx = offdiag_vector_index(i, j)
        return MOI.convert(MOI.ScalarAffineFunction{T}, bridge.variables[idx][2])
    end
end
function MOIB.Variable.unbridged_map(
    bridge::ScaledDiagonallyDominantBridge{T},
    vis::Vector{MOI.VariableIndex}) where T

    SAF = MOI.ScalarAffineFunction{T}
    umap = Pair{MOI.VariableIndex, SAF}[]
    k = 0
    z = zero(SAF)
    saf(i) = convert(SAF, vis[i])
    # vis[trimap(j, j)] is replaced by a sum of several variables.
    # The strategy is to replace all of them by zero except one.
    for j in 1:bridge.side_dimension
        for i in 1:(j-1)
            k += 1
            if i == 1 && j == 2
                push!(umap, bridge.variables[k][1] => saf(trimap(1, 1)))
            else
                push!(umap, bridge.variables[k][1] => z)
            end
            push!(umap, bridge.variables[k][2] => saf(trimap(i, j)))
            if i == 1
                push!(umap, bridge.variables[k][3] => saf(trimap(j, j)))
            else
                push!(umap, bridge.variables[k][3] => z)
            end
        end
    end
    return umap
end

max_stars_repo_name: blegat/SumOfSquares

max_issues_repo_name: blegat/SumOfSquares

max_forks_repo_name: blegat/SumOfSquares


text:
using Revise
using MccDaqHats
using Infiltrator
includet(joinpath(@__DIR__, "scan_utils.jl"))

#Constants
const READ_ALL_AVAILABLE = -1
const CURSOR_BACK_2 = "\x1b[2D"
const ERASE_TO_END_OF_LINE = "\x1b[0K"

"""
    function continuous_scan()
Example for continuous scan
MCC 172 Functions Demonstrated:
    mcc172_iepe_config_write
    mcc172_a_in_clock_config_write
    mcc172_a_in_clock_config_read
    mcc172_a_in_scan_start
    mcc172_a_in_scan_read
    mcc172_a_in_scan_stop
    mcc172_a_in_scan_cleanup

Purpose:
    Performa a continuous acquisition on 1 or more channels.

Description:
    Continuously acquires blocks of analog input data for a
    user-specified group of channels until the acquisition is
    stopped by the user.  The RMS voltage for each channel
    is displayed for each block of data received from the device.
"""
function continuous_scan()

    # Channel 1 is Force and Channel 2 is acceleration
    # can be passed as a parameter to the MCC 172 functions.
    channels = [0, 1]
    sensitivity = [0.16,  1.007489892]     # mV/N, mV/m/s^2
    unit = ["mV/N", "mV/m/s^2"]
    num_channels = length(channels)

    samples_per_channel = 0

    options = Set([:OPTS_CONTINUOUS])

    scan_rate = 10240.0
    # Turn on IEPE supply?
    iepe_enable = true

    try
        # Select an MCC 172 HAT device to use.
        hat = select_hat_devices(:MCC_172, 1)
        address = hat[1].address

        println("\nSelected MCC 172 HAT device at address $address")

        mcc172_open(address)
        for c in channels
            mcc172_iepe_config_write(address, c, iepe_enable)
            mcc172_a_in_sensitivity_write(address, c, sensitivity[c])
        end

        # Configure the clock and wait for sync to complete.
        mcc172_a_in_clock_config_write(address, :SOURCE_LOCAL, scan_rate)

        synced = false
        actual_scan_rate = 0.0  # initialize to Float64
        while !synced
            (_source_type, actual_scan_rate, synced) = mcc172_a_in_clock_config_read(address)
            if !synced
                sleep(0.005)
            end
        end

        println("    IEPE power: $(iepe_enable ? "on" : "off")")
        println("    Channels: $(join(channels, ", "))")
        println("    Sensitivity: $(sensitivity[1]) $(unit[1]) and $(sensitivity[1]) $(unit[1])")
        println("    Requested scan rate: $scan_rate")
        println("    Actual scan rate: $(round(actual_scan_rate, digits=3))")
        println("    Samples per channel $(round(samples_per_channel, digits=3))")
        println("    Options: $(options)")

        try
            println("\nPress 'Enter' to continue")
            readline()
        catch
            error("^C to end program")
        end

        # Configure and start the scan.
        # Since the continuous option is being used, the samples_per_channel
        # parameter is ignored if the value is less than the default internal
        # buffer size (10000 * num_channels in this case). If a larger internal
        # buffer size is desired, set the value of this parameter accordingly.
        # Configure and start the scan.
        channel_mask = chan_list_to_mask(channels)
        mcc172_a_in_scan_start(address, channel_mask, samples_per_channel, options)

        # Display the header row for the data table.
        println("Samples Read    Scan Count")
        for (chan, item) in enumerate(channels)
            print("       Channel $item")
        end
        println()

        try
            read_and_store_data(address, samples_per_channel, num_channels)

        catch e # KeyboardInterrupt
            if isa(e, InterruptException)
                # Clear the "^C" from the display.
                println("$CURSOR_BACK_2, $ERASE_TO_END_OF_LINE, Aborted\n")
            else
                println("\n $e")
            end
            mcc172_a_in_scan_stop(address)
        end

    catch err  # (HatError, ValueError) as err
        println("\n $err")
    finally
        # Turn off IEPE supply
        for channel in chans[i]
            mcc172_iepe_config_write(hat.address, channel, false)
        end

        mcc172_a_in_scan_cleanup(address)
end
    return nothing
end

"""
    function read_and_display_data(address, samples_per_channel, num_channels):
Reads data from the specified channels on the specified DAQ HAT devices
and updates the data on the terminal display.  The reads are executed in a
loop that continues until the user stops the scan or an overrun error is
detected.

Args:
hat (mcc172): The mcc172 HAT device object.
num_channels (int): The number of channels to display.

Returns:
Nothing
"""
function read_and_display_data(address::Integer, samples_per_channel::Integer, num_channels::Integer)
    total_samples_read = 0
    read_request_size = READ_ALL_AVAILABLE

    # When doing a continuous scan, the timeout value will be ignored in the
    # call to a_in_scan_read because we will be requesting that all available
    # samples (up to the default buffer size) be returned.
    timeout = 5.0

    # Read all of the available samples (up to the size of the read_buffer which
    # is specified by the user_buffer_size).  Since the read_request_size is set
    # to -1 (READ_ALL_AVAILABLE), this function returns immediately with
    # whatever samples are available (up to user_buffer_size) and the timeout
    # parameter is ignored.
    while true
        #sleep(0.001)
        #read_request_size = 10000
        #@show(address, read_request_size, num_channels, timeout)
        #println("Here I Am before")
        # @infiltrate
        resultcode, statuscode, result, samples_read =
        mcc172_a_in_scan_read(address, Int32(read_request_size), num_channels, timeout)
        #print("Here I Am after")

        # Check for an overrun error
        status = mcc172_status_decode(statuscode)
        if status.hardwareoverrun
            println("\n\nHardware overrun\n")
            break
        elseif status.bufferoverrun
            println("\n\nBuffer overrun\n")
            break
        end

        samples_read_per_channel = length(result) ÷ num_channels
        total_samples_read += samples_read_per_channel

        print("\r$samples_read_per_channel   $total_samples_read")

        # Display the RMS voltage for each channel.
        if samples_read_per_channel > 0
            for i in 1:num_channels
                #@show(length(result), i, num_channels, samples_read_per_channel)
                value = calc_rms(result, i, num_channels, samples_read_per_channel)
                print(" chan $i $(round(value, digits=5)) Vrms ")
            end
            #stdout.flush()

            sleep(0.1)
        end

        println("\n")
    end
end

max_stars_repo_name: Spectrum-Tec/MccDaqHats

max_issues_repo_name: Spectrum-Tec/MccDaqHats

max_forks_repo_name: Spectrum-Tec/MccDaqHats


text:
using Bobby
using Test
using Printf

Bb = Bobby
@testset "Bobby.jl" begin
    include("test_cvt.jl")
    include("test_bitboard.jl")
    include("test_king.jl")
    include("test_knight.jl")
    include("test_rook.jl")
    include("test_bishop.jl")
    include("test_magic.jl")
    include("test_move.jl")

    @printf("\n")
end


text:
const PAR_MAGIC = "PAR1"
const SZ_PAR_MAGIC = length(PAR_MAGIC)
const SZ_FOOTER = 4
const SZ_VALID_PAR = 2*SZ_PAR_MAGIC + SZ_FOOTER

# page is the unit of compression
mutable struct Page
    colchunk::ColumnChunk
    hdr::PageHeader
    pos::Int
    data::Vector{UInt8}
end

mutable struct PageLRU
    refs::Dict{ColumnChunk,DRef}
    function PageLRU()
        new(Dict{ColumnChunk,DRef}())
    end
end

function cacheget(lru::PageLRU, chunk::ColumnChunk, nf)
    if chunk in keys(lru.refs)
        poolget(lru.refs[chunk])
    else
        data = nf(chunk)
        lru.refs[chunk] = poolset(data)
        data
    end
end

# parquet file.
# Keeps a handle to the open file and the file metadata.
# Holds a LRU cache of raw bytes of the pages read.
mutable struct ParFile
    path::AbstractString
    handle::IOStream
    meta::FileMetaData
    schema::Schema
    page_cache::PageLRU
end

function ParFile(path::AbstractString)
    f = open(path)
    try
        return ParFile(path, f)
    catch ex
        close(f)
        rethrow(ex)
    end
end

function ParFile(path::AbstractString, handle::IOStream; maxcache::Integer=10)
    # TODO: maxcache should become a parameter to MemPool
    is_par_file(handle) || error("Not a parquet format file: $path")
    meta_len = metadata_length(handle)
    meta = metadata(handle, path, meta_len)
    ParFile(path, handle, meta, Schema(meta.schema), PageLRU())
end

##
# layer 1 access
# can access raw (uncompressed) bytes from pages

schema(par::ParFile) = par.schema
schema(conv::T, par::ParFile, schema_name::Symbol) where {T<:SchemaConverter} = schema(conv, par.schema, schema_name)

colname(col::ColumnChunk) = colname(col.meta_data)
colname(col::ColumnMetaData) = join(col.path_in_schema, '.')
colnames(rowgroup::RowGroup) = [colname(col) for col in rowgroup.columns]
function colnames(par::ParFile)
    s = Set{AbstractString}()
    for rg in rowgroups(par)
        push!(s, colnames(rg)...)
    end
    collect(s)
end

ncols(par::ParFile) = length(colnames(par))
nrows(par::ParFile) = par.meta.num_rows

coltype(col::ColumnChunk) = coltype(col.meta_data)
coltype(col::ColumnMetaData) = col._type

# return all rowgroups in the par file
rowgroups(par::ParFile) = par.meta.row_groups

# Return rowgroups that stores all the columns mentioned in `cnames`.
# Returned row groups can be further queried to get the range of rows.
rowgroups(par::ParFile, colname::AbstractString, rowrange::UnitRange=1:typemax(Int64)) = rowgroups(par, [colname], rowrange)
function rowgroups(par::ParFile, cnames, rowrange::UnitRange=1:typemax(Int64))
    R = RowGroup[]
    L = length(cnames)
    beginrow = 1
    for rowgrp in rowgroups(par)
        cnamesrg = colnames(rowgrp)
        found = length(intersect(cnames, cnamesrg))
        endrow = beginrow + rowgrp.num_rows - 1
        (found == L) && (length(intersect(beginrow:endrow)) > 0) && push!(R, rowgrp)
        beginrow = endrow + 1
    end
    R
end

columns(par::ParFile, rowgroupidx::Integer) = columns(par, rowgroups(par)[rowgroupidx])
columns(par::ParFile, rowgroup::RowGroup) = rowgroup.columns
columns(par::ParFile, rowgroup::RowGroup, colname::AbstractString) = columns(par, rowgroup, [colname])
function columns(par::ParFile, rowgroup::RowGroup, cnames)
    R = ColumnChunk[]
    for col in columns(par, rowgroup)
        (colname(col) in cnames) && push!(R, col)
    end
    R
end

function _pagevec(par::ParFile, col::ColumnChunk)
    # read pages from the file
    pos = page_offset(col)
    endpos = end_offset(col)
    io = par.handle
    pagevec = Page[]
    while pos < endpos
        seek(io, pos)
        pagehdr = read_thrift(io, PageHeader)

        buff = Array{UInt8}(undef, page_size(pagehdr))
        page_data_pos = position(io)
        data = read!(io, buff)
        page = Page(col, pagehdr, page_data_pos, data)
        push!(pagevec, page)
        pos = position(io)
    end
    pagevec
end
pages(par::ParFile, rowgroupidx::Integer, colidx::Integer) = pages(par, columns(par, rowgroupidx), colidx)
pages(par::ParFile, cols::Vector{ColumnChunk}, colidx::Integer) = pages(par, cols[colidx])
pages(par::ParFile, col::ColumnChunk) = cacheget(par.page_cache, col, col->_pagevec(par,col))

function bytes(page::Page, uncompressed::Bool=true)
    data = page.data
    codec = page.colchunk.meta_data.codec
    if uncompressed && (codec != CompressionCodec.UNCOMPRESSED)
        uncompressed_sz = page.hdr.uncompressed_page_size
        if codec == CompressionCodec.SNAPPY
            data = Snappy.uncompress(data)
        elseif codec == CompressionCodec.GZIP
            data = transcode(GzipDecompressor, data)
        else
            error("Unknown compression codec for column chunk: $codec")
        end
        (length(data) == uncompressed_sz) || error("failed to uncompress page. expected $(uncompressed_sz), got $(length(data)) bytes")
    end
    data
end

##
# layer 2 access
# can access decoded values from pages

map_dict_vals(valdict::Vector{T1}, vals::Vector{T2}) where {T1, T2} = isempty(valdict) ? vals : [valdict[v+1] for v in vals]

values(par::ParFile, rowgroupidx::Integer, colidx::Integer) = values(par, columns(par, rowgroupidx), colidx)
values(par::ParFile, cols::Vector{ColumnChunk}, colidx::Integer) = values(par, cols[colidx])
function values(par::ParFile, col::ColumnChunk)
    ctype = coltype(col)
    pgs = pages(par, col)
    valdict = Int[]
    jtype = PLAIN_JTYPES[ctype+1]
    if (ctype == _Type.BYTE_ARRAY) || (ctype == _Type.FIXED_LEN_BYTE_ARRAY)
        jtype = Vector{jtype}
    end
    vals = Array{jtype}(undef, 0)
    defn_levels = Int[]
    repn_levels = Int[]

    for pg in pgs
        typ = pg.hdr._type
        valtup = values(par, pg)
        if (typ == PageType.DATA_PAGE) || (typ == PageType.DATA_PAGE_V2)
            @debug("reading a data page for columnchunk")
            _vals, _defn_levels, _repn_levels = valtup
            enc, defn_enc, rep_enc = page_encodings(pg)
            if enc == Encoding.PLAIN_DICTIONARY || enc == Encoding.RLE_DICTIONARY
                append!(vals, map_dict_vals(valdict, _vals))
            else
                append!(vals, _vals)
            end
            append!(defn_levels, _defn_levels)
            append!(repn_levels, _repn_levels)
        elseif typ == PageType.DICTIONARY_PAGE
            _vals = valtup[1]
            @debug("read a dictionary page for columnchunk with $(length(_vals)) values")
            valdict = isempty(valdict) ? _vals : append!(valdict, _vals)
        end
    end
    vals, defn_levels, repn_levels
end

function read_levels(io::IO, max_val::Integer, enc::Int32, num_values::Integer)
    bw = bitwidth(max_val)
    (bw == 0) && (return Int[])
    @debug("reading levels. enc:$enc ($(Thrift.enumstr(Encoding,enc))), max_val:$max_val, num_values:$num_values")

    if enc == Encoding.RLE
        read_hybrid(io, num_values, bw)
    elseif enc == Encoding.BIT_PACKED
        read_bitpacked_run_old(io, num_values, bw)
    elseif enc == Encoding.PLAIN
        # levels should never be of this type though
        read_plain_values(io, _Type.INT32, num_values)
    else
        error("unsupported encoding $enc ($(Thrift.enumstr(Encoding,enc))) for levels")
    end
end

function read_values(io::IO, enc::Int32, typ::Int32, num_values::Integer)
    @debug("reading values. enc:$enc ($(Thrift.enumstr(Encoding,enc))), num_values:$num_values")

    if enc == Encoding.PLAIN
        read_plain_values(io, num_values, typ)
    elseif enc == Encoding.PLAIN_DICTIONARY || enc == Encoding.RLE_DICTIONARY
        read_rle_dict(io, num_values)
    else
        error("unsupported encoding $enc for pages")
        #@debug("unsupported encoding $enc ($(Thrift.enumstr(Encoding,enc))) for pages")
        #return Int[]
    end
end

function values(par::ParFile, page::Page)
    ctype = coltype(page.colchunk)
    rawbytes = bytes(page)
    io = IOBuffer(rawbytes)
    encs = page_encodings(page)
    num_values = page_num_values(page)
    typ = page.hdr._type

    if (typ == PageType.DATA_PAGE) || (typ == PageType.DATA_PAGE_V2)
        read_levels_and_values(io, encs, ctype, num_values, par, page)
    elseif typ == PageType.DICTIONARY_PAGE
        (read_plain_values(io, num_values, ctype),)
    else
        ()
    end
end

function read_levels_and_values(io::IO, encs::Tuple, ctype::Int32, num_values::Integer, par::ParFile, page::Page)
    cname = colname(page.colchunk)
    enc, defn_enc, rep_enc = encs

    #@debug("before reading defn levels bytesavailable in page: $(bytesavailable(io))")
    # read definition levels. skipped if column is required
    defn_levels = isrequired(par.schema, cname) ? Int[] : read_levels(io, max_definition_level(par.schema, cname), defn_enc, num_values)

    #@debug("before reading repn levels bytesavailable in page: $(bytesavailable(io))")
    # read repetition levels. skipped if all columns are at 1st level
    repn_levels = ('.' in cname) ? read_levels(io, max_repetition_level(par.schema, cname), rep_enc, num_values) : Int[]

    #@debug("before reading values bytesavailable in page: $(bytesavailable(io))")
    # read values
    vals = read_values(io, enc, ctype, num_values)

    vals, defn_levels, repn_levels
end

# column and page metadata
open(par::ParFile, col::ColumnChunk) = open(par.handle, par.path, col)
close(par::ParFile, col::ColumnChunk, io) = (par.handle == io) || close(io)
function open(io, path::AbstractString, col::ColumnChunk)
    if isfilled(col, :file_path)
        @debug("opening file $(col.file_path) to read column metadata at $(col.file_offset)")
        open(col.file_path)
    else
        if io == nothing
            @debug("opening file $path to read column metadata at $(col.file_offset)")
            open(path)
        else
            @debug("reading column metadata at $(col.file_offset)")
            io
        end
    end
end

function metadata(io, path::AbstractString, col::ColumnChunk)
    fio = open(io, path, col)
    seek(fio, col.file_offset)
    meta = read_thrift(fio, ColumnMetaData)
    (fio !== io) && close(fio)
    meta
end

function page_offset(col::ColumnChunk)
    colmeta = col.meta_data
    offset = colmeta.data_page_offset
    Thrift.isfilled(colmeta, :index_page_offset) && (offset = min(offset, colmeta.index_page_offset))
    Thrift.isfilled(colmeta, :dictionary_page_offset) && (offset = min(offset, colmeta.dictionary_page_offset))
    offset
end
end_offset(col::ColumnChunk) = page_offset(col) + col.meta_data.total_compressed_size

page_size(page::PageHeader) = Thrift.isfilled(page, :compressed_page_size) ? page.compressed_page_size : page.uncompressed_page_size

page_encodings(page::Page) = page_encodings(page.hdr)
function page_encodings(page::PageHeader)
    Thrift.isfilled(page, :data_page_header) ? page_encodings(page.data_page_header) :
    Thrift.isfilled(page, :data_page_header_v2) ? page_encodings(page.data_page_header_v2) :
    Thrift.isfilled(page, :dictionary_page_header) ? page_encodings(page.dictionary_page_header) : ()
end
page_encodings(page::DictionaryPageHeader) = (page.encoding,)
page_encodings(page::DataPageHeader) = (page.encoding, page.definition_level_encoding, page.repetition_level_encoding)
page_encodings(page::DataPageHeaderV2) = (page.encoding, Encoding.RLE, Encoding.RLE)

page_num_values(page::Page) = page_num_values(page.hdr)
function page_num_values(page::PageHeader)
    Thrift.isfilled(page, :data_page_header) ? page_num_values(page.data_page_header) :
    Thrift.isfilled(page, :data_page_header_v2) ? page_num_values(page.data_page_header_v2) :
    Thrift.isfilled(page, :dictionary_page_header) ? page_num_values(page.dictionary_page_header) : 0
end
page_num_values(page::Union{DataPageHeader,DataPageHeaderV2,DictionaryPageHeader}) = page.num_values

# file metadata
read_thrift(buff::Array{UInt8}, ::Type{T}) where {T} = read(TCompactProtocol(TMemoryTransport(buff)), T)
read_thrift(io::IO, ::Type{T}) where {T} = read(TCompactProtocol(TFileTransport(io)), T)
read_thrift(t::TR, ::Type{T}) where {TR<:TTransport,T} = read(TCompactProtocol(t), T)

function metadata_length(io)
    sz = filesize(io)
    seek(io, sz - SZ_PAR_MAGIC - SZ_FOOTER)

    # read footer size as little endian signed Int32
    ProtoBuf.read_fixed(io, Int32)
end

function metadata(io, path::AbstractString, len::Integer=metadata_length(io))
    @debug("metadata len = $len")
    sz = filesize(io)
    seek(io, sz - SZ_PAR_MAGIC - SZ_FOOTER - len)
    meta = read_thrift(io, FileMetaData)
    # go through all column chunks and read metadata from file offsets if required
    for grp in meta.row_groups
        for col in grp.columns
            if !isfilled(col, :meta_data)
                # need to read metadata from an offset
                col.meta_data = metadata(io, path, col)
            end
        end
    end
    meta
end

metadata(par::ParFile) = par.meta

# file format verification
function is_par_file(fname::AbstractString)
    open(fname) do io
        return is_par_file(io)
    end
end

function is_par_file(io)
    sz = filesize(io)
    (sz > SZ_VALID_PAR) || return false

    seekstart(io)
    magic = Array{UInt8}(undef, 4)
    read!(io, magic)
    (String(magic) == PAR_MAGIC) || return false

    seek(io, sz - SZ_PAR_MAGIC)
    magic = Array{UInt8}(undef, 4)
    read!(io, magic)
    (String(magic) == PAR_MAGIC) || return false

    true
end


text:
for c = 1:40, b = 1:c-1, a = 1:b-1
  if a * a + b * b == c * c
    println("$a, $b, $c")
  end
end

max_stars_repo_name: idrougge/ple

max_issues_repo_name: idrougge/ple

max_forks_repo_name: idrougge/ple


text:
abstract type Field end
abstract type Vortex end
abstract type CoreShape end
abstract type VortexGroup end

"""
    Torus <: Field

Specify toroidal boundary conditions: doubly periodic in 2D.

# Arguments
- `ψ::Array{Complex{Float64},2}`: wavefunction.
- `x::Vector{Float64}`: `x` coordinates.
- `y::Vector{Float64}`: `y` coordinates.

See also: [`Sphere`](@def), [`Field`](@ref)
"""
mutable struct Torus <: Field
    ψ::Array{Complex{Float64},2}
    x::Vector{Float64}
    y::Vector{Float64}
end

"""
    Sphere <: Field

Specify spherical boundary conditions in 2D.

# Arguments
- `ψ::Array{Complex{Float64},2}`: wavefunction.
- `x::Vector{Float64}`: `x` coordinates.
- `y::Vector{Float64}`: `y` coordinates.

See also: [`Torus`](@def), [`Field`](@ref)
"""
mutable struct Sphere <: Field
    ψ::Array{Complex{Float64},2}
    x::Vector{Float64}
    y::Vector{Float64}
end

mutable struct PointVortex <: Vortex
    xv::Float64
    yv::Float64
    qv::Int64
end

struct Ansatz <: CoreShape
    f::Interpolations.GriddedInterpolation{Float64,1,Float64,Gridded{Linear},Tuple{Array{Float64,1}}}
    ξ::Float64
    Λ::Float64
end

struct Exact <: CoreShape
    f::Interpolations.GriddedInterpolation{Float64,1,Float64,Gridded{Linear},Tuple{Array{Float64,1}}}
    ξ::Float64
end

mutable struct ScalarVortex{T<:CoreShape} <: Vortex
    core::T
    vort::PointVortex
end

mutable struct Dipole <: VortexGroup
    vp::PointVortex
    vn::PointVortex
    Dipole(v1,v2) = (charge(v1) > 0 && charge(v2) < 0) ? new(v1,v2) : new(v2,v1)
end

mutable struct Cluster <: VortexGroup
    vortices::Array{PointVortex,1}
    tree::Array{LightGraphs.SimpleGraphs.SimpleEdge{Int64},1}
end
Cluster(vort::Array{PointVortex,1}) = Cluster(vort,spanning_tree(xpos(vort),ypos(vort)))


text:
using SIMDPoly
using Test

@testset "complex poly" begin
    @test evalpolydiv(2.0+0.0im, (1.,2.,3.,4.,5.), (2.,4.,6.,8.,10.)) == evalpoly(2.0+0.0im, (1.,2.,3.,4.,5.))/evalpoly(2.0+0.0im+0.0im+0.0im, (2.,4.,6.,8.,10.))
end


text:
module Runtests

using Test

@testset "Test TinyModia with simulation" begin

    include("TestVariables.jl")
    include("TestFirstOrder.jl")
    include("TestFirstOrder2.jl")
    include("TestStateSelection.jl")
    include("TestFilterCircuit.jl")

    include("TestUnits.jl")
    include("TestUncertainties.jl")
    include("TestUnitsAndUncertainties.jl")

    include("TestTwoInertiasAndIdealGear.jl")
    include("TestTwoInertiasAndIdealGearWithUnits.jl")
    include("TestTwoInertiasAndIdealGearWithUnitsAndUncertainties.jl")
    include("TestTwoInertiasAndIdealGearWithMonteCarlo.jl")
    #include("TestTwoInertiasAndIdealGearWithUnitsAndMonteCarlo.jl")  # MonteCarlo and Unitful not yet supported

    include("TestSingularLRRL.jl")

    include("TestStateSpace.jl")
    include("TestHeatTransfer.jl")

    include("../examples/runexamples.jl")

end

end


text:
function radon_linear_robust(In, Par, itype,lp,imax);
h = Par[1];
p = Par[2];
dt = Par[3];
np = length(p);
nh = length(h);

if itype==1;
    m=In;
    (nt,np)=size(m);
    d=zeros(nt,nh);
end;

if itype == -1;
    d = In;
    (nt, nh) = size(d);
    m = zeros(nt,np);
end;
aux=zeros(np*2,1);
for itau=1:nt;
    for ip=1:np;
        k=1;
        aux=0*aux;
        for ih=1:nh;
            t=(itau-1)*dt+p[ip]*h[ih];
            it=t/dt+1;
            it1=floor(Int,it);
            it2=it1+1;
            a=it-it1;
            if it2<nt && it1>1
                if itype==1;
                    d[it1,ih]=d[it1,ih]+(1-a)*m[itau,ip];
                    d[it2,ih]=d[it2,ih]+(a)*m[itau,ip];
                else
                    aux[k,1]=(1-a)*d[it1,ih]+(a)*d[it2,ih];
                    k=k+1;
                end
            end
        end
        m[itau,ip]=LP_norm(ones(k-1,1),aux[1:k-1,1],lp,imax,1);
    end
end

if itype==1;
    Out=d;
end;
if itype == -1;
    Out = m;
end;
return Out
end

max_stars_repo_name: Ji-seismic/SeisTest

max_issues_repo_name: Ji-seismic/SeisTest

max_forks_repo_name: Ji-seismic/SeisTest


text:
module ModuleEstimation
using MultipleViewGeometry.ModuleTypes, MultipleViewGeometry.ModuleDataNormalization, MultipleViewGeometry.ModuleOperators, MultipleViewGeometry.ModuleMathAliases
using MultipleViewGeometry.ModuleMoments, MultipleViewGeometry.ModuleCarriers, MultipleViewGeometry.ModuleCostFunction
using MultipleViewGeometry.ModuleTransform, MultipleViewGeometry
using StaticArrays, LinearAlgebra
using LsqFit
export estimate
include("estimate_twoview.jl")
end


text:
using MultivariateStats
using LinearAlgebra
using Test

## testing data

function pwdists(X)
    S = sum(abs2, X, dims=1)
    D2 = S .+ S' - 2 * (X'X)
    D2[diagind(D2)] .= 0.0
    return sqrt.(D2)
end

X0 = randn(3, 5)
G0 = X0'X0
D0 = pwdists(X0)

## conversion between dmat and gram

@assert issymmetric(D0)
@assert issymmetric(G0)

D = gram2dmat(G0)
@test issymmetric(D)
@test D ≈ D0

G = dmat2gram(D0)
@test issymmetric(G)
@test gram2dmat(G) ≈ D0

## classical MDS

X = classical_mds(D0, 3)
@test pwdists(X) ≈ D0

#Test MDS embeddings in dimensions >= number of points
@test classical_mds([0 1; 1 0], 2, dowarn=false) == [-0.5 0.5; 0 0]
@test classical_mds([0 1; 1 0], 3, dowarn=false) == [-0.5 0.5; 0 0; 0 0]

#10 - dmat2gram produces negative definite matrix
let D = [1.0 0.5 0.3181818181818182 0.38095238095238093 0.6111111111111112 0.36363636363636365 0.3333333333333333 0.4444444444444444 0.45454545454545453 0.38095238095238093
        0.5 1.0 0.3333333333333333 0.4166666666666667 0.32 0.37037037037037035 0.35 0.38461538461538464 0.38461538461538464 0.4
        0.3181818181818182 0.3333333333333333 1.0 0.38095238095238093 0.125 0.22727272727272727 0.47058823529411764 0.35294117647058826 0.2608695652173913 0.21052631578947367
        0.38095238095238093 0.4166666666666667 0.38095238095238093 1.0 0.17391304347826086 0.2962962962962963 0.38095238095238093 0.45454545454545453 0.3076923076923077 0.5454545454545454
        0.6111111111111112 0.32 0.125 0.17391304347826086 1.0 0.5 0.15 0.391304347826087 0.38095238095238093 0.3333333333333333
        0.36363636363636365 0.37037037037037035 0.22727272727272727 0.2962962962962963 0.5 1.0 0.2857142857142857 0.44 0.4 0.3181818181818182
        0.3333333333333333 0.35 0.47058823529411764 0.38095238095238093 0.15 0.2857142857142857 1.0 0.42105263157894735 0.45 0.23529411764705882
        0.4444444444444444 0.38461538461538464 0.35294117647058826 0.45454545454545453 0.391304347826087 0.44 0.42105263157894735 1.0 0.5416666666666666 0.5555555555555556
        0.45454545454545453 0.38461538461538464 0.2608695652173913 0.3076923076923077 0.38095238095238093 0.4 0.45 0.5416666666666666 1.0 0.4
        0.38095238095238093 0.4 0.21052631578947367 0.5454545454545454 0.3333333333333333 0.3181818181818182 0.23529411764705882 0.5555555555555556 0.4 1.0],
    Xt = [-0.27529104101488666 0.006134513718202863 0.33298809606740326 0.2608994458893664 -0.46185275796909575 -0.23734315039370618 0.29972782027671513 0.03827901455843394 -0.04096713097883363 0.07742518984640051
        -0.08177061420820278 -0.0044504235228030225 -0.3271919093638943 0.28206254638779243 -0.0954706915166714 -0.07137742126520012 -0.30754933764853587 0.18582658369448027 -0.03715307349750036 0.45707434094053534]',
    M = classical_mds(D, 2)'
    @test M ≈ Xt .* cis.(angle.(sum(conj.(Xt) .* M, dims=1)))
end

#10 - test degenerate problem
@test classical_mds(zeros(10, 10), 3, dowarn=false) == zeros(3, 10)