text:
__precompile__()
module SeisDownload

include("utils.jl")
include("downloadfunc.jl")

using .Utils
using .DownloadFunc

using SeisIO, Dates, Printf, JLD2, FileIO, Distributed

#------------------------------------------------------------------#
#For the time being, we need remove_response function from obspy
#This will be replaced by SeisIO modules in the near future.
#Please activate obspy enviroment before launching Julia.
# include("remove_response_obspy.jl")
# using .Remove_response_obspy
#------------------------------------------------------------------#

export seisdownload

"""
    seisdownload(InputDict::Dict)

    Request seismic data and save into jld2.
# Arguments
- `InputDict`    : dictionary which contains request information
"""
function seisdownload(InputDict::Dict)

    Utils.initlogo()

	printparams(InputDict)

	DownloadType    = InputDict["DownloadType"]

	fodir = ""
	sp = splitpath(InputDict["fopath"])
	for i = 1:length(sp)-1
		fodir = joinpath(fodir, sp[i])
	end
	tmppath = joinpath(fodir, "./seisdownload_tmp")
	InputDict["tmppath"] = tmppath
	mkpath(tmppath)

    if DownloadType == "Noise" || DownloadType == "noise"

		#stationlist
		stationlist     = InputDict["stationinfo"]["stationlist"]
		starttime       = InputDict["starttime"]
		endtime         = InputDict["endtime"]
		DL_time_unit    = InputDict["DL_time_unit"]
		DownloadType    = InputDict["DownloadType"]
		fopath          = InputDict["fopath"]

		if mod((endtime - starttime).value,  DL_time_unit) != 0 || (endtime - starttime).value < DL_time_unit
			error("Total download time cannot be devided by Download Time unit; this may cause unexpected result. Abort.")
		end

		# calculate start time list (starttimelist) with each Donwload_time_unit
		starttimelist = Utils.get_starttimelist(starttime, endtime, DL_time_unit)
		# generate DLtimestamplist and ststationlist
		DLtimestamplist = Utils.get_timestamplist(starttimelist)

		InputDict["starttimelist"] = starttimelist
		InputDict["DLtimestamplist"] = DLtimestamplist

		#----Restrict number of processors------#
		#NEVER CHANGE THIS THRESHOLD OTHERWISE IT OVERLOADS THE DATA SERVER
		np = nprocs()
		if np > 100 throw(DomainError(np, "np must be smaller than 100.")) end
		#---------------------------------------#

        # Test download to evaluate use of memory and estimate download time.
		InputDict_test = deepcopy(InputDict) # to avoid overwriting InputDict; unknown bug while deepcopying in testdownload function
		testdownload(InputDict_test, length(starttimelist))

		# Start downloading data
		t_download = @elapsed pmap(x -> seisdownload_NOISE(x, InputDict), 1:length(starttimelist))

		# convert intermediate file to prescibed file format (JLD2, ASDF, ...)
		t_convert = @elapsed convert_tmpfile(InputDict)

		println("---Summary of computational time---")
		println(@sprintf("Total download time:%8.4f[s]", t_download))
		println(@sprintf("Total convert time:%8.4f[s]", t_convert))

		if !InputDict["Istmpfilepreserved"]
			rm(tmppath, recursive=true, force=true)
		end

    elseif  DownloadType == "Earthquake" || DownloadType == "earthquake"

		method		    = InputDict["method"]
		event		    = InputDict["event"]
		reg			    = InputDict["reg"]
		fopath          = InputDict["fopath"]

		#save info into jld2
		jldopen(fopath, "w") do file
			file["info/method"]  = method;
			file["info/event"]   = event;
			file["info/reg"]     = reg
			file["info/fopath"]  = fopath
		end

		#Test download to evaluate use of memory and estimate download time.
		max_num_of_processes_per_parallelcycle = testdownload(InputDict, length(event), MAX_MEM_PER_CPU)

		if max_num_of_processes_per_parallelcycle < 1
			error("Memory allocation is not enought (currently $MAX_MEM_PER_CPU [GB]). Please inclease MAX_MEM_PER_CPU or decrease number of stations")
		end

		if max_num_of_processes_per_parallelcycle >= length(event)

			S = pmap(x -> seisdownload_EARTHQUAKE(x, InputDict), 1:length(event))

			# save data to jld2
			file = jldopen(fopath, "r+")
			unavalilablefile = jldopen(join([fopath[1:end-5], "_unavailablestations.jld2"]), "w+")

			for ii = 1:size(S)[1] #loop at each starttime
				varname = joinpath("event",InputDict["event"][ii]["origin"]["time"][1:end-4])
				file[varname] = S[ii]
			end

			JLD2.close(file)
			JLD2.close(unavalilablefile)

		else

			#parallelization by time
			pitr = 1

			while pitr <= length(event)

				startid1 = pitr
				startid2 = pitr + max_num_of_processes_per_parallelcycle - 1

				if startid1 == length(event)
					#use one
					S = pmap(x -> seisdownload_EARTHQUAKE(x, InputDict), startid1:startid1)

				elseif startid2 <= length(event)
					#use all processors
					S = pmap(x -> seisdownload_EARTHQUAKE(x, InputDict), startid1:startid2)

				else
					#use part of processors
					startid2 = startid1 + mod(length(event), max_num_of_processes_per_parallelcycle) - 1
					S = pmap(x -> seisdownload_EARTHQUAKE(x, InputDict), startid1:startid2)
				end

				# save data to jld2
				file = jldopen(fopath, "r+")
				unavalilablefile = jldopen(join([fopath[1:end-5], "_unavailablestations.jld2"]), "w+")

				for ii = 1:size(S)[1] #loop at each starttime
					varname = joinpath("event",InputDict["event"][startid1+ii-1]["origin"]["time"][1:end-4])
					file[varname] = S[ii]
				end

				JLD2.close(file)
				JLD2.close(unavalilablefile)

				pitr += max_num_of_processes_per_parallelcycle

				#println("pitr: $pitr")
			end
		end

    else
        println("Download type is not known (chose Noise or Earthquake).")
    end

    println("Downloading and Saving data is successfully done.\njob ended at "*string(now()))
    return 0

end

end


text:
const googlefontstart = "<link ref=\"https://fonts.googleapis.com/css?family="
const googlefontend   = "\", rel=\"stylesheet\">"

function googlefontlink(family::String; weight::Vector{Int}=Int[], style::Vector{String} = [""], display::String="")
    str = googlefont(family, weight, style, display)
    return string(googlefontstart, str, googlefontend)
end

function googlefontlink(families::Vector{String};
                        weights::Vector{Vector{Int}}=fill(Int[], length(families)),
                        styles::Vector{Vector{String}}=fill([""], length(families)), display::String="")
    str = googlefont(families, weights, styles, display)
    return string(googlefontstart, str, googlefontend)
end

function googlefont(family::String, weight::Vector{Int}=Int[], style::Vector{String} = [""], display::String="")
    str = join(googlefont_family(family), "+")
    if isempty(weight) && !isempty(style) && !all(isempty.(style))
        weight = [400]
    end
    if isempty(weight)
        weight_style = ""
    else
        weight_style = string(":", googlefont_weight_style(weight, style))
    end
    if !isempty(display)
        display = string("&display=", display)
    end
    return string(str, weight_style, display)
end

function googlefont(families::Vector{String},
                    weights::Vector{Vector{Int}}=fill(Int[], length(families)),
                    styles::Vector{Vector{String}}=fill([""], length(families)), display::String="")
   (length(families) == length(weights) == length(styles)) ||
       throw(ErrorException("inputs ($(length(families)) $(length(weights)) $(length(styles))) must be of the same length"))
   fonts = []
   for (family,weight,style) in zip(families, weights, styles)
        push!(fonts, googlefont(family, weight, style))
   end
   return join(fonts, "|")
end

function googlefont_family(family::String)
    return uppercasefirst.(splitstring(family, " "))
end

function googlefont_weight_style(weight::Vector{Int}, style::Vector{String})
    strs = String[]
    for w in weight
        push!(strs, googlefont_weight_style(w, style))
    end
    return join(strs, ",")
end

function googlefont_weight_style(weight::Int, style::Vector{String})
    nstyles = length(style)
    weightstr = string(weight)
    strs = String[]
    for s in style
        push!(strs, string(weightstr, s))
    end
    return join(strs, ",")
end

ext: jl


text:
using .MathematicalSystems

_throw_not_noisy(sys) = isnoisy(sys) ? sys :  throw(ArgumentError("system is not noisy"))
_throw_not_affine(sys) = isaffine(sys) ? sys :  throw(ArgumentError("system is not affine"))
_throw_not_constrained(sys) = isconstrained(sys) ? sys :  throw(ArgumentError("system is not constrained"))
_throw_not_autonomous(sys) = !iscontrolled(sys) ? sys :  throw(ArgumentError("system is not autonomous"))
_throw_not_controlled(sys) = iscontrolled(sys) ? sys :  throw(ArgumentError("system is not controlled"))

# =============
# robust invariance
# =============
# function minimal_RPI_set(sys::AbstractDiscreteSystem, K;
#                          algorithm=RakovicLazy(), kwargs...)
#     sys |> _throw_not_affine |>
#            _throw_not_noisy
#
#     W = noiseset(sys)
#     S = sys.A - sys.B*K
#     return minimal_RPI_set(S, W)
# end

# function maximal_RPI_set(sys::AbstractDiscreteSystem, K;
#                          algorithm=MaxRPIIterative(), kwargs...)
#     sys |> _throw_not_affine |>
#            _throw_not_noisy |>
#            _throw_not_constrained |>
#            _throw_not_controlled
#
#     X = stateset(sys)
#     W = noiseset(sys)
#     S = sys.A - sys.B*K
#     return maximal_RPI_set(algorithm, S, W, X; kwargs...)
# end

# ===========
# Invariant Sets
# ===========

function state_constraints(sys::AbstractDiscreteSystem, K)
    sys |> _throw_not_constrained |>
           _throw_not_controlled

    X = stateset(sys)
    U = inputset(sys)
    return state_constraints(X, U, K)
end

function preset(sys::AbstractDiscreteSystem, P::LazySet)
    sys |> _throw_not_affine

    A = sys.A
    if !iscontrolled(sys)
        return preset(A, P)
    else
        sys |> _throw_not_constrained
        U = inputset(sys)
        B = sys.B
        return preset(A, B, U, P)
    end
end

function maximum_invariant_set(sys::AbstractDiscreteSystem, max_iter::Integer=100)
    sys |> _throw_not_affine |>
           _throw_not_constrained |>
           _throw_not_autonomous

    A = sys.A
    X = stateset(sys)
    return maximum_invariant_set(A, X, max_iter)
end
function maximum_invariant_set(sys::AbstractDiscreteSystem, K, max_iter::Integer=100)
    sys |> _throw_not_affine |>
           _throw_not_constrained |>
           _throw_not_controlled

    X = stateset(sys)
    U = inputset(sys)
    A = state_matrix(sys)
    B = input_matrix(sys)
    return maximum_invariant_set(A, B, X, U, K, max_iter)
end

function maximum_control_invariant_set(sys::AbstractDiscreteSystem, max_iter=100;
                                       preset=_preset_vrep, kwargs...)
   sys |> _throw_not_affine |>
          _throw_not_constrained |>
          _throw_not_controlled

   X = stateset(sys)
   U = inputset(sys)
   A = state_matrix(sys)
   B = input_matrix(sys)
   return maximum_control_invariant_set(A, B, X, U, max_iter; preset=preset, kwargs...)
end

# =============
# Control Systems
# =============

# function tightened_constraints(sys::AbstractDiscreteSystem, Ω::LazySet, K;
#                                algorithm=RakovicLazy(), kwargs...)
#     sys |> _throw_not_affine |>
#            _throw_not_constrained
#
#     X = stateset(sys)
#     U = inputset(sys)
#     return tightened_constraints(X, U, Ω, K)
# end

function terminal_set(sys::AbstractDiscreteSystem, K; kwargs...)
    sys |> _throw_not_affine |>
          _throw_not_constrained |>
          _throw_not_controlled

    X = stateset(sys)
    U = inputset(sys)
    A = state_matrix(sys)
    B = input_matrix(sys)
    return terminal_set(A, B, X, U, K; kwargs...)
end

function feasible_set(sys::AbstractDiscreteSystem, 𝕏f::LazySet, N::Integer; kwargs...)
    sys |> _throw_not_affine |>
           _throw_not_constrained |>
           _throw_not_controlled

    X = stateset(sys)
    U = inputset(sys)
    A = state_matrix(sys)
    B = input_matrix(sys)
    return feasible_set(A, B, X, U, 𝕏f, N; kwargs...)
end

max_stars_repo_name: mforets/InvariantSets

max_issues_repo_name: mforets/InvariantSets


text:
function VImintqs_q0(p0, ps,q0, k, T, Npi,IRScale,UVScale, mfun, lamfun)
    -hquadrature(
        kprim -> FRGRealTime.dkVImintqs(
            p0,
            ps,
            q0,
            k,
            kprim,
            mfun(kprim),
            T,
            Npi,
            lamfun(kprim),
        ),
        k,
        UVScale,
        rtol = 1e-2,
        atol = 1e-2,
    )[1]
end

using HCubature,FiniteDifferences
plot(
    q0 -> VImintqs_q0(
        0.01,
        2.0,
        q0,
        10.0,
        145.0,
        4.0,
        1.0,
        800.0,
        m2fun,
        lamdafun,
    ),-50.0,50.0,title="V(q0)",xaxis="q0"
)

plot(
    q0 -> central_fdm(5, 1)(x->VImintqs_q0(
        10.0,
        2.0,
        x,
        10.0,
        145.0,
        4.0,
        1.0,
        800.0,
        m2fun,
        lamdafun,
    ),q0),-50.0,50.0,title="V'(q0)",xaxis="q0"
)


text:
#=
function serial_job(args...)
    cmd = ""
    for (i,arg) in enumerate(args)
       cmd *= " "*arg
    end
    run(`scuff-neq $cmd`)
end
=#

"""
    scuff_job(frequency :: Float64, geometry_file :: String ,  output_file :: String )

Serial job for a single frequency given as a float
# Arguments
- `frequency :: Float64` : frequency in radHz. Note that to get the real frequncies, one has to multiply by 3e14.
- `geometry_file :: String` : the *.scuffgeo file
- `output_file :: String` : string to be added to the *.EMTPFT file
"""
function scuff_job(frequency :: Float64, geometry_file :: String ,  output_file :: String )
    freq = string(frequency)
    output = output_file*"freq="*freq*"_radHz"
    run(`scuff-neq --Geometry $geometry_file --Omega $frequency --EMTPFT --FileBase $output`)
end

"""
    scuff_job(frequency :: String, geometry_file :: String ,  output_file :: String )

when using a file containing one or many frequencies
# Arguments
- `frequency :: String` : the name of the file containing the frequencies.
"""
function scuff_job(frequencyfile :: String, geometry_file :: String ,  output_file :: String )
    lines  = readlines(frequencyfile)
    output = output_file*"freq="*lines[1]*"-"*lines[end]
    run(`scuff-neq --Geometry $geometry_file --OmegaFile $frequencyfile --EMTPFT --FileBase $output`)
end

"""
    scuff_parallel(frequencies, geometry_file :: String ,  output_file :: String )

Launch parallel computation for a range of frequencies
# Arguments
- `frequencies` : range object for the frequencies.

# Examples

```jldoctest
julia> frequencies = 1e2:1e2:1e3
100.0:100.0:1000.0
julia> frequencies = linspace(1e2,1e3,10)
100.0:100.0:1000.0
julia> frequencies = logspace(2,3,10)
10-element Array{Float64,1}:
  100.0
  129.155
  166.81
  215.443
  278.256
  359.381
  464.159
  599.484
  774.264
 1000.0
```
"""
function scuff_parallel(frequencies :: AbstractArray, geometry_file :: String , output_file :: String)
    freqv = collect(frequencies)
    scuff_par(frequency) = scuff_job(frequency, geometry_file, output_file)
    pmap(scuff_par, freqv)
end

"""
    scuff_parallel(frequencies :: AbstractString , geometry_file :: AbstractString , output_file :: AbstractString)
when frequencies is a filename, reads the data
"""

function scuff_parallel(frequencies :: AbstractString , geometry_file :: AbstractString , output_file :: AbstractString)
    freqv = readdlm(frequencies)
    scuff_par(frequency) = scuff_job(frequency, geometry_file, output_file)
    pmap(scuff_par, freqv)
end


text:
julia_version = "1.4"
org = "charleskawczynski"
pkg_name = "MyGeneratedPkg.jl"
directory = joinpath(@__DIR__, "..", "..", "GeneratedPackages")

ci = ContinuousIntegration(
    azure=Azure(;org=org, pkg_name=pkg_name, julia_version=julia_version),
    travis=Travis(;julia_version=julia_version),
    appveyor=AppVeyor(;julia_version=julia_version),
)

t = Template(;
    julia_version=julia_version,
    directory=directory,
    pkg_name=pkg_name,
    ci=ci,
)

generate(t)

max_stars_repo_name: charleskawczynski/PkgTemplates

max_issues_repo_name: charleskawczynski/PkgTemplates

max_forks_repo_name: charleskawczynski/PkgTemplates


text:
#: Suite of differentiable functions for TV optimization

export LeastSquaresMisfit, leastsquares_misfit

# Least-squares misfit

"""
Least-squares misfit function associated to the problem
min_x 0.5*||A*x-y||^2
"""
struct LeastSquaresMisfit{T,N1,N2}<:DifferentiableFunction{T,N1}
    A::AbstractLinearOperator{T,N1,N2}
    y::AbstractArray{T,N2}
end

leastsquares_misfit(A::AbstractLinearOperator{T,N1,N2}, y::AbstractArray{T,N2}) where {T,N1,N2} = LeastSquaresMisfit{T,N1,N2}(A, y)

function grad!(f::LeastSquaresMisfit{T,N1,N2}, x::AbstractArray{T,N1}, g::Union{Nothing,AbstractArray{T,N1}}; eval::Bool=true) where {T,N1,N2}
    r = f.A*x-f.y
    eval ? (fval = T(0.5)*norm(r)^2) : (fval = nothing)
    g !== nothing && (g .= adjoint(f.A)*r)
    return fval
end

max_stars_repo_name: grizzuti/FastSolversForWeightedTV

max_issues_repo_name: grizzuti/FastSolversForWeightedTV

max_forks_repo_name: grizzuti/FastSolversForWeightedTV


text:
module DeepDiffs

export deepdiff, added, removed, changed, before, after
export SimpleDiff, VectorDiff, StringDiff, DictDiff

# Helper function for comparing two instances of a type for equality by field
function fieldequal(x::T, y::T) where T
    for f in fieldnames(T)
        isequal(getfield(x, f), getfield(y, f)) || return false
    end
    true
end

hascolor(io::IO) = get(IOContext(io), :color, false)

"""
diff = deepdiff(obj1, obj2)

deepdiff computes the structural difference between two objects and returns
a diff representing "edits" needed to transform obj1 into obj2. This diff
supports the `added`, `removed`, and `modified` functions that return `Set`s of
dictionary keys or array indices.
"""
function deepdiff end

abstract type DeepDiff end

# fallback diff that just stores two values
struct SimpleDiff{T1, T2} <: DeepDiff
    before::T1
    after::T2
end

Base.:(==)(lhs::SimpleDiff, rhs::SimpleDiff) = fieldequal(lhs, rhs)

before(d::SimpleDiff) = d.before
after(d::SimpleDiff) = d.after

deepdiff(x, y) = SimpleDiff(x, y)

include("arrays.jl")
include("dicts.jl")
include("strings.jl")

end # module

max_stars_repo_name: ssfrr/StructureDiffs

max_issues_repo_name: ssfrr/StructureDiffs

max_forks_repo_name: ssfrr/StructureDiffs


text:
# Translating
# https://proplot.readthedocs.io/en/latest/colorbars_legends.html#Axes-colorbars-and-legends
# to Julia

# Import proplot using PyCall
using PyCall
proplot = pyimport("proplot")

proplot.rc.abc=true
f, axs = proplot.subplots(ncols=2, share=0)

# Colorbars
ax = axs[1]
m = ax.heatmap(rand(10, 10), colorbar="t", cmap="dusk")
ax.colorbar(m, loc="r")
ax.colorbar(m, loc="ll", label="colorbar label")
ax.format(title="Axes colorbars", suptitle="Axes colorbars and legends demo")

# Legends
ax = axs[2]
ax.format(title="Axes legends", titlepad="0em")
hs = ax.plot(
    cumsum(rand(10, 5) .- 0.5, dims=1), linewidth=3,
    cycle="ggplot", legend="t",
    labels=["a","b","c","d","e"], legend_kw=Dict("ncols" => 5, "frame" => false)
)
ax.legend(hs, loc="r", ncols=1, frame=false)
ax.legend(hs, loc="ll", label="legend label")
[ax.format(xlabel="xlabel", ylabel="ylabel") for ax in axs]

f.savefig("figs/Axes-colorbars-and-legends-1.svg", transparent=false)

max_stars_repo_name: natgeo-wong/JuliaProPlotExamples

max_issues_repo_name: natgeo-wong/JuliaProPlotExamples

max_forks_repo_name: natgeo-wong/JuliaProPlotExamples


text:
module Maps

export Map, Building
export Circular, Polygon

abstract Map

abstract Building

immutable Circular <: Building

    center::(Float64, Float64)
    radius::Float64
end

immutable Polygon <: Building

    points::Matrix    # Nx3 matrix representing building points
    vertices::Vector  # vertices representing the point order
end

end # module

max_stars_repo_name: etotheipluspi/ParticleFilters

max_issues_repo_name: etotheipluspi/ParticleFilters

max_forks_repo_name: etotheipluspi/ParticleFilters


text:
using Glob
using Weave
cd(Base.source_dir())
weave("images_for_carl.jmd"; doctype="md2html", out_path=:pwd)
rm.(glob("jl_*"), recursive=true)


text:
type ObjectiveFuncEvaluator
  f::Function
  archive::Archive

  ObjectiveFuncEvaluator(objectiveFunc, archive = false) = begin
    archive = archive || TopListArchive()
    new(objectiveFunc, archive)
  end
end

function evaluate(candidate, evaluator::Evaluator)
  fitness = evaluator.f(candidate)
  add!(fitness, candidate, evaluator.archive)
  fitness
end


text:
# Use original distribution when resampling from SampleStore
_get_dist(rv::RandomVariable) = (rv.dist isa SampleStore ? rv.dist.dist : rv.dist)


text:
using FTCTests  # reexport FaulTolerantControl
using Transducers
using Random
using Test

"""
This is for the training of cost function.
"""
function run_multiple_sim_training_cost(N=1;
        manoeuvres=[:hovering],
        methods=[:adaptive],
        t0=0.0, tf=20.0,
        h_threshold=nothing,  # m (nothing: no constraint)
        actual_time_limit=nothing,  # s (nothing: no constraint)
        N_thread=Threads.nthreads(),
        collector=Transducers.tcollect, will_plot=false, seed=2021,
    )
    fault_time = 0.0
    # single fault
    # rotor index ∈ {1, 2, 3, 4, 5, 6}
    # effectiveness ∈ [0, 1]
    _faults = 1:N |> Map(i -> [LoE(fault_time, rand(1:6), rand(1)[1])]) |> collect  # randomly sampled N faults
    # for trajectory parameters (terminal position constraint)
    multicopter = LeeHexacopter()  # dummy
    min_nt, max_nt = FTCTests.distribution_info(:hovering)
    # Control points of Bezier curve; randomly generated constant position tracking
    θs_array = 1:N |> Map(i -> [FTCTests.sample(multicopter, min_nt, max_nt)[1]]) |> collect  # sample(...)[1] corresponds to the position vector
    run_multiple_sim(N,
                     manoeuvres,
                     methods,
                     _faults,
                     θs_array,
                     t0, tf,
                     h_threshold, actual_time_limit, N_thread, collector,
                     will_plot, seed
                    )
end


text:
export sequential_contract

function _contraction_code(fg::FGraph, i::Int, idn::Int)
    DIM = ndims(fg.tensors[i])
    ids = zeros(Int, DIM)
    for ie=1:ne(fg)
        leg = fg.legmap[i,ie]
        if leg != 0
            @inbounds ids[leg] = ie
        end
    end
    for i=1:DIM
        if ids[i] == 0
            idn += 1
            ids[i] = idn
        end
    end
    ids, idn
end

"""get the contraction code from a factor graph"""
function contraction_code(fg::FGraph)
    codes = Any[]
    idn = ne(fg)
    for i in 1:nv(fg)
        code, idn = _contraction_code(fg, i, idn)
        push!(codes, code)
    end
    IC = FactorGraph.leg_analysis(codes)[3]
    codes, IC
end

using TensorOperations: BinaryTreeNode
order2tree(order) = length(order) == 1 ? order[1] : BinaryTreeNode(order2tree(order[1:end-1]), order[end])

@adjoint contraction_code(fg) = contraction_code(fg), adjy->nothing
@adjoint order2tree(order) = order2tree(order), adjy->nothing

"""
    sequential_contract(fg::FGraph, order)

Contract tensors in a factor graph by `order` of tensors.
"""
function sequential_contract(fg::FGraph, order)
    ixs, iy = contraction_code(fg)
    tss = Tuple(fg.tensors)
    tree = order2tree(order)
    treecontract(tree, Tuple.(ixs), tss, Tuple(iy))
end


text:
@testset "execution" begin

############################################################################################

@eval exec_dummy() = return nothing

@testset "cufunction" begin
    @test_throws UndefVarError cufunction(exec_undefined_kernel, ())

    cufunction(dev, exec_dummy, ())
    cufunction(dev, exec_dummy, Tuple{})

    # NOTE: other cases are going to be covered by tests below,
    #       as @cuda internally uses cufunction
end

############################################################################################

@testset "@cuda" begin

@test_throws UndefVarError @cuda (1,1) exec_undefined_kernel()

@testset "reflection" begin
    @cuda (1,1) exec_dummy()

    @grab_output CUDAnative.@code_lowered @cuda (1,1) exec_dummy()
    @grab_output CUDAnative.@code_typed @cuda (1,1) exec_dummy()
    @grab_output CUDAnative.@code_warntype @cuda (1,1) exec_dummy()
    @grab_output CUDAnative.@code_llvm @cuda (1,1) exec_dummy()
    @grab_output CUDAnative.@code_ptx @cuda (1,1) exec_dummy()
    @grab_output CUDAnative.@code_sass @cuda (1,1) exec_dummy()

    @grab_output CUDAnative.@code_lowered exec_dummy()
    @grab_output CUDAnative.@code_typed exec_dummy()
    @grab_output CUDAnative.@code_warntype exec_dummy()
    @grab_output CUDAnative.@code_llvm exec_dummy()
    @grab_output CUDAnative.@code_ptx exec_dummy()
    @grab_output CUDAnative.@code_sass exec_dummy()
end

@testset "reflection with argument conversion" begin
    @eval exec_dummy_array(a, i) = (a[1] = i; return nothing)

    a = CuArray{Float32}(1)

    @grab_output CUDAnative.@code_llvm @cuda (1,1) exec_dummy_array(a, 1)
    @grab_output CUDAnative.@code_ptx @cuda (1,1) exec_dummy_array(a, 1)
    @grab_output CUDAnative.@code_sass @cuda (1,1) exec_dummy_array(a, 1)

    @grab_output CUDAnative.@code_llvm exec_dummy_array(a, 1)
    @grab_output CUDAnative.@code_ptx exec_dummy_array(a, 1)
    @grab_output CUDAnative.@code_sass exec_dummy_array(a, 1)
end

@testset "dimensions" begin
    @cuda (1,1) exec_dummy()
    @test_throws ArgumentError @cuda (0,0) exec_dummy()
end

@testset "shared memory" begin
    @cuda (1,1,1) exec_dummy()
end

@testset "streams" begin
    s = CuStream()
    @cuda (1,1,1,s) exec_dummy()
end

@testset "external kernels" begin
    @eval module KernelModule
        export exec_external_dummy
        exec_external_dummy() = return nothing
    end
    @cuda (1,1) KernelModule.exec_external_dummy()
    @eval begin
        using KernelModule
        @cuda (1,1) exec_external_dummy()
    end

    @eval module WrapperModule
        using CUDAnative
        @eval exec_dummy() = return nothing
        wrapper() = @cuda (1,1) exec_dummy()
    end
    WrapperModule.wrapper()
end

@testset "return values" begin
    @eval exec_return_int_inner() = return 1
    @test_throws ErrorException @cuda (1,1) exec_return_int_inner()

    @eval function exec_return_int_outer()
        exec_return_int_inner()
        return nothing
    end
    @cuda (1,1) exec_return_int_outer()
end

@testset "calling device function" begin
    @eval @noinline exec_devfun_child(i) = return i+1
    @eval exec_devfun_parent() = (exec_devfun_child(1); return nothing)

    @cuda (1,1) exec_devfun_parent()
end

end

############################################################################################

@testset "argument passing" begin

dims = (16, 16)
len = prod(dims)

@testset "manually allocated" begin
    @eval function exec_pass_ptr(input, output)
        i = (blockIdx().x-1) * blockDim().x + threadIdx().x

        val = unsafe_load(input, i)
        unsafe_store!(output, val, i)

        return nothing
    end

    input = round.(rand(Float32, dims) * 100)

    input_dev = CuArray(input)
    output_dev = similar(input_dev)

    @cuda (1,len) exec_pass_ptr(input_dev.devptr, output_dev.devptr)
    output = Array(output_dev)
    @test input ≈ output
end

@testset "scalar through single-value array" begin
    @eval function exec_pass_scalar(a, x)
        i = (blockIdx().x-1) * blockDim().x + threadIdx().x
        max = gridDim().x * blockDim().x
        if i == max
            val = unsafe_load(a, i)
            unsafe_store!(x, val)
        end

        return nothing
    end

    arr = round.(rand(Float32, dims) * 100)
    val = Float32[0]

    arr_dev = CuArray(arr)
    val_dev = CuArray(val)

    @cuda (1,len) exec_pass_scalar(arr_dev.devptr, val_dev.devptr)
    @test arr[dims...] ≈ Array(val_dev)[1]
end

@testset "scalar through single-value array, using device function" begin
    @eval @noinline function exec_pass_scalar_devfun(a, x)
        i = (blockIdx().x-1) * blockDim().x + threadIdx().x
        max = gridDim().x * blockDim().x
        if i == max
            val = exec_pass_scalar_devfun_child(a, i)
            unsafe_store!(x, val)
        end

        return nothing
    end
    @eval function exec_pass_scalar_devfun_child(a, i)
        return unsafe_load(a, i)
    end

    arr = round.(rand(Float32, dims) * 100)
    val = Float32[0]

    arr_dev = CuArray(arr)
    val_dev = CuArray(val)

    @cuda (1,len) exec_pass_scalar_devfun(arr_dev.devptr, val_dev.devptr)
    @test arr[dims...] ≈ Array(val_dev)[1]
end

@testset "ghost function parameters" begin
    # bug: ghost type function parameters are elided by the compiler

    len = 60
    a = rand(Float32, len)
    b = rand(Float32, len)

    d_a = CuArray(a)
    d_b = CuArray(b)
    d_c = similar(d_a)

    @eval struct ExecGhost end

    @eval function exec_pass_ghost(ghost, a, b, c)
        i = (blockIdx().x-1) * blockDim().x + threadIdx().x
        unsafe_store!(c, unsafe_load(a,i)+unsafe_load(b,i), i)

        return nothing
    end
    @cuda (1,len) exec_pass_ghost(ExecGhost(), d_a.devptr, d_b.devptr, d_c.devptr)

    c = Array(d_c)
    @test a+b == c
end

@testset "tuples" begin
    # issue #7: tuples not passed by pointer

    @eval function exec_pass_tuples(keeps, out)
        if keeps[1]
            unsafe_store!(out, 1)
        else
            unsafe_store!(out, 2)
        end
        nothing
    end

    keeps = (true,)
    d_out = CuArray{Int}(1)

    @cuda (1,1) exec_pass_tuples(keeps, d_out.devptr)
    @test Array(d_out) == [1]
end

@testset "immutables" begin
    # issue #15: immutables not passed by pointer

    @eval function exec_pass_immutables(A, b)
        unsafe_store!(A, imag(b))
        nothing
    end

    A = CuArray(zeros(Float32, (1,)))
    x = Complex64(2,2)

    @cuda (1, 1) exec_pass_immutables(A.devptr, x)
    @test Array(A) == Float32[imag(x)]
end

end

############################################################################################

@testset "profile" begin

CUDAnative.@profile @cuda (1,1) exec_dummy()

end

############################################################################################

end


text:
export @pcyan,@pblue,@pred,@pgreen

for (fn,clr) in ((:pcyan,   :cyan),
                 (:pblue,   :blue),
                 (:pred,    :red),
                 (:pgreen,  :green))

    @eval macro $fn(msg...)
        Base.print_with_color($(Expr(:quote, clr)), msg... )
    end

end


text:
using Test
using MPI: mpiexec

# Do the MPI tests first so we do not have mpi running inside MPI
mpi_tests = ("parallelstencil_mpi.jl", "petsc_mpi.jl")

# Perform the MPI tests
for file in mpi_tests
    @testset "MPI test: $file" begin
        @test mpiexec() do mpi_cmd
            cmd =
                `$mpi_cmd -n 4 $(Base.julia_cmd()) --startup-file=no --project $file`
            success(pipeline(cmd, stderr = stderr))
        end
    end
end

# Perform non-MPI tests
include("./parallelstencil_singlecore.jl")
include("./petsc_singlecore.jl")


text:
module BLB

# Test that the BLB functions work properly
using MixedModelsBLB, Random, Distributions, DataFrames, StatsModels, Ipopt, Test, StableRNGs, WiSER

# Simulate a dataset
# βtrue  = ones(3); σ²true = 1; Σtrue  = [1.0 0; 0 1.0]
rng = StableRNG(123)
# Random.seed!(rng, 123)
N = 1000; reps = 10
x1 = rand(rng, Normal(0, 1), reps * N)
x2 = rand(rng, Normal(0, 1), reps * N)
x3 = repeat(["M", "F"], inner = reps * N >> 1)
rand_slope = zeros(reps * N)
for j in 1:N
    rand_slope[(reps * (j-1) + 1) : reps * j] = x1[(reps * (j-1) + 1) : reps * j] .* rand(rng, Normal(0, 1), 1)
end
y = 1 .+ x1 + x2 + # fixed effects
    repeat(rand(rng, Normal(0, 1), N), inner = reps) + # random intercept, standard normal
    rand_slope +
    rand(rng, Normal(0, 1), reps * N) # error, standard normal
id = repeat(1:N, inner = reps)
dat = DataFrame(y = y, x1 = x1, x2 = x2, x3 = x3, id = id)

# Parametric BLB with Maximum Likelihood
result0 = blb_full_data(
        StableRNG(123),
        dat;
        feformula = @formula(y ~ 1 + x1 + x2 + x3),
        reformula = @formula(y ~ 1 + x1),
        id_name = "id",
        cat_names = ["x3"], #Array{String,1}(),
        subset_size = 200,
        n_subsets = 2,
        n_boots = 100,
        method = :ML,
        solver = Ipopt.IpoptSolver(print_level=0, max_iter=100, mehrotra_algorithm = "yes", warm_start_init_point = "yes", warm_start_bound_push = 1e-9),
        verbose = false,
        nonparametric_boot = false
)
mean_β = fixef(result0)
mean_Σ, mean_σ² = vc(result0)
ci_β, ci_Σ, ci_σ² = MixedModelsBLB.confint(result0)

@testset "Parametric BLB Result for Method = :ML" begin
    # Since RNG is not stable between Julia versions, the tests below can fail on any version.
    # One solution is to use StableRNG https://github.com/JuliaRandom/StableRNGs.jl
    # Otherwise we just don't test the actual numbers. Instead, just test n_boots, n_subsets etc.

    @test mean_β[1] ≈ 0.9523749104428877 atol = 1e-3
    @test mean_Σ[1, 1] ≈ 1.094385387720442 atol = 1e-3
    @test mean_σ² ≈ 1.0150958802746968 atol = 1e-3

    @test ci_β[1, 1] ≈ 0.7406733280497191 atol = 1e-3
    @test ci_Σ[1, 1] ≈ 0.8490824216956047 atol = 1e-3
    @test ci_σ²[1, 1] ≈ 0.9453143174741494 atol = 1e-3
end

# Nonparametric BLB with Maximum Likelihood
result1 = blb_full_data(
        StableRNG(123),
        dat;
        feformula = @formula(y ~ 1 + x1 + x2 + x3),
        reformula = @formula(y ~ 1 + x1),
        id_name = "id",
        cat_names = ["x3"], #Array{String,1}(),
        subset_size = 200,
        n_subsets = 2,
        n_boots = 100,
        method = :ML,
        solver = Ipopt.IpoptSolver(print_level=0, max_iter=100, mehrotra_algorithm = "yes", warm_start_init_point = "yes", warm_start_bound_push = 1e-9),
        verbose = false,
        nonparametric_boot = true
)
mean_β = fixef(result1)
mean_Σ, mean_σ² = vc(result1)
ci_β, ci_Σ, ci_σ² = MixedModelsBLB.confint(result1)

@testset "Nonparametric BLB Result for Method = :ML" begin
    # Since RNG is not stable between Julia versions, the tests below can fail on any version.
    # One solution is to use StableRNG https://github.com/JuliaRandom/StableRNGs.jl
    # Otherwise we just don't test the actual numbers. Instead, just test n_boots, n_subsets etc.

    @test mean_β[1] ≈ 0.9681988820743825 atol = 1e-3
    @test mean_Σ[1, 1] ≈ 1.0331613190126188 atol = 1e-3
    @test mean_σ² ≈ 1.011932405830439 atol = 1e-3

    @test ci_β[1, 1] ≈ 0.8811621102951823 atol = 1e-3
    @test ci_Σ[1, 1] ≈ 0.9379266586770029 atol = 1e-3
    @test ci_σ²[1, 1] ≈ 0.9776523913354277 atol = 1e-3
end

# Nonparametric BLB with WiSER
result2 = blb_full_data(
        StableRNG(123),
        dat;
        feformula    = @formula(y ~ 1 + x1 + x2 + x3),
        reformula    = @formula(y ~ 1 + x1),
        wsvarformula = @formula(y ~ 1),
        id_name      = "id",
        cat_names    = ["x3"], #Array{String,1}(),
        subset_size  = 200,
        n_subsets    = 2,
        method       = :WiSER,
        n_boots      = 100,
        solver       = Ipopt.IpoptSolver(print_level=0),
        verbose      = false,
        nonparametric_boot = true
)
mean_β = fixef(result2)
mean_Σ, mean_σ² = vc(result2)
ci_β, ci_Σ, ci_σ² = MixedModelsBLB.confint(result2)

@testset "Nonparametric BLB Result for Method = :WiSER" begin
    @test mean_β[1] ≈ 0.9682015976735611 atol = 1e-3
    @test mean_Σ[1, 1] ≈ 1.0331762861674212 atol = 1e-3
    @test mean_σ² ≈ 1.0119282472432958 atol = 1e-3

    @test ci_β[1, 1] ≈ 0.8811874356556482 atol = 1e-3
    @test ci_Σ[1, 1] ≈ 0.9379372642984771 atol = 1e-3
    @test ci_σ²[1, 1] ≈ 0.9776457919670338 atol = 1e-3
end

end


text:
"""
gives the ordering for nedelec2 (divergence)
"""
function relorientation_gen(face::SArray{Tuple{3},T,1,3}, tet::SArray{Tuple{4},T,1,4}) where {T}
    # v = setdiff(tet,face)
    # length(v) == 1 || return 0

    # a = something(findfirst(isequal(v[1]), tet),0)

    a = findfirst(x -> !(x in face), tet)
    a == nothing && return 0

    @assert 1 <= a <= 4
    q = findnext(x ->!(x in face), tet, a+1)
    q == nothing || return 0

    v = tet[a]

    Face = deleteat(tet,a)
    # @show a tet face Face indexin(face,Face)
    p = parity(Vector{Int}(indexin(face,Face)))
    return a * (-1)^a * (-1)^p

    # face = [v[1],face[1],face[2],face[3]]
    # w = sortperm(face)
    # b = parity(w)
    # return -a*(-1)^b
end

const RelOp34 = cat(
    [0 0 0 0; 0 0 4 -3; 0 -4 0 2; 0 3 -2 0],
    [0 0 -4 3; 0 0 0 0; 4 0 0 -1; -3 0 1 0],
    [0 4 0 -2; -4 0 0 1; 0 0 0 0; 2 -1 0 0],
    [0 -3 2 0; 3 0 -1 0; -2 1 0 0; 0 0 0 0],
    dims=3)

function relorientation(face::SArray{Tuple{3},T,1,3}, tet::SArray{Tuple{4},T,1,4}) where {T}

    i = findfirst(==(face[1]),tet); i == nothing && return 0
    j = findfirst(==(face[2]),tet); j == nothing && return 0
    k = findfirst(==(face[3]),tet); k == nothing && return 0
    return RelOp34[i,j,k]

end

"""
gives the ordering for nedelec1 (curl)
if a tetrahedron is given by simplex(a,b,c,d), the order is:
    a-b
    a-c
    a-d
    b-c
    b-d
    c-d
"""

const _relorient24 = @SMatrix[
    0 1 2 3
    -1 0 4 5
    -2 -4 0 6
    -3 -5 -6 0]

# function relorientation(edge::SArray{Tuple{2},Int64,1,2}, tet::SArray{Tuple{4},Int64,1,4})
@inline function relorientation(edge::SVector{2,Int}, tet::SVector{4,Int})

    # relorient24 = [
    # 0 1 2 3
    # -1 0 4 5
    # -2 -4 0 6
    # -3 -5 -6 0]

    # # println("relorientation 2-4")
    # e1 = edge[1]
    # a = 0
    # for v in tet
    #     a = a + 1
    #     if e1 == v
    #         break
    #     end
    # end
    # a == 5 && return 0

    # e2 = edge[2]
    # b = 0
    # for v in tet
    #     b = b + 1
    #     if e2 == v
    #         break
    #     end
    # end
    # b == 5 && return 0
    # return relorient24[a,b]

    a = findfirst(==(edge[1]), tet)
    a == nothing && return 0

    b = findfirst(==(edge[2]), tet)
    b == nothing && return 0

    return _relorient24[a,b]

end

# function relorientation(edge::SArray{Tuple{2},Int64,1,2}, tet::SArray{Tuple{4},Int64,1,4})
#     # v = setdiff(tet, edge)
#     # length(v) == 2 || return 0

#     count(x -> !(x in edge), tet) != 2 && return 0

#     w1 = findfirst(isequal(edge[1]), tet)
#     w2 = findfirst(isequal(edge[2]), tet)

#     @assert w1 != nothing
#     @assert w2 != nothing

#     t = tetschoice(w1,w2)
#     s = sign(w2-w1)

#     return s*t
# end

const RelOp23 = [0 3 -2; -3 0 1; 2 -1 0]
function relorientation(face::SVector{2,Int}, cell::SVector{3,Int})
    i = findfirst(==(face[1]), cell); i == nothing && return 0
    j = findfirst(==(face[2]), cell); j == nothing && return 0
    return RelOp23[i,j]
end

function tetschoice(w1::Int64,w2::Int64)
    a = min(w1,w2)
    b = max(w1,w2)
    if a == 1
        return b-1
    elseif a == 2
        return b+1
    elseif a == 3
        return b+2
    end
end

function relorientation(node::SArray{Tuple{1},Int64,1,1}, face::SArray{Tuple{N},Int64,1,N} where {N})
    idx = something(findfirst(isequal(node[1]), face), 0)
    sgn = (-1)^(idx-1)
    return sgn * idx
end

function relorientation(face, simplex)

    v = setdiff(simplex, face)
    length(v) == 1 || return 0

    # find the position of the missing vertex
    v = v[1]
    i = something(findfirst(isequal(v), simplex),0)
    s = (-1)^(i-1)

    # remove that vertex from the simplex
    face2 = Array{Int}(undef,length(simplex)-1)
    for j in 1 : i-1
        face2[j] = simplex[j]
    end
    for j in i : length(simplex)-1
        face2[j] = simplex[j+1]
    end

    # get the permutation that maps face to face2
    p = [something(findfirst(isequal(v),face2),0) for v in face]

    return s * levicivita(p) * i
end

relorientation(p::SVector{4}, q::SVector{3}) = relorientation(q,p)
relorientation(p::SVector{4}, q::SVector{2}) = relorientation(q,p)
relorientation(p::SVector{4}, q::SVector{1}) = relorientation(q,p)
relorientation(p::SVector{3}, q::SVector{2}) = relorientation(q,p)
relorientation(p::SVector{3}, q::SVector{1}) = relorientation(q,p)
relorientation(p::SVector{2}, q::SVector{1}) = relorientation(q,p)

max_stars_repo_name: krcools/CompScienceMeshes

max_issues_repo_name: krcools/CompScienceMeshes

max_forks_repo_name: krcools/CompScienceMeshes


text:
dipparams = DipoleParameters(theta_rad = deg2rad(15),
    phi_rad = deg2rad(30),
    speed_m_s = 1e7,
    t_k = 10.0,
)

@test dipole_temperature([0, 0, 1], params = dipparams) ≈ 0.32717489089672647
@test dipole_temperature([0, 0, 1], 1e12, params = dipparams) ≈ 0.3475227551239395
@test dipole_temperature([0, 0, 1], 1e13, params = dipparams) ≈ 0.5713065804441115
@test dipole_temperature([0, 1, 0], params = dipparams) ≈ 0.03776458768801616
@test dipole_temperature([0, 1, 0], 1e12, params = dipparams) ≈ 0.0436209260747376
@test dipole_temperature([0, 1, 0], 1e13, params = dipparams) ≈ 0.047637665018382125


text:
using Plots, Statistics, MLDatasets, LinearAlgebra, Clustering, Printf, Dates, Random
path1=realpath(dirname(@__FILE__)*"/..")
include(string(path1,"/src/ANNBN.jl"))

noise_perc=0.1;vars=5;
ff(x)=-x[1]+(1/2)x[2].^2-(1/3)x[3].^3+(1/4)x[4].^4-(1/5)x[5].^5
i_train=1000; xx_train=0.8*rand(i_train,vars).+0.1
yy_train=zeros(i_train); for i=1:i_train yy_train[i]=ff(xx_train[i,:]) end
yy_train.+=noise_perc.*(1/2 .-rand(i_train))
mi=minimum(yy_train);yy_train.-=mi;ma=maximum(yy_train);yy_train./=ma;yy_train.*=0.8;yy_train.+=0.1
i_test=1000;xx_test=0.8*rand(i_test,vars).+0.1
yy_test=zeros(i_test); for i=1:i_test yy_test[i]=ff(xx_test[i,:]) end
# yy_test.+=noise_perc.*rand(i_test)
yy_test.-=mi;yy_test./=ma;yy_test.*=0.8;yy_test.+=0.1

using DecisionTree, Statistics
num_trees=100;
regr_3 = RandomForestRegressor(n_trees=num_trees);
tmp_xx_train=copy(xx_train);tmp_yy_train=copy(yy_train);
tmp_xx_test=copy(xx_test);tmp_yy_test=copy(yy_test);
DecisionTree.fit!(regr_3, tmp_xx_train, tmp_yy_train[:,1]);
predl1 = DecisionTree.predict(regr_3, tmp_xx_train);
predl1_test = DecisionTree.predict(regr_3, tmp_xx_test);
mean(abs.(predl1-yy_train))
mae=mean(abs.(predl1_test-yy_test))
scatter(yy_test,predl1_test,label="Random Forests",markersize=4,legend=:topleft)

using PyCall,ScikitLearn
@sk_import ensemble: AdaBoostRegressor
sci_model = fit!(AdaBoostRegressor(n_estimators=100),xx_train,yy_train)
predl1=predict(sci_model,xx_train)
predl1_test=predict(sci_model,xx_test)
maetr=mean(abs.(yy_train-predl1))
maete=mean(abs.(yy_test-predl1_test))
scatter!(yy_test,predl1_test,label="AdaBoost",markersize=4,legend=:topleft)

using XGBoost
num_round = 10
bst = xgboost(copy(xx_train), num_round, label = copy(yy_train), eta = 1, max_depth = 7)
predl1 = XGBoost.predict(bst, copy(xx_train))
predl1_test = XGBoost.predict(bst, copy(xx_test))
mean(abs.(predl1-yy_train))
mae=mean(abs.(predl1_test-yy_test))
scatter!(yy_test,predl1_test,label="XGBoost",markersize=4,legend=:topleft)

##### ANNBN RBF
neurons=Int64(floor(i_train/(vars+1))) # this automatically generates number of neurons== number of clusters (§3.1)
inds_all,n_per_part=ANNBN.___clustering(neurons,xx_train,200)
cc1=1.0
a_all,a_layer1,layer1=ANNBN.train_layer_1_rbf(neurons,vars,i_train,n_per_part,inds_all,xx_train,yy_train,cc1)
predl1,layer1_train=ANNBN.predict_new_rbf(a_all,a_layer1,xx_train,i_train,neurons,n_per_part,inds_all,xx_train,vars,cc1)
predl1_test,layer1_test=ANNBN.predict_new_rbf(a_all,a_layer1,xx_test,i_test,neurons,n_per_part,inds_all,xx_train,vars,cc1)
maetr=mean(abs.(yy_train-predl1))
maete=mean(abs.(yy_test-predl1_test))
scatter!(yy_test,predl1_test,label="ANNBN-RBF",markersize=4,legend=:topleft)
# savefig("sort-err.pdf")

# RBF FOLDS
nof_folds=10 # you may change these
i_fold=Int64(floor(0.99i_train))
maetrs,a_all_all,a_layer1_all,n_per_part_all,inds_all_all,xx_fold_all,neurons_all=
    ANNBN.fit_nfolds!(xx_train,yy_train,nof_folds,vars,i_train,i_fold,cc1,neurons);
predl1_all=ANNBN.predict_nfolds(a_all_all,a_layer1_all,n_per_part_all,inds_all_all,xx_train,i_train,vars,cc1,xx_fold_all,nof_folds,maetrs,neurons_all)
predl1_train=sum(predl1_all)./nof_folds# ./sum(1 ./maetrs)
maetr=mean(abs.(yy_train-predl1_train))
predl1_all=ANNBN.predict_nfolds(a_all_all,a_layer1_all,n_per_part_all,inds_all_all,xx_test,i_test,vars,cc1,xx_fold_all,nof_folds,maetrs,neurons_all)
predl1_test=sum(predl1_all)./nof_folds#./sum(1 ./maetrs)
maete=mean(abs.(yy_test-predl1_test))
scatter!(yy_test,predl1_test,label="ANNBN-RBF-FOLDS",markersize=4,legend=:topleft)
# check folds history
maes_all=Vector{Float64}()
for i=1:nof_folds
    predl1_i=sum(predl1_all[1:i])./i#./sum(1 ./maetrs[1:i])
    maete_i=mean(abs.(yy_test-predl1_i))
    push!(maes_all,maete_i)
end
# scatter(maes_all)

# Sigmoid
nof_folds=100 # you may change these
i_fold=Int64(floor(0.95i_train))
maetrs,a_all_all,a_layer1_all,n_per_part_all,inds_all_all,xx_fold_all,neurons_all=
    ANNBN.fit_nfolds_sigmoid(xx_train,yy_train,nof_folds,vars,i_train,i_fold,neurons);
predl1_all=ANNBN.predict_nfolds_sigmoid(a_all_all,a_layer1_all,n_per_part_all,inds_all_all,xx_train,i_train,vars,xx_fold_all,nof_folds,maetrs,neurons_all)
predl1_train=sum(predl1_all)./nof_folds# ./sum(1 ./maetrs)
maetr=mean(abs.(yy_train-predl1_train))
predl1_all=ANNBN.predict_nfolds_sigmoid(a_all_all,a_layer1_all,n_per_part_all,inds_all_all,xx_test,i_test,vars,xx_fold_all,nof_folds,maetrs,neurons_all)
predl1_test=sum(predl1_all)./nof_folds#./sum(1 ./maetrs)
maete=mean(abs.(yy_test-predl1_test))
scatter!(yy_test,predl1_test,label="ANNBN-SIGMOID-FOLDS",markersize=4,legend=:topleft)
# check folds history
maes_all=Vector{Float64}()
for i=1:nof_folds
    predl1_i=sum(predl1_all[1:i])./i#./sum(1 ./maetrs[1:i])
    maete_i=mean(abs.(yy_test-predl1_i))
    push!(maes_all,maete_i)
end
# scatter(maes_all)


text:
function __init__ptblm()
    DEPNAME = "PTBLM"
    TRAINFILE = "ptb.train.txt"
    TESTFILE  = "ptb.test.txt"

    register(DataDep(
        DEPNAME,
        """
        Dataset: Penn Treebank sentences for language modeling
        Website: https://github.com/tomsercu/lstm

        The files are available for download at the github
        repository linked above. Note that using the data
        responsibly and respecting copyright remains your
        responsibility.
        """,
        "https://raw.githubusercontent.com/tomsercu/lstm/master/data/" .* [TRAINFILE, TESTFILE],
        "218f4e6c7288bb5efeb03cc4cb8ae9c04ecd8462ebfba8e13e3549fab69dc25f",
    ))
end

"""
    PTBLM(; split=:train, dir=nothing)
    PTBLM(split; [dir])

The PTBLM dataset consists of Penn Treebank sentences
for language modeling, available from https://github.com/tomsercu/lstm.
The unknown words are replaced with <unk> so that the
total vocaburary size becomes 10000.
"""
struct PTBLM <: SupervisedDataset
    metadata::Dict{String, Any}
    split::Symbol
    features::Vector{Vector{String}}
    targets::Vector{Vector{String}}
end

PTBLM(; split=:train, dir=nothing) = PTBLM(split; dir)

function PTBLM(split::Symbol; dir=nothing)
    DEPNAME = "PTBLM"
    @assert split ∈ [:train, :test]
    FILE = split == :train ? "ptb.train.txt" : "ptb.test.txt"
    path = datafile(DEPNAME, FILE, dir)

    lines = open(readlines, path)
    @assert all(x -> x isa String, lines)
    features = map(l -> Vector{String}(Base.split(chomp(l))), lines)

    targets = map(features) do x
        y = copy(x)
        popfirst!(y)
        push!(y, "<eos>")
    end

    metadata = Dict{String,Any}("n_observations" => length(features))
    return PTBLM(metadata, split, features, targets)
end

# DEPRECATED INTERFACE, REMOVE IN v0.7 (or 0.6.x)
function Base.getproperty(::Type{PTBLM}, s::Symbol)
    if s === :traindata
        @warn "PTBLM.traindata() is deprecated, use `PTBLM(split=:train)[]` instead." maxlog=2
        function traindata(; dir=nothing)
            PTBLM(; split=:train, dir)[]
        end
        return traindata
    elseif s === :testdata
        @warn "PTBLM.testdata() is deprecated, use `PTBLM(split=:test)[]` instead."  maxlog=2
        function testdata(; dir=nothing)
            PTBLM(; split=:test, dir)[]
        end
        return testdata
    else
        return getfield(PTBLM, s)
    end
end


text:
import Base: *, sort
using Random, Statistics, DataFrames, CSV
using BioSequences, XLSX,ArgParse

# Overload the concatenation operator to combine a sequence and
# a single base, i.e. dna"AGCGTGC" * DNA_T
function *(seq::LongDNASeq, base::DNA)
    seq * LongDNASeq([base])
end

# Overload sort to order DNA codes by degeneracy:
#   A,G,C,T,M,R,W,S,Y,K,V,H,D,B,N
function sort(seq::LongSequence{DNAAlphabet{4}}; rev::Bool=false)
    sort(convert(Vector{DNA}, seq), by=degeneracy, rev=rev)
end

"""
    isvalid(seq, sites)

Check if any of the sites appear in the sequence.
"""
function isvalid(seq, sites)
    for site in sites
        if occursin(site, seq)
            return false
        end
    end
    return true
end

"""
    find_valid_randomer(prefix, bases, sites)

Randomly select bases to add the prefix, ensuring no site
appears in the randomer.
"""
function find_valid_randomer(prefix, bases, sites)
    n = length(bases)
    pre_n = length(prefix)
    randomer = prefix * (dna"-" ^ n)
    for i = 1:n
        candidates = shuffle(bases[i])
        for cand in candidates
            randomer[pre_n+i] = cand
            if isvalid(randomer[1:pre_n+i], sites)
                break
            else
                randomer[pre_n+i] = DNA_Gap
            end
        end
        if randomer[pre_n+i] == DNA_Gap
            break
        end
    end
    return randomer
end

"""
    degeneracy(base::DNA[, uselog2=true])

Calculate the degeneracy of a single base. The degeneracy
is the number of bases in the code, so
    degeneracy(DNA_N) == 4
and
    degeneracy(DNA_A) == 1

If `uselog2`, the degeneracy is log2 transformed. Because degeneracy(DNA_Gap)=0, we approximate
the log2 transformation as 2^-100.
"""
function degeneracy(base::DNA; uselog2=true)
    deg = 0.0
    if base == DNA_N
        deg = 4.0
    elseif base == DNA_A || base == DNA_C || base == DNA_G || base == DNA_T
        deg = 1.0
    elseif base == DNA_B || base == DNA_D || base == DNA_H || base == DNA_V
        deg = 3.0
    elseif (base == DNA_K || base == DNA_Y || base == DNA_S || base == DNA_W
            || base == DNA_R || base == DNA_M)
        deg = 2.0
    else # DNA_Gap
        deg = 0.0
    end

    if uselog2
        if base == DNA_Gap
            deg = 2^-100
        end
        return log2(deg)
    else
        return deg
    end
end

"""
    degeneracy(seq::LongSequence{DNAAlphabet{4}}; uselog2=true)

Calculate the degeneracy of a DNA sequence.
"""
function degeneracy(seq::LongSequence{DNAAlphabet{4}}; uselog2=true)
    if uselog2
        return sum(degeneracy.(seq, uselog2=true))
    else
        return prod(degeneracy.(seq ;uselog2=false))
    end
end

# Run Monte Carlo simulations using a random base policy over a fixed horizon.
function simulate_random(prefix, bases, sites; nsims=1000, horizon=length(bases))
    degs = zeros(nsims)
    true_horizon = min(horizon, length(bases)) # don't go beyond the randomer length
    for i = 1:nsims
        degs[i] = degeneracy(find_valid_randomer(prefix, bases[1:true_horizon], sites))
    end
    return mean(degs)
end

# Run a 1-step greedy lookahead policy.
function simulate_greedy(prefix, bases, sites; horizon=length(bases))
    true_horizon = min(horizon, length(bases)) # don't go beyond the randomer length
    for h = 1:true_horizon
        # Iterate through the bases by decending degeneracy, so the greedy approach
        # is to stop as soon as we've found one.
        found = false # define out here for scoping rules
        for base in sort(bases[h], rev=true)
            found = false
            if isvalid(prefix * base, sites)
                prefix *= base
                found = true
                break
            end
        end
        if ~found
            # None of the bases were valid, so stop and leave the DNA_Gap
            break
        end
    end
    return degeneracy(prefix)
end

"""
    cutfree_rollout(bases::Array{LongSequence{DNAAlphabet{4}},1}, sites::Array{LongSequence{DNAAlphabet{4}},1}; simulate=simulate_random, kwargs...)

Run the rollout algorithm to solve the CutFree MDP.

# Arguments
- `bases`: An array allowed bases at each postion, usually all 15 degenerate bases. Option is given because some companies restrict which degernate bases are allowed.
- `sites`: An array of restriciton enzyme recognition sequences to be blocked in the random barcode.
- `simulate`: The choice of policy for rollout simulations. simulate_random will use a random rollout policy. simulate_greedy will use a greedy 1 step lookahead policy.
#Optional Keyword Arugments
- `nsims`: The number of rollout simulations per action
"""
function cutfree_rollout(bases::Array{LongSequence{DNAAlphabet{4}},1}, sites::Array{LongSequence{DNAAlphabet{4}},1}; simulate=simulate_random, kwargs...)
    n = length(bases)
    randomer = dna"-" ^ n

    # If the restriction site is non-palindromic, we also need to block
    # the reverse complement.
    for i = 1:length(sites)
        if ~ispalindromic(sites[i])
            push!(sites, reverse_complement(sites[i]))
        end
    end

    # Find the longest restriction site and use this as the horizon.
    max_len = map(length, sites) |> maximum
    horizon = max_len

    for i = 1:n
        best_base = DNA_Gap
        best_deg = -1.0
        start = max(i - max_len, 1) # only need to look back max_len when checking the oligos
        for base in bases[i]
            if ~isvalid(randomer[start:i-1] * base, sites)
                # adding this bases creates a restriction site; skip it
                continue
            end
            deg = simulate(randomer[start:i-1] * base, bases[i+1:end], sites; horizon=horizon, kwargs...)

            # Update use this base if:
            #   1. the mean randomer degeneracy is higher than the previous best, OR
            #   2. the randomer degeneracy is tied with the previous best, but the new base has
            #      higher individual degeneracy
            if deg > best_deg || (deg == best_deg && degeneracy(base) > degeneracy(best_base))
                best_deg = deg
                best_base = base
            end
        end

        if best_base == DNA_Gap
            # The sequence has terminated; there are no bases that can be added without making
            # a restriction site.
            break
        else
            randomer[i] = best_base
        end
    end

    return randomer
end

# Test code
#all_bases = dna"AGCTMRWSYKVHDBN"
#all_ns = [all_bases for _ in 1:40] # make a length 10 randomer
#sites = [dna"GGTCTC"] # BsaI

#randomer = cutfree_rollout(all_ns, sites, simulate=simulate_random, nsims=1000)
#randomer = cutfree_rollout(all_ns, sites, simulate=simulate_greedy)

#deg = degeneracy(randomer)
#natdeg = deg / log2(exp(1))
#println("Final randomer: $randomer")
#println("Final degeneracy: $deg ($natdeg)")

# Run the benchmark CutFree set
function runtimes(data)
    data[:random_objval] = 0.0
    data[:random_runtime] = 0.0
    data[:random_code] = ""

    data[:greedy_objval] = 0.0
    data[:greedy_runtime] = 0.0
    data[:greedy_code] = ""

    #data.objval = data.objval ./ log(2)

    #bases = [all_bases for _ in 1:20] #Used for running sites data
    sites = map(x -> split(x, ',') .|> LongDNASeq, data.sites)
    for i = 1:nrow(data)
        bases = [all_bases for _ in 1:data.oligo_lengths[i]] #used for running lengths data
        stats = @timed cutfree_rollout(bases, sites[i], nsims=1000)
        data[i,:random_objval] = degeneracy(stats.value; uselog2=false)
        data[i,:random_runtime] = stats.time
        data[i,:random_code] = convert(String, stats.value)

        stats = @timed cutfree_rollout(bases, sites[i], simulate=simulate_greedy)
        data[i,:greedy_objval] = degeneracy(stats.value;uselog2=false)
        data[i,:greedy_runtime] = stats.time
        data[i,:greedy_code] = convert(String, stats.value)
        println(i)
    end
    return data
end

#data = runtimes(CSV.read("./results.csv")[1:10,:])
#data = runtimes(CSV.read("./results_sites_4_1_21.csv"))
#CSV.write("rollout_results_sites_4_1_21.csv",data)
#data = runtimes(CSV.read("./results_length_4_15_21_combined.csv"))
#CSV.write("rollout_results_length_4_15_21.csv",data)
function read_restriction_sites(file="cutfree_rebase_data.xlsx") #File taken from the rebase.R code. Filtered for the 56 commerially available, palindromic, 6 bp, restriction sites
    xf=XLSX.readxlsx(file);
    x=xf[XLSX.sheetnames(xf)[1]]
    blocking_sites=[]
    dims=size(x[:]);
    n=dims[1]
 for i in 2:n #start from second row because excel uses row 1 as the label
        seq=LongDNASeq(string(x[i,7])) #the sequences are the 7th column of this dataframe
        push!(blocking_sites,seq)
    return blocking_sites
end
end

function random_subset(data;size=10)
    shuffled=shuffle(data)
    if size >length(data)
        return shuffled
    else
        return shuffled[1:size]
    end
end

function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table s begin
        "--sequence","-s"
            help="Starting DNA sequence that should be blocked from containing restriction sites. To generate a set of barcodes with the highest diversity, start with a string of N's the length of your oligo. "
            required=true
        "--restrictionsites","-r"
            help = "Sequences to block from the oligo pools. Separate multiple sequences by commas. Do not include spaces."
            required=true
        "--nsims","-n"
            help = "Number of simulations per rollout"
            arg_type= Int
            default=100
    end

    return parse_args(s)
end

function main()
    parsed_args = parse_commandline()
    println("Parsed args:")
    for (arg,val) in parsed_args
        println("  $arg  =>  $val")
    end

    allowedbasedict=Dict(DNA_A => dna"A", DNA_T => dna"T", DNA_C => dna"C", DNA_G => dna"G",DNA_R =>dna"AG",DNA_Y=>dna"CT",DNA_S=>dna"GC",
    DNA_W=>dna"AT",DNA_K=>dna"GT",DNA_M=>dna"AC",DNA_B=>dna"CGTYSKB",DNA_D=>dna"AGTRWKD",DNA_H=>dna"ACTYWMH",DNA_V=>dna"ACGRSMV",DNA_N=>dna"ACGTRYSWKMBDHVN",DNA_Gap=>dna"-")
    sequence=parsed_args["sequence"]
    sites=parsed_args["restrictionsites"]
    nsims=parsed_args["nsims"]
    sequence=LongDNASeq.(sequence)
    bases=[dna"-" for i=1:length(sequence)]
    for i=1:length(sequence)
        bases[i]=allowedbasedict[sequence[i]]
    end
    sites=LongDNASeq.(split(sites,","))
    randomer=cutfree_rollout(bases,sites;nsims=nsims)
    randomer=String(randomer)
    show(stdout,"text/plain",randomer)
    println("\n")
end

main()


text:
using DeformableRegistration: Regularizer, Examples, ImageProcessing, Transformation, Interpolation, Distance
 using DeformableRegistration.regOptions
 using Base.Test

#using DeformableRegistration.Visualization # does not work on travis-ci

# setup logging

@testset "artificial data SSD" begin
    # create test images
    data = zeros(120,120); data[31:90,21:60] = 1
    temImg = createImage(data)
    dataR = deepcopy(data)
    dataR[41:80,41:90] = 1
    refImg = createImage(dataR)
    options = DeformableRegistration.regOptions()
    options.levels = [5,4,3,2]
    options.matrixFree = true;
    options.interpolateToReferenceImage = true
    options.regularizerWeight = 1

    ## register images nonparametric
    displacement = registerImagesNonParametric(refImg, temImg, options) # regularizerOperator=createDiffusiveOperatorCentered
    ssdvalue = ssdDistance(refImg, temImg, getCellCenteredGrid(refImg).data + displacement.data, options=options)[1]
    ngfvalue = ngfDistance(refImg, temImg, getCellCenteredGrid(refImg).data + displacement.data, options=options)[1]
    #visualizeResults(refImg, temImg, displacement=displacement)

    @test ssdvalue ≈ 42.63 atol=1e-1
    @test ngfvalue ≈ 365.01 atol=1e-1

end

@testset "artificial data NGF" begin

	# NGF, smooth data first to generate gradients
	data = zeros(120,120); data[31:90,21:60] = 1
	dataT = copy(data)
	data = conv2([0.2 0.5 0.2
				 0.5 1.0 0.5
				 0.2 0.5 0.2],data)[1:120,1:120]
	temImg = createImage(data)
	dataT[41:80,41:90] = 1
	dataT = conv2([0.2 0.5 0.2
				 0.5 1.0 0.5
				 0.2 0.5 0.2],dataT)[1:120,1:120]
	refImg = createImage(dataT)

	options = regOptions()
	options.levels = [4,3,2]
	options.edgeParameterR = 1
	options.edgeParameterT = options.edgeParameterR
	options.regularizerWeight = 100
	options.interpolateToReferenceImage = true
	@time displacement = registerImagesNonParametric(refImg, temImg, options, measureDistance=ngfDistance)

    ssdvalue = ssdDistance(refImg,temImg,getCellCenteredGrid(refImg).data+displacement.data)[1]
	ngfvalue = ngfDistance(refImg,temImg,getCellCenteredGrid(refImg).data+displacement.data)[1]

    Logging.debug(ssdvalue)
    Logging.debug(ngfvalue)
    # These results seem to depend on the current phase of the moon, your platform or something similar.
	# This seems to be a bug in ngfDistance or at least in there and is notet in an issue in GitHub. Please
	# make the test more specific once this is solved.
    # visualizeResults(refImg, temImg, displacement=displacement)

	@test_skip ssdvalue ≈ 498.07 atol=2e2
	@test_skip ngfvalue ≈ 399.403 atol=2e2
##
end

##
@testset "SSD/NGF shifted"  begin
    for D=[ssdDistance]#, ngfDistance]
##
	testimage = dirname(Base.source_path()) * "/testdata/luebeck.jpg"
	referenceImage = loadImage(testimage)
	# define a cell centered grid, transform it and create a template image
	centeredGrid = getCellCenteredGrid(referenceImage)
	affineParametersInitial = [1.0,0.0,-50,0.0,1.0,-50]
	transformedGrid = transformGridAffine(centeredGrid,affineParametersInitial)
	temImg = interpolateImage(referenceImage,transformedGrid,interpolationScheme=InterpLinearFast)
	templateImage = createImage(temImg)
	options = regOptions()
	options.levels = [6,5,4,3]
	options.parametricOnly = false;
	options.regularizerWeight = 100
	options.matrixFree = true;
	options.interpolateToReferenceImage = false
	options.useEdgeParameterInNumerator = true

	displacement = registerImagesNonParametric(referenceImage, templateImage, options, measureDistance=D)
	@test mean(displacement.data) ≈ -affineParametersInitial[3] atol = 10

    end
end


text:
@testset "interface" begin
    problem = MultivariateProblems.UnconstrainedProblems.examples["Exponential"]
    f = MVP.objective(problem)
    g! = MVP.gradient(problem)
    h! = MVP.hessian(problem)
    nd = NonDifferentiable(f, fill!(similar(problem.initial_x), 0))
    od = OnceDifferentiable(f, g!, fill!(similar(problem.initial_x), 0))
    td = TwiceDifferentiable(f, g!, h!, fill!(similar(problem.initial_x), 0))
    tdref = TwiceDifferentiable(f, g!, h!, fill!(similar(problem.initial_x), 0))
    ref = optimize(tdref, problem.initial_x, Newton(), Optim.Options())
    # test AbstractObjective interface
    for obj in (nd, od, td)
        res = []
        push!(res, optimize(obj, problem.initial_x))

        push!(res, optimize(obj, problem.initial_x, Optim.Options()))

        for r in res
            @test norm(Optim.minimum(ref)-Optim.minimum(r)) < 1e-6
        end
    end
    ad_res = optimize(od, problem.initial_x, Newton())
    @test norm(Optim.minimum(ref)-Optim.minimum(ad_res)) < 1e-6
    ad_res2 = optimize(od, problem.initial_x, Newton())
    @test norm(Optim.minimum(ref)-Optim.minimum(ad_res2)) < 1e-6
    # test f, g!, h! interface
    for tup in ((f,), (f, g!), (f, g!, h!))
        fgh_res = []
        push!(fgh_res, optimize(tup..., problem.initial_x))
        for m in (NelderMead(), LBFGS(), Newton())
            push!(fgh_res, optimize(tup..., problem.initial_x; f_tol = 1e-8))
            push!(fgh_res, optimize(tup..., problem.initial_x, m))
            push!(fgh_res, optimize(tup..., problem.initial_x, m, Optim.Options()))
        end
        for r in fgh_res
            @test norm(Optim.minimum(ref)-Optim.minimum(r)) < 1e-6
        end
    end
end

@testset "only_fg!, only_fgh!" begin
    f(x) = (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2
    function g!(G, x)
      G[1] = -2.0 * (1.0 - x[1]) - 400.0 * (x[2] - x[1]^2) * x[1]
      G[2] = 200.0 * (x[2] - x[1]^2)
      G
    end
    function h!(H, x)
      H[1, 1] = 2.0 - 400.0 * x[2] + 1200.0 * x[1]^2
      H[1, 2] = -400.0 * x[1]
      H[2, 1] = -400.0 * x[1]
      H[2, 2] = 200.0
      H
    end
    function fg!(F,G,x)
      G == nothing || g!(G,x)
      F == nothing || return f(x)
      nothing
    end
    function fgh!(F,G,H,x)
      G == nothing || g!(G,x)
      H == nothing || h!(H,x)
      F == nothing || return f(x)
      nothing
    end

    result_fg! = Optim.optimize(Optim.only_fg!(fg!), [0., 0.], Optim.LBFGS()) # works fine
    @test result_fg!.minimizer ≈ [1,1]
    result_fgh! = Optim.optimize(Optim.only_fgh!(fgh!), [0., 0.], Optim.Newton())
    @test result_fgh!.minimizer ≈ [1,1]
end


text:
#  MIPLearn: Extensible Framework for Learning-Enhanced Mixed-Integer Optimization
#  Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
#  Released under the modified BSD license. See COPYING.md for more details.

using DataStructures
abstract type VariableBranchingRule end

struct Variable
    index::Any
end

mutable struct MIP
    constructor::Any
    optimizers::Vector
    binary_variables::Vector{Variable}
    sense::Float64
    lp_iterations::Int64
end

struct Node
    mip::MIP
    index::Int
    depth::Int
    obj::Float64
    status::Symbol
    branch_variables::Array{Variable}
    branch_values::Array{Float64}
    fractional_variables::Array{Variable}
    fractional_values::Array{Float64}
    parent::Union{Nothing,Node}
end

Base.@kwdef mutable struct History
    avg_pseudocost_up::Float64 = 1.0
    avg_pseudocost_down::Float64 = 1.0
end

mutable struct VariableHistory
    fractional_values::Array{Float64}
    obj_change_up::Array{Float64}
    obj_change_down::Array{Float64}
    obj_ratio_up::Array{Float64}
    obj_ratio_down::Array{Float64}
    pseudocost_up::Float64
    pseudocost_down::Float64

    VariableHistory() = new(
        Float64[], # fractional_values
        Float64[], # obj_change_up
        Float64[], # obj_change_down
        Float64[], # obj_ratio_up
        Float64[], # obj_ratio_up
        0.0, # pseudocost_up
        0.0, # pseudocost_down
    )
end

Base.@kwdef mutable struct NodePool
    mip::MIP
    pending::PriorityQueue{Node,Float64} = PriorityQueue{Node,Float64}()
    processing::PriorityQueue{Node,Float64} = PriorityQueue{Node,Float64}()
    processed::Int = 0
    next_index::Int = 1
    lock::ReentrantLock = ReentrantLock()
    primal_bound::Float64 = Inf
    dual_bound::Float64 = Inf
    gap::Float64 = Inf
    history::History = History()
    var_history::Dict{Variable,VariableHistory} = Dict()
end


text:
export proxTV!, normTV

"""
    proxTV!(x::Vector{T}, λ::Float64; kargs...) where T

proximal map for ansitropic TV regularization using the Fast Gradient Projection algorithm.

# Arguments
* `x::Array{T}`             - Vector to apply proximal map to
* `λ::Float64`              - regularization paramter
* `shape::NTuple=[]`        - size of the underlying image
* `iterationsTV::Int64=20`  - number of FGP iterations
* `weights::Array=[]`       - weights to apply to the image gradients
"""
function proxTV!(x::Vector{T}, λ::Float64; shape=[], iterationsTV::Int64=20, weights::Array=[], kargs...) where T
  m,n = shape

  # initialize dual variables
  p = zeros(eltype(x), m-1,n)
  q = zeros(eltype(x), m,n-1)
  r = zeros(eltype(x), m-1,n)
  s = zeros(eltype(x), m,n-1)
  pOld = zeros(eltype(x), m-1,n)
  qOld = zeros(eltype(x), m,n-1)
  weights!=[] ? weights=reshape(weights,shape) : weights=ones(shape)

  t = 1
  for i=1:iterationsTV
    pOld[:] = p[:]
    qOld[:] = q[:]

    # gradient projection step for dual variables
    p,q = [r,s] + collect( Φ_hermitian( x-λ*Φ(r,s), shape ) ) ./(8*λ)
    p = restrictMagnitude(p, weights[1:m-1,:])
    q = restrictMagnitude(q, weights[:,1:n-1])

    # form linear combinaion of old and new estimates
    tOld = t
    t = ( 1+sqrt(1+4*tOld^2) )/2.0
    r = p + (tOld-1)/t*(p-pOld)
    s = q + (tOld-1)/t*(q-qOld)
  end

  x[:] = x[:]-λ*vec(Φ(p,q))
end

#
# primal-dual mapping
#
function Φ(p::Matrix{T}, q::Matrix{T}) where T
  m = size(q,1)
  n = size(p,2)
  x = zeros(T, m, n)

  # points in the corners
  x[1,1] = p[1,1]+q[1,1]
  x[m,1] = q[m,1]-p[m-1,1]
  x[1,n] = p[1,n]-q[1,n-1]
  x[m,n] = -p[m-1,n]-q[m,n-1]
  # remaining points with i=1 (first row)
  x[1,2:n-1] = p[1,2:n-1]+q[1,2:n-1]-q[1,1:n-2]
  # remaingin points with j=1 (first column)
  x[2:m-1,1] = p[2:m-1,1]+q[2:m-1,1]-p[1:m-2,1]
  #remaining points with i=m (last row)
  x[m,2:n-1] = q[m,2:n-1]-p[m-1,2:n-1]-q[m,1:n-2]
  #remaingin points with j=n (last column)
  x[2:m-1,n] = p[2:m-1,n]-p[1:m-2,n]-q[2:m-1,n-1]
  #all remaining (inner) matrix elements
  x[2:m-1,2:n-1] = p[2:m-1,2:n-1]+q[2:m-1,2:n-1]-p[1:m-2,2:n-1]-q[2:m-1,1:n-2]

  return vec(x)
end

#
# hermitian conjugate of the primal-dual mapping
#
function Φ_hermitian(x::Vector{T}, shape::NTuple) where T
  m,n = shape
  x = reshape(x,m,n)
  p = zeros(T,m-1,n)
  q = zeros(T,m,n-1)
  p[:,:] = x[1:m-1,1:n]-x[2:m,1:n]
  q[:,:] = x[1:m,1:n-1]-x[1:m,2:n]

  return p,q
end

# restrict x to a number smaller then one
# restrictMagnitude(x::Vector) = x/max(1, abs(x))
function restrictMagnitude(x::Array, w::Array=[])
  w != [] ? maxval = w : maxval = ones(size(x))
  return x./max.(1.0, abs.(x)./maxval)
end

"""
    normTV(x::Vector{T},λ::Float64;shape::NTuple=[],kargs...) where T

returns the value of the ansisotropic TV-regularization term.
Arguments are the same as in `proxTV!`
"""
function normTV(x::Vector{T},λ::Float64;shape::NTuple=[],kargs...) where T
  x = reshape(x,shape)
  tv = norm(vec(x[1:end-1,1:end-1]-x[2:end,1:end-1]),1)
  tv += norm(vec(x[1:end-1,1:end-1]-x[1:end-1,2:end]),1)
  tv += norm(vec(x[1:end-1,end]-x[2:end,end]),1)
  tv += norm(vec(x[end,1:end-1]-x[end,2:end]),1)
  return λ*tv
end


text:
@inline tuplejoin(x) = x
@inline tuplejoin(x, y) = (x..., y...)
@inline tuplejoin(x, y, z...) = (x..., tuplejoin(y, z...)...)

function exponentofRterm(T::RExpr)
    return v->parse(Algebra.deg(T,v))
end

function listofRterms(P::RExpr)
    return RExpr.(split(String(P),r"\+|\-"))
end

function exponentofRpoly(P::RExpr, vars)
    return [exponentofRterm(t).(vars) for t in listofRterms(P)]
end

function Rcoeffs(P::RExpr, var)
    return collect(RExpr.(parse(Algebra.coeff(P,var))))
end

# function Rcoeffs(P::RExpr, var)
#     return RExpr.(split(String(Algebra.coeff(q,var))[2:end-1], ","))
# end

function degOfPoly(P::RExpr)
    return v -> parse(Algebra.deg(P,v))
end

function RpolyToTensor(P::RExpr,vars)
    dims = Tuple(1 .+ degOfPoly(P).(vars))
    N = size(dims,1)-1
    out = zeros(RExpr,dims)
    out[:,fill(1,N)...] = Rcoeffs(P,vars[1])
    for n in 1:N
        for I in CartesianIndices(dims[1:n])
            # l = Rcoeffs(out[(Tuple(I)..., fill(1, N-n+1)...)...],vars[n+1])
            l = Rcoeffs(out[I, fill(1, N-n+1)...],vars[n+1])
            out[I, 1:size(l,1), fill(1, N-n)...] = l
        end
    end
    return out
end

function RpolyToNemo(RP::RExpr,Rvars,NR::MPolyRing,Nvars=Nemo.gens(R))
    exp = exponentofRpoly(RP, Rvars)
    coeffs = parse.(RpolyToTensor(RP, Rvars))
    newP = AbstractAlgebra.MPolyBuildCtx(NR)
    for e in exp
        AbstractAlgebra.push_term!(newP, base_ring(NR)(coeffs[(1 .+ e)...]), e)
    end
    return AbstractAlgebra.finish(newP)
end

function NemoToRpoly(NP)
    return RExpr(NP)
end


text:
# Julia1_2.4.jl Operators and Elementary Math Functions
#
#Operators + - * / \ ^
e=Base.MathConstants.e
a=2//3 + 3.4*π/e
b=(3//4)^2 - 55\5  # x=a\b is x=b/a

V=[1, 2, 3]    # vector
3V

M=[1 2; 3 4] # 2x2 matrix
Y=3M         # * implied, 3* each entry

Z1=M/3  # matrix M divided by 3
Z2=3\M  # divide into -- same result

W1=(3+3im)\M
W2=(3+3im)W1
# op= Operators
i=0
i+=1
M*=2 # Stays integer, all antries times 2
M.-=1 # Stays integer, -1 from each entry
M*=2.0 # now floating, all entries * 2.0
# Special functions for complex numbers
real(1+2im)
imag(1+2im)
conj(1+2im)

# Basic function library
d=3+3im
a=sqrt(d)
a^2

# notice sqrt(-1) gives an error
sqrt(-1)
# but sqrt(-1 + 0im) works

a1=abs(1+2im)
a2=abs2(1+2im)
a3=sqrt(-1+0im)
a4=sqrt(-1+2im)
a5=sin(1+2im)
a6=cosh(1+2im)
a7=exp(1+2im)
a8=log(1+2im)

# Miscellaneous functions for Ints and Floats
# also round(x), floor(x), ceil(x), trunc(x),
# rem(x,y), mod(x,y), gcd(x,y), lcm(x,y) and many variations
# also hypot(x,y) hypotenuse, cbrt(x) cube root,
# various log  and exp functions
println(ceil(4.3))
println(round(4.3))
println(floor(4.3))
println(hypot(3,4))
println(cbrt(27))
#


text:
type GPUArray{T<:FloatingPoint)
	id = ccall((:new_array, "libGPUArray"), )

	finalizer(gpu, gpu -> ccall((:gpu_free, "libGPUArray"), Void, ()))
end

immutable type GPUArray{T<:FloatingPoint, N}
	id::Int
	dim::NTuple{N, Int}

	GPUArray(d::NTuple{N, Int}) = ( id = ccall((:allocate_array, "libGPUArray"), Int, (Any,) , (&d)); new(id, d) )
end

GPUArray{T<:FloatingPoint, N}(n::Int...) = GPUArray{T,N}(n)

max_stars_repo_name: magania/JuliaNV

max_issues_repo_name: magania/JuliaNV

max_forks_repo_name: magania/JuliaNV


text:
module ptgibbs

include("constrained_distributions.jl")
include("helpers.jl")
include("processors.jl")
include("plotters.jl")

import StatsBase: rle, pweights
import RLEVectors: rep
import DataFrames: DataFrame, colwise

using Statistics
using Distributions
using LinearAlgebra
using Lazy
using ProgressMeter

export make_gibbs_update
function make_gibbs_update(dat::DataFrame, param::Array, hyp::Tuple, alpha::Array, labels::Array)
    NIW = map(x -> x[1], param)
    prop = map(x -> x[2], param)
    z = map(x -> x[3], param)

    kappa0, mu0, Psi0 = hyp

    nw = size(prop, 1)
    nm = size(prop[1],1)
    n, dm = size(dat)

    # Sample new mean estimates from constrained MVN distribution in each cluster
    NIW_out = copy(NIW)
    for i in 1:nw
        rles = @> begin
                @inbounds z[i]
                @>sort()
                @>rle()
            end
        nz = zeros(nm)
        @inbounds nz[rles[1]] = rles[2]

        # xbar is an array of d dictionaries
        xbar = @views colwise(x -> tapply_mean(z[i], x), dat)
        for m in 1:nm
            if nz[m] >= dm
                # Compute this only once because it gets reused a lot here
                xbarmap = map(x -> x[m], xbar)

                # Extract the current estimate for Sigma for this cluster
                Sigma_hat = get(NIW[i,1][m], "Sigma", 0)

                # Pre-compute some inverses we use multiple times
                iS0 = inv(Sigma_hat ./ max(kappa0[m], dm))
                iS = inv(Sigma_hat)

                # Check if we are ok for mu
                likelihood_check = rep(false, times = dm)
                rcMVN_in = inv(iS0 + nz[m] .* iS) * (iS0 * mu0[m,:] + nz[m] .* iS * xbarmap)

                for j=1:dm
                    if sign(rcMVN_in[j]) != sign(labels[m][j]) & sign(labels[m][j]) != 0
                        likelihood_check[j] = true
                    end
                end

                if any(likelihood_check)
                    # Sample from the prior for mu
                    @inbounds mu =
                    rand_constrained_MVN(
                        round.(
                            Matrix(Hermitian(Sigma_hat)); digits = 8),
                        mu0[m,:],
                        labels[m]
                    )
                else
                    @inbounds mu =
                        rand_constrained_MVN(
                            inv(iS0 + nz[m] .* iS),
                            rcMVN_in,
                            labels[m]
                        )
                end
                @inbounds NIW_out[i,1][m]["mu"] = mu
            else
                # Draw from the prior
                Sigma_hat = get(NIW[i,1][m], "Sigma", 0)
                @inbounds mu = rand_constrained_MVN(
                    round.(
                        Matrix(Hermitian(Sigma_hat)); digits = 8),
                    mu0[m,:],
                    labels[m]
                )
                @inbounds NIW_out[i,1][m]["mu"] = mu
            end
        end
    end

    # Sample new cluster labels
    zout = copy(z)
    for i in 1:nw
        distns = map(x -> MvNormal(x["mu"], Matrix(Hermitian(x["Sigma"]))), NIW_out[i,1])
        p = Array{Float64,2}(undef,n,nm)
        for m in 1:nm
            @inbounds p[:,m] = pdf(distns[m], Matrix(dat)') * prop[i][m]
        end
        @inbounds zout[i] = mapslices(x -> sample(1:nm, pweights(x)), p, dims = 2)[:,1]
    end

    # Sample new cluster weights
    propout = copy(prop)
    for i in 1:nw
        rles = @> begin
            zout[i]
            @>sort()
            @>rle()
        end
        nz = zeros(nm)
        nz[rles[1]] = rles[2]
        propout[i] = rand(Dirichlet(alpha .+ nz))
    end

    (zout, NIW_out, propout)
end

export log_likelihood
function log_likelihood(dat::DataFrame,
                        param::Array,
                        nz::Float64,
                        walker_num::Int64,
                        cluster_num::Int64,
                        curr_Sigma::Array,
                        Sigma_star::Array,
                        mu_hat::Array)
    NIW = map(x -> x[1], param)
    prop = map(x -> x[2], param)
    z = map(x -> x[3], param)

    n, dm = size(dat)
    nm = size(prop[walker_num], 1)
    subdat = @inbounds Matrix(dat[z[walker_num] .== cluster_num,:])

    # If there are enough observations in the cluster, compute the log difference of the likelihoods
    if nz >= dm
        mvn_distn_star = MvNormal(mu_hat, Sigma_star)
        mvn_distn_curr = MvNormal(mu_hat, curr_Sigma)
        ll_normal = sum(logpdf(mvn_distn_star, subdat')) - sum(logpdf(mvn_distn_curr, subdat'))
    else
        ll_normal = 0.0
    end

    # Establish the multivariate distributions that describe the mixture
    distns_curr = map(x -> MvNormal(get(x, "mu", 0), Matrix(Hermitian(get(x, "Sigma", 0)))), NIW[walker_num,1])
    distns_star = deepcopy(distns_curr)
    distns_star[cluster_num] = MvNormal(get(NIW[walker_num,1][cluster_num], "mu", 0), Sigma_star)

    # Compute density of the data in each mixture component
    p_curr = Array{Float64,2}(undef,n,nm)
    for m in 1:nm
        @inbounds p_curr[:,m] = pdf(distns_curr[m], Matrix(dat)') * prop[walker_num][m]
    end

    p_star = deepcopy(p_curr)
    @inbounds p_star[:,cluster_num] = pdf(distns_star[cluster_num], Matrix(dat)') * prop[walker_num][cluster_num]

    # Convert density to probabilities
    @views prob_curr = @> begin
            p_curr
            @>> mapslices(x -> all(x .== 0.0) ? rep(1, times = nm) : x, dims = 2)
            @>> mapslices(x -> x/sum(x), dims = 2)
    end

    @views prob_star = @> begin
            p_star
            @>> mapslices(x -> all(x .== 0.0) ? rep(1, times = nm) : x, dims = 2)
            @>> mapslices(x -> x/sum(x), dims = 2)
    end

    # Plug probabilities into the multinomial likelihood
    mult_distns_curr = mapslices(x -> Multinomial(1,x), prob_curr, dims = 2)
    ll_mult_curr = zeros(n)
    for i in 1:n
        occ = Int64.(zeros(nm))
        @inbounds occ[z[walker_num][i]] = 1
        @inbounds ll_mult_curr[i] = logpdf(mult_distns_curr[i], occ)
    end

    mult_distns_star = mapslices(x -> Multinomial(1,x), prob_star, dims = 2)
    ll_mult_star = zeros(n)
    for i in 1:n
        occ = Int64.(zeros(nm))
        @inbounds occ[z[walker_num][i]] = 1
        @inbounds ll_mult_star[i] = logpdf(mult_distns_star[i], occ)
    end

    sum1 = sum(ll_mult_star)
    sum2 = sum(ll_mult_curr)
    if isnan(sum1-sum2)
        return ll_normal
    else
        return sum1 - sum2 + ll_normal
    end
end

export logprior
function logprior(curr_Sigma::Array, Sigma_star::Array, mu_hat::Array, hyp::Tuple, cluster_num::Int64)
    """
    The prior is based on two main parts
        1. the density of curr_Sigma and Sigma_star given the prior on Sigma
        2. the density of the current estimate of mu given the MVNormal
            parameterized by either curr_Sigma or Sigma_star
    """
    kappa0, mu0, Psi0 = hyp
    dm = size(curr_Sigma, 1)
    mnz = @inbounds max(kappa0[cluster_num], dm + 2)

    #norm_distn_curr = @inbounds MvNormal(mu0[cluster_num,:], curr_Sigma .* old_tune_df ./ mnz)
    #norm_distn_star = @inbounds MvNormal(mu0[cluster_num,:], Sigma_star .* tune_df ./ mnz)
    norm_distn_curr = @inbounds MvNormal(mu0[cluster_num,:], curr_Sigma)
    norm_distn_star = @inbounds MvNormal(mu0[cluster_num,:], Sigma_star)
    invwish_distn = @inbounds InverseWishart(mnz, Psi0[:,:,cluster_num])

    linvwish_ratio = logpdf(invwish_distn, Sigma_star ./ (mnz - dm - 1)) - logpdf(invwish_distn, curr_Sigma ./ (mnz - dm - 1))
    lnorm_ratio = logpdf(norm_distn_star, mu_hat) - logpdf(norm_distn_curr, mu_hat)

    linvwish_ratio + lnorm_ratio
end

export get_lhastings
function get_lhastings(curr_Sigma::Array, Sigma_star::Array, tune_df::Float64)
    """
    The ratio of two Wishart distributions (all indicators cancel out)
    """
    dm = size(Sigma_star, 1)

    term1 = ((2*tune_df-dm-1) / 2) * (logdet(curr_Sigma) - logdet(Sigma_star))
    term2 = (tr(inv(curr_Sigma) * Sigma_star) - tr(inv(Sigma_star) * curr_Sigma)) .* tune_df ./ 2

    term1 + term2
end

export propose_Sigma
function propose_Sigma(curr_Sigma::Array, lab::Array, tune_df::Float64)
    """
      Draw a covariance from the constrained Wishart distribution
    """
    rand_constrained_Wish(curr_Sigma, tune_df, lab) ./ tune_df
end

export make_mcmc_move
function make_mcmc_move(dat::DataFrame,
                        param::Array,
                        hyp::Tuple,
                        alpha::Array{Float64,1},
                        labels::Array,
                        tune_df::Array{Float64,1})
    """
      Make proposals cluster-at-a-time
      For current parameter estimates and current cluster,
        propose a new covariance according to the constrained Wishart
      Return acpt, which tracks acceptances in each cluster
        (1 => accept, 0 => reject)
    """
    NIW = map(x -> x[1], param)
    prop = map(x -> x[2], param)
    z = map(x -> x[3], param)
    kappa0, mu0, Psi0 = hyp

    nm = size(prop[1], 1)
    nw = size(prop, 1)

    # Tracking the acceptance rate
    acpt = zeros(nm)

    # For each walker, for each cluster, sample a covariance from the cWISH
    for i in 1:nw
            rles = @> begin
                    @inbounds z[i]
                    @>sort()
                    @>rle()
                end
            nz = zeros(nm)
            @inbounds nz[rles[1]] = rles[2]

            mu_hats = map(x -> get(x, "mu", 0), NIW[i,1])
            curr_Sigmas = map(x -> get(x, "Sigma", 0), NIW[i,1])
        for m in 1:nm
            # Draw Sigma star from constrained Wishart distribution
            Sigma_star = @inbounds propose_Sigma(curr_Sigmas[m], labels[m], tune_df[m])
            lhaste = @inbounds get_lhastings(curr_Sigmas[m], Sigma_star, tune_df[m])

            # Compute the log prior for this proposal - log prior for the current estimate
            lp = @inbounds logprior(Matrix(Hermitian(curr_Sigmas[m])), Matrix(Hermitian(Sigma_star)), mu_hats[m], hyp, m)

            # Compute the log likelihood for this proposal - log likelihood for the current estimate
            ll = @inbounds log_likelihood(dat, param, max(nz[m], kappa0[m]), i, m, Matrix(Hermitian(curr_Sigmas[m])), Matrix(Hermitian(Sigma_star)), mu_hats[m])

            # If random uniform small enough, update Sigma to Sigma_star
            if log(rand(Uniform(0,1))[1]) < (ll + lp + lhaste)
                NIW[i,1][m]["Sigma"] = Sigma_star
                acpt[m] = 1
            end
        end
    end

    (NIW, acpt)
end

export run_mcmc
function run_mcmc(dat::DataFrame,
                    param::Array,
                    hyp::Tuple,
                    alpha::Array{Float64,1},
                    nstep::Int64,
                    labels::Array{String,1},
                    tune_df::Array{Float64,1};
                    opt_rate::Float64 = 0.3)
    """
      The function will return `(chain, acpt_chain, tune_df_chain)` where
        each returned chain value has an extra dimension appended counting steps of the
        chain (so `chain` is of shape `(ndim, nwalkers, nstep)`, for example).
      * acpt_chain tracks the acceptance rate for each cluster across the chain
      * tune_df_chain tracks the tuning degrees of freedom in the wishart proposal across the chain
      * dat is an n x nd array of observations
      * alpha is an nm array of hyperparameters for the mixing proportions
      * nstep = the number of steps of the already tuned mcmc
      * nd = number of dimensions
      * nw = number of walkers
      * param contains current parameter estimates for each dimension
            across walkers
    """

    labels = hcat([(parse.(Int64, split(x, "")) .-1) for x in labels])

    nw = size(param, 1)
    n, nd = size(dat)
    nm = size(param[1,1][2],1)

    chain = Array{
                    Tuple{Array{Dict{String,Array{Float64,N} where N},1},
                    Array{Float64,1},
                    Array{Int64,1}}
                }(undef, nw, nstep + 1)

    for j in 1:nw
        chain[j,1] = deepcopy(param[j,1])
    end
    acpt_tracker = zeros(nm)

    # For tracking the acceptance rate, per cluster, for the adaptive tuning variance
    win_len = min(nstep, 50)
    acpt_win = zeros(nm, win_len)
    acpt_chain = zeros(nm, nstep)

    old_tune_df = deepcopy(tune_df)
    tune_df_chain = zeros(nm, nstep)

    @showprogress 1 "Running the MCMC..."  for i in 1:nstep
        # Proposes IW cluster at a time, accepts/rejects/returns new IW
        NIW, acpt = make_mcmc_move(dat, chain[:,i], hyp, alpha, labels, tune_df)
        acpt_tracker = acpt_tracker .+ acpt
        @inbounds acpt_win[:, (i-1) % win_len + 1] = acpt

        # Makes Gibbs updates of means, mixing weights, and class labels
        newz, newNIW, newprop = make_gibbs_update(dat, chain[:,i], hyp, alpha, labels)

        if i > 50
            # Update tuning parameter per cluster
            gamma1 = 10 / (i ^ 0.8)
            old_tune_df = copy(tune_df)
            @inbounds [tune_df[m] = update_tune_df(tune_df[m], mean(acpt_win[m,:]), opt_rate, gamma1) for m in 1:1:nm]
        end

        for j in 1:nw
            chain_link = deepcopy(param)

            @inbounds map!(x -> x, chain_link[j,1][1], newNIW[j,1])
            @inbounds map!(x -> x, chain_link[j,1][2], newprop[j,1])
            @inbounds map!(x -> x, chain_link[j,1][3], newz[j,1])

            @inbounds chain[j,i+1] = deepcopy(chain_link[j,1])
            @inbounds acpt_chain[:,i] = acpt_tracker ./ i
            @inbounds tune_df_chain[:,i] = tune_df
        end
    end
    (chain, acpt_chain, tune_df_chain)
end

end # end the module


text:
# This file was generated, do not modify it. # hide
pred_rfr = predict(rfr_m, rows=test);
rms(pred_rfr, y[test])

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*--------------------------------------------------------------------
#* Author: Gerson Balbueno Bicca
#* $Id$
#*--------------------------------------------------------------------
type Tube_Heat_Transfer
	Tube_Heat_Transfer()=begin
		new(
			positive(Dict{Symbol,Any}(
				:Brief=>"Tube Side Reynolds Number",
				:Default=>1000,
				:Lower=>1
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Nusselt Number",
				:Default=>0.5,
				:Lower=>1e-8
			)),
			heat_trans_coeff(Dict{Symbol,Any}(
				:Brief=>"Tube Side Film Coefficient",
				:Default=>1,
				:Lower=>1e-12,
				:Upper=>1e6
			)),
			fricfactor(Dict{Symbol,Any}(
				:Brief=>"Friction Factor",
				:Default=>0.05,
				:Lower=>1e-10,
				:Upper=>2000
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Tube Side Prandtl Number",
				:Default=>0.5,
				:Lower=>1e-8
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Phi Correction",
				:Default=>1,
				:Lower=>1e-3
			)),
			velocity(Dict{Symbol,Any}(
				:Brief=>"Tube Side Velocity",
				:Lower=>1e-8
			)),
			[:Re,:Nu,:htube,:fi,:PR,:Phi,:Vtube,]
		)
	end
	Re::positive
	Nu::positive
	htube::heat_trans_coeff
	fi::fricfactor
	PR::positive
	Phi::positive
	Vtube::velocity
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export Tube_Heat_Transfer
function atributes(in::Tube_Heat_Transfer,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=false
	fields[:Brief]="to be documented"
	fields[:Info]="to be documented"
	drive!(fields,_)
	return fields
end
Tube_Heat_Transfer(_::Dict{Symbol,Any})=begin
	newModel=Tube_Heat_Transfer()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
using Plots;
using Printf;
gr();
anim = Animation();
f(x) = 1/2*cos(8pi*(x-1/2));

function forward(u_old::Array, i::Int64)
	return -(u[ir2[i]] - 2u[ir[i]] + 9u[i] - 10u[il[i]] + 2u[il2[i]])C / 6;
end

function backward(u_old::Array, i::Int64)
	return -(-2u[ir2[i]] + 10u[ir[i]] - 9u[i] + 2[il[i]] - u[il2[i]])C / 6;
end

function rk4()
	##First Substep
	phi1 = zeros(Float64, nx);
	@inbounds for i = 1:nx
		phi1[i] = u[i] + dt * forward(u, i);
	end

	##Second Substep
	phi2 = zeros(Float64, nx);
	@inbounds for i = 1:nx
		phi2[i] = u[i] + dt * forward((u + phi1) ./ 2, i);
	end

	##Third Substep
	phi3 = zeros(Float64, nx);
	@inbounds for i = 1:nx
		phi3[i] = u[i] + dt * forward((u + phi2) ./ 2, i);
	end

	##Forth Substep
	phi4 = zeros(Float64, nx);
	@inbounds for i = 1:nx
		phi4[i] = u[i] +  dt * forward(phi3, i);
	end

	@inbounds for i = 1:nx
		global u[i] = (phi1[i] + 2phi2[i] + 2phi3[i] + phi4[i]) / 6.0;
	end
end

function main()
#=********Numerical Set Up*********=#
	global c = 1.0;
	global nx = 101;
	global lx = 100;
	global x = range(0.0, stop=lx, length=nx);
	global dx = lx /(nx-1);
	global u  = zeros(Float64, nx);
	@inbounds for i = 1:nx
		if 40.0 < x[i] < 60.0;
			u[i] = f(x[i]);
		end
	end
	global t = 0.0;
	global dt = 5e-1;
	global tlims = 100;
	global C = dt * c / dx;
	global ir = zeros(Int64, nx);
	global ir2 = zeros(Int64, nx);
	global il = zeros(Int64, nx);
	global il2 = zeros(Int64, nx);
	@inbounds for i = 1:nx
		ir[i] = i + 1;
		ir2[i] = i + 2;
		il[i] = i - 1;
		il2[i] = i - 2;
	end
	ir[end] = 1;
	ir2[end - 1] = 1;
	ir2[end] = 2;
	il[begin] = nx;
	il2[begin + 1] = nx;
	il2[begin] = nx - 1;
#=******Finish Numerical Set Up*************=#

	while t < tlims
		rk4();
		plt = plot(x, u, xlims=(0, lx), xticks=0:lx/5:lx, ylims=(-0.2, 1.2), yticks=-0.2:0.2:1.2,
		xlabel="x", ylabel="u", color="blue", title="KK scheme by using RK4", legend=false);
		frame(anim, plt);
		t += dt;
		@printf("t = %.4f\n", t);
	end
	gif(anim, "KK_rk4.gif", fps=5);
end

main();


text:
# Must increase the number of tests and add some that actually do numeric comparisons

@testset "PROJ4" begin
	geod([0, 0], 90, 111000)
	geod([0, 0], [30 60], 111000)
	geod([0, 0], 30, [50, 111], unit=:k)
	geod([0, 0], [30 60], [50, 111], unit=:k)
	geod([0., 0], [15., 45], [[0, 10000, 50000, 111000.], [0., 50000]])[1]
	invgeod([0., 0], [1., 0])
	@test_throws ErrorException("'azimuth' MUST be either a scalar or a 1-dim array, and 'distance' may also be a Vector{Vector}") geod([0, 0], [30 8; 1 1], [50 111], unit=:k)
	GMT.proj_info()

	pj = GMT.proj_create_crs_to_crs("EPSG:4326", "+proj=utm +zone=32 +datum=WGS84", C_NULL)	# target, also EPSG:32632
	#@test GMT.is_latlong(GMT.proj_create("+proj=longlat +datum=WGS84 +no_defs"))
	@test circgeo(0.,0, radius=50, dataset=true, unit=:k).data[1] == 0.0
	@test_throws ErrorException("Bad shape name (a)") circgeo(0.,0, radius=50, shape=:a)
	circgeo([0 0; 0 2], radius=50, unit=:k)
	buffergeo(mat2ds([0 0; 5 5])[1], width=100, unit=:k, tol=0)
	buffergeo(mat2ds([178 73; -175 74]), width=100, unit=:k)
	wkt = epsg2wkt(4326)
	prj = epsg2proj(4326)
	proj2wkt(prj)
	wkt2proj(wkt)
	loxo = loxodrome_direct(0,0,45, 10000)
	loxo = loxodrome([0 0; 30 50], step=500, unit=:k);
	loxo = loxodrome(mat2ds([0 0; 30 50])[1], step=500, unit=:k);
	loxo = loxodrome(0, 0, 30, 50, step=500, unit=:k);
	orto = orthodrome(mat2ds([0 0; 30 50]), step=500, unit=:k);
	orthodrome(mat2ds([0 0; 15 25; 30 50])[1], step=500, unit=:k)
	orthodrome(0, 0, 30, 50, step=500, unit=:k);
	dist, azim = GMT.loxodrome_inverse(0,0,5,5)
end


text:
# FrontendOptions
function getModulesEmbedFilesNum(x::FrontendOptions)
    @check_ptrs x
    n = clang_FrontendOptions_getModulesEmbedFilesNum(x)
    files = Vector{Ptr{Cuchar}}(undef, n)
    clang_FrontendOptions_getModulesEmbedFiles(x, files, n)
    return unsafe_string.(files)
end

function PrintStats(x::FrontendOptions)
    @check_ptrs x
    return clang_FrontendOptions_PrintStats(x)
end