text:
δmulticast(f, x::AbstractFloat) = f(x)

function δmulticast(f, x::AbstractArray)
    iter = Base.CartesianRange(size(x))
    temp = next(iter, start(iter))
    I = temp[1]
    st = temp[2]
    val = f(x[I])
    v = similar(Array{typeof(val[1])}, size(x))
    ∇v = similar(Array{typeof(val[2])}, size(x))
    v[I] = val[1]
    ∇v[I] = val[2]
    while !done(iter, st)
        I, st = next(iter, st)
        v[I], ∇v[I] = f(x[I])
    end
    v, u -> ∇multicast(∇v, u)
end

function ∇multicast(∇v, x)
    iter = Base.CartesianRange(size(x))
    st = start(iter)
    I, st = next(iter, st)
    val = ∇v[I](x[I])
    v = similar(Array{typeof(val)}, size(x))
    v[I] = val
    while !done(iter, st)
        I, st = next(iter, st)
        v[I] = ∇v[I](x[I])
    end
    v
end

δmulticast(f, x1::AbstractFloat, x2::AbstractFloat) = f(x1, x2)
δmulticast_const1(f, x1::Number, x2::AbstractFloat) = f(x1, x2)
δmulticast_const2(f, x1::AbstractFloat, x2::Number) = f(x1, x2)

function δmulticast{T<:AbstractArray}(f, x1::T, x2::T)
    iter = Base.CartesianRange(size(x1))
    temp = next(iter, start(iter))
    I = temp[1]
    st = temp[2]
    val = f(x1[I], x2[I])
    v = similar(Array{typeof(val[1])}, size(x1))
    ∇v = similar(Array{typeof(val[2])}, size(x1))
    v[I] = val[1]
    ∇v[I] = val[2]
    while !done(iter, st)
        I, st = next(iter, st)
        v[I], ∇v[I] = f(x1[I], x2[I])
    end
    v, u -> ∇multicast2(∇v, u)
end

function ∇multicast2(∇v, x)
    iter = Base.CartesianRange(size(x))
    st = start(iter)
    I, st = next(iter, st)
    val = ∇v[I](x[I])
    v1 = similar(Array{typeof(val[1])}, size(x))
    v2 = similar(Array{typeof(val[2])}, size(x))
    v1[I] = val[1]
    v2[I] = val[2]
    while !done(iter, st)
        I, st = next(iter, st)
        v1[I], v2[I] = ∇v[I](x[I])
    end
    v1, v2
end

function δmulticast_const1(f, x1::Number, x2::AbstractArray)
    iter = Base.CartesianRange(size(x2))
    temp = next(iter, start(iter))
    I = temp[1]
    st = temp[2]
    val = f(x1, x2[I])
    v = similar(Array{typeof(val[1])}, size(x2))
    ∇v = similar(Array{typeof(val[2])}, size(x2))
    v[I] = val[1]
    ∇v[I] = val[2]
    while !done(iter, st)
        I, st = next(iter, st)
        v[I], ∇v[I] = f(x1, x2[I])
    end
    v, u -> ∇multicast(∇v, u)
end

function δmulticast_const2(f, x1::AbstractArray, x2::Number)
    iter = Base.CartesianRange(size(x1))
    temp = next(iter, start(iter))
    I = temp[1]
    st = temp[2]
    val = f(x1[I], x2)
    v = similar(Array{typeof(val[1])}, size(x1))
    ∇v = similar(Array{typeof(val[2])}, size(x1))
    v[I] = val[1]
    ∇v[I] = val[2]
    while !done(iter, st)
        I, st = next(iter, st)
        v[I], ∇v[I] = f(x1[I], x2)
    end
    v, u -> ∇multicast(∇v, u)
end

function δmulticast(f, x1::AbstractVector, x2::AbstractMatrix)
    idx = length(x1) == size(x2, 1) ? 1 : 2
    iter = Base.CartesianRange(size(x2))
    temp = next(iter, start(iter))
    I = temp[1]
    st = temp[2]
    val = f(x1[I.I[idx]], x2[I])
    v = similar(Array{typeof(val[1])}, size(x2))
    ∇v = similar(Array{typeof(val[2])}, size(x2))
    v[I] = val[1]
    ∇v[I] = val[2]
    while !done(iter, st)
        I, st = next(iter, st)
        v[I], ∇v[I] = f(x1[I.I[idx]], x2[I])
    end
    v, u -> ∇multicast_vector_matrix(∇v, u, idx)
end

function ∇multicast_vector_matrix(∇v, x, idx)
    iter = Base.CartesianRange(size(x))
    st = start(iter)
    I, st = next(iter, st)
    val = ∇v[I](x[I])
    v1 = zeros(typeof(val[1]), size(x, idx))
    v2 = similar(Array{typeof(val[2])}, size(x))
    @simd for I in iter
        vv1, v2[I] = ∇v[I](x[I])
        v1[I.I[idx]] += vv1
    end
    v1, v2
end

function δmulticast(f, x1::AbstractMatrix, x2::AbstractVector)
    idx = length(x2) == size(x1, 1) ? 1 : 2
    iter = Base.CartesianRange(size(x1))
    temp = next(iter, start(iter))
    I = temp[1]
    st = temp[2]
    val = f(x1[I], x2[I.I[idx]])
    v = similar(Array{typeof(val[1])}, size(x1))
    ∇v = similar(Array{typeof(val[2])}, size(x1))
    v[I] = val[1]
    ∇v[I] = val[2]
    while !done(iter, st)
        I, st = next(iter, st)
        v[I], ∇v[I] = f(x1[I], x2[I.I[idx]])
    end
    v, u -> ∇multicast_matrix_vector(∇v, u, idx)
end

function ∇multicast_matrix_vector(∇v, x, idx)
    iter = Base.CartesianRange(size(x))
    st = start(iter)
    I, st = next(iter, st)
    val = ∇v[I](x[I])
    v1 = similar(Array{typeof(val[1])}, size(x))
    v2 = zeros(typeof(val[2]), size(x, idx))
    @simd for I in iter
        v1[I], vv2 = ∇v[I](x[I])
        v2[I.I[idx]] += vv2
    end
    v1, v2
end

function δmulticast(f, x1::AbstractFloat, x2::AbstractArray)
    iter = Base.CartesianRange(size(x2))
    temp = next(iter, start(iter))
    I = temp[1]
    st = temp[2]
    val = f(x1, x2[I])
    v = similar(Array{typeof(val[1])}, size(x2))
    ∇v = similar(Array{typeof(val[2])}, size(x2))
    v[I] = val[1]
    ∇v[I] = val[2]
    while !done(iter, st)
        I, st = next(iter, st)
        v[I], ∇v[I] = f(x1, x2[I])
    end
    v, u -> ∇multicast_float_array(∇v, u)
end

function ∇multicast_float_array(∇v, x)
    iter = Base.CartesianRange(size(x))
    st = start(iter)
    I, st = next(iter, st)
    val = ∇v[I](x[I])
    v2 = similar(Array{typeof(val[2])}, size(x))
    v1 = val[1]
    v2[I] = val[2]
    @simd for I in 2:length(x)
        vv1, v2[I] = ∇v[I](x[I])
        v1 += vv1
    end
    v1, v2
end

function δmulticast(f, x1::AbstractArray, x2::AbstractFloat)
    iter = Base.CartesianRange(size(x1))
    temp = next(iter, start(iter))
    I = temp[1]
    st = temp[2]
    val = f(x1[I], x2)
    v = similar(Array{typeof(val[1])}, size(x1))
    ∇v = similar(Array{typeof(val[2])}, size(x1))
    v[I] = val[1]
    ∇v[I] = val[2]
    while !done(iter, st)
        I, st = next(iter, st)
        v[I], ∇v[I] = f(x1[I], x2)
    end
    v, u -> ∇multicast_array_float(∇v, u)
end

function ∇multicast_array_float(∇v, x)
    iter = Base.CartesianRange(size(x))
    st = start(iter)
    I, st = next(iter, st)
    val = ∇v[I](x[I])
    v1 = similar(Array{typeof(val[1])}, size(x))
    v1[I] = val[1]
    v2 = val[2]
    @simd for I in 2:length(x)
        v1[I], vv2 = ∇v[I](x[I])
        v2 += vv2
    end
    v1, v2
end


text:
module Graphons

abstract Graphon
export Graphon

include("utilities.jl")
include("sbms.jl")
include("statistics.jl")
include("simple_graphons.jl")

end # module

ext: jl


text:
function send_inducing_points_to_global(xi_low::Matrix{Float64}, yi_low::Vector{Float64}, root, comm)

    smsg = MPI.Send(xi_low, root, 0, comm)
    # println("xi_low sent to root")

    smsg = MPI.Send(yi_low, root, 1, comm)
    # println("mean_low sent to root")

end #send_inducing_points_to_global

function receive_inducing_points_in_global(ni_low::Int64, p_low::Int64, recv_id, comm)

    xi_low = zeros(ni_low, p_low)
    rmsg = MPI.Recv!(xi_low, recv_id, 0, comm)
    # println("xi_low received in root\n")

    yi_low = zeros(ni_low)
    rmsg = MPI.Recv!(yi_low, recv_id, 1, comm)
    # println("mean_low received in root\n")

    return xi_low, yi_low

end #receive_inducing_points_in_global

function send_prior_to_local(pred_data::Tuple{Vector{Float64},Vector{Float64}}, recv_id, comm)
    prior_mean, prior_sd = pred_data

    smsg = MPI.Send(prior_mean, recv_id, 2, comm)
    # println("prior_mean sent to local")

    smsg = MPI.Send(prior_sd, recv_id, 3, comm)
    # println("prior_cov sent to local")

end # send_prior_to_local

function receive_prior_in_local(ni_low, root, comm)

    prior_mean = zeros(ni_low)
    prior_sd   = zeros(ni_low)

    rmsg = MPI.Recv!(prior_mean, root, 2, comm)
    # println("prior_mean received in local")
    rmsg = MPI.Recv!(prior_sd,  root, 3, comm)
    # println("prior_cov received in local")

    return prior_mean, prior_sd

end # receive_prior_in_local

function send_trace(trace_l, root ,comm)

    smsg = MPI.Send(trace_l, root, 4, comm)
    # println("trace_l sent to root \n")

end # send_trace

function receive_trace(n_local_iters, recv_id, comm)

    trace_l = zeros(n_local_iters)

    rmsg = MPI.Recv!(trace_l, recv_id, 4, comm)
    # println("trace_l received in root \n")

    return trace_l

end # receive_trace

function send_svgp_params(svgp_params, root ,comm)

    smsg = MPI.Send(svgp_params, root, 5, comm)
    # println("trace_l sent to root \n")

end # send_svgp_params

function receive_svgp_params(recv_id, comm)

    recv_params = SVGP_params(SVGP_data("dummy"), 2)

    rmsg = MPI.Recv!(recv_params, recv_id, 5, comm)
    # println("trace_l received in root \n")

    return recv_params

end # receive_svgp_params

function check_for_query(comm)
    is_message, status  = MPI.Iprobe(MPI.MPI_ANY_SOURCE, MPI.MPI_ANY_TAG, comm)
    if is_message
        recv_id = status.source
        tag_ind = status.tag
        return true, recv_id, tag_ind
    else
        return false, -1, -1
    end
end

function send_data_to_neighbor(send_x, send_y, send_ind, comm)
    smsg = MPI.Send(send_x, send_ind, 0, comm)
    smsg = MPI.Send(send_y, send_ind, 1, comm)
end

function receive_neighbor_data(n::Int64, p::Int64, recv_id, comm)
    x = zeros(n, p)
    rmsg = MPI.Recv!(x, recv_id, 0, comm)

    y = zeros(n)
    rmsg = MPI.Recv!(y, recv_id, 1, comm)

    return x, y
end

function query_for_neighbor_data(nbor_ind, part_ind, comm)
    smsg = MPI.Send([0], nbor_ind, part_ind, comm)
end

max_stars_repo_name: ayanbiswas/PRISM

max_issues_repo_name: ayanbiswas/PRISM

max_forks_repo_name: ayanbiswas/PRISM


text:
println("Starting...")
@info "Starting..."
time_start = time()
@info "...at time " time_start
# using Revise
# using ExSup
using ArgParse

function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table s begin
        "--block_size"
            help = "Block size given by scalar Int"
			arg_type = Int
			default = 64
        "--input_path", "-i"
            help = "Input path, The default shape is generated if 'truncated_sphere' string is given."
			default = nothing
        "--output_path", "-o"
            help = "output path"
			default = "."
        "--input_path_in_datasets", "-d"
            help = "Input path relative to Pio3d.jl dataset path"
			default = nothing
        "--threshold"
            help = "another option with an argument"
            arg_type = Int
            default = 7000
        "--stepz"
            help = "Every stepz-th voxel in z-axis is readed"
            arg_type = Int
			default = 1
        "--stepxy"
            help = "Every stepxy-th voxel in xy-axis is readed"
            arg_type = Int
			default = 1
        "--crop"
            help = "Integer number describing the size of crop of input volume for all axes."
            arg_type = Int
			default = nothing
        "--output_csv_file"
            help = "Path to outpu CSV file"
            default = "exp_surface_extraction_times.csv"
        "--label"
            help = "label used in output filename"
            default = "data"
        "--taubin_lambda"
            help = "Taubin smoothing parameter. lambda=0.33, mu=-0.34"
            arg_type = Float64
			default = 0.33
        "--taubin_mu"
            help = "Taubin smoothing parameter. lambda=0.33, mu=-0.34"
            arg_type = Float64
			default = -0.34
        "--taubin_n"
            help = "Taubin smoothing parameter. Number of iterations "
            arg_type = Int
			default = 10
        "--n_procs"
            help = "Number of required CPU-cores"
            arg_type = Int
			default = 4
        "--show"
            help = "Show 3D visualization"
            action = :store_true
        "--skip_smoothing"
            help = "Skip smoothing procedure"
            action = :store_true
        "--color"
			nargs = 4
            help = "Visualization color, RGBA"
			arg_type = Float64
			default = [1, 0, 1, 0.1]
            # action = :store_true
        # "arg1"
        #     help = "a positional argument"
        #     required = true
    end

    return parse_args(s)
end
args = parse_commandline()

# println("args: $(args)")
# println("color")
# println(args["color"])
# exit()

using Test
using Logging
using SparseArrays
using ExSup
using JLD2
using ViewerGL
@info "Distributed init..."
using Distributed
if nprocs() == 1
    addprocs(args["n_procs"]-1)
end
block_size_scalar = args["block_size"]

using LarSurf
# @everywhere using LarSurf

# fn = "exp_surface_extraction_ircad_times.csv"
output_csv_file = args["output_csv_file"]

@info "before everywhere using"
@info "time from start: $(time()-time_start) [s]"
# using Plasm
# @everywhere using LarSurf
# @everywhere using Distributed

@info "after everywhere using, time from start: $(time()-time_start) [s]"

# data_id = 1
show = false
taubin = true
# taubin_n = 5
# taubin_lambda = 0.4
# taubin_mu = -0.2

crop_px = args["crop"]
if crop_px == nothing
	do_crop = false
else
	do_crop = true
end
cropx = crop_px
cropy = crop_px
cropz = crop_px
# stepxy = 4
# block_size = [128, 128, 128]
# block_size = [128, 128, 128]
# block_size = [32, 32, 32]
# data_size1 = 128
# data_size1 = 256
# data_size1 = 512

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

data = LarSurf.Experiments.report_init_row(@__FILE__)
LarSurf.set_time_data(data)

# data["nprocs"] = nprocs()
# data["fcn"] = String(Symbol(fcni))

data["using done"] = time()-time_start
# segmentation = LarSurf.data234()
@info "Generate data..."
@info "time from start: $(time() - time_start) [s]"
# pth = Pio3d.datasets_join_path("medical/orig/3Dircadb1.$data_id/MASKS_DICOM/liver")

# for mask_label in mask_labels

	if args["input_path"] == nothing
		using Pio3d
		if args["input_path_in_datasets"] == nothing
			# pth = Pio3d.datasets_join_path("medical/orig/jatra_mikro_data/Nejlepsi_rozliseni_nevycistene")
			pth = Pio3d.datasets_join_path("medical/processed/corrosion_cast/nrn10.pklz")
		else
			pth = Pio3d.datasets_join_path(args["input_path_in_datasets"])
		end

	else
		pth = args["input_path"]
	end

	data["input_path"] = pth
	if pth[end-4:end] == ".jld2"
		# @info "Surface model given in .jld file. Skipping surface extraction"
		# uu = @JLD2.load pth
		# if :V in uu
		# 	datap_readed = false
		# 	@JLD2.load pth V FV
		# 	FVtri = triangulate_quads(FV)
		# elseif :datap in uu
		@JLD2.load pth datap
		# else
		# 	@error "Expected F and FV or datap in jld2 file"
		# end
		# uu = nothing
		data3d_full = datap["data3d"]
		voxelsize_mm = datap["voxelsize_mm"]
	elseif args["input_path"] == "truncated_sphere"
		data3d_full = LarSurf.generate_truncated_sphere(30)
		voxelsize_mm = [1., 1., 1.]
		args["threshold"] = 0
	else
		using Pio3d
		datap = Pio3d.read3d(pth)
		data3d_full = datap["data3d"]
		voxelsize_mm = datap["voxelsize_mm"]
	end

# V1 is V or Vs accoring to smoothing parameter
V1, FVtri = LarSurf.Experiments.experiment_make_surf_extraction_and_smoothing(
	data3d_full, voxelsize_mm;
	output_path=args["output_path"],
	threshold=args["threshold"],
	mask_label = args["label"],
	stepz = args["stepz"],
	stepxy = args["stepxy"],
	do_crop=do_crop, cropx=cropx, cropy=cropy, cropz=cropz,
	block_size_scalar=block_size_scalar, data=data, time_start=time_start,
	output_csv_file = output_csv_file,
	taubin=taubin,
	taubin_lambda = args["taubin_lambda"],
	taubin_mu = args["taubin_mu"],
	taubin_n = args["taubin_n"],
	smoothing= !args["skip_smoothing"]
)
show=args["show"]
if show
	@info "ViewerGL init ..."
	using ViewerGL
	c = args["color"]
	# print("color: $c")
	EVtri = LarSurf.Smoothing.get_EV_quads(FVtri)
	ViewerGL.VIEW([
	    ViewerGL.GLGrid(V1, FVtri, ViewerGL.Point4d(c[1], c[2], c[3], c[4]))
	    # ViewerGL.GLGrid(V1, FVtri, ViewerGL.Point4d(1,0,1,0))
		ViewerGL.GLGrid(V1, EVtri,ViewerGL.Point4d(0.,0.,0.,0.9))
		# ViewerGL.GLGrid(V1, EVtri, ViewerGL.Point4d(0,0,0,1))
		# ViewerGL.GLGrid(V1, FVtri)
		# ViewerGL.GLAxis(ViewerGL.Point3d(-1, -1, -1),ViewerGL.Point3d(1, 1, 1))
	])
end


text:
function entropy(partitionSet, N)
    return mapreduce(x -> -(length(x) / N) * log(length(x) / N), +, partitionSet)
end

function getPartitionSet(vector, nrOfClusters)
    d = Dict(zip(1:nrOfClusters,[Int[] for _ in 1:nrOfClusters]));
    for i in 1:length(vector)
        push!(d[vector[i]], i)
    end
    return collect(values(d))
end

function mutualInformation(partitionSetA, partitionSetB, N)

    anonFunc = (x) -> (
                intersection = size(intersect(partitionSetA[x[1]], partitionSetB[x[2]]),1);
                return intersection > 0 ?
                        (intersection/N) * log((intersection * N) / (size(partitionSetA[x[1]],1) * size(partitionSetB[x[2]],1))) : 0;
    )

    mapreduce(anonFunc, + , product(1:length(partitionSetA),1:length(partitionSetB)))
end

function normalizedMutualInformation(partitionSetA, partitionSetB, N)
    mutualInformation(partitionSetA, partitionSetB, N) / ((entropy(partitionSetA, N)  + entropy(partitionSetB, N)) / 2)
end


text:
__precompile__(true)

module preGS

using DataFrames
using UnicodePlots
using ProgressMeter
using CSV
using Printf
using Distributed
using Statistics
using DelimitedFiles
using LinearAlgebra

include("preGS_current.jl")

end


text:
import CausalityToolsBase: RectangularBinning
using StaticArrays

t1 = TransferOperatorGridTest(binning = RectangularBinning(5), ηs = [-4, -3, 3, 5])
t2 = TransferOperatorGridTest(binning = RectangularBinning(5), ηs = 1:5)
@test_throws ArgumentError PredictiveAsymmetryTest(t1)
@test_throws ArgumentError PredictiveAsymmetryTest(t2)

x, y = rand(1000), rand(1000)
single_binning = RectangularBinning(5)
multiple_binnings = [RectangularBinning(n) for n = 3:5]

test1 = TransferOperatorGridTest(binning = single_binning, ηs = -10:10)
test2 = TransferOperatorGridTest(binning = single_binning, ηs = 5)
test3 = TransferOperatorGridTest(binning = multiple_binnings, ηs = -10:10)
test4 = TransferOperatorGridTest(binning = multiple_binnings, ηs = 5)

test5 = VisitationFrequencyTest(binning = single_binning, ηs = -10:10)
test6 = VisitationFrequencyTest(binning = single_binning, ηs = 5)
test7 = VisitationFrequencyTest(binning = multiple_binnings, ηs = -10:10)
test8 = VisitationFrequencyTest(binning = multiple_binnings, ηs = 5)

f = 1.0
@test causality(x, y, NormalisedPredictiveAsymmetryTest(test1, f = f)) isa Vector{T} where T
@test causality(x, y, NormalisedPredictiveAsymmetryTest(test2, f = f)) |> typeof <: Real
@test causality(x, y, NormalisedPredictiveAsymmetryTest(test3, f = f)) isa Vector{T} where T
@test causality(x, y, NormalisedPredictiveAsymmetryTest(test4, f = f)) |> typeof<: Real

@test causality(x, y, NormalisedPredictiveAsymmetryTest(test5, f = f)) isa Vector{T} where T
@test causality(x, y, NormalisedPredictiveAsymmetryTest(test6, f = f)) |> typeof<: Real
@test causality(x, y, NormalisedPredictiveAsymmetryTest(test7, f = f)) isa Vector{T} where T
@test causality(x, y, NormalisedPredictiveAsymmetryTest(test8, f = f)) |> typeof<: Real


text:
module OteraEngine

using TOML

export Template

function parse_config(filename::String)
    if filename[end-3:end] != "toml"
        throw(ArgumentError("Suffix of config file must be `toml`! Now, it is `$(filename[end-3:end])`."))
    end
    config = ""
    open(filename, "r") do f
        config = read(f, String)
    end
    return TOML.parse(config)
end

"""
    Template(html::String; path::Bool=true, config_path::String="",
        config::Dict{String, T} = Dict(
            "code_block_start"=>"```",
            "code_block_stop"=>"```"
        )
    )
This is the only structure and function of this package.
This structure has 2 parameter,
- `html` is the path to the HTML file or HTML of String type.
- `path` determines whether the parameter `html` represents the file path. The default value is `true`.
- `config_path` is path to config file. The suffix of config file must be `toml`.
- `config` is configuration of Template. It is type of `Dict`, and now there are two settings bellow.
    - `code_block_start` : The string at the start of code blocks.
    - `code_block_stop` : The string at the end of code blocks.

# Config File
A config file must be written in TOML format. like this:
```
code_block_start = "{{"
code_block_stop = "}}"
```
The item is the same as the argiment `config`.

# HTML
You can write the code of Template in JuliaLang, and just write the variables you want to output to a HTML at the end of the code.
The code needs to be enclosed by ```(This can be changed by `config` variable).

For exmaple, this HTML work:
```
<html>
    <head><title>OteraEngine Test</title></head>
    <body>
        Hello, ```usr```!
    </body>
</html>
```

# Rendering
After you create a Template, you just have to execute the codes! For this, you use the Function-like Object of Template structure.`tmp(; init::Dict{String, T}) where T <: Any` variables are initialized by `init`(`init` is the parameter for Function-like Object). `init` must be `Dict`type. If you don't pass the `init`, the initialization won't be done.
Please see the example below.

# Example
```@repl
tmp = Template("./test1.html") #The last HTML code
init = Dict("usr"=>"OteraEngine")
result = tmp(init)
println(result)
```
"""
struct Template
    splitted_html::Array{String, 1}
    codes::Array{String, 1}
end

function Template(html::String; path::Bool=true, config_path::String="",
        config::Dict{String, T} = Dict(
            "code_block_start"=>"```",
            "code_block_stop"=>"```"
        )
    ) where T<:Any
    if path
        open(html, "r") do f
            html = read(f, String)
        end
    end
    if config_path!=""
        conf_file = parse_config(config_path)
        for v in keys(conf_file)
            config[v] = conf_file[v]
        end
    end
    code = ""
    codes = Array{String}(undef, 0)
    splitted_html = Array{String}(undef, 0)
    i = 1
    code_block_start, code_block_stop = config["code_block_start"], config["code_block_stop"]
    start_len, stop_len = length(code_block_start), length(code_block_stop)
    regex = Regex(code_block_start*"[\\s\\S]*?"*code_block_stop)
    codes_indices = findall(regex, html)
    for index in codes_indices
        s, e = index.start, index.stop
        push!(codes, replace(html[s+start_len:e-stop_len], "\n"=>"; ")) #Does it work without replace?
        push!(splitted_html, html[i:s-1])
        i = e+1
    end
    push!(splitted_html, html[i:end])
    return Template(splitted_html, codes)
end

function (tmp::Template)(init::Dict{String, T}) where T <: Any
    html = tmp.splitted_html[1]
    arg_string = ""
    for v in keys(init)
        arg_string*=(v*",")
    end
    for (part_of_html, code) in zip(tmp.splitted_html[2:end], tmp.codes)
        eval(Meta.parse("function f("*arg_string*"); "*code*";end"))
        html*=string(Base.invokelatest(f, values(init)...))*part_of_html
    end
    return html
end

function (tmp::Template)()
    html = tmp.splitted_html[1]
    for (part_of_html, code) in zip(tmp.splitted_html[2:end], tmp.codes)
        eval(Meta.parse("function f();"*code*";end"))
        html*=string(Base.invokelatest(f))*part_of_html
    end
    return html
end

end


text:
module Molly

using StaticArrays
using ProgressMeter
using BioStructures
using Distributions

using LinearAlgebra: norm, normalize, dot, ×

include("setup.jl")
include("md.jl")
include("analysis.jl")

end


text:
include("prelude.jl")

@testset "println" begin
    interact(`cat`, 1) do proc
        println(proc, "hello world")
        @require expect!(proc, "\n") == "hello world"
        @require process_running(proc)
    end
end

@testset "readstring" begin
    interact(`printf 'hello\nworld\n'`, 1) do proc
        @require read(proc, String) == "hello\nworld\n"
        @require eof(proc)
    end
end

@testset "readuntil" begin
    interact(`printf 'hello\nworld\n'`, 1) do proc
        @require readuntil(proc, "\n", keep=true) == "hello\n"
        @require readuntil(proc, '\n', keep=true) == "world\n"
        @require eof(proc)
    end
end

@testset "readline" begin
    interact(`printf 'hello\nworld\n'`, 1) do proc
        @require readline(proc) == "hello"
        @require !eof(proc)
    end
end


text:
"""
    ArbFloat(x)
    ArbFloat(x; [bits])
    ArbFloat(x; [digits [, base=10]])

- ArbFloat(x) == ArbFloat(x, bits=workingprecision(ArbFloat))
- ArbFloat(x, digits=digits, base=2) == ArbFloat(x, bits=digits)
- ArbFloat(x, digits=digits, base=10) == ArbFloat(x, digits=digits)

ArbFloats are arbitrary precision binary floating-point numbers with some specific precision established at construction.

The precision of an ArbFloat must be >= 24 bits or >= 8 digits.

Negative zero, unsigned Infinity, NaNs with distinct payloads are not supported.

See also: [`ArbReal`](@ref), [`ArbComplex`](@ref), [`workingprecision`](@ref), [`displayprecision`](@ref)
""" ArbFloat

mutable struct ArbFloat{P} <: AbstractFloat    # P is the precision in bits
    exp::Int        # fmpz         exponent of 2 (2^exp)
    size::UInt      # mp_size_t    nwords and sign (lsb holds sign of significand)
    d1::UInt        # significand  unsigned, immediate value or the initial span
    d2::UInt        #   (d1, d2)   the final part indicating the significand, or 0

    function ArbFloat{P}() where {P}
        z = new{P}(0,0,0,0)
        ccall(@libarb(arf_init), Cvoid, (Ref{ArbFloat},), z)
        finalizer(arf_clear, z)
        return z
    end
end

# for use within structs, e.g ArbFloatMatrix
const PtrToArbFloat = Ptr{ArbFloat} # arf_ptr
const PtrToPtrToArbFloat = Ptr{Ptr{ArbFloat}} # arf_ptr*

arf_clear(x::ArbFloat{P}) where {P} = ccall(@libarb(arf_clear), Cvoid, (Ref{ArbFloat},), x)

ArbFloat{P}(x::ArbFloat{P}) where {P} = x
ArbFloat(x::ArbFloat{P}) where {P} = x

ArbFloat{P}(x::Missing) where {P} = missing
ArbFloat(x::Missing) = missing

@inline sign_bit(x::ArbFloat{P}) where {P} = isodd(x.size)

function ArbFloat(x::T; bits::Int=0, digits::Int=0, base::Int=iszero(bits) ? 10 : 2) where {T<:Number}
    bits > 0 && bits < MINIMUM_PRECISION_BASE2 && throw(DomainError("bit precision $bits < $MINIMUM_PRECISION_BASE2"))
    digits > 0 && digits < MINIMUM_PRECISION_BASE10 && throw(DomainError("digit precision $digits < $MINIMUM_PRECISION_BASE10"))
    if base === 10
        bits = digits > 0 ? bits4digits(digits)+extrabits() : (bits > 0 ? bits+extrabits() : DEFAULT_PRECISION.x)
    elseif base === 2
        bits = bits > 0 ? bits+extrabits() : (digits > 0 ? digits+extrabits() : DEFAULT_PRECISION.x)
    else
        throw(ErrorException("base expects 2 or 10"))
    end
    ArbFloat{bits}(x)
end

swap(x::ArbFloat{P}, y::ArbFloat{P}) where {P} = ccall(@libarb(arf_swap), Cvoid, (Ref{ArbFloat}, Ref{ArbFloat}), x, y)

function copy(x::ArbFloat{P}) where {P}
    z = ArbFloat{P}()
    ccall(@libarb(arf_set), Cvoid, (Ref{ArbFloat}, Ref{ArbFloat}), z, x)
    return z
end

function copy(x::ArbFloat{P}, bitprecision::Int, roundingmode::RoundingMode) where {P}
    z = ArbFloat{P}()
    rounding = match_rounding_mode(roundingmode)
    rounddir = ccall(@libarb(arf_set_round), Cint, (Ref{ArbFloat}, Ref{ArbFloat}, Clong, Cint), z, x, bitprecision, rounding)
    return z
end

copy(x::ArbFloat{P}, roundingmode::RoundingMode) where {P} = copy(x, P, roundingmode)
copy(x::ArbFloat{P}, bitprecision::Int) where {P} = copy(x, bitprecision, RoundNearest)

function ArbFloat{P}(x::Int32) where {P}
    z = ArbFloat{P}()
    ccall(@libarb(arf_set_si), Cvoid, (Ref{ArbFloat}, Clong), z, x)
    return z
end
ArbFloat{P}(x::T) where {P, T<:Union{Int8, Int16}} = ArbFloat{P}(Int32(x))
ArbFloat{P}(x::T) where {P, T<:Union{Int64, Int128}} = ArbFloat{P}(BigInt(x))

function ArbFloat{P}(x::UInt32) where {P}
    z = ArbFloat{P}()
    ccall(@libarb(arf_set_ui), Cvoid, (Ref{ArbFloat}, Culong), z, x)
    return z
end
ArbFloat{P}(x::T) where {P, T<:Union{UInt8, UInt16}} = ArbFloat{P}(UInt32(x))
ArbFloat{P}(x::T) where {P, T<:Union{UInt64, UInt128}} = ArbFloat{P}(BigInt(x))

function ArbFloat{P}(x::Float64) where {P}
    z = ArbFloat{P}()
    ccall(@libarb(arf_set_d), Cvoid, (Ref{ArbFloat}, Cdouble), z, x)
    return z
end
ArbFloat{P}(x::T) where {P, T<:Union{Float16, Float32}} = ArbFloat{P}(Float64(x))

function ArbFloat{P}(x::BigFloat) where {P}
    z = ArbFloat{P}()
    ccall(@libarb(arf_set_mpfr), Cvoid, (Ref{ArbFloat}, Ref{BigFloat}), z, x)
    return z
end
ArbFloat{P}(x::BigInt) where {P} = ArbFloat{P}(BigFloat(x))
ArbFloat{P}(x::Rational{T}) where {P, T<:Signed} = ArbFloat{P}(BigFloat(x))

function ArbFloat{P}(x::Irrational{S}) where {P,S}
    prec = precision(BigFloat)
    newprec = max(prec, P + 32)
    setprecision(BigFloat, newprec)
    y = BigFloat(x)
    z = ArbFloat{P}(y)
    setprecision(BigFloat, prec)
    return z
end

function Int64(x::ArbFloat{P}, roundingmode::RoundingMode) where {P}
    rounding = match_rounding_mode(roundingmode)
    z = ccall(@libarb(arf_get_si), Clong, (Ref{ArbFloat}, Cint), x, rounding)
    return z
end
Int32(x::ArbFloat{P}, roundingmode::RoundingMode) where {P} = Int32(Int64(x), roundingmode)
Int16(x::ArbFloat{P}, roundingmode::RoundingMode) where {P} = Int16(Int64(x), roundingmode)

BigFloat(x::ArbFloat{P}) where {P} = BigFloat(x, RoundNearest)
function BigFloat(x::ArbFloat{P}, roundingmode::RoundingMode) where {P}
    rounding = match_rounding_mode(roundingmode)
    z = BigFloat(0, workingprecision(x))
    roundingdir = ccall(@libarb(arf_get_mpfr), Cint, (Ref{BigFloat}, Ref{ArbFloat}, Cint), z, x, rounding)
    return z
end
BigFloat(x::ArbFloat{P}, bitprecision::Int) where {P} = BigFloat(x, bitprecision, RoundNearest)
function BigFloat(x::ArbFloat{P}, bitprecision::Int, roundingmode::RoundingMode) where {P}
    rounding = match_rounding_mode(roundingmode)
    z = BigFloat(0, bitprecision)
    roundingdir = ccall(@libarb(arf_get_mpfr), Cint, (Ref{BigFloat}, Ref{ArbFloat}, Cint), z, x, rounding)
    return z
end

BigInt(x::ArbFloat{P}) where {P} = BigInt(trunc(BigFloat(x)))

function Base.Integer(x::ArbFloat{P}) where {P}
    if isinteger(x)
       abs(x) <= typemax(Int64) ? Int64(x) : BigInt(x)
    else
       throw(InexactError("$x"))
    end
end

for (F,A) in ((:floor, :arf_floor), (:ceil, :arf_ceil))
    @eval begin
        function $F(x::ArbFloat{P}) where {P}
            z = ArbFloat{P}()
            ccall(@libarb($A), Cvoid, (Ref{ArbFloat}, Ref{ArbFloat}), z, x)
            return z
        end
        $F(::Type{T}, x::ArbFloat{P}) where {P, T<:Integer} = T($F(x))
    end
end

trunc(x::ArbFloat{P}) where {P} = signbit(x) ? ceil(x) : floor(x)
trunc(::Type{T}, x::ArbFloat{P}) where {P, T<:Integer} = T(trunc(x))

midpoint(x::ArbFloat{P}) where {P} = x
radius(x::ArbFloat{P}) where {P} = zero(ArbFloat{P})

function modf(x::ArbFloat{P}) where {P}
    ipart = trunc(x)
    fpart = x - ipart
    return (fpart, ipart)
end

fmod(fpartipart::Tuple{ArbFloat{P}, ArbFloat{P}}) where {P} =
    fpartipart[1] + fpartipart[2]
fmod(fpart::ArbFloat{P}, ipart::ArbFloat{P}) where {P} =
    fpart + ipart

div(x::ArbFloat{P}, y::ArbFloat{P}) where {P} =
    trunc(x / y)

rem(x::ArbFloat{P}, y::ArbFloat{P}) where {P} =
    x - (div(x,y) * y)

function divrem(x::ArbFloat{P}, y::ArbFloat{P}) where {P}
    dv = div(x,y)
    rm = x - (dv * y)
    return (dv, rm)
end

fld(x::ArbFloat{P}, y::ArbFloat{P}) where {P} =
    floor(x / y)

mod(x::ArbFloat{P}, y::ArbFloat{P}) where {P} =
    x - (fld(x,y) * y)

function fldmod(x::ArbFloat{P}, y::ArbFloat{P}) where {P}
    fd = fld(x,y)
    md = x - (fd * y)
    return (fd, md)
end

cld(x::ArbFloat{P}, y::ArbFloat{P}) where {P} =
    ceil(x / y)

# a type specific hash function helps the type to 'just work'
const hash_arbfloat_lo = (UInt === UInt64) ? 0x37e642589da3416a : 0x5d46a6b4
const hash_0_arbfloat_lo = hash(zero(UInt), hash_arbfloat_lo)
hash(z::ArbFloat{P}, h::UInt) where {P} =
    hash(reinterpret(UInt,z.d1) ⊻ z.exp,
         (h ⊻ hash(z.d2 ⊻ (~reinterpret(UInt,P)), hash_arbfloat_lo) ⊻ hash_0_arbfloat_lo))


text:
function csr_to_sparse_test()
    i = [1;2;3]
    j = [3;4;4]
    v = [8;9;10]
    (rp,ci,ai,m) = sparse_to_csr(i,j,v)
    (nzi,nzj,nzv) = csr_to_sparse(rp,ci,ai)
    A = sparse(nzi,nzj,nzv,length(rp)-1,maximum(ci))

    # more tests added here
    # clique to sparse test
    rp = collect(1:5:26)
    ci = vec(reshape(repmat(1:5,5,1)',25,1))
    ai = ones(Int64,25)
    A = csr_to_sparse_matrix(rp,ci,ai,5,5)
    if !isequal(full(A),ones(5,5))
        error("csr_to_sparse_test failed")
    end

    # 100 random trials
    for t = 1:100
        A = sprand(100,80,0.01)
        (rp,ci,ai) = sparse_to_csr(A)
        A2 = csr_to_sparse_matrix(rp,ci,ai,100,80)
        if ~isequal(A,A2)
            error("csr_to_sparse_test failed")
        end
    end
    return true
end


text:
module ssa_varimax

export varimax_rotate!

using LinearAlgebra

"""
Varimax rotation as per
   Groth, A. and M. Ghil, 2011: Multivariate singular spectrum analysis and the
      road to phase synchronization, Physical Review E, 84, 036206.
   Groth, A. and M. Ghil, 2015: Monte Carlo Singular Spectrum Analysis (SSA)
      revisited: Detecting oscillator clusters in multivariate datasets, Journal
      of Climate, 28, 7873-7893.

Code based on Matlab implementation by Andreas Groth at
https://www.mathworks.com/matlabcentral/fileexchange/65939-multichannel-singular-spectrum-analysis-varimax-tutorial

Copyright (c) 2018, Andreas Groth
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

* Redistributions of source code must retain the above copyright notice, this
  list of conditions and the following disclaimer.

* Redistributions in binary form must reproduce the above copyright notice,
  this list of conditions and the following disclaimer in the documentation
  and/or other materials provided with the distribution

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
"""
function varimax(A::Array{float_type, 3}, reltol=sqrt(eps(float_type)),
                 maxit=1000, normalize=true, G=[]) where {float_type<:AbstractFloat}
   M, D, S = size(A)
   B = A

   B = reshape(B, D*M, S, 1)[:, :, 1]
   D = D*M

   T = diagm(0=>ones(S))

   for iter=1:maxit
      maxTheta = 0
      for i = 1:(S-1)
         for j = (i+1):S
            u = B[:, i].^2 - B[:,j].^2
            v = 2*B[:,i].*B[:,j]
            if M>1
               u = sum(reshape(u, M, :), dims=1)'
               v = sum(reshape(v, M, :), dims=1)'
            end

            usum = sum(u, dims=1);
            vsum = sum(v, dims=1);
            numer = 2*u'*v - 2*usum*vsum/D;
            denom = u'*u - v'*v - (usum^2 - vsum^2)/D;
            theta = atan(numer[1, 1], denom[1, 1])/4;
            maxTheta = max(maxTheta, abs(theta));
            Tij = [cos(theta) -sin(theta); sin(theta) cos(theta)];
            B[:, [i, j]] = B[:, [i, j]] * Tij;
            T[:, [i, j]] = T[:, [i, j]] * Tij;
         end
      end
      if (maxTheta < reltol)
         break
      end
   end

   return B, T
end

function varimax_rotate!(EW, EV, M, D, S)
   EVscaled=EV[:,1:S]*diagm(0=>sqrt.(EW[1:S]))
   EVreshape=reshape(EVscaled,M,D,S)

   _, T = varimax(EVreshape)

   EV[:,1:S]=EV[:,1:S]*T
   EW[1:S]=diag(T'*diagm(0=>EW[1:S])*T)

   EV = EV[:, sortperm(EW, rev=true)]
   EW = sort(EW, rev=true)

   return EW, EV
end

end


text:
# This file was generated, do not modify it. # hide
X = source(Xraw)
y = source(yraw; kind=:target)

max_stars_repo_name: ven-k/MLJTutorials

max_issues_repo_name: ven-k/MLJTutorials

max_forks_repo_name: ven-k/MLJTutorials


text:
module Day7
export parse_input, build_dag, part_1, part_2

pattern = r"Step (\w) must be finished before step (\w) can begin"

function parse_input(input::String)::Array{Array{String, 1}, 1}
    return map(
        l -> match(pattern, l).captures,
        split(strip(input), "\n")
    )
end

function build_dag(parsed_input::Array{Array{String, 1}, 1})
    dag = Dict{String, Array{String}}()
    for (dependency, target) in parsed_input
        if !in(target, keys(dag))
            dag[target] = [dependency]
        else
            push!(dag[target], dependency)
        end

        if !in(dependency, keys(dag))
            dag[dependency] = []
        end
    end
    return dag
end

function part_1(dag)
    dag = deepcopy(dag)
    order = ""
    while length(dag) > 0
        possible = collect(filter(k -> length(dag[k]) == 0, keys(dag)))
        sort!(possible)

        next = possible[1]
        pop!(dag, next)
        order *= next

        for dependencies in values(dag)
            filter!(l -> l != next, dependencies)
        end
    end
    return order
end

function part_2(dag; workers=2, per_task=0)
    dag = deepcopy(dag)
    done_at = Dict{String, Int}()
    free_workers = workers

    t = -1
    while (length(dag) > 0) | (length(done_at) > 0)
        t += 1

        for task_done in filter(k -> done_at[k] == t, keys(done_at))
            pop!(done_at, task_done)
            free_workers += 1

            for dependencies in values(dag)
                filter!(l -> l != task_done, dependencies)
            end

        end

        if free_workers == 0
            continue
        end

        possible = collect(filter(k -> length(dag[k]) == 0, keys(dag)))

        if length(possible) == 0
            continue
        end

        sort!(possible)
        n_possible = min(length(possible), free_workers)

        for i in 1:n_possible
            next = possible[i]
            pop!(dag, next)
            free_workers -= 1
            done_at[next] = t + per_task + next[1] - 'A' + 1
        end
    end

    return t
end

end


text:
"""
    _colon(I:J)

`_colon(I, J)` works equivelently to `I:J`, it's used to backward support julia v"1.0".
"""
_colon(I, J) = I:J
if v"1.0" <= VERSION < v"1.1"
    _colon(I::CartesianIndex{N}, J::CartesianIndex{N}) where N =
        CartesianIndices(map((i,j) -> i:j, Tuple(I), Tuple(J)))
end

# patch for ColorVectorSpace 0.9
# for CVS < 0.9, we can just use the fallback solution in Distances
if isdefined(ImageCore.ColorVectorSpace, :⊙)
    # Because how abs2 calculated in color vector space is ambiguious, abs2(::RGB) is un-defined
    # since ColorVectorSpace 0.9
    # https://github.com/JuliaGraphics/ColorVectorSpace.jl/pull/131
    @inline _abs2(c::Colorant) = mapreducec(v->float(v)^2, +, 0, c)
    @inline _abs2(x) = abs2(x)
    @inline _abs(c::Colorant) = mapreducec(v->abs(float(v)), +, 0, c)
    @inline _abs(x) = abs(x)
else
    const _abs2 = abs2
    const _abs = abs
end


text:
# Support Vector Machines Module
module SVM
export WeightsLinearSVM, train_linear
export WeightsSVM, train
export cost, predict, predict_proba

import Random
import Statistics

include("../tools.jl")
import .JuTools

"""
The weight data struct for training linear SVM
"""
mutable struct WeightsLinearSVM
    C::AbstractFloat
    w::Array{T} where T<:AbstractFloat
    b::AbstractFloat
end

"""
Cost function for linear SVM\\
If `X_data` has shape (M,N)\\
`Y_data` has shape (M,)\\
`weights.w` has shape (N,)\\
Returns the computed loss
"""
function cost(X_data::Array{T} where T<:Number, Y_data::Array{T} where T<:Number, weights::WeightsLinearSVM)::AbstractFloat
    @assert ndims(Y_data) == ndims(weights.w) == 1
    @assert size(X_data) == (size(Y_data)[1], size(weights.w)[1])
    loss_w = 0.5 * (weights.w' * weights.w)
    loss_inner = 1.0 .- Y_data .* vec(X_data * weights.w .+ weights.b)
    loss_inner .= map(m->max(0.0,m), loss_inner)
    loss = loss_w + weights.C * sum(loss_inner) / size(X_data)[1]
    return loss
end

"""
Learning function for linear SVM, using Gradient Descent\\
If `X_data` has shape (M,N)\\
`Y_data` has shape (M,)\\
`weights.w` has shape (N,)\\
Update the weights in place
"""
function learn!(X_data::Array{T} where T<:Number, Y_data::Array{T} where T<:Number, weights::WeightsLinearSVM, momentum::WeightsLinearSVM, alpha::AbstractFloat)
    @assert ndims(Y_data) == ndims(weights.w) == 1
    @assert size(X_data) == (size(Y_data)[1], size(weights.w)[1])
    # compute deciding feature
    decide = (Y_data .* (X_data * weights.w .+ weights.b)) .< 1 # (? < 1) will be 1, otherwise 0
    # update w
    gradient_w = weights.w .+ (weights.C / size(X_data)[1]) .* vec(-(Y_data .* decide)' * X_data)
    gradient_w .= gradient_w .* alpha
    momentum.w .= gradient_w .+ (0.9 .* momentum.w)
    weights.w .= weights.w .- momentum.w
    # update b
    gradient_b = (weights.C / size(X_data)[1]) * sum(-(Y_data .* decide))
    gradient_b *= alpha
    momentum.b = gradient_b + (0.9 * momentum.b)
    weights.b = weights.b - momentum.b
    return nothing
end

"""
Prediction function for linear SVM\\
If `X_data` has shape (M,N)\\
`weights.w` has shape (N,)\\
Returns the actual computed values
"""
function predict_proba(X_predict::Array{T} where T<:Number, weights::WeightsLinearSVM)::Array
    @assert 1 <= ndims(X_predict) <= 2
    if ndims(X_predict) == 1
        X_predict = reshape(X_predict, (1, size(X_predict)[1]))
    end
    @assert size(X_predict)[2] == size(weights.w)[1]
    prediction = vec(X_predict * weights.w .+ weights.b)
    return prediction
end

"""
Prediction function for linear SVM\\
If `X_data` has shape (M,N)\\
`weights.w` has shape (N,)\\
Returns the converted predictions, in {-1, 1}
"""
function predict(X_predict::Array{T} where T<:Number, weights::WeightsLinearSVM)::Array
    @assert 1 <= ndims(X_predict) <= 2
    if ndims(X_predict) == 1
        X_predict = reshape(X_predict, (1, size(X_predict)[1]))
    end
    @assert size(X_predict)[2] == size(weights.w)[1]
    prediction = vec(X_predict * weights.w .+ weights.b)
    prediction .= map(m -> m >= 0 ? 1.0 : -1.0, prediction)
    return prediction
end

"""
Training function for linear SVM\\
If `X_data` has shape (M,N)\\
`Y_data` has shape (M,)\\
`C` is the constraint for SVM decision margin\\
`max_iter` should be >= 0
------
Returns trained weights as `WeightsLinearSVM` object

------
Set `early_stop` to `false`, to force run maximum iteractions\\
Set `random_weights` to `false` to initialize weights of 0.0
"""
function train_linear(X_data::Array{T} where T<:Number, Y_data::Array{T} where T<:Number, C::AbstractFloat;
        learning_rate::AbstractFloat=0.1, max_iter::Integer=1000, verbose::Bool=false, random_weights::Bool=true)::WeightsLinearSVM
    @assert ndims(X_data) == ndims(Y_data) + 1 == 2
    @assert size(X_data)[1] == size(Y_data)[1]
    @assert max_iter >= 0
    X_data = Float64.(X_data)
    Y_data = Float64.(Y_data)
    # is it better to use zero weights than normal weights ?
    weights = nothing
    if random_weights
        weights = WeightsLinearSVM(C, Random.randn(size(X_data)[2]), Random.randn())
    else
        weights = WeightsLinearSVM(C, zeros(size(X_data)[2]), 0.0)
    end
    momentum = WeightsLinearSVM(C, zeros(size(X_data)[2]), 0.0)
    best_weights = WeightsLinearSVM(C, copy(weights.w), weights.b)
    best_cost = nothing
    for i in 1:max_iter
        JuTools.shuffle_data!(X_data, Y_data)
        learn!(X_data, Y_data, weights, momentum, learning_rate)
        new_cost = cost(X_data, Y_data, weights)
        if verbose
            acc = JuTools.compute_accuracy(predict(X_data, weights), Y_data)
            println("Iter: $i")
            println("Cost = $new_cost")
            println("Accuracy = $acc")
            println()
        end
        if best_cost === nothing
            best_cost = new_cost
        else
            if new_cost <= best_cost
                best_cost = min(new_cost, best_cost)
                best_weights.w = copy(weights.w)
                best_weights.b = weights.b
            end
        end
    end
    return best_weights
end

"""
The weight data struct for training common SVM (with kernels)
"""
mutable struct WeightsSVM
    C::AbstractFloat                # constraint
    b::AbstractFloat                # threshold
    gamma::AbstractFloat            # parameter used for polynomial, rbf, and sigmoid kernels
    r::AbstractFloat                # parameter used for polynomial, and sigmoid kernels
    d::AbstractFloat                # parameter used for polynomial kernel
    tol_alpha::AbstractFloat        # tolerance for alpha
    tol_error::AbstractFloat        # tolerance for error
    alpha::Array{T} where T<:Number # alpha array
    error::Array{T} where T<:Number # array for error cache
    kernel::String                  # kernel function name
end

"""
Kernel Function (linear)
"""
function kernel_linear(X1::Array{T} where T<:Number, X2::Array{T} where T<:Number)::Array
    @assert ndims(X1) == ndims(X2) == 2
    @assert size(X1)[2] == size(X2)[2]
    result = X1 * X2'
    return result
end
"""
Kernel Function (polynomial)\\
`d` is `degree` parameter in training function\\
`r` is `coef` parameter in training function\\
`gamma` is computed from `gamma` parameter in training function
"""
function kernel_polynomial(X1::Array{T} where T<:Number, X2::Array{T} where T<:Number;
        d::AbstractFloat=1.0, r::AbstractFloat=0.0, gamma::AbstractFloat=1.0)::Array
    @assert ndims(X1) == ndims(X2) == 2
    @assert size(X1)[2] == size(X2)[2]
    result = (gamma .* (X1 * X2') .+ r) .^ d
    return result
end
"""
Kernel Function (rbf)\\
`gamma` is computed from `gamma` parameter in training function
"""
function kernel_rbf(X1::Array{T} where T<:Number, X2::Array{T} where T<:Number; gamma::AbstractFloat=1.0)::Array
    @assert ndims(X1) == ndims(X2) == 2
    @assert size(X1)[2] == size(X2)[2]
    result = (sum(X1 .^ 2, dims=2) * ones(size(X2)[1])') .+ (ones(size(X1)[1]) * sum(X2 .^ 2, dims=2)') .- 2.0 .* (X1 * X2')
    result .= broadcast(m->max(0.0, m), result) # ignore very small negative outputs, due to precision
    result .= sqrt.(result)
    result .= (-gamma) .* result
    result .= exp.(result)
    return result
end
"""
Kernel Function (sigmoid)\\
`r` is `coef` parameter in training function\\
`gamma` is computed from `gamma` parameter in training function
"""
function kernel_sigmoid(X1::Array{T} where T<:Number, X2::Array{T} where T<:Number; gamma::AbstractFloat=1.0, r::AbstractFloat=0.0)::Array
    @assert ndims(X1) == ndims(X2) == 2
    @assert size(X1)[2] == size(X2)[2]
    result = gamma .* (X1 * X2') .+ r
    result .= tanh.(result)
    return result
end

"""
Cost function for SVM\\
If `X_data` has shape (M,N)\\
`weights.alpha` has shape (M,)\\
Returns the computed cost\\
__Note__:\\
Cost is the opposite of objective value
"""
function cost(X_data::Array{T} where T<:Number, weights::WeightsSVM)::AbstractFloat
    @assert ndims(X_data) == ndims(weights.alpha) + 1 == 2
    @assert size(X_data)[1] == size(weights.alpha)[1]
    result = nothing
    if weights.kernel == "linear"
        result = kernel_linear(X_data, X_data)
    elseif weights.kernel == "polynomial"
        result = kernel_polynomial(X_data, X_data, d=weights.d, r=weights.r, gamma=weights.gamma)
    elseif weights.kernel == "rbf"
        result = kernel_rbf(X_data, X_data, gamma=weights.gamma)
    elseif weights.kernel == "sigmoid"
        result = kernel_sigmoid(X_data, X_data, gamma=weights.gamma, r=weights.r)
    else
        throw(ArgumentError("Error: kernel function $weights.kernel is not recognized"))
    end
    result = 0.5 * (weights.alpha' * result * weights.alpha) - sum(weights.alpha)
    return result
end

"""
Step function for learning SVM\\
If `X_data` has shape (M,N)\\
`Y_data` has shape (M,)\\
`weights.alpha` has shape (M,)\\
`weights.error` has shape (M,)\\
`id1` and `id2` are index in `weights.alpha`\\
Tries to update `weights.alpha` in place\\
Returns `1` if step further, else `0`
"""
function learn_step!(X_data::Array{T} where T<:Number, Y_data::Array{T} where T<:Number,
        weights::WeightsSVM, id1::Integer, id2::Integer)::Integer
    @assert ndims(X_data) == ndims(Y_data) + 1 == 2
    @assert size(X_data)[1] == size(Y_data)[1]
    @assert size(X_data)[1] == size(weights.alpha)[1]
    @assert size(weights.error) == size(weights.alpha)
    @assert id1 >= 1
    @assert id2 >= 1
    # if choosing same alpha, skip
    if id1 == id2
        return 0
    end
    # prepare data
    alpha1 = weights.alpha[id1]
    alpha2 = weights.alpha[id2]
    Y1 = Y_data[id1]
    Y2 = Y_data[id2]
    error1 = weights.error[id1]
    error2 = weights.error[id2]
    # compute L & H
    L = nothing
    H = nothing
    if Y1 != Y2
        L = max(0.0, alpha2 - alpha1)
        H = min(weights.C, weights.C + alpha2 - alpha1)
    else
        L = max(0.0, alpha1 + alpha2 - weights.C)
        H = min(weights.C, alpha1 + alpha2)
    end
    if L == H
        return 0
    end
    # compute kernel results and 2nd derivative eta
    k11 = nothing
    k12 = nothing
    k22 = nothing
    n_features = size(X_data)[2]
    X_id1 = reshape(X_data[id1, :], (1, n_features))
    X_id2 = reshape(X_data[id2, :], (1, n_features))
    if weights.kernel == "linear"
        k11 = kernel_linear(X_id1, X_id1)[1]
        k12 = kernel_linear(X_id1, X_id2)[1]
        k22 = kernel_linear(X_id2, X_id2)[1]
    elseif weights.kernel == "polynomial"
        k11 = kernel_polynomial(X_id1, X_id1, d=weights.d, r=weights.r, gamma=weights.gamma)[1]
        k12 = kernel_polynomial(X_id1, X_id2, d=weights.d, r=weights.r, gamma=weights.gamma)[1]
        k22 = kernel_polynomial(X_id2, X_id2, d=weights.d, r=weights.r, gamma=weights.gamma)[1]
    elseif weights.kernel == "rbf"
        k11 = kernel_rbf(X_id1, X_id1, gamma=weights.gamma)[1]
        k12 = kernel_rbf(X_id1, X_id2, gamma=weights.gamma)[1]
        k22 = kernel_rbf(X_id2, X_id2, gamma=weights.gamma)[1]
    elseif weights.kernel == "sigmoid"
        k11 = kernel_sigmoid(X_id1, X_id1, gamma=weights.gamma, r=weights.r)[1]
        k12 = kernel_sigmoid(X_id1, X_id2, gamma=weights.gamma, r=weights.r)[1]
        k22 = kernel_sigmoid(X_id2, X_id2, gamma=weights.gamma, r=weights.r)[1]
    else
        throw(ArgumentError("Error: kernel function $weights.kernel is not recognized"))
    end
    eta = 2.0 * k12 - k11 - k22
    # compute new alpha2 (a2)
    a2 = nothing
    if eta < 0.0
        a2 = alpha2 - Y2 * (error1 - error2) / eta
        a2 = min(a2, H)
        a2 = max(a2, L)
    else
        c1 = eta / 2.0
        c2 = Y2 * (error1 - error2) - eta * alpha2
        Lobj = c1 * L * L + c2 * L
        Hobj = c1 * H * H + c2 * H
        if Lobj > (Hobj + weights.tol_alpha)
            a2 = L
        elseif Lobj < (Hobj - weights.tol_alpha)
            a2 = H
        else
            a2 = alpha2
        end
    end
    # push to 0 or C
    if a2 < 1e-8
        a2 = 0.0
    elseif a2 > (weights.C - 1e-8)
        a2 = weights.C
    end
    # skip if cannot be optimized
    if abs(a2 - alpha2) < weights.tol_alpha * (a2 + alpha2 + weights.tol_alpha)
        return 0
    end
    # compute new alpha1 (a1)
    a1 = alpha1 + (Y1 * Y2) * (alpha2 - a2)
    if a1 < 0.0
        a2 += (Y1 * Y2) * a1
        a1 = 0.0
    elseif a1 > weights.C
        a2 += (Y1 * Y2) * (a1 - weights.C)
        a1 = weights.C
    end
    # update threshold
    b1 = error1 + Y1 * (a1 - alpha1) * k11 + Y2 * (a2 - alpha2) * k12 + weights.b
    b2 = error2 + Y1 * (a1 - alpha1) * k12 + Y2 * (a2 - alpha2) * k22 + weights.b
    b_new = nothing
    if 0 < a1 < weights.C
        b_new = b1
    elseif 0 < a2 < weights.C
        b_new = b2
    else
        b_new = (b1 + b2) * 0.5
    end
    # update error cache
    non_optimized_ids = [i for i in 1:size(X_data)[1] if (0 < weights.alpha[i] < weights.C)]
    kerr1 = nothing
    kerr2 = nothing
    if weights.kernel == "linear"
        kerr1 = vec(kernel_linear(X_id1, X_data[non_optimized_ids, :]))
        kerr2 = vec(kernel_linear(X_id2, X_data[non_optimized_ids, :]))
    elseif weights.kernel == "polynomial"
        kerr1 = vec(kernel_polynomial(X_id1, X_data[non_optimized_ids, :], d=weights.d, r=weights.r, gamma=weights.gamma))
        kerr2 = vec(kernel_polynomial(X_id2, X_data[non_optimized_ids, :], d=weights.d, r=weights.r, gamma=weights.gamma))
    elseif weights.kernel == "rbf"
        kerr1 = vec(kernel_rbf(X_id1, X_data[non_optimized_ids, :], gamma=weights.gamma))
        kerr2 = vec(kernel_rbf(X_id2, X_data[non_optimized_ids, :], gamma=weights.gamma))
    elseif weights.kernel == "sigmoid"
        kerr1 = vec(kernel_sigmoid(X_id1, X_data[non_optimized_ids, :], gamma=weights.gamma, r=weights.r))
        kerr2 = vec(kernel_sigmoid(X_id2, X_data[non_optimized_ids, :], gamma=weights.gamma, r=weights.r))
    end
    weights.error[non_optimized_ids] .= weights.error[non_optimized_ids] .+
        ((Y1*(a1-alpha1)) .* kerr1) .+ ((Y2*(a2-alpha2)) .* kerr2) .- (b_new - weights.b)
    weights.error[id1] = 0.0
    weights.error[id2] = 0.0
    # update alpha and b
    weights.b = b_new
    weights.alpha[id1] = a1
    weights.alpha[id2] = a2
    return 1
end

"""
Learning function for SVM, using SMO\\
If `X_data` has shape (M,N)\\
`Y_data` has shape (M,)\\
`weights.alpha` has shape (M,)\\
Tries to update `weights.alpha` with step function\\
Returns `1` if step further, else `0`
"""
function learn!(X_data::Array{T} where T<:Number, Y_data::Array{T} where T<:Number,
        weights::WeightsSVM, id::Integer; verbose::Bool=false)::Integer
    @assert ndims(X_data) == ndims(Y_data) + 1 == 2
    @assert size(X_data)[1] == size(Y_data)[1]
    @assert size(X_data)[1] == size(weights.alpha)[1]
    @assert size(weights.error) == size(weights.alpha)
    @assert id >= 1
    Y = Y_data[id]
    alpha = weights.alpha[id]
    error = weights.error[id]
    r = error * Y
    if ((r < -weights.tol_error) && (alpha < weights.C)) || ((r > weights.tol_error) && (alpha > 0))
        alpha_target = [i for (i, m) in enumerate(weights.alpha) if (0.0 < m < weights.C)]
        # try argmax E1 - E2
        if verbose
            println("Trying argmax(abs(E1 - E2))")
        end
        new_id = 0
        tmax = 0
        for i in alpha_target
            tmp = abs(error - weights.error[i])
            if(tmp > tmax)
                tmax = tmp
                new_id = i
            end
        end
        if new_id >= 1
            step = learn_step!(X_data, Y_data, weights, new_id, id)
            if step > 0
                return step
            end
        end
        # loop non-bound alphas, randomly
        if verbose
            println("Trying random non-bound alphas")
        end
        for new_id in alpha_target[Random.randperm(length(alpha_target))]
            step = learn_step!(X_data, Y_data, weights, new_id, id)
            if step > 0
                return step
            end
        end
        # else loop all alphas, randomly
        if verbose
            println("Trying random remaining alphas")
        end
        for new_id in Random.randperm(length(weights.alpha))
            if new_id in alpha_target
                continue # skip the alpha ids that already looked at
            end
            step = learn_step!(X_data, Y_data, weights, new_id, id)
            if step > 0
                return step
            end
        end
    end
    return 0
end

"""
Prediction function for SVM\\
If `X_data` has shape (M,N)\\
`Y_data` has shape (M,)\\
`weights.alpha` has shape (M,)\\
`X_predict` has shape (k,N)\\
Returns the actual predictions
"""
function predict_proba(X_predict::Array{T} where T<:Number, X_data::Array{T} where T<:Number,
        Y_data::Array{T} where T<:Number, weights::WeightsSVM)::Array
    if ndims(X_predict) == 1
        X_predict = reshape(X_predict, (1, size(X_predict)[1]))
    end
    @assert ndims(X_predict) == ndims(X_data) == ndims(Y_data) + 1 == 2
    @assert size(X_predict)[2] == size(X_data)[2]
    result = nothing
    if weights.kernel == "linear"
        result = kernel_linear(X_data, X_predict)
    elseif weights.kernel == "polynomial"
        result = kernel_polynomial(X_data, X_predict, d=weights.d, r=weights.r, gamma=weights.gamma)
    elseif weights.kernel == "rbf"
        result = kernel_rbf(X_data, X_predict, gamma=weights.gamma)
    elseif weights.kernel == "sigmoid"
        result = kernel_sigmoid(X_data, X_predict, gamma=weights.gamma, r=weights.r)
    else
        throw(ArgumentError("Error: kernel function $weights.kernel is not recognized"))
    end
    prediction = vec((weights.alpha .* Y_data)' * result)
    return prediction
end

"""
Prediction function for SVM\\
If `X_data` has shape (M,N)\\
`Y_data` has shape (M,)\\
`weights.alpha` has shape (M,)\\
`X_predict` has shape (k,N)\\
Returns the converted predictions, in {-1, 1}
"""
function predict(X_predict::Array{T} where T<:Number, X_data::Array{T} where T<:Number,
        Y_data::Array{T} where T<:Number, weights::WeightsSVM)::Array
    if ndims(X_predict) == 1
        X_predict = reshape(X_predict, (1, size(X_predict)[1]))
    end
    @assert ndims(X_predict) == ndims(X_data) == ndims(Y_data) + 1 == 2
    @assert size(X_predict)[2] == size(X_data)[2]
    result = nothing
    if weights.kernel == "linear"
        result = kernel_linear(X_data, X_predict)
    elseif weights.kernel == "polynomial"
        result = kernel_polynomial(X_data, X_predict, d=weights.d, r=weights.r, gamma=weights.gamma)
    elseif weights.kernel == "rbf"
        result = kernel_rbf(X_data, X_predict, gamma=weights.gamma)
    elseif weights.kernel == "sigmoid"
        result = kernel_sigmoid(X_data, X_predict, gamma=weights.gamma, r=weights.r)
    else
        throw(ArgumentError("Error: kernel function $weights.kernel is not recognized"))
    end
    prediction = vec((weights.alpha .* Y_data)' * result)
    prediction .= map(m -> m >= 0 ? 1.0 : -1.0, prediction)
    return prediction
end

"""
Training function for linear SVM\\
If `X_data` has shape (M,N)\\
`Y_data` has shape (M,)\\
`C` is the constraint for SVM decision margin\\
`tol_alpha` sets tolerance for alpha array\\
`tol_error` sets tolerance for error cache\\
`kernel` selects a kernel function, possible values are:
* "rbf"
* "linear"
* "sigmoid"
* "polynomial"

`gamma` selects a way to compute gamma, possible values are:
* "scale"
* "auto"

`degree` defines the degree in `polynomial` kernel function\\
`coef` defines the `r` in `polynomial` and `sigmoid` kernel functions
------
Returns trained alphas as `WeightsSVM` object
"""
function train(X_data::Array{T} where T<:Number, Y_data::Array{T} where T<:Number, C::AbstractFloat;
        tol_alpha::AbstractFloat=0.01, tol_error::AbstractFloat=0.01, kernel::String="rbf", gamma::String="scale",
        degree::AbstractFloat=1.0, coef::AbstractFloat=0.0, verbose::Bool=false)::WeightsSVM
    @assert ndims(X_data) == ndims(Y_data) + 1 == 2
    @assert size(X_data)[1] == size(Y_data)[1]
    X_data = Float64.(X_data)
    Y_data = Float64.(Y_data)
    # gamma is computed the same way sklearn does
    gamma_num = nothing
    if gamma == "scale"
        gamma_num = 1.0 / (size(X_data)[2] * Statistics.var(X_data))
    elseif gamma == "auto"
        gamma_num = 1.0 / size(X_data)[2]
    else
        throw(ArgumentError("Error: gamma $gamma is not recognized, possible values are 'scale' and 'auto'"))
    end
    weights = WeightsSVM(C, 0.0, gamma_num, coef, degree, tol_alpha, tol_error, Float64.(zeros(size(X_data)[1])), -copy(Y_data), kernel)
    num_changed = 0
    examine_all = true
    total_steps = 0
    while (num_changed > 0) || examine_all
        num_changed = 0
        if examine_all
            if verbose
                println("Scanning all training data")
            end
            for i in 1:size(X_data)[1]
                step = learn!(X_data, Y_data, weights, i, verbose=verbose)
                num_changed += step
                if step > 0 && verbose
                    obj = -cost(X_data, weights)
                    println("1 step further, objective = $obj")
                end
            end
        else
            if verbose
                println("Scanning data whose alpha is not at limit")
            end
            alpha_target = [i for (i, m) in enumerate(weights.alpha) if (m != 0.0 && m != weights.C)]
            for i in alpha_target
                step = learn!(X_data, Y_data, weights, i, verbose=verbose)
                num_changed += step
                if step > 0 && verbose
                    obj = -cost(X_data, weights)
                    println("1 step further, objective = $obj")
                end
            end
        end
        if examine_all
            examine_all = false
        elseif num_changed <= 0
            examine_all = true
        end
        total_steps += num_changed
    end
    if verbose
        println("Training Complete\nTotal steps: $total_steps")
    end
    return weights
end

end

max_stars_repo_name: teamclouday/JuliaLearn

max_issues_repo_name: teamclouday/JuliaLearn

max_forks_repo_name: teamclouday/JuliaLearn


text:
using IterativeSolvers
using FactCheck

facts("rsvd_fnkz") do

context("Rank 2 matrix") do
    A = reshape(1:64,8,8)
    S = rsvd_fnkz(A, 2, ϵ=1e-9)
    @fact vecnorm(A - S[:U]*Diagonal(S[:S])*S[:Vt]) --> less_than(1e-9)
end

context("Linearly dependent matrix") do
    A = [collect(1:8) collect(1:8) collect(1:8) collect(1:8)]
    S = rsvd_fnkz(A, 4)
    @fact vecnorm(A - S[:U]*Diagonal(S[:S])*S[:Vt]) --> less_than(1e-9)
end

#SVD of a matrix of zeros should be empty
context("Zero matrix") do
    S = rsvd_fnkz(zeros(10,10), 1)
    @fact S[:U] --> zeros(10, 0)
    @fact S[:S] --> []
    @fact S[:Vt] --> zeros(0, 10)
end

end


text:
using Pkg
Pkg.activate("c:/git/JLBoost/jlboost-test")
@time using JLBoost
@time using DataFrames, JDF

###############################################################################
# fitting GiveMeSomeCredit
###############################################################################
if !isdir("c:/data/GiveMeSomeCredit/cs-training.jdf")
	using CSV
	a = CSV.read("c:/data/GiveMeSomeCredit/cs-training.csv", missingstrings=["NA"])
	rename!(a, Symbol("")=>:column1)
	create_missing!(a, :MonthlyIncome)
	create_missing!(a, :NumberOfDependents)
	type_compress!(a, compress_float=true)
	savejdf("c:/data/GiveMeSomeCredit/cs-training.jdf", a)
end

###############################################################################
# testing find_best_split
###############################################################################
@time a = loadjdf("c:/data/GiveMeSomeCredit/cs-training.jdf")

# this won't work until missing won't work
@time jlboost(a, :SeriousDlqin2yrs)

@time jlboost(a, :SeriousDlqin2yrs, max_leaves = 32)


text:
using FiniteElements
using Test

@testset "FiniteElements.jl" begin
    a = Abc(1.2)
    @test a.hallo == 1.2
end
@testset "FiniteElements2" begin
    a = Abc(1.2)
    @test_broken a.hallo == 1.3
end
@testset "PlattenConstructor" begin
    @test Platte(1.0, 1.0, 10, 10, 7874.0, 80.0, 449.0, 273.15) != nothing
end

@testset "Tempertaturiteration" begin
    p = Platte(1.0, 1.0, 10, 10, 7874.0, 80.0, 449.0, 273.15)
    pneu = tempiter(p, 1.0)
    @test pneu.temp ≈ p.temp
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "Pango"
version = v"1.44.5"

# Collection of sources required to build Pango
sources = [
    "http://ftp.gnome.org/pub/GNOME/sources/pango/1.44/pango-1.44.5.tar.xz" =>
    "8527dfcbeedb4390149b6f94620c0fa64e26046ab85042c2a7556438847d7fc1",
    "https://github.com/mesonbuild/meson/releases/download/0.51.1/meson-0.51.1.tar.gz" =>
    "f27b7a60f339ba66fe4b8f81f0d1072e090a08eabbd6aa287683b2c2b9dd2d82",
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/pango-*/
MESON="$(readlink -f ../meson-0.*/meson.py)"
mkdir build
cd build

cat > cross_compile.txt << EOF
[binaries]
c = '${CC}'
cpp = '${CXX}'
ar = '${AR}'
strip = '${STRIP}'
pkgconfig = '/usr/bin/pkg-config'

[host_machine]
system = 'linux'
cpu_family = 'x86_64'
cpu = 'i686'
endian = 'little'

[target_machine]
system = 'linux'
cpu_family = 'x86'
cpu = 'i686'
endian = 'little'

[paths]
prefix = '${prefix}'
libdir = 'lib'
bindir = 'bin'
EOF

$MESON .. -Dintrospection=false --cross-file cross_compile.txt
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:x86_64, libc=:glibc)
]

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libpango", :libpango)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    "https://github.com/staticfloat/GlibBuilder/releases/download/v2.54.2-2/build.jl"
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: giordano/PangoBuilder

max_issues_repo_name: giordano/PangoBuilder

max_forks_repo_name: giordano/PangoBuilder


text:
using MBSinBLG
using Test

@testset "MBSinBLG.jl" begin
    # Write your tests here.
end

max_stars_repo_name: fernandopenaranda/MBSinBLG

max_issues_repo_name: fernandopenaranda/MBSinBLG

max_forks_repo_name: fernandopenaranda/MBSinBLG


text:
include("automata.jl")

"""
A Deterministic Finite Automaton (DFA)

Has only 3 components
 - the current state the DFA is in
 - the start state (same as `.aut.q₀`)
 - all other components are encapsulated in Automaton
   since they are static and need not have to be mutable
"""
mutable struct DFA <: AbstractFiniteAutomaton
    cur_state::State
    aut::Automaton
end

function DFA(aut::Automaton)
    DFA(aut.q₀, aut)
end

DFA(δ::AbstractDict,
    q₀::State,
    F::AbstractSet{State}) = DFA(Automaton(Dict(k => [v,] for (k, v) in δ),
                                           q₀, F))

"""
    transition!(dfa::DFA, char::Char)::Bool

Applies a transition from current state using the symbol
char. Returns `false` when transition fails (Note: does
not throw an exception for performance reasons) and true
when the transition succeeds while updating the current
state of the DFA.
"""
function transition!(M::DFA, a)::Bool
    key = (M.cur_state, a)

    if !haskey(M.aut.δ, key)
        return false
    end

    M.cur_state = M.aut.δ[key][1]
    return true
end

"""
    reset!(dfa::DFA)

Resets the DFA to it's start state.
"""
function reset!(M::DFA)
    M.cur_state = M.aut.q₀
end

function accepting(M::DFA)::Bool
    M.cur_state ∈ M.aut.F
end

max_issues_repo_name: MadhavJivrajani/DFA-Sim

max_forks_repo_name: MadhavJivrajani/DFA-Sim


text:
using LinearAlgebra

let
	rFTvals = [
		rFFT(nᵢ=(7,7),pᵢ=(1.0,2.0)),
		rFFT(nᵢ=(7,8),pᵢ=(1.0,2.0)),
		rFFT(nᵢ=(8,7),pᵢ=(1.0,2.0)),
		rFFT(nᵢ=(8,8),pᵢ=(1.0,2.0)),
		rFFT(nᵢ=(128,),pᵢ=(2.0,)),
		rFFT(nᵢ=(127,),pᵢ=(2.0,)),
	]
	cFTvals = [
		cFFT(nᵢ=(7,7),pᵢ=(1.0,2.0)),
		cFFT(nᵢ=(7,8),pᵢ=(1.0,2.0)),
		cFFT(nᵢ=(8,7),pᵢ=(1.0,2.0)),
		cFFT(nᵢ=(8,8),pᵢ=(1.0,2.0)),
		cFFT(nᵢ=(128,),pᵢ=(2.0,)),
		cFFT(nᵢ=(127,),pᵢ=(2.0,)),
	]

	for (rFT,cFT) ∈ zip(rFTvals, cFTvals)

		rFTplan, rg = rFT |> x->(plan(x),Grid(x))
		cFTplan, cg = cFT |> x->(plan(x),Grid(x))
		rFTplanᴴ = adjoint(rFTplan)
		cFTplanᴴ = adjoint(cFTplan)

		cF = zeros(Complex{Float64},prod(rg.nxi),prod(rg.nxi))
		for i = 1:prod(rg.nxi)
		    imls = zeros(Complex{Float64}, prod(rg.nxi))
		    imls[i] = 1
		    cF[:,i] = vec(cFTplan * reshape(imls,rg.nxi))
		end
		cFᴴ  = adjoint(cF)
		cFᵀ  = transpose(cF)
		cF⁻¹ = inv(cF)
		cF⁻ᴴ = inv(cFᴴ)

		##
		# Need to check on subset of full fourier matrix
		# These are the tuples (i,j) of the real fft frequencies
		rCi = CartesianIndices(Base.OneTo.(rg.nki))
		# This is the matrix of linear indices into complex fft frequencies
		cLi = LinearIndices(Base.OneTo.(cg.nki))
		# This gets the linear indices of real fft frequencies in the full matrix
		rsub_cLi = cLi[rCi]

		##
		# First covariance modeling in pixel space
		nrm = zeros(prod(rg.nxi), prod(rg.nxi))
		for xi ∈ pixels(rFT)
		    nrm .+= (vec(xi) .- vec(xi)').^2
		end
		nrm = sqrt.(nrm)
		Σ1   = @. exp(-nrm*2)
		M   = Diagonal(vec(1 .+ sin.(pixels(rFT)[1])))
		Σ1  = M * Σ1 * M
		cΓ1 = cF * Σ1 * cFᴴ
		cC1 = cF * Σ1 * cFᵀ
		rΓ1 = cΓ1[vec(rsub_cLi), vec(rsub_cLi)]
		rC1 = cC1[vec(rsub_cLi), vec(rsub_cLi)]

		##
		# Second covariance model in pixel space
		kw = wavenumber(cFT)
		Ck = 1 ./ (1 .+ abs.(kw) ./ 10)
		# wnx = randn(rg.nxi...) ./ sqrt(rg.Ωx) # white noise
		# wnk = rFTplan * wnx
		# fk = sqrt.(Ck) .* wnk
		δk₀ = 1/cg.Ωk
		Σ2  = real.(cF⁻¹ * Diagonal(vec(Ck * δk₀)) * cF⁻ᴴ)
		cΓ2 = cF * Σ2 * cFᴴ
		cC2 = cF * Σ2 * cFᵀ
		rΓ2 = cΓ2[vec(rsub_cLi), vec(rsub_cLi)]
		rC2 = cC2[vec(rsub_cLi), vec(rsub_cLi)]

		##
		#now test
		rFFTimpulses, CI, LI, get_dual_ci = get_rFFTimpulses(rFT)
		ve = f -> vec(f)
		rex = vf -> reshape(vf, rg.nxi)

		EZXᵀ1test  = fill(0.0im, prod(rg.nki), prod(rg.nki))
		EZYᵀ1test  = fill(0.0im, prod(rg.nki), prod(rg.nki))
		EZXᵀ2test  = fill(0.0im, prod(rg.nki), prod(rg.nki))
		EZYᵀ2test  = fill(0.0im, prod(rg.nki), prod(rg.nki))
		for (li,ci) ∈ zip(LI, CI)
		    φ, iφ = rFFTimpulses(ci)
		    EZXᵀ1test[:,li] = ve(rFTplan * rex(Σ1 * ve(rFTplanᴴ * φ)))
		    EZYᵀ1test[:,li] = ve(rFTplan * rex(Σ1 * ve(rFTplanᴴ * iφ)))
		    EZXᵀ2test[:,li] = ve(rFTplan * rex(Σ2 * ve(rFTplanᴴ * φ)))
		    EZYᵀ2test[:,li] = ve(rFTplan * rex(Σ2 * ve(rFTplanᴴ * iφ)))
		end

		Σ1XX = real.(EZXᵀ1test)
		Σ1YY = imag.(EZYᵀ1test)
		Σ1YX = imag.(EZXᵀ1test)
		Σ1XY = real.(EZYᵀ1test)

		Σ2XX = real.(EZXᵀ2test)
		Σ2YY = imag.(EZYᵀ2test)
		Σ2YX = imag.(EZXᵀ2test)
		Σ2XY = real.(EZYᵀ2test)

		rΓ1test2 = complex.(Σ1XX .+ Σ1YY, Σ1YX .- Σ1XY)
		rC1test2 = complex.(Σ1XX .- Σ1YY, Σ1YX .+ Σ1XY)
		rΓ2test2 = complex.(Σ2XX .+ Σ2YY, Σ2YX .- Σ2XY)
		rC2test2 = complex.(Σ2XX .- Σ2YY, Σ2YX .+ Σ2XY)

		for (A,B) ∈ zip((rΓ1test2, rC1test2, rΓ2test2, rC2test2), (rΓ1, rC1, rΓ2, rC2))
	    	@test sum(abs2, A.- B) ≈ 0.0 atol=1e-15
		end
	end # for
end # let


text:
#=

    abc.jl

    Contains code to conduct approximate Bayesian computation (ABC)

    Author:     Alexander P. Browning
                ======================
                School of Mathematical Sciences
                Queensland University of Technology
                ======================
                ap.browning@icloud.com
                alexbrowning.me

=#

##############################################################
## ABC Sequential Monte-Carlo (SMC)
##############################################################
"""
    abc_smc(dist,prior;n=100,α=0.5,S₀=100,c=0.01,ptarget=0.005,maxsteps=20,f=2.0)

Perform ABC SMC.
"""
function abc_smc(dist,prior;n=100,kwargs...)
    P = sample(prior,dist,n)    # Initialise using prior
    print("Initialised $n particles. $(round(minimum(P).d,sigdigits=4)) ≤ d ≤ $(round(maximum(P).d,sigdigits=4)).\n")
    return abc_smc(P,dist,prior;kwargs...)
end
function abc_smc(P,dist,prior;step=1,α=0.5,ptarget=0.01,maxsteps=30,f=2.0)

    # Loop until target reached
    p = 1.0
    while p > ptarget && step ≤ maxsteps
        @printf("Step %-3s (ε = %-8s\b...",
            "$step",
            "$(round(sort(distances(P))[(Int ∘ floor)(α * length(P))],sigdigits=4))",
        )

        p = abc_smc_step!(P,dist,prior,α,f)
        step += 1

        @printf("\b\b\b dmin = %-8s p = %-8s ESS = %-8s)\n",
            "$(round(minimum(P).d,sigdigits=4))",
            "$(round(p,sigdigits=2))",
            "$(round(1 / sum(weights(P).^2),sigdigits=5))"
        )

    end

    # Finish
    return P
end

"""
    abc_smc_step!(P,dist,prior,α,f=2.0)

Perform a single step of ABC SMC, where α * n particles are dicarded.

Proposal is estimated as f * Σ, where Σ is the (weighted) covariance matrix of the non-dicarded particles.

"""
function abc_smc_step!(P,dist,prior,α,f=2.0)
    N = length(P); Nα = (Int ∘ floor)(α * N)

    sort!(P)                # Sort particles
    ε = P[N-Nα].d           # Current threshold
    Σ = f*cov(P[1:N-Nα])    # Proposal covariance
    prop = MvNormal(Σ)      # Proposal

    # Resample and update particles
    tries = 0
    @threads for j = (N - Nα + 1) : N
        θs,d = similar(P[1].θ),Inf
        while d > ε
            θs = copy(sample(P[1:N-Nα],weights(P[1:N-Nα])).θ) + rand(prop)
            if insupport(prior,θs)
                d = dist(θs)
                tries += 1
            end
        end
        P[j] = Particle(θs,d,pdf(prior,θs) / sum(p.w * pdf(prop,θs - p.θ) for p in P[1:N-Nα]))
    end

    # Reweight
    w = weights(P)
    w[1:N-Nα] /= sum(w[1:N-Nα])         # Normalise
    w[N-Nα+1:N] /= sum(w[N-Nα+1:N])     # Normalise
    w[1:N-Nα] *= (N - Nα) / N           # Re-weight
    w[N-Nα+1:N] *= Nα / N               # Re-weight
    [p.w = wᵢ for (p,wᵢ) in zip(P,w)]

    # Report acceptance probability (exclude trials where we propose a value outside the prior)
    return Nα / tries
end

##############################################################
## ABC Markov-chain Monte-Carlo (MCMC)
##############################################################
"""
    abc_mcmc(dist,prior,θ₀,Σ,ε;n=1000,nchains=nthreads(),...)

Perform ABC MCMC using the Metropolis-Hastings algorithm (implemented through `AdvancedMH.jl`).
"""
function abc_mcmc(dist,prior,θ₀,Σ,K::Function;n=1000,nchains=Threads.nthreads(),kwargs...)
    parallel_args = nchains == 1 ? n : (MCMCThreads(),n,nchains)
    sample(
        DensityModel(θ -> insupport(prior,θ) ? log(K(dist(θ))) : -Inf),
        RWMH(MvNormal(Σ)),
        parallel_args...;
        init_params=θ₀,
        chain_type=Chains,
        kwargs...
    )
end
abc_mcmc(dist,prior,θ₀,Σ,ε::Number;kwargs...) = abc_mcmc(dist,prior,θ₀,Σ,d -> d ≤ ε ? 1.0 : 0.0;kwargs...)

##############################################################
## ABC MCMC Discrepancy Metrics
##############################################################

logistic(d,ε,k) = 1 / (1 + exp(k * (d - ε)))
logit(p,ε,k) = ε + log((1 - p) / p) / k

##############################################################
## MCMCChains.Chains helper functions
##############################################################

# Convert MCMC chain output to a matrix
Base.Matrix(C::Chains) = permutedims(hcat([Matrix(C[:,i,:])[:] for i = 1:size(C,2)-1]...))

acceptance_rate(C::Chains) = mean(hcat([cumsum(C[1:end-1,1,i] .!= C[2:end,1,i],dims=1)[:] ./ (1:size(C,1)-1) for i = 1:size(C,3)]...),dims=2)[:]

MCMCChains.missing_datetime(T::Type) = missing

StatsBase.sample(C::Chains,n) = (CM = Matrix(C); [CM[:,i] for i in sample(1:size(CM,2),n)])
StatsBase.sample(C::Chains) = sample(C,1)[1]

max_stars_repo_name: ap-browning/internalisation

max_issues_repo_name: ap-browning/internalisation

max_forks_repo_name: ap-browning/internalisation


text:
@reexport module _DataStructure

import Base: getindex, push!, pop!
import Nemo: cols, rows

export
    flatten

# https://groups.google.com/d/msg/julia-users/1QrIhbRA8hs/9PcNeO2N9wQJ
flatten(a::Array{T,1}) where {T} = any(x -> Base.isa(x, Array), a) ? flatten(vcat(map(flatten, a)...)) : a
flatten(a::Array{T}) where {T} = reshape(a, prod(size(a)))
flatten(a) = a

include("iterator.jl")
include("vector.jl")
include("matrix.jl")
include("intset.jl")
include("stack.jl")

end


text:
"""

    lossy_fit(xs,ys;Rdc=1e0,Rs=1e-6,L0=1e-9,epsilon=1e-3,C0=1e-12,show_trace=false,maxIter=1000)

lossy_fit accepts two positional arguments

xs::Float{64,1} - Array of sampled frequencies
ys::Float{64,1} - Array of frequency response values

and 7 keyword arguments

Rdc::Float64
Rs::Float64
L0::Float64
epsilon::Float64
C0::Float64
show_trace::Bool
maxIter::Int32

where the first 5 keyword arguments are the sizes (in SI units) of the guesses
of each of the fitted parameter values and the last two control the behavior of
LsqFit.jl.

lossy_fit computes a fit of the parameters Rdc (DC resistive loss), Rs (skin effect loss),
L0 (series inductance), loss tangent, C0 (shunt capacitance) provided the
frequency-domain response of a device under test assuming the model described in
[[RINEHART2017]].
"""
function lossy_fit(xs,ys;Rdc=1e0,Rs=1e-6,L0=1e-9,epsilon=1e-3,C0=1e-12,show_trace=false,maxIter=1000)
    z(p) = f -> tline_z_fd(p...)(f)
    # re-define the model to conform to LsqFit expectations
    model(f,p) = vcat(real(z(p).(f)), imag(z(p).(f)))

    # ForwardDiff.jacobian takes a function that only depends on the parameters to
    # be solved for.
    function model_jac(w,p)
        return ForwardDiff.jacobian(pj -> model(w,pj), p)
    end
    # guess
    p0 = [Rdc, Rs, L0, epsilon, C0]
    fit = LsqFit.curve_fit(model, model_jac, xs, ys, p0,
                           show_trace=show_trace,
                           lower=zeros(5)+eps(),
                           maxIter=maxIter,
                        )
    fit_errors = LsqFit.estimate_errors(fit, .95)
    return Dict("vals" => fit.param, "errors" => fit_errors)
end

"""
    impedance(ZR,ZL,ZG,ZC)

returns the impedance function to be fit.
"""
function impedance(ZR,ZL,ZG,ZC)
    frac(p) = f -> (ZR(p[1],p[2])(f)+ZL(p[3])(f)) / (ZG(p[4],p[5])(f) + ZC(p[5])(f))
    function Z(p)
        return function(f)
            f==0?sqrt(p[1]):abs(frac(p)(f))*cis(angle(frac(p)(f)))
        end
    end
    return Z
end


text:
is_optional(member::SpecStructMember) = member.name == :pNext || member.requirement ∈ [OPTIONAL, POINTER_OPTIONAL] || is_default_count(member)
is_optional(param::SpecFuncParam) = param.requirement ∈ [OPTIONAL, POINTER_OPTIONAL]

"""
Represent an integer that gives the start of a C pointer.
"""
function is_pointer_start(spec::Spec)
    params = children(parent_spec(spec))
    any(params) do param
        !isempty(param.arglen) &&
            spec.type == :UInt32 &&
            string(spec.name) == string("first", uppercasefirst(replace(string(param.name), r"Count$" => "")))
    end
end

is_semantic_ptr(type) = is_ptr(type) || type == :Cstring
needs_deps(spec::SpecStruct) = any(is_semantic_ptr, spec.members.type) && !spec.is_returnedonly
must_return_success_code(spec::SpecFunc) = length(spec.success_codes) > 1 && :VK_INCOMPLETE ∉ spec.success_codes
must_repeat_while_incomplete(spec::SpecFunc) = !must_return_success_code(spec) && :VK_INCOMPLETE ∈ spec.success_codes
is_data_with_retrievable_size(spec::SpecFuncParam) = is_data(spec) && len(spec).requirement == POINTER_REQUIRED
is_opaque_data(spec) = is_data(spec) && len(spec).requirement ≠ POINTER_REQUIRED
is_implicit_return(spec::SpecFuncParam) =
    !is_opaque_data(spec) &&
    !spec.is_constant &&
    is_ptr(spec.type) &&
    !is_length(spec) &&
    spec.type ∉ extension_types &&
    ptr_type(spec.type) ∉ extension_types
has_implicit_return_parameters(spec::SpecFunc) = any(is_implicit_return, children(spec))

function is_default_count(spec::Spec)
    is_length(spec) && @match parent(spec) begin
        # see `descriptorCount` at https://www.khronos.org/registry/vulkan/specs/1.1-extensions/html/vkspec.html#VkWriteDescriptorSet
        :VkWriteDescriptorSet => true
        _ => false
    end
end

function default_count(spec::Spec)
    @match parent(spec) begin
        :VkWriteDescriptorSet => :(max($((:(pointer_length($(wrap_identifier(arg)))) for arg in arglen(spec))...)))
    end
end

function is_specific_count(spec::Spec)
    @match spec begin
        GuardBy(is_length) => @match parent(spec) begin
            :VkDescriptorSetLayoutBinding => true
            _ => false
        end
        _ => false
    end
end


text:
using BinaryProvider # requires BinaryProvider 0.3.0 or later

dependencies = [
    "build_Zlib.v1.2.11.jl",
    "build_libpng.v1.0.0.jl",
    "build_libjpeg.v9.0.0-b.jl",
    "build_libtiff.v4.0.9.jl"
]

for elem in dependencies
    # it's a bit faster to run the build in an anonymous module instead of
    # starting a new julia process
    m = Module(:__anon__)
    Core.eval(m, :(Main.include($(joinpath(@__DIR__, elem)))))
end

# Parse some basic command-line arguments
const verbose = "--verbose" in ARGS
const prefix = Prefix(get([a for a in ARGS if a != "--verbose"], 1, joinpath(@__DIR__, "usr")))
products = [
    LibraryProduct(prefix, String["libMagickWand"], :libwand),
]
# Download binaries from hosted location
version = v"6.9.10-4"

bin_prefix = "https://github.com/JuliaIO/ImageMagickBuilder/releases/download/v$version"

# Listing of files generated by BinaryBuilder:
download_info = Dict(
    Linux(:aarch64, :glibc) => ("$bin_prefix/ImageMagick.v$version.aarch64-linux-gnu.tar.gz", "8780c6eec7f2c34dc35355fdb332e5c0fcd021105293e8e4f6b6e75032e35b97"),
    Linux(:aarch64, :musl) => ("$bin_prefix/ImageMagick.v$version.aarch64-linux-musl.tar.gz", "629a373cf7990232b8272024aa649bd44443930706c43c060b1b384913b840f7"),
    Linux(:armv7l, :glibc, :eabihf) => ("$bin_prefix/ImageMagick.v$version.arm-linux-gnueabihf.tar.gz", "4b5500f31a8fcf89137f9ca40ebd60bdb3499cdc8eab4e321a1429c61b3ecd45"),
    Linux(:armv7l, :musl, :eabihf) => ("$bin_prefix/ImageMagick.v$version.arm-linux-musleabihf.tar.gz", "2232234af05d29c85378d589e04f29b38c35e05d25d6c732f7a01a9ff394a0f2"),
    Linux(:i686, :glibc) => ("$bin_prefix/ImageMagick.v$version.i686-linux-gnu.tar.gz", "34db87e26bff290dedec73394e8266aaaaca816984513b57766c933cf6ec1a26"),
    Linux(:i686, :musl) => ("$bin_prefix/ImageMagick.v$version.i686-linux-musl.tar.gz", "64d9c1867ed370407adfae04ed7ef48bfe5f29ba843069b5c76299bd5a5da155"),
    Windows(:i686) => ("$bin_prefix/ImageMagick.v$version.i686-w64-mingw32.tar.gz", "382bc4d3c9e709711eafa2aa795eb3ba1c88dacefb6932ef32c44675a4b11132"),
    Linux(:powerpc64le, :glibc) => ("$bin_prefix/ImageMagick.v$version.powerpc64le-linux-gnu.tar.gz", "4ba60630e955f6c7fe952c21d5ea403bfe3725cdd3e7730937dc604f6c33ceb7"),
    MacOS(:x86_64) => ("$bin_prefix/ImageMagick.v$version.x86_64-apple-darwin14.tar.gz", "b4e6b9a32e1ebaf314295fb569022014f64de0f25518f1fbea1f2588a94cbe28"),
    Linux(:x86_64, :glibc) => ("$bin_prefix/ImageMagick.v$version.x86_64-linux-gnu.tar.gz", "908fd88a0e52f038f3f19d192036df3c325d8560a8290a5abb9bd2d9307d8aa8"),
    Linux(:x86_64, :musl) => ("$bin_prefix/ImageMagick.v$version.x86_64-linux-musl.tar.gz", "d7c083bf912cedff59cfac4740c9a15977882c92896f181dd7114f3d5508a126"),
    FreeBSD(:x86_64) => ("$bin_prefix/ImageMagick.v$version.x86_64-unknown-freebsd11.1.tar.gz", "ca8f099f7c3639e4063213ac00a17a72b01f45a6713ec2306cd5af5d9398f989"),
    Windows(:x86_64) => ("$bin_prefix/ImageMagick.v$version.x86_64-w64-mingw32.tar.gz", "c3d0d85f0c44ae3d10b42a80514a62ea566b84d30feaf9dbe574bb73aa41b085"),
)

# Install unsatisfied or updated dependencies:
unsatisfied = any(!satisfied(p; verbose=verbose) for p in products)
if haskey(download_info, platform_key())
    url, tarball_hash = download_info[platform_key()]
    if unsatisfied || !isinstalled(url, tarball_hash; prefix=prefix)
        # Download and install binaries
        install(url, tarball_hash; prefix=prefix, force=true, verbose=verbose)
    end
elseif unsatisfied
    # If we don't have a BinaryProvider-compatible .tar.gz to download, complain.
    # Alternatively, you could attempt to install from a separate provider,
    # build from source or something even more ambitious here.
    error("Your platform $(triplet(platform_key())) is not supported by this package!")
end

# Write out a deps.jl file that will contain mappings for our products
write_deps_file(joinpath(@__DIR__, "deps_im.jl"), products)

function include_deps(name)
    """
    module $name
        path = joinpath(@__DIR__, $(repr(string("deps_", name, ".jl"))))
        isfile(path) || error("$name wasn't build correctly. Please run Pkg.build(\\\"ImageMagick\\\")")
        Main.include(path)
    end
    using .$name
    """
end

open("deps.jl", "w") do io
    for dep in (:zlib, :png, :jpeg, :tiff, :im)
        println(io, include_deps(dep))
    end
    println(io, """
    libversion() = $(repr(version))
    const libwand = im.libwand

    function check_deps()
        zlib.check_deps()
        png.check_deps()
        jpeg.check_deps()
        tiff.check_deps()
        im.check_deps()
    end
    """)
end


text:
#=

Modela diferentes tipos de anemômetros

=#

"""
    HWSensor(I, E, R, T)

Condições de operação do sensor de um anemômetro de fio quente.
"""
struct HWSensor
    "Corrente passando no sensor"
    I::Float64
    "Tenão no sensor"
    E::Float64
    "Resistência do sensor"
    R::Float64
    "Temmperatura do sensor"
    T::Float64
    "Potência dissipada no sensor"
    P::Float64
    HWSensor(I, E, R, T) = new(I, E, R, T, I*E)
end

"""
    const_current(mA, R::Thermistor, U=5.0, d=3.0, Ta=20.0, Pa=101.3, dterr=1e-6, nmax=20000)

Modela um anemômetro trabalhando a corrente constante. Não simula circuitos elétricos.
Utiliza um esquema iterativo.

Argumentos:
 * `mA` Corrente em mA
 * `R` Objeto que modela o termistor
 * `U` Velocidade em m/s
 * `d` Diâmetro em mm
 * `Ta` Temperatura ao longe em °C
 * `Pa` Pressão atmosférica em kPa
 * `dterr` Erro admissível para convergência da temperatura.
 * `nmax` Número máximo de iterações.

```jldoctest
julia> R = Thermistor(5e3, 3200, 20)
THW.Thermistor(5000.0, 3200.0, 293.15)

julia> const_current(20, R, 5.0, 2.0)
THW.HWSensor(0.02, 9.395209322307162, 469.7604661153581, 101.07772779650266, 0.18790418644614323)

julia> const_current(20, R, 1.0, 2.0)
THW.HWSensor(0.02, 5.872996120670518, 293.64980603352586, 122.8352988426391, 0.11745992241341036)

julia> const_current(20, R, 10.0, 2.0)
THW.HWSensor(0.02, 11.70131737284239, 585.0658686421194, 91.71177677466494, 0.2340263474568478)

julia> const_current(20, R, 10.0, 10.0)
THW.HWSensor(0.02, 44.971641701710844, 2248.582085085542, 43.15632417996514, 0.8994328340342169)
```
"""
function const_current(mA, R::Thermistor, U=5.0, d=3.0, Ta=20.0, Pa=101.3, dterr=1e-6, nmax=20000)

    I = mA/1000

    D = d/1000
    r = D/2
    A = 4π*r^2

    Tw = Ta + 0.01

    i = 0
    P = 0.0
    for i = 1:nmax
        h = hconvect(U, D, Tw, Ta, Pa)
        #println(h)
        P = R(Tw)*I^2
        Tw1 = P/(h*A) + Ta
        dt = Tw1-Tw
        if abs(dt) < dterr
            Tw = Tw1
            break
        end
        Tw = Tw + 0.4*dt
    end

    res = R(Tw)
    E = res*I
    return HWSensor(I, E, res, Tw)
end

"""
    const_temperature(Tw, R::Thermistor, U=5.0, d=3.0, Ta=20.0, Pa=101.3)

Modela um anemômetro trabalhando a temperatura constante. Não simula circuitos elétricos.

Argumentos:
 * `Tw` Temperatura de operação do sensor
 * `R` Objeto que modela o termistor
 * `U` Velocidade em m/s
 * `d` Diâmetro em mm
 * `Ta` Temperatura ao longe em °C
 * `Pa` Pressão atmosférica em kPa

```jldoctest
julia> R = Thermistor(5e3, 3200, 20)
THW.Thermistor(5000.0, 3200.0, 293.15)

julia> const_temperature(100.0, R, 5.0, 2.0)
THW.HWSensor(0.019627341883449854, 9.45069202339691, 481.50646580247894, 100.0, 0.18549196337840362)

julia> const_temperature(100.0, R, 1.0, 2.0)
THW.HWSensor(0.01383013241982859, 6.65929818305195, 481.50646580247894, 100.0, 0.09209897569473241)

julia> const_temperature(100.0, R, 10.0, 2.0)
THW.HWSensor(0.023241402160834233, 11.190885414757389, 481.50646580247894, 100.0, 0.26009186846019067)

julia> const_temperature(100.0, R, 20.0, 2.0)
THW.HWSensor(0.027741857167247384, 13.357883599398457, 481.50646580247894, 100.0, 0.3705724988712284)

julia> const_temperature(100.0, R, 20.0, 3.0)
THW.HWSensor(0.037797691428912955, 18.19983281542853, 481.50646580247894, 100.0, 0.6879116648153717)
```
"""
function const_temperature(Tw, R::Thermistor, U=5.0, d=3.0, Ta=20.0, Pa=101.3)

    D = d/1000
    r = D/2
    A = 4π*r^2

    Rw = R(Tw)
    h = hconvect(U, D, Tw, Ta, Pa)
    i = sqrt(h*A*(Tw-Ta)/Rw)

    return HWSensor(i, Rw*i, Rw, Tw)
end

function const_temperature_cyl(Tw, R::AbstractResistor, U=5.0, d=5.0, l=3, Ta=20.0, Pa=101.325)

    D = d/1e6
    L = l / 1000

    r = D/2

    A = L * π * D

    Rw = R(Tw)
    h = hconvectcyl(U, D, Tw, Ta, Pa)
    i = sqrt(h*A*(Tw-Ta)/Rw)

    return i*1000
end


text:
using Test
using LinearAlgebra
using FlightMechanics

@testset "quaternion <-> euler" begin
    quat = (0.8660254037844387, 0.0, 0.5, 0.0)
    euler_exp = [0.0, 1.04719755, 0.0]

    euler = quaternion2euler(quat...)
    @test isapprox(euler, euler_exp)

    quat_exp = [quat...]
    quat = euler2quaternion(euler_exp...)
    @test isapprox(quat, quat_exp)

    quat = [0.5, 0.5, 0.0, 0.0]
    euler_exp = [0.0, 0.0, pi/2.0]

    euler = quaternion2euler(quat...)
    @test isapprox(euler, euler_exp)

    quat_exp = [quat...]
    quat = euler2quaternion(euler_exp...)
    @test isapprox(quat, quat_exp / norm(quat_exp))

    psi, theta, phi = pi/4.0, pi/6.0, pi/12.0
    quat = euler2quaternion(psi, theta, phi)
    euler = quaternion2euler(quat...)
    @test isapprox([psi, theta, phi], euler)

    quat = [0.5, 0.1, 0.2, 0.7]
    quat = quat / norm(quat)
    euler = quaternion2euler(quat...)
    quat3 = euler2quaternion(euler...)
    @test isapprox(quat, quat3)
end

@testset "body <-> horizon" begin
    ones_ = [1.0, 1.0, 1.0]

    # body2hor Euler
    # no rotation
    @test ones_ ≈ body2hor(ones_..., 0., 0., 0.)

    angles1 = [0., 45*pi/180., 0.]
    angles2 = [0., 0., 45*pi/180]
    angles3 = [45*pi/180, 0., 0.]
    exp_b2h_1 = [2*0.70710678118654757, 1, 0]
    exp_b2h_2 = [1, 0, 2*0.70710678118654757]
    exp_b2h_3 = [0, 2*0.70710678118654757, 1]

    @test exp_b2h_1 ≈ body2hor(ones_..., angles1...)
    @test exp_b2h_2 ≈ body2hor(ones_..., angles2...)
    @test exp_b2h_3 ≈ body2hor(ones_..., angles3...)

    # hor2body Euler
    @test ones_ ≈ hor2body(ones_..., 0., 0., 0.)
    @test ones_ ≈ hor2body(exp_b2h_1..., angles1...)
    @test ones_ ≈ hor2body(exp_b2h_2..., angles2...)
    @test ones_ ≈ hor2body(exp_b2h_3..., angles3...)

    # body2hor quaternion
    @test ones_ ≈ body2hor(ones_..., 0., 0., 0.)

    quat1 = euler2quaternion(angles1...)
    quat2 = euler2quaternion(angles2...)
    quat3 = euler2quaternion(angles3...)

    @test exp_b2h_1 ≈ body2hor(ones_..., quat1...)
    @test exp_b2h_2 ≈ body2hor(ones_..., quat2...)
    @test exp_b2h_3 ≈ body2hor(ones_..., quat3...)

    # hor2body quaternionr
    @test ones_ ≈ hor2body(ones_..., 0., 0., 0.)
    @test ones_ ≈ hor2body(exp_b2h_1...,  quat1...)
    @test ones_ ≈ hor2body(exp_b2h_2...,  quat2...)
    @test ones_ ≈ hor2body(exp_b2h_3...,  quat3...)

    # rotation matrix body to hor (euler)
    r_hb1 = rot_matrix_body2hor(angles1...)
    @test exp_b2h_1 ≈ r_hb1 * ones_
    r_hb2 = rot_matrix_body2hor(angles2...)
    @test exp_b2h_2 ≈ r_hb2 * ones_
    r_hb3 = rot_matrix_body2hor(angles3...)
    @test exp_b2h_3 ≈ r_hb3 * ones_

    # rotation matrix body <-> hor (quaternion)
    r_hb1q = rot_matrix_body2hor(quat1...)
    @test exp_b2h_1 ≈ r_hb1q * ones_
    r_hb2q = rot_matrix_body2hor(quat2...)
    @test exp_b2h_2 ≈ r_hb2q * ones_
    r_hb3q = rot_matrix_body2hor(quat3...)
    @test exp_b2h_3 ≈ r_hb3q * ones_

    # r_hb equivalent with euler and quaternion
    @test isapprox(r_hb1, r_hb1q)
    @test isapprox(r_hb2, r_hb2q)
    @test isapprox(r_hb3, r_hb3q)

    # rotation matrix hor to body (euler)
    r_bh1 = rot_matrix_hor2body(angles1...)
    @test ones_ ≈ r_bh1 * exp_b2h_1
    r_bh2 = rot_matrix_hor2body(angles2...)
    @test ones_ ≈ r_bh2 * exp_b2h_2
    r_bh3 = rot_matrix_hor2body(angles3...)
    @test ones_ ≈ r_bh3 * exp_b2h_3

    # rotation matrix hor2body (quaternion)
    r_bh1q = rot_matrix_hor2body(quat1...)
    @test ones_ ≈ r_bh1q * exp_b2h_1
    r_bh2q = rot_matrix_hor2body(quat2...)
    @test ones_ ≈ r_bh2q * exp_b2h_2
    r_bh3q = rot_matrix_hor2body(quat3...)
    @test ones_ ≈ r_bh3q * exp_b2h_3

    # r_bh equivalent with euler and quaternion
    @test isapprox(r_bh1, r_bh1q)
    @test isapprox(r_bh2, r_bh2q)
    @test isapprox(r_bh3, r_bh3q)

    # Test that quaternion and euler produce the same transformation
    psi, theta, phi = pi/4.0, pi/6.0, pi/12.0
    quat = euler2quaternion(psi, theta, phi)
    xh, yh, zh = 100.0, 10.0, -1.0

    xyz_b_e = hor2body(xh, yh, zh, psi, theta, phi)
    xyz_b_q = hor2body(xh, yh, zh, quat...)

    @test isapprox(xyz_b_e, xyz_b_q)

    xyz_h_e = body2hor(xyz_b_e..., psi, theta, phi)
    xyz_h_q = body2hor(xyz_b_e..., quat...)

    @test isapprox(xyz_h_e, xyz_h_q)
end

@testset "wind <-> hor/body" begin
    ones_ = [1.0, 1.0, 1.0]
    angles1 = [0., 45*pi/180., 0.]
    angles2 = [0., 0., 45*pi/180]
    angles3 = [45*pi/180, 0., 0.]
    exp_b2h_1 = [2*0.70710678118654757, 1, 0]
    exp_b2h_2 = [1, 0, 2*0.70710678118654757]
    exp_b2h_3 = [0, 2*0.70710678118654757, 1]
    #wind2hor
    @test ones_ ≈ wind2hor(ones_..., 0., 0., 0.)
    @test exp_b2h_1 ≈ wind2hor(ones_..., angles1...)
    @test exp_b2h_2 ≈ wind2hor(ones_..., angles2...)
    @test exp_b2h_3 ≈ wind2hor(ones_..., angles3...)
    #hor2wind
    @test ones_ ≈ hor2wind(ones_..., 0., 0., 0.)
    @test ones_ ≈ hor2wind(exp_b2h_1..., angles1...)
    @test ones_ ≈ hor2wind(exp_b2h_2..., angles2...)
    @test ones_ ≈ hor2wind(exp_b2h_3..., angles3...)
    #wind2body
    @test ones_ ≈ wind2body(ones_..., 0., 0.)
    @test [0, 1, 2*0.70710678118654757] ≈ wind2body(ones_..., 45*pi/180., 0.)
    @test exp_b2h_3 ≈ wind2body(ones_..., 0., 45*pi/180.)
    #body2wind
    @test ones_ ≈ body2wind(ones_..., 0., 0.)
    @test ones_ ≈ body2wind(0, 1, 2*0.70710678118654757, 45*pi/180., 0.)
    @test ones_ ≈ body2wind(exp_b2h_3..., 0., 45*pi/180.)
end

@testset "hor <-> ECEF" begin
    xecef, yecef, zecef = 1.0, 10.0, 100.0
    lat, lon = 0.0, 0.0
    exp_xyz_hor = [100.0, 10.0 ,-1.0]

    xyz_hor =  ecef2hor(xecef, yecef, zecef, lat, lon)
    @test isapprox(xyz_hor, exp_xyz_hor)

    r_hecef = rot_matrix_ecef2hor(lat, lon)
    @test isapprox(exp_xyz_hor, r_hecef * [xecef, yecef, zecef])

    exp_xyz_ecef = [xecef, yecef, zecef]
    xyz_ecef = hor2ecef(exp_xyz_hor..., lat, lon)
    @test isapprox(xyz_ecef, exp_xyz_ecef)

    r_ecefh = rot_matrix_hor2ecef(lat, lon)
    @test isapprox(exp_xyz_ecef, r_ecefh * exp_xyz_hor)

    lat, lon = pi/2.0, 0.0
    exp_xyz_hor = [-1.0, 10.0 ,-100.0]

    xyz_hor =  ecef2hor(xecef, yecef, zecef, lat, lon)
    @test isapprox(xyz_hor, exp_xyz_hor)

    r_hecef = rot_matrix_ecef2hor(lat, lon)
    @test isapprox(exp_xyz_hor, r_hecef * [xecef, yecef, zecef])

    exp_xyz_ecef = [xecef, yecef, zecef]
    xyz_ecef = hor2ecef(exp_xyz_hor..., lat, lon)
    @test isapprox(xyz_ecef, exp_xyz_ecef)

    r_ecefh = rot_matrix_hor2ecef(lat, lon)
    @test isapprox(exp_xyz_ecef, r_ecefh * exp_xyz_hor)

    lat, lon = 0.0, pi/2.0
    exp_xyz_hor = [100.0, -1.0 ,-10.0]

    xyz_hor =  ecef2hor(xecef, yecef, zecef, lat, lon)
    @test isapprox(xyz_hor, exp_xyz_hor)

    r_hecef = rot_matrix_ecef2hor(lat, lon)
    @test isapprox(exp_xyz_hor, r_hecef * [xecef, yecef, zecef])

    exp_xyz_ecef = [xecef, yecef, zecef]
    xyz_ecef = hor2ecef(exp_xyz_hor..., lat, lon)
    @test isapprox(xyz_ecef, exp_xyz_ecef)

    r_ecefh = rot_matrix_hor2ecef(lat, lon)
    @test isapprox(exp_xyz_ecef, r_ecefh * exp_xyz_hor)

    lat, lon = pi/2.0, pi/2.0
    exp_xyz_hor = [-10.0, -1.0 ,-100.0]

    xyz_hor =  ecef2hor(xecef, yecef, zecef, lat, lon)
    @test isapprox(xyz_hor, exp_xyz_hor)

    r_hecef = rot_matrix_ecef2hor(lat, lon)
    @test isapprox(exp_xyz_hor, r_hecef * [xecef, yecef, zecef])

    exp_xyz_ecef = [xecef, yecef, zecef]
    xyz_ecef = hor2ecef(exp_xyz_hor..., lat, lon)
    @test isapprox(xyz_ecef, exp_xyz_ecef)

    r_ecefh = rot_matrix_hor2ecef(lat, lon)
    @test isapprox(exp_xyz_ecef, r_ecefh * exp_xyz_hor)
end

@testset "body <-> ECEF" begin
    # uses ecef2hor and hor2body so not intensive testing necessary
    xecef, yecef, zecef = 1.0, 10.0, 100.0
    psi, theta, phi = pi/4.0, pi/6.0, pi/12.0
    quat = euler2quaternion(psi, theta, phi)
    lat, lon = pi/3, pi/6
    xh, yh, zh = ecef2hor(xecef, yecef, zecef, lat, lon)
    exp_xyz_b = hor2body(xh, yh, zh, psi, theta, phi)

    xyz_body = ecef2body(xecef, yecef, zecef, lat, lon , psi, theta, phi)
    @test isapprox(xyz_body, exp_xyz_b)

    r_becef = rot_matrix_ecef2body(lat, lon , psi, theta, phi)
    @test isapprox(r_becef * [xecef, yecef, zecef], exp_xyz_b)

    xyz_body = ecef2body(xecef, yecef, zecef, lat, lon , quat...)
    @test isapprox(xyz_body, exp_xyz_b)

    r_becef_q = rot_matrix_ecef2body(lat, lon , quat...)
    @test isapprox(r_becef_q * [xecef, yecef, zecef], exp_xyz_b)

    exp_xyz_ecef = [xecef, yecef, zecef]

    xyz_ecef = body2ecef(exp_xyz_b..., lat, lon, psi, theta, phi)
    @test isapprox(xyz_ecef, exp_xyz_ecef)

    r_ecefb = rot_matrix_body2ecef(lat, lon , psi, theta, phi)
    @test isapprox(r_ecefb * exp_xyz_b, exp_xyz_ecef)

    xyz_ecef = body2ecef(exp_xyz_b..., lat, lon, quat...)
    @test isapprox(xyz_ecef, exp_xyz_ecef)

    r_ecefb_q = rot_matrix_body2ecef(lat, lon , quat...)
    @test isapprox(r_ecefb_q * exp_xyz_b, exp_xyz_ecef)
end

@testset "llh <-> ECEF" begin
    # llh ECEF (using data from
    # - [1] Bowring, B. R. (1976). Transformation from spatial to geographical
    # coordinates. Survey review, 23(181), 323-327.)
    using FlightMechanics: Ellipsoid

    # From [1] Example 1 (page 325)
    x = 4114496.258  # m
    y = 0.0          # m
    z = 4870157.031  # m
    # Ellipsoid does not match any known standard
    a_bowring1976 = 6378249.145  # m
    b_bowring1976 = 6356514.870  # m
    finv_bowring1976 = a_bowring1976 / (a_bowring1976 - b_bowring1976)
    Bowring1976 = Ellipsoid(a_bowring1976, finv_bowring1976)

    exp_llh = [deg2rad(50.0), 0, 10000.0]
    llh = ecef2llh(x, y, z; ellipsoid=Bowring1976)
    @test isapprox(llh, exp_llh)

    exp_xyz = [4114496.258, 0.0, 4870157.031]  # m
    xyz = llh2ecef(deg2rad(50.0), 0, 10000.0; ellipsoid=Bowring1976)
    @test isapprox(xyz, exp_xyz)

    # http://www.sysense.com/products/ecef_lla_converter/index.html
    x, y, z = (4114496.258, 0.0, 4870157.031)  # m
    exp_llh = [deg2rad(49.996908), 0.000000, 9907.31]
    llh = ecef2llh(x, y, z)
    # Longitude should be exact, latitude max error: 0.0018" according to [1]
    # Altitude 0.17 m
    # The implemented function must be better, but I don't know which one is using
    # this web page. The ellipsoid is assumed to be WGS84
    @test isapprox(llh[:2], exp_llh[:2], atol=deg2rad(0.0018/3600.0))
    @test isapprox(llh[3], exp_llh[3], atol=0.17)

    exp_xyz = [4114291.97, 0.00, 4870449.48]  # m
    xyz = llh2ecef(deg2rad(50.0), 0, 10000.0)
    @test isapprox(xyz, exp_xyz)

    # non-zero longitude
    x, y, z = (3814496.258, 1514496.258, 4870157.031)  # m
    exp_llh = [deg2rad(50.068095), deg2rad(21.654910), 3263.97]
    llh = ecef2llh(x, y, z)
    # Longitude should be exact, latitude max error: 0.0018" according to [1]
    # Altitude 0.17 m
    # The implemented function must be better, but I don't know which one is using
    # this web page. The ellipsoid is assumed to be WGS84
    @test isapprox(llh[:2], exp_llh[:2], atol=deg2rad(0.0018/3600.0))
    @test isapprox(llh[3], exp_llh[3], atol=0.17)

    exp_xyz = [3814496.258, 1514496.258, 4870157.031]  # m
    xyz = llh2ecef(deg2rad(50.068095), deg2rad(21.654910), 3263.97)
    @test isapprox(xyz, exp_xyz)
end


text:
function othersum(si::Vector{Int64},s::Vector{Float64},m::Int64)
  # OTHERSUM Compute the sum of each column of the matrix, without each
  # individual entry.  This corresponds to a matrix where each entry is the
  # sum of each column but with each entry subtracted.
  # rowsum = accumarray(si,s,[m,1])
  rowsum = zeros(Float64,m)
  for i = 1:length(si)
    rowsum[si[i]] += s[i]
  end
  os = rowsum[si] - s
  return os
end