text:
# initalize an array
function Init!(X,L)
    size=length(X)
    h=L/size
    for i=0:size-1
        if i>floor(size//8) && i<floor(size//2)+floor(size//8)
            X[i+1]=1.-2*(i-floor(size/8))*h/L;
        else
            X[i+1]=0.0
        end
    end
end

# proc[1-3]! functions make the same computation, with different programing
# styles.

function proc1!(In::Array{Float64,1},Out::Array{Float64,1},niter::Int64)
    const size=length(In)
    h2::Float64= (1./size)^2
    for it=1:niter
        Out[2:size-1]= h2*(In[1:size-2]- 2.0*In[2:size-1]+ In[3:size])
        In,Out=Out,In
    end
end
function proc2!(In::Array{Float64,1},Out::Array{Float64,1},niter::Int64)
    const size=length(In)
    h2::Float64= (1./size)^2
    for it=1:niter
        @.  Out[2:size-1]= h2*(In[1:size-2]- 2.0*In[2:size-1]+ In[3:size])
        In,Out=Out,In
    end
end
function proc3!(In::Array{Float64,1},Out::Array{Float64,1},niter::Int64)
    const size=length(In)
    h2::Float64= (1./size)^2
    for it=1:niter
        @simd for i=2:size-1
            Out[i]= h2*(In[i-1]- 2.0*In[i]+ In[i+1])
        end
        In,Out=Out,In
    end
end
# measure as precisely as possible the execution time for a given
# p= proc[1-5]!
function test(p,In::Array{Float64,1},Out::Array{Float64,1},nit::Int)

    niter=nit
    # be sure to run once before actually running the benchmark!
    Init!(In,1.)
    Init!(Out,1.)
    p(In,Out,niter)
    #
    T=0.
    while true
        Init!(In,1.)
        Init!(Out,1.)
        t1 = time_ns()
        p(In,Out,niter)
        t = (time_ns() -t1)/niter
        if abs(t-T)/t<0.025
            break
    else
            T=t
            if niter> 150000
                error("too many iterations",niter)
            else
                niter*=2
            end
    end

end
    T,niter
end
DD=Dict("proc1!"=>"Vectorisation with braces         ",
        "proc2!"=>"Vectorisation with @. and braces  ",
        "proc3!"=>"Unrolled loop                     "
        )

fw=open("RunningOn"*gethostname()*"_lapl_1","w")

# computation starts here:
size=32
sizemax=10^5
const niter=2
while size<sizemax
    println("size: ",size)
    In=Array{Float64}(size)
    Out=Array{Float64}(size)
    tbest=10.^29
    best=0
    t=0.0
    for p= [proc1!,proc2!,proc3!]
        t,it=test(p,In,Out,niter)
        if t<tbest
            tbest=t
            best=p
        end
        t*=10.0^(-9)
        println(DD[string(p)]," : t= ",t," seconds ")
    end

    nflops=size*4
    flops=nflops/tbest
    println("\nbest: ",DD[string(best)])

    write(fw,string(size)," ",string(tbest*10.0^(-9)),"\n")

    println("nb. flops: ",nflops, ", Gflops/s (best): ",flops)
    println("-------")

    size*=2
    println()
end
close(fw)

max_stars_repo_name: wolfv/BenchmarksPythonJuliaAndCo

max_issues_repo_name: wolfv/BenchmarksPythonJuliaAndCo

max_forks_repo_name: wolfv/BenchmarksPythonJuliaAndCo


text:
using Documenter
using POMDPSimulators

makedocs(
    format = :html,
    sitename = "POMDPSimulators.jl"
)

deploydocs(
    repo = "github.com/JuliaPOMDP/POMDPSimulators.jl.git",
    julia = "1.0",
    target = "build",
    deps = nothing,
    make = nothing
)

ext: jl


text:
using PCLCommon
using PCLIO
using PCLFilters
using PCLVisualization
using Cxx

pcd_file = Pkg.dir("PCLIO", "test", "data", "table_scene_lms400.pcd")
cloud = PointCloud{PointXYZ}(pcd_file)
cloud_filtered = PointCloud{PointXYZ}()

sor = StatisticalOutlierRemoval{PointXYZ}()
setInputCloud(sor, cloud)
setMeanK(sor, 50)
setStddevMulThresh(sor, 1.0)
filter(sor, cloud_filtered)

viewer = PCLVisualizer("pcl visualizeer", create_interactor=false)
@assert !hasInteractor(viewer)
setOffScreenRendering(viewer, true)

red_handler = PointCloudColorHandlerCustom(cloud, 255, 0, 0)
green_handler = PointCloudColorHandlerCustom(cloud, 0, 255, 0)
addPointCloud(viewer, cloud, red_handler, id="cloud")
addPointCloud(viewer, cloud_filtered, green_handler, id="cloud_filtered")

dst_filename = string(tempname() |> x -> split(x, '/')[end], ".png")
PCLVisualization.renderToPng(viewer, dst_filename)
@assert isfile(dst_filename)

println("renderToPng: $dst_filename was created")
rm(dst_filename)
println("$dst_filename was sucessfully removed")


text:
"""
    gegenbauer(x, p::Integer)

Evaluate the Gegenbauer polynomial with parameter `α` of degree `p` at `x` using
the three term recursion.
"""
function gegenbauer(x, p::Integer, α)
    T = typeof( 2α*x )
    p₀ = one(T)
    p₁ = 2α*x

    if p <= 0
        return p₀
    elseif p == 1
        return p₁
    end

    for n in 2:p
        p₀, p₁ = p₁, ( 2x*(n+α-1)*p₁ - (2α+n-2)*p₀ ) / n
    end

    p₁
end


text:
include("Util_single.jl")
include("Util_RKHS.jl")

@everywhere module MultipleParticle
using StatsBase
using KernelDensity
using SingleParticle
using KernelTrick

export generate_init_particles, generate_init_particles2, resample, get_values, get_weights, get_obj_func_vals, get_stepsizes,
       normalize_weights, reset_weights, drop_zero_weights, update_weights, calc_ESS, calcPA,
       MCstep_all, MCstep_all2, MCstep_all_parallel, MCstep_all_parallel2, update_stepsize_inRKHS_all

function generate_init_particles(N, lbs, ubs, obj_func)
    """
    N: number of samplers
    """
    particles = [init_particle(lbs, ubs, obj_func, 1 / N) for i in 1:N]
    return particles
end

function resample(particles)
    N = length(particles)
    wv = [particle.weight for particle in particles]
    wv = Weights(wv)
    samples = similar(particles)
    for i in 1:N
        samples[i] = deepcopy(sample(particles, wv))
    end
    return samples
end

function resample(particles, N)
    wv = [particle.weight for particle in particles]
    wv = Weights(wv)
    samples = Array{eltype(particles), 1}(N)
    for i in 1:N
        samples[i] = deepcopy(sample(particles, wv))
    end
    return samples
end

function get_values(particles)
    return [particle.x_present for particle in particles]
end

function get_weights(particles)
    return [particle.weight for particle in particles]
end

function get_obj_func_vals(particles)
    return [particle.obj_func_val for particle in particles]
end

function get_stepsizes(particles)
    return [particle.stepsize for particle in particles]
end

function normalize_weights(particles)
    candidates = deepcopy(particles)
    wv = [candidate.weight for candidate in candidates]
    weight_sum = sum(wv)
    for candidate in candidates
        candidate.weight = candidate.weight / weight_sum
    end
    return candidates
end

function reset_weights(particles)
    N = length(particles)
    for particle in particles
        particle.weight = 1 / N
    end
    return particles
end

function drop_zero_weights(particles)
    wv = get_weights(particles)
    particles_new = eltype(particles)[]
    for (i, w) in enumerate(wv)
        if w > 0
            push!(particles_new, particles[i])
        end
    end
    return particles_new
end

function update_weights(particles,
                        obj_func_val_old::Array{Float64,1},
                        ϵ_old::Float64,
                        ϵ_new::Float64)
    for (i, particle) in enumerate(particles)
        W = particle.weight
        F_new = particle.obj_func_val < ϵ_new
        F_old = obj_func_val_old[i] < ϵ_old
        particle.weight = W * (F_new / F_old)
    end
    return particles
end

function calc_ESS(particles)
    """
    NOTE: weights must be normalized before.
    """
    wv = get_weights(particles)
    wv = wv / sum(wv)
    ESS = sum(wv.^2)^(-1)
    return ESS
end

function calcPA(particles, α)
    F = get_obj_func_vals(particles)
    return quantile(F, α)
end

function MCstep_all(particles, ϵ)

    # kernel density estimation of prior
    params = hcat(get_values(particles)...)'
    prior = [kde(params[:, i], npoints=16) for i in 1:size(params)[2]]

    N = length(particles)
    flags = zeros(length(particles))
    for i in 1:N
        particles[i], flags[i] = MCstep(particles[i], prior, ϵ)
    end
    return particles, flags
end

function MCstep_all_parallel(particles, ϵ)
    """
    to fetch results of @parallel for macro from several worker process,
    SharedArray must be required.
    """
    N = length(particles)

    # kernel density estimation of prior
    params = hcat(get_values(particles)...)'
    proposed = SharedArray{Float64}(size(params))
    obj_func_vals = SharedArray{Float64}(N)
    flags = SharedArray{Float64}(N)
    indecies = SharedArray{Int}(N)

    prior = [kde(params[:, i], npoints=16) for i in 1:size(params)[2]]

    @sync @parallel for i in 1:N
        indecies[i] = i
        proposed[i,:], obj_func_vals[i], flags[i] = MCstep_parallel(particles[i], prior, ϵ)
    end

    for i in indecies
        particles[i].x_present = proposed[i,:]
        particles[i].obj_func_val = obj_func_vals[i]
    end

    return particles, flags
end

function update_stepsize_inRKHS_all(particles, γ, nu, acceptance_target, acceptance; parallel=false, N=100, adaptive=false)
    if length(particles) < 3
        nothing
    else
        N = minimum([length(particles), N])
        sampled = sample(particles, N, replace=false)
        σ = bandwidth(sampled)
        if σ < 0
            nothing
        else
            num_param = length(particles[1].x_present)
            C = eye(N)-1/N
            for (i, particle) in enumerate(particles)
                sampled = sample(particles, N, replace=false)
                M = GaussianKernelInRKHS_M(sampled, particle, σ)
                Σ = κ_KAMH(sampled, particle, C, M, γ=γ, nu=nu)
                particle.cov = Σ
            end
            if adaptive==true
                γ = γ * (1. + rand()*(acceptance - acceptance_target))
                nu = nu * (1. + rand()*(acceptance - acceptance_target))
            end
        end
    end
    return particles, γ, nu
end

end


text:
# Bring a few packages into scope:
using BesselK, ForwardDiff, StaticArrays
import BesselK: gamma, adbesselkxv

# Matern covariance function, using the rescaled (x^v)*besselk(v, x):
function matern(x, y, params)
  (sg, rho, nu) = params
  dist = norm(x-y)
  iszero(dist) && return sg*sg
  arg = sqrt(2*nu)*dist/rho
  (sg*sg*(2^(1-nu))/gamma(nu))*adbesselkxv(nu, arg)
end

# Slightly fancy matrix assembly and factorization:
function assemble_matrix(points, params)
  buf = Matrix{eltype(params)}(undef, length(points), length(points))
  Threads.@threads for k in eachindex(points) # nice multi-threading
    ptk = points[k]
    buf[k,k] = matern(ptk, ptk, params)
    @inbounds for j in 1:(k-1) # turns  off array bounds checking
      buf[j,k] = matern(points[j], ptk, params)
    end
  end
  cholesky!(Symmetric(buf)) # in-place Cholesky factorization to avoid heap allocs.
end

# Negative log-likelihood, with a slight trick for the quadratic form to just
# have to solve with the triangular factor once:
function nll(points, data, params)
  Sigma = assemble_matrix(points, params)
  (logdet(Sigma) + sum(abs2, Sigma.U'\data))/2
end

# Sample locations and data, using stack-allocated arrays via StaticArrays.jl to
# make sure that the autodiff derivatives don't make any heap allocations.
# This is just a random example to demonstrate how to create the single-arg
# closure that you can pass to ForwardDiff.
const LOCS = [@SVector rand(2) for _ in 1:1000]
const DAT  = randn(length(LOCS))

# create single-argument closure for the log-likelihood:
objective(p) = nll(LOCS, DAT, p)

# AD-generated gradient and hessian. It's that easy! Plug in to your favorite
# optimizer and you're good to go.
objective_grad(p) = ForwardDiff.gradient(objective, p)
objective_hess(p) = ForwardDiff.hessian(objective, p)


text:
module RunTests

using JuMP

import Ipopt
import PowerModels
import SymbolicAD
import Test

function runtests()
    for name in names(@__MODULE__; all = true)
        if startswith("$(name)", "test_")
            Test.@testset "$(name)" begin
                getfield(@__MODULE__, name)()
            end
        end
    end
    return
end

"""
    _run_clnlbeam_benchmark(; N::Int)

This is a good benchmark for the perils of Symbolic AD, because it has a
nonlinear objective function, and the nonlinear objective function has thousands
of terms.
"""
function _run_clnlbeam_benchmark(; N::Int)
    h = 1 / N
    model = Model()
    @variable(model, -1 <= t[1:(N+1)] <= 1)
    @variable(model, -0.05 <= x[1:(N+1)] <= 0.05)
    @variable(model, u[1:(N+1)])
    @NLobjective(
        model,
        Min,
        sum(
            h / 2 * (u[i+1]^2 + u[i]^2) +
            350 * h / 2 * (cos(t[i+1]) + cos(t[i])) for i in 1:N
        ),
    )
    for i in 1:N
        @NLconstraint(model, x[i+1] - x[i] == h / 2 * (sin(t[i+1]) + sin(t[i])))
        @constraint(model, t[i+1] - t[i] == h / 2 * (u[i+1] - u[i]))
    end
    # Different rtol because the points we pick can be quite bad
    SymbolicAD.run_unit_benchmark(model; direct_model = true, rtol = 1e-6)
    SymbolicAD.run_solution_benchmark(model, Ipopt.Optimizer)
    return
end

function power_model(case::String)
    pm = PowerModels.instantiate_model(
        joinpath(@__DIR__, "data", case),
        PowerModels.ACPPowerModel,
        PowerModels.build_opf,
    )
    return pm.model
end

function test_case5_pjm_unit()
    model = power_model("pglib_opf_case5_pjm.m")
    SymbolicAD.run_unit_benchmark(model; direct_model = false)
    SymbolicAD.run_unit_benchmark(model; direct_model = true)
    return
end

function test_case5_pjm_solution()
    model = power_model("pglib_opf_case5_pjm.m")
    SymbolicAD.run_solution_benchmark(model, Ipopt.Optimizer)
    return
end

function test_clnlbeam()
    _run_clnlbeam_benchmark(; N = 10)
    _run_clnlbeam_benchmark(; N = 100)
    return
end

function test_optimizer_clnlbeam(; N::Int = 10)
    h = 1 / N
    model = Model(() -> SymbolicAD.Optimizer(Ipopt.Optimizer))
    @variable(model, -1 <= t[1:(N+1)] <= 1)
    @variable(model, -0.05 <= x[1:(N+1)] <= 0.05)
    @variable(model, u[1:(N+1)])
    @NLobjective(
        model,
        Min,
        sum(
            h / 2 * (u[i+1]^2 + u[i]^2) +
            350 * h / 2 * (cos(t[i+1]) + cos(t[i])) for i in 1:N
        ),
    )
    for i in 1:N
        @NLconstraint(model, x[i+1] - x[i] == h / 2 * (sin(t[i+1]) + sin(t[i])))
        @constraint(model, t[i+1] - t[i] == h / 2 * (u[i+1] - u[i]))
    end
    optimize!(model)
    Test.@test isapprox(objective_value(model), 350; atol = 1e-6)
    t_sol = value.(t)
    u_sol = value.(u)
    Test.@test isapprox(
        sum(
            h / 2 * (u_sol[i+1]^2 + u_sol[i]^2) +
            350 * h / 2 * (cos(t_sol[i+1]) + cos(t_sol[i])) for i in 1:N
        ),
        350.0;
        atol = 1e-6,
    )
    return
end

function test_optimizer_case5_pjm()
    model = power_model("pglib_opf_case5_pjm.m")
    set_optimizer(model, () -> SymbolicAD.Optimizer(Ipopt.Optimizer))
    optimize!(model)
    symbolic_obj = objective_value(model)
    set_optimizer(model, Ipopt.Optimizer)
    optimize!(model)
    Test.@test isapprox(objective_value(model), symbolic_obj, atol = 1e-6)
    return
end

function test_optimizer_case5_pjm_optimize_hook()
    model = power_model("pglib_opf_case5_pjm.m")
    set_optimizer(model, Ipopt.Optimizer)
    set_optimize_hook(model, SymbolicAD.optimize_hook)
    optimize!(model)
    symbolic_obj = objective_value(model)
    set_optimize_hook(model, nothing)
    optimize!(model)
    Test.@test isapprox(objective_value(model), symbolic_obj, atol = 1e-6)
    return
end

"""
    test_optimize_hook()

This model is chosen to have a number of unique features that cover the range of
node types in JuMP's nonlinear expression data structure.
"""
function test_optimize_hook()
    model = Model(Ipopt.Optimizer)
    @variable(model, 0 <= x <= 2π, start = π)
    @NLexpression(model, y, sin(x))
    @NLparameter(model, p == 1)
    @NLconstraint(model, sqrt(y + 2) <= p + 1)
    @NLobjective(model, Min, p * y)
    set_optimize_hook(model, SymbolicAD.optimize_hook)
    optimize!(model)
    Test.@test isapprox(objective_value(model), -1; atol = 1e-6)
    Test.@test isapprox(value(x), 1.5 * π; atol = 1e-6)
    return
end

end  # module

RunTests.runtests()


text:
module DataTypes

using DataStructures
import Base.convert, Base.show, Base.showerror, Base.log, Base.log10,
       Base.exp, Base.isless, Base.ifelse, Base.sqrt, Base.abs,
       Base.==, Base.min, Base.max

export Parameter, Forcing, State, Variable, Observed
export Species, Compartment, Reactant, Component
export MEReaction, Reaction
export MESource, ReactionSource, OdeSource, OdeSorted, OdeModel
export Dimension, Unit

###################################################
###  Dimensional analysis and unit conversion  ####
###################################################

immutable Dimension
  K::Rational{Int64}
  s::Rational{Int64}
  m::Rational{Int64}
  g::Rational{Int64}
  cd::Rational{Int64}
  mol::Rational{Int64}
  A::Rational{Int64}
end

immutable Unit
  d::Dimension
  f::Number
end

# Unit error (when we try to add or substract different units)
type UnitError <: Exception
end
showerror(io::IO, e::UnitError) = print(io, "There was an error with the units");

# Algebra with dimensions

# +
+(d::Dimension, s::Number) = d != none.d ? (throw(UnitError())) : d
+(s::Number, d::Dimension) = d + s
+(u1::Dimension, u2::Dimension) = u1 == u2 ? u1 : throw(UnitError())
+(d::Dimension) = d

# -
-(d::Dimension, s::Number) = d != none.d ? (throw(UnitError())) : d
-(s::Number, d::Dimension) = d - s
-(u1::Dimension, u2::Dimension) = u1 == u2 ? u1 : throw(UnitError())
-(d::Dimension) = d

# *
*(u1::Dimension, u2::Dimension) =
  Dimension(u1.K + u2.K, u1.s + u2.s, u1.m + u2.m,
            u1.g + u2.g, u1.cd + u2.cd, u1.mol + u2.mol, u1.A + u2.A)
*(u1::Dimension, s::Number) = u1
*(s::Number, u1::Dimension) = u1

# /
/(u1::Dimension, u2::Dimension) =
  Dimension(u1.K - u2.K, u1.s - u2.s, u1.m - u2.m,
            u1.g - u2.g, u1.cd - u2.cd, u1.mol - u2.mol, u1.A - u2.A)
/(u1::Dimension, s::Number) = u1
/(s::Number, u1::Dimension) = u1^(-s)
exp(d::Dimension) = d != none.d ? (throw(UnitError())) : none.d

# ^
^(d1::Dimension, d2::Dimension) = d2 != none.d ? (throw(UnitError())) : d1
^(u1::Dimension, s::Integer) =
  Dimension(u1.K*s, u1.s*s, u1.m*s, u1.g*s, u1.cd*s, u1.mol*s, u1.A*s)
^(u1::Dimension, s::Rational) =
  Dimension(u1.K*s, u1.s*s, u1.m*s, u1.g*s, u1.cd*s, u1.mol*s, u1.A*s)
^(u1::Dimension, s::Float64) =
  Dimension(u1.K*s, u1.s*s, u1.m*s, u1.g*s, u1.cd*s, u1.mol*s, u1.A*s)
^(d1::Number, d2::Dimension) = d2 != none.d ? (throw(UnitError())) : none.d

# exp
exp(d::Dimension) = d != none.d ? (throw(UnitError())) : d

# log
log(d::Dimension) = d != none.d ? (throw(UnitError())) : d

# log10
log10(d::Dimension) = d != none.d ? (throw(UnitError())) : d

# sqrt
sqrt(d::Dimension) = d^(1//2)

# isless
isless(d1::Dimension, d2::Dimension) = d1 != d2 && (throw(UnitError()))
isless(d1::Dimension, s::Number) = d1 != none.d ? (throw(UnitError())) : d1
isless(s::Number, d1::Dimension) = isless(d1, s)

# ifelse comparison
ifelse(d1::Dimension, d2::Dimension, s::Number) = exp(d2::Dimension)
ifelse(d1::Dimension, s::Number, d2::Dimension) = d2 != none.d ? (throw(UnitError())) : none.d
ifelse(d1::Dimension, d2::Dimension, d3::Dimension) = d2 != d3 ? (throw(UnitError())) : d2
ifelse(d1::Dimension, s1::Number, s2::Number) = none.d

# abs
abs(d::Dimension) = d

# <
<(d1::Dimension, d2::Dimension) = d1 != d2 ? (throw(UnitError())) : none.d
<(d1::Dimension, n::Number) = d1 != none.d ? (throw(UnitError())) : none.d
<(n::Number, d1::Dimension) = d1 < n

# >
>(d1::Dimension, d2::Dimension) = d1 != d2 ? (throw(UnitError())) : none.d
>(d1::Dimension, n::Number) = d1 != none.d ? (throw(UnitError())) : none.d
>(n::Number, d1::Dimension) = d1 < n

# max function
max(d1::Dimension, d2::Dimension) = d1 != d2 ? (throw(UnitError())) : d1
max(d1::Dimension, n::Number) = d1 != none.d ? (throw(UnitError())) : none.d
max(n::Number, d1::Dimension) = max(d1,n)

# min function
min(d1::Dimension, d2::Dimension) = d1 != d2 ? (throw(UnitError())) : d1
min(d1::Dimension, n::Number) = d1 != none.d ? (throw(UnitError())) : none.d
min(n::Number, d1::Dimension) = min(d1,n)

# Basic algebra with units (for composite and derived units)
*(u1::Unit, u2::Unit) = Unit(u1.d * u2.d, u1.f*u2.f)
/(u1::Unit, u2::Unit) = Unit(u1.d / u2.d, u1.f/u2.f)
+(u1::Unit, u2::Unit) = u1.d == u2.d ? u1 : throw(UnitError())
-(u1::Unit, u2::Unit) = u1.d == u2.d ? u1 : throw(UnitError())
^(u1::Unit, s::Integer) = Unit(u1.d^s, u1.f^s)
^(u1::Unit, s::Rational) = Unit(u1.d^s, u1.f^s)
*(u1::Unit, s::Number) = Unit(u1.d, u1.f*s)
*(s::Number, u1::Unit) = u1*s
/(u1::Unit, s::Number) = Unit(u1.d, u1.f/s)
/(s::Number, u1::Unit) = u1^(-s)

# Show pretty representation (still need to improve more...)
function show(io::IO, u::Unit)
  unit = u.d.mol != 0 ? ( u.d.mol - 1 < eps(1.0) ? "*mol" : ( u.d.mol < 0 ? "*1/mol^$(-float(u.d.mol))" : "*mol^$(float(u.d.mol))") ) : ""
  unit *= u.d.g != 0 ? ( u.d.g == 1 ? "*g" : ( u.d.g < 0 ? "*1/g^$(-float(u.d.g))" : "*g^$(float(u.d.g))") ) : ""
  unit *= u.d.m != 0 ? ( u.d.m == 1 ? "*m" : ( u.d.m < 0 ? "*1/m^$(-float(u.d.m))" : "*m^$(float(u.d.m))") ) : ""
  unit *= u.d.s != 0 ? ( u.d.s == 1 ? "*s" : ( u.d.s < 0 ? "*1/s^$(-float(u.d.s))" : "*s^$(float(u.d.s))") ) : ""
  unit *= u.d.K != 0 ? ( u.d.K == 1 ? "*K" : ( u.d.K < 0 ? "*1/K^$(-float(u.d.K))" : "*K^$(float(u.d.K))") ) : ""
  unit *= u.d.cd != 0 ? ( u.d.cd == 1 ? "*cd" : ( u.d.cd < 0 ? "*1/cd^$(-float(u.d.cd))" : "*cd^$(float(u.d.cd))") ) : ""
  unit *= u.d.A != 0 ? ( u.d.A == 1 ? "*A" : ( u.d.A < 0 ? "*1/A^$(-float(u.d.A))" : "*A^$(float(u.d.A))") ) : ""
  # Remove trailing "* "
  unit = unit[2:end]
  print(io, "$(float(u.f)) $unit")
end

function show(io::IO, d::Dimension)
  unit = d.mol != 0 ? ( d.mol == 1 ? "*mol" : ( d.mol < 0 ? "*1/mol^$(-float(d.mol))" : "*mol^$(float(d.mol))") ) : ""
  unit *= d.g != 0 ? ( d.g == 1 ? "*g" : ( d.g < 0 ? "*1/g^$(-float(d.g))" : "*g^$(float(d.g))") ) : ""
  unit *= d.m != 0 ? ( d.m == 1 ? "*m" : ( d.m < 0 ? "*1/m^$(-float(d.m))" : "*m^$(float(d.m))") ) : ""
  unit *= d.s != 0 ? ( d.s == 1 ? "*s" : ( d.s < 0 ? "*1/s^$(-float(d.s))" : "*s^$(float(d.s))") ) : ""
  unit *= d.K != 0 ? ( d.K == 1 ? "*K" : ( d.K < 0 ? "*1/K^$(-float(d.K))" : "*K^$(float(d.K))") ) : ""
  unit *= d.cd != 0 ? ( d.cd == 1 ? "*cd" : ( d.cd < 0 ? "*1/cd^$(-float(d.cd))" : "*cd^$(float(d.cd))") ) : ""
  unit *= d.A != 0 ? ( d.A == 1 ? "*A" : ( d.A < 0 ? "*1/A^$(-float(d.A))" : "*A^$(float(d.A))") ) : ""
  # Remove trailing "* "
  unit = unit[2:end]
  print(io, "$unit")
end

# Prefixes
const prefixes = {"" => 1//1,
                  "E" => int(1e18)//1, "P" => int(1e15)//1,
                  "T" => int(1e12)//1, "G" => int(1e9)//1, "M" => int(1e6), "k" => int(1e3),
                  "h" => int(1e2)//1, "da" => int(1e1)//1, "d" => 1//int(1e1), "c" => 1//int(1e2),
                  "m" => 1//int(1e3), "u" => 1//int(1e6), "n" => 1//int(1e9), "p" => 1//int(1e12),
                  "f" => 1//int(1e15), "a" => 1//int(1e18)}

# SI units
const SI_units = {"K" => Unit(Dimension(1//1,0,0,0,0,0,0),1//1),
                  "s" => Unit(Dimension(0,1//1,0,0,0,0,0),1//1),
                  "m" => Unit(Dimension(0, 0, 1//1, 0, 0, 0, 0), 1//1),
                  "g" => Unit(Dimension(0,0,0,1//1,0,0,0),1//1),
                  "cd" => Unit(Dimension(0,0,0,0,1//1,0,0),1//1),
                  "mol" => Unit(Dimension(0,0,0,0,0,1//1,0),1//1),
                  "A" => Unit(Dimension(0,0,0,0,0,0,1//1),1//1),
                  "none" => Unit(Dimension(0,0,0,0,0,0,0),1//1)}

# Automagically generate the SI units
for (pr,prval) = prefixes
    for (key,val) = SI_units
        unit_name = parse(pr*key)
        unit = Unit(val.d, val.f*prval)
        @eval const $unit_name = $unit
    end
end

# Derived units
const derived_units = {"M" => mol/dm^3,
                       "L" => dm^3,
                      "ha" => 100//1*m^2,
                      "Hz" => 1/s,
                      "N" => m*kg/s^2,
                      "Pa" => kg/m/s^2,
                      "J" => m^2*kg/s^2,
                      "W" => m^2*kg/s^3,
                      "C" => s*A,
                      "V" => m^2*kg/s^3/A,
                      "F" => s^4*A^2/kg/m^2,
                      "O" => m^2*kg/s^3/A^2,
                      "S" => s^3*A^2/kg/m^2}

# Automagically generate the derived units
for (pr,prval) = prefixes
    for (key,val) = derived_units
        unit_name = parse(pr*key)
        unit = Unit(val.d, val.f*prval)
        @eval const $unit_name = $unit
    end
end

# Create a unit object from a string that is parsed with these constants
Unit(s::ASCIIString) = s == "" ? none : eval(parse(s))
Dimension(u::Unit) = Dimension(u.d.K, u.d.s, u.d.m, u.d.g, u.d.cd, u.d.mol, u.d.A)

######################################
###  Some convenient conversions  ####
######################################
# Facilitate translating dicts into ordered dicts
function convert{T <: Any}(::Type{OrderedDict{String, T}}, x::Dict{String, T})
    out =  OrderedDict{String, T}()
    for (key,val) in x
        out[key] = val
    end
end

function convert{T <: Any}(::Type{OrderedDict{String, T}}, x::Nothing)
    return OrderedDict{String, T}()
end

##################################################
###  Datatypes for type description of model  ####
##################################################

# Basic datatypes
type Parameter
    Value::Float64
    Units::Unit
end
==(p1::Parameter, p2::Parameter) = p1.Value == p2.Value &&
                                   p1.Units == p2.Units

type Forcing
    Time::Array{Float64, 1}
    Value::Array{Float64, 1}
    Units::Unit
end
==(f1::Forcing, f2::Forcing) = f1.Time == f2.Time &&
                               f1.Value == f2.Value &&
                               f1.Units == f2.Units

type Variable
    Expr::ASCIIString
    Units::Unit
end
==(v1::Variable, v2::Variable) = v1.Expr == v2.Expr &&
                                 v1.Units == v2.Units

type Derivative
    Expr::ASCIIString
    State::ASCIIString
    Units::Unit
end
==(d1::Derivative, d2::Derivative) = d1.Expr == d2.Expr &&
                                     d1.State == d2.State &&
                                     d1.Units == d2.Units

type Observed
    Expr::ASCIIString
    Units::Unit
end
==(o1::Observed, o2::Observed) = o1.Expr == o2.Expr &&
                                 o1.Units == o2.Units

type State
    Value::Float64
    Units::Unit
    QSS::Bool
end
==(s1::State, s2::State) = s1.Value == s2.Value &&
                           s1.Units == s2.Units &&
                           s1.QSS == s2.QSS

type Reactant
    Name::ASCIIString
    Stoichiometry::Float64
end
==(r1::Reactant, r2::Reactant) = r1.Name == r2.Name &&
                                 r1.Stoichiometry == r2.Stoichiometry

type Component
    Name::ASCIIString
    Species::ASCIIString
    Values::Array{Float64,1}
    Forms::Array{ASCIIString,1}
end
==(c1::Component, c2::Component) = c1.Species == c2.Species &&
                                   c1.Values == c2.Values &&
                                   c1.Forms == c2.Forms

type Species
    Value::Float64
    Compartment::ASCIIString
    Units::Unit
    Components::Vector{Component}
    QSS::Bool
end
==(s1::Species, s2::Species) = s1.Value == s2.Value &&
                               s1.Compartment == s2.Compartment &&
                               s1.Units == s2.Units &&
                               s1.Components == s2.Components &&
                               s1.QSS == s2.QSS

type Reaction
    Substrates::Array{Reactant, 1}
    Products::Array{Reactant, 1}
    Expr::ASCIIString
    Compartment::ASCIIString
    Units::Unit
end
==(r1::Reaction, r2::Reaction) = r1.Substrates == r2.Substrates &&
                                 r1.Products == r2.Products &&
                                 r1.Expr == r2.Expr &&
                                 r1.Compartment == r2.Compartment &&
                                 r1.Units == r2.Units

type Equation
    Expr::Union(Symbol, Float64, Expr)
    Units::Unit
end

# Container for the structured description of the simulation language
type ReactionSource
    DynamicType::ASCIIString
    Constants::OrderedDict{ASCIIString, Parameter}
    Parameters::OrderedDict{ASCIIString, Parameter}
    Forcings::OrderedDict{ASCIIString, Forcing}
    Species::OrderedDict{ASCIIString, Species}
    Reactions::OrderedDict{ASCIIString, Reaction}
    States::OrderedDict{ASCIIString, State}
    Derivatives::OrderedDict{ASCIIString, Derivative}
    Variables::OrderedDict{ASCIIString, Variable}
    Observed::OrderedDict{ASCIIString, Observed}
end
type OdeSource
    DynamicType::ASCIIString
    Constants::OrderedDict{ASCIIString, Parameter}
    Parameters::OrderedDict{ASCIIString, Parameter}
    Forcings::OrderedDict{ASCIIString, Forcing}
    States::OrderedDict{ASCIIString, State}
    Derivatives::OrderedDict{ASCIIString, Derivative}
    Variables::OrderedDict{ASCIIString, Variable}
    Observed::OrderedDict{ASCIIString, Observed}
end
# Container where the different equations are parsed and stored in a single dict (we still conserve their meaning)
type OdeModel
    DynamicType::ASCIIString
    Constants::OrderedDict{ASCIIString, Parameter}
    Parameters::OrderedDict{ASCIIString, Parameter}
    Forcings::OrderedDict{ASCIIString, Forcing}
    States::OrderedDict{ASCIIString, State}
    Equations::OrderedDict{ASCIIString, Equation}
    NamesDerivatives::Vector{ASCIIString}
    NamesObserved::Vector{ASCIIString}
end
# Container where all rules have been transformed into an ODE system and the equations have been sorted
type OdeSorted
    DynamicType::ASCIIString
    Constants::OrderedDict{ASCIIString, Parameter}
    Parameters::OrderedDict{ASCIIString, Parameter}
    Forcings::OrderedDict{ASCIIString, Forcing}
    States::OrderedDict{ASCIIString, State}
    SortedEquations::Vector{OrderedDict{ASCIIString, Equation}}
    NamesDerivatives::Vector{ASCIIString}
    NamesObserved::Vector{ASCIIString}
end

function deepcopy(model::OdeSorted)

  new_eqns = OrderedDict{ASCIIString, dt.Equation}[]
  for i in 1:length(sorted_model.SortedEquations)
    push!(new_eqns, OrderedDict{ASCIIString, dt.Equation}())
    for (key,val) in sorted_model.SortedEquations[i]
      new_eqns[i][key] = dt.Equation(val.Expr, val.Units)
    end
  end

  OdeSorted(model.DynamicType,
            copy(model.Constants),
            copy(model.Parameters),
            copy(model.Forcings),
            copy(model.States),
            neq_eqns,
            copy(NamesDerivatives),
            copy(NamesObserved))

end

##################################################
#######  Methods to print the datatypes  ########
##################################################

function paste(sep, s...)
    output = ""
    for i in s
        if(isa(i,ASCIIString))
            output *= i * sep
            elseif isa(i,Array{ASCIIString,1})
            output *= paste(sep,i) * sep
        end
    end
    return output
end
function paste{T <: String}(sep, s::Array{T,1})
    output = ""
    for i in s
        output *= i * sep
    end
    return output
end

function show(io::IO, model::OdeSource)
  print(io,
"""
Constants:
----------
$(paste("\n", ["$(i[1]) = $(i[2])" for i in model.Constants]))
Parameters:
-----------
$(paste("\n", ["$(i[1]) = $(i[2])" for i in model.Parameters]))
Forcings:
---------
$(paste("\n", ["$(i[1]) = $(i[2])" for i in model.Forcings]))
States:
-------
$(paste("\n", ["$(i[1]) = $(i[2])" for i in model.States]))
Variables:
----------
$(paste("\n\n", ["$(i[1]) = $(parse(replace(i[2].Expr,"1.0", ""))) " for i in model.Variables]))
Observed:
----------
$(paste("\n\n", ["$(i[1]) = $(parse(replace(i[2].Expr,"1.0", ""))) " for i in model.Observed]))
Derivatives:
----------
$(paste("\n\n", ["$(i[1]) = $(parse(replace(i[2].Expr,"1.0", ""))) " for i in model.Derivatives]))
""")
end

function show(io::IO, model::OdeSorted)
  print(io,
"""
Constants:
----------
$(paste("\n", ["$(i[1]) = $(i[2])" for i in model.Constants]))
Parameters:
-----------
$(paste("\n", ["$(i[1]) = $(i[2])" for i in model.Parameters]))
Forcings:
---------
$(paste("\n", ["$(i[1]) = $(i[2])" for i in model.Forcings]))
States:
-------
$(paste("\n", ["$(i[1]) = $(i[2])" for i in model.States]))
Equations:
----------
$(prod([paste("\n\n", "$(j[2]) ($(j[2].Units))") for j in model.SortedEquations[1]]))
$(prod([paste("\n\n", ["$(i[1]) ($(i[2].Units.d)) = $(i[2])" for i in j]) for j in model.SortedEquations[2:end]]))
""")
end
show(io::IO, par::Parameter) = print(io,"$(par.Value*par.Units)")
show(io::IO, eq::Variable) = print(io,"$(string(eq.Expr)) $(string(eq.Units))")
show(io::IO, eq::Derivative) = print(io,"$(string(eq.Expr))")
show(io::IO, eq::Observed) = print(io,"$(string(eq.Expr))")
show(io::IO, f::Forcing) = print(io,"$(f.Value*f.Units.f) $(f.Units.d) at $(f.Time)")

end


text:
println([reduce(op, 1:5) for op in [+, -, *]])
println([foldl(op, 1:5) for op in [+, -, *]])
println([foldr(op, 1:5) for op in [+, -, *]])

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using Plots; pyplot()

movies = ["Annie Hall", "Ben-Hur", "Casablanca", "Gandhi", "West S1de Story"]
num_oscars = [5, 11, 3, 8, 10]

years = ["pf-5", "pf-6", "pf-7", "pf-8", "pf-9"]
gdp = [8014, 7957, 7929, 8041, 8048]

plot(years, gdp, seriestype= :bar)
plot!(ylabel = "Counts out of 40026", ylims = (0,8500),legend=false )
plot!(title= "Alogrithms based on xid" )


text:
#=

This file contains the entry point for running the RT simulation, rt_run.

There are two implementations: one that accepts the raw parameters, and one that accepts
the model. The latter should generally be used by users.

=#

"""
    $(FUNCTIONNAME)(pol_type, obs_geom::ObsGeometry, τ_rayl, τ_aer, quad_points::QuadPoints, max_m, aerosol_optics, greek_rayleigh, τ_abs, brdf, architecture::AbstractArchitecture)

Perform Radiative Transfer calculations using given parameters

"""
function rt_run(pol_type::AbstractPolarizationType,   # Polarization type (IQUV)
                obs_geom::ObsGeometry,                # Solar Zenith, Viewing Zenith, Viewing Azimuthal
                τ_rayl,                               # Rayleigh optical depth
                τ_aer,                                # Aerosol optical depth and single-scattering albedo
                quad_points::QuadPoints,              # Quadrature points and weights
                max_m,                                # Max Fourier terms
                aerosol_optics,                       # AerosolOptics (greek_coefs, ω̃, k, fᵗ)
                greek_rayleigh::GreekCoefs,           # Greek coefficients of Rayleigh Phase Function
                τ_abs,                                # nSpec x Nz matrix of absorption
                brdf,                                 # BRDF surface type
                architecture::AbstractArchitecture)   # Whether to use CPU / GPU

    @unpack obs_alt, sza, vza, vaz = obs_geom   # Observational geometry properties
    @unpack qp_μ, wt_μ, qp_μN, wt_μN, iμ₀Nstart, μ₀, iμ₀,Nquad = quad_points # All quadrature points
    FT = eltype(sza)                    # Get the float-type to use
    Nz = length(τ_rayl)                 # Number of vertical slices
    nSpec = size(τ_abs, 1)              # Number of spectral points
    arr_type = array_type(architecture) # Type of array to use
    SFI = true                          # SFI flag
    NquadN = Nquad * pol_type.n         # Nquad (multiplied by Stokes n)
    dims = (NquadN,NquadN)              # nxn dims
    nAer  = length(aerosol_optics)      # Number of aerosols

    # Need to check this a bit better in the future!
    FT_dual = length(τ_aer) > 0 ? typeof(τ_aer[1]) : Float64

    # Output variables: Reflected and transmitted solar irradiation at TOA and BOA respectively # Might need Dual later!!
    R = zeros(FT_dual, length(vza), pol_type.n, nSpec)
    T = zeros(FT_dual, length(vza), pol_type.n, nSpec)
    R_SFI = zeros(FT_dual, length(vza), pol_type.n, nSpec)
    T_SFI = zeros(FT_dual, length(vza), pol_type.n, nSpec)

    # Notify user of processing parameters
    msg =
    """
    Processing on: $(architecture)
    With FT: $(FT)
    Source Function Integration: $(SFI)
    Dimensions: $((NquadN, NquadN, nSpec))
    """
    @info msg

    # Create arrays
    @timeit "Creating layers" added_layer         = make_added_layer(FT_dual, arr_type, dims, nSpec)
    @timeit "Creating layers" added_layer_surface = make_added_layer(FT_dual, arr_type, dims, nSpec)
    @timeit "Creating layers" composite_layer     = make_composite_layer(FT_dual, arr_type, dims, nSpec)
    @timeit "Creating arrays" Aer𝐙⁺⁺ = arr_type(zeros(FT_dual, (dims[1], dims[2], nAer)))
    @timeit "Creating arrays" Aer𝐙⁻⁺ = similar(Aer𝐙⁺⁺)
    @timeit "Creating arrays" I_static = Diagonal(arr_type(Diagonal{FT}(ones(dims[1]))));

    println("Finished initializing arrays")

    # Loop over fourier moments
    for m = 0:max_m - 1

        println("Fourier Moment: ", m, "/", max_m-1)

        # Azimuthal weighting
        weight = m == 0 ? FT(0.5) : FT(1.0)

        # Compute Z-moments of the Rayleigh phase matrix
        # For m>=3, Rayleigh matrices will be 0, can catch with if statement if wanted
        @timeit "Z moments" Rayl𝐙⁺⁺, Rayl𝐙⁻⁺ = Scattering.compute_Z_moments(pol_type, Array(qp_μ), greek_rayleigh, m, arr_type = arr_type);

        # Need to make sure arrays are 0:
        # TBD here

        # Compute aerosol Z-matrices for all aerosols
        for i = 1:nAer
            @timeit "Z moments"  Aer𝐙⁺⁺[:,:,i], Aer𝐙⁻⁺[:,:,i] = Scattering.compute_Z_moments(pol_type, Array(qp_μ), aerosol_optics[i].greek_coefs, m, arr_type = arr_type)
        end

        @timeit "Creating arrays" τ_sum_old = arr_type(zeros(FT, nSpec)) # Suniti: declaring τ_sum to be of length nSpec

        # Loop over all layers and pre-compute all properties before performing core RT
        @timeit "Computing Layer Properties" computed_atmosphere_properties = construct_all_atm_layers(FT, nSpec, Nz, NquadN, τ_rayl, τ_aer, aerosol_optics, Rayl𝐙⁺⁺, Rayl𝐙⁻⁺, Aer𝐙⁺⁺, Aer𝐙⁻⁺, τ_abs, arr_type, qp_μ, μ₀, m)

        # Loop over vertical layers:
        @showprogress 1 "Looping over layers ..." for iz = 1:Nz  # Count from TOA to BOA

            # Construct the atmospheric layer
            # From Rayleigh and aerosol τ, ϖ, compute overall layer τ, ϖ
            computed_layer_properties = get_layer_properties(computed_atmosphere_properties, iz, arr_type)

            # Perform Core RT (doubling/elemental/interaction)
            rt_kernel!(pol_type, SFI, added_layer, composite_layer, computed_layer_properties, m, quad_points, I_static, architecture, qp_μN, iz)
        end

        # Create surface matrices:
        create_surface_layer!(brdf, added_layer, SFI, m, pol_type, quad_points, arr_type(computed_atmosphere_properties.τ_sum_all[:,end]), architecture);

        # One last interaction with surface:
        @timeit "interaction" interaction!(computed_atmosphere_properties.scattering_interfaces_all[end], SFI, composite_layer, added_layer, I_static)

        # Postprocess and weight according to vza
        postprocessing_vza!(iμ₀, pol_type, composite_layer, vza, qp_μ, m, vaz, μ₀, weight, nSpec, SFI, R, R_SFI, T, T_SFI)
    end

    # Show timing statistics
    print_timer()
    reset_timer!()

    # Return R_SFI or R, depending on the flag
    return SFI ? R_SFI : R
end

"""
    $(FUNCTIONNAME)(model::vSmartMOM_Model, i_band::Integer = -1)

Perform Radiative Transfer calculations using parameters passed in through the
vSmartMOM_Model struct

"""
function rt_run(model::vSmartMOM_Model; i_band::Integer = -1)

    # Number of bands total
    n_bands = length(model.params.spec_bands)

    # Check that i_band is valid
    @assert (i_band == -1 || i_band in collect(1:n_bands)) "i_band is $(i_band) but there are only $(n_bands) bands"

    # User wants a specific band
    if i_band != -1
        return rt_run(model.params.polarization_type,
                      model.obs_geom::ObsGeometry,
                      model.τ_rayl[i_band],
                      model.τ_aer[i_band],
                      model.quad_points,
                      model.params.max_m,
                      model.aerosol_optics[i_band],
                      model.greek_rayleigh,
                      model.τ_abs[i_band],
                      model.params.brdf[i_band],
                      model.params.architecture)

    # User doesn't specify band, but there's only one
    elseif n_bands == 1

        return rt_run(model.params.polarization_type,
                      model.obs_geom::ObsGeometry,
                      model.τ_rayl[1],
                      model.τ_aer[1],
                      model.quad_points,
                      model.params.max_m,
                      model.aerosol_optics[1],
                      model.greek_rayleigh,
                      model.τ_abs[1],
                      model.params.brdf[1],
                      model.params.architecture)

    # User wants all bands
    else

        Rs = []

        for i in 1:n_bands

            println("------------------------------")
            println("Computing R for band #$(i)")
            println("------------------------------")

            R = rt_run(model.params.polarization_type,
                       model.obs_geom::ObsGeometry,
                       model.τ_rayl[i],
                       model.τ_aer[i],
                       model.quad_points,
                       model.params.max_m,
                       model.aerosol_optics[i],
                       model.greek_rayleigh,
                       model.τ_abs[i],
                       model.params.brdf[i],
                       model.params.architecture)
            push!(Rs, R);
        end

        return Rs
    end

end


text:
function filterdelay(fobj::Vector)
    return (length(fobj) - 1) ÷ 2
end

function default_fir_filterorder(responsetype::FilterType, samplingrate::Number)
    # filter settings are the same as firfilt eeglab plugin (Andreas Widmann) and MNE Python.
    # filter order is set to 3.3 times the reciprocal of the shortest transition band
    # transition band is set to either
    # min(max(l_freq * 0.25, 2), l_freq)
    # or
    # min(max(h_freq * 0.25, 2.), nyquist - h_freq)
    #
    # That is, 0.25 times the frequency, but maximally 2Hz
    %

    transwidthratio = 0.25 # magic number from firfilt eeglab plugin
    fNyquist = samplingrate ./ 2
    cutOff = responsetype.w * samplingrate
    # what is the maximal filter width we can have
    if typeof(responsetype) <: Highpass
        maxDf = cutOff

        df = minimum([maximum([maxDf * transwidthratio, 2]), maxDf])

    elseif typeof(responsetype) <: Lowpass
        #for lowpass we have to look back from nyquist
        maxDf = fNyquist - cutOff
        df = minimum([maximum([cutOff * transwidthratio, 2]), maxDf])

    end

    filterorder = 3.3 ./ (df ./ samplingrate)
    filterorder = Int(filterorder ÷ 2 * 2) # we need even filter order
    return filterorder
end

#######################################
#
# Filter compensation
#
#######################################

function compensate_for_filter(d::Dict, spectrum::AbstractArray, fs::Real)
    frequencies = range(0, stop = 1, length = Int(size(spectrum, 1))) * fs / 2

    key_name = "filter"
    key_numb = 1
    key = string(key_name, key_numb)

    while haskey(d, key)

        spectrum = compensate_for_filter(d[key], spectrum, frequencies, fs)

        @debug("Accounted for $key response in spectrum estimation")

        key_numb += 1
        key = string(key_name, key_numb)
    end

    return spectrum
end

"""
    compensate_for_filter(filter::FilterCoefficients, spectrum::AbstractArray, frequencies::AbstractArray, fs::Real)

Recover the spectrum of signal by compensating for filtering done.

# Arguments

* `filter`: The filter used on the spectrum
* `spectrum`: Spectrum of signal
* `frequencies`: Array of frequencies you want to apply the compensation to
* `fs`: Sampling rate

# Returns

Spectrum of the signal after comensating for the filter
"""
function compensate_for_filter(
    filter::FilterCoefficients,
    spectrum::AbstractArray,
    frequencies::AbstractArray,
    fs::Real,
)
    filter_response = [freqresp(filter, f * ((2pi) / fs)) for f in frequencies]

    for f = 1:length(filter_response)
        spectrum[f, :, :] = spectrum[f, :, :] ./ abs.(filter_response[f])^2
    end

    return spectrum
end


text:
using Konig
using Test


text:
using Documenter, ModiaMath

makedocs(
  modules  = [ModiaMath],
  format   = :html,
  sitename = "ModiaMath",
  authors  = "Martin Otter (DLR-SR)",
  pages    = [
     "Home"   => "index.md",
     "Manual" => [
        "man/Overview.md",
        "man/Plans.md"
        ],
     "Library" => [
        "lib/SimulationEngine.md"
        "lib/Result.md"
        "lib/DAE.md"
        "lib/Variables.md"
        "lib/Logging.md"
        "lib/Frames.md"
        "lib/Utilities.md"
        "lib/ModiaMath.md"
     ]
  ]
)


text:
module TextSearch
    import Base: broadcastable
    import StatsBase: fit, predict
    include("textconfig.jl")
    include("bow.jl")
    include("svec.jl")
    include("io.jl")
    include("basicmodels.jl")
    include("distmodel.jl")
    include("entmodel.jl")
    include("invindex.jl")
    include("rocchio.jl")
    include("neardup.jl")
end


text:
if basename(pwd()) == "aoc"
    cd("2015/23")
end

mutable struct Registers
    a::BigInt
    b::BigInt
end
Registers() = Registers(0, 0)

function interpret(filename, rs = Registers())
    program = readlines(filename)
    pc = 1
    while checkbounds(Bool, program, pc)
        line = program[pc]
        if startswith(line, "j")
            if startswith(line, "jmp")
                pc += parse(Int, line[4:end])
            else
                cond = startswith(line, "jie") ? iseven : isone
                x = parse(Int, line[7:end])
                pc += cond(getfield(rs, Symbol(line[5]))) ? x : 1
            end
        else
            r = Symbol(line[5])
            if startswith(line, "hlf")
                setfield!(rs, r, getfield(rs, r) ÷ 2)
            elseif startswith(line, "tpl")
                setfield!(rs, r, getfield(rs, r) * 3)
            elseif startswith(line, "inc")
                setfield!(rs, r, getfield(rs, r) + 1)
            end
            pc += 1
        end
    end
    rs
end

@assert interpret("example.txt").a == 2
interpret("input.txt").b
interpret("input.txt", Registers(1, 0)).b

max_stars_repo_name: stellartux/advent-of-code

max_issues_repo_name: stellartux/advent-of-code

max_forks_repo_name: stellartux/advent-of-code


text:
# plot probability distribution of sensor data
module ProbDist
    using DataFrames, CSV, Plots, FreqTables, NamedArrays
    pyplot()

    function main(is_test=false)
        data_path = joinpath(split(@__FILE__, "src")[1], "data/sensor_data_200.txt")
        df_200_mm = CSV.read(data_path, DataFrame,
                             header=["date", "time", "ir", "lidar"],
                             delim=' ')

        tbl = freqtable(df_200_mm.lidar)
        freqs = [value for(name, value) in enumerate(tbl)]
        observs = [Int(obsv) for obsv in names(tbl, 1)]
        probs = freqs ./ length(df_200_mm.lidar)

        println("FreqTable: $(tbl)")
        println("Frequency: $(freqs)")
        println("Observations: $(observs)")
        println("Probability: $(probs)")

        bar(probs, xticks=observs, label="prob dist")

        if is_test == false
            save_path = joinpath(split(@__FILE__, "src")[1], "src/prob_stats/prob_dist/prob_dist.png")
            savefig(save_path)
        end
    end
end

max_stars_repo_name: ShisatoYano/JuliaAutonomy

max_issues_repo_name: ShisatoYano/JuliaAutonomy

max_forks_repo_name: ShisatoYano/JuliaAutonomy


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule utf8proc_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("utf8proc")
JLLWrappers.@generate_main_file("utf8proc", UUID("00992c89-a35c-5347-9984-e6609dacc59a"))
end  # module utf8proc_jll


text:
function collect_nodes( subsys )  #=

collects independent nodes in a given circuit. It is assumed that the information of the circuit has been stored in the variable "subsys" given as the input argument of this function. This function will return the collected independent nodes as the array "nodes" containing the node names.

Author: Taku Noda
Started on Nov. 11, 2019
Update on June 29, 2021

Copyright 2021 Central Research Institute of Electric Power Industry

Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

=#

nodes = String[]

gnd_exists = false
for elm in eachelement( subsys )
  if elm.name == "branch"
    nodes_str = ""
    if !haskey( elm, "nph" )  # if single-phase,
      nodes_str = ( elm[ "lhs" ] )*';'*( elm[ "rhs" ] )
    else  # if multiphase,
      for elm1 in eachelement( elm )
        if elm1.name == "lhs"
          nodes_str = ( elm1.content )*';'
          break
        end
      end
      for elm1 in eachelement( elm )
        if elm1.name == "rhs"
          nodes_str = nodes_str*( elm1.content )
          break
        end
      end
    end
    nodes1 = split( nodes_str, ";" )
    for node1 in nodes1
      if ( !gnd_exists )&&( node1 == "GND" )
        gnd_exists = true
      elseif node1 != "GND"
        already_exists = false
        for node in nodes
          if node == node1
            already_exists = true
            break
          end
        end
        if !already_exists
          push!( nodes, node1 )
        end
      end
    end
  end
end

if !gnd_exists
  omegasim_error( "In a circuit description, at least one node must be named \"GND\" that is used as the zero-voltage reference. However, the given circuit has no node named \"GND\"." )
end

return nodes

end

max_stars_repo_name: takunoda/omegasim

max_issues_repo_name: takunoda/omegasim

max_forks_repo_name: takunoda/omegasim


text:
# This file was generated, do not modify it. # hide
for model in models(matching(X, y))
       print("Model Name: " , model.name , " , Package: " , model.package_name , "\n")
end


text:
abstract ChromoType
type Recurrent   <: ChromoType end
type FeedForward <: ChromoType end

type Chromosome
    id::Int64
    inputCnt::Int64
    outputCnt::Int64
    node_gene_type::ChromoType
    conn_gene_type::Symbol
    connection_genes::Dict{(Int64,Int64),ConnectionGene}
    node_genes::Vector{NodeGene}
    fitness::Float64
    species_id::Int64
    parent1_id::Int64
    parent2_id::Int64
    node_order::Vector{Int64}
    # a chromosome for general recurrent neural networks.
    function Chromosome(g::Global, parent1_id::Int64, parent2_id::Int64, node_gene_type::ChromoType, conn_gene_type::Symbol)

        new(incChromeId!(g),
            g.cf.input_nodes, g.cf.output_nodes,

            # the type of NodeGene and ConnectionGene the chromosome carries
            node_gene_type, conn_gene_type,

            # how many genes of the previous type the chromosome has
            Dict{(Int64,Int64),ConnectionGene}(), # dictionary of connection genes
            [], # empty array of node_genes
            0., # stub for fitness function
            0, # species_id

            # parents ids: helps in tracking chromosome's genealogy
            parent1_id, parent2_id,
            [] # node order (only used by feedforward networks)
            )
    end
end

function incChromeId!(g::Global)
    g.chromosomeCnt += 1
    return g.chromosomeCnt
end

function mutate(ch::Chromosome, g::Global)
    # Mutates the chromosome
    if rand() < g.cf.prob_addnode
        mutate_add_node!(ch, g, ch.node_gene_type)
    elseif rand() < g.cf.prob_addconn
        mutate_add_connection!(ch, g, ch.node_gene_type)
    else
        map(cg -> mutate!(cg[2], g.cf), ch.connection_genes) # mutate weights
        map(ng -> mutate!(ng, g.cf), ch.node_genes[ch.inputCnt+1:end]) # mutate bias, response, and etc...
    end
    return ch
end

function crossover(g::Global, self::Chromosome, other::Chromosome)
    # Crosses over parents' chromosomes and returns a child

    # This can't happen! Parents must belong to the same species.
    @assert self.species_id == other.species_id

    # TODO: if they're of equal fitnesses, choose the shortest
    parent1, parent2 = self.fitness >= other.fitness? (self,other):(other,self)

    # create a new child
    child = Chromosome(g, self.id, other.id, self.node_gene_type, self.conn_gene_type)
    inherit_genes!(g, child ,parent1, parent2)
    child.species_id = parent1.species_id
    child.fitness = parent1.fitness
    #child._input_nodes = parent1._input_nodes

    return child
end

function inherit_genes!(g::Global, child::Chromosome, parent1::Chromosome, parent2::Chromosome)
    # Applies the crossover operator.
    @assert parent1.fitness >= parent2.fitness

    # Crossover connection genes
    for (key, cg1) in parent1.connection_genes
        if haskey(parent2.connection_genes, key)
            cg2 = parent2.connection_genes[cg1.key]
            gene = is_same_innov(cg1, cg2)? deepcopy(get_child(cg1, cg2)) : deepcopy(cg1)
            child.connection_genes[cg1.key] = gene
        else # Copy excess or disjoint genes from the fittest parent
            child.connection_genes[cg1.key] = deepcopy(cg1)
        end
    end

    # Crossover node genes
    for i = 1:length(parent1.node_genes)
        ng1 = parent1.node_genes[i]
        # matching node genes: randomly selects the neuron's bias and response
        if length(parent2.node_genes) >= i && parent2.node_genes[i].id == ng1.id
            push!(child.node_genes, get_child(ng1, parent2.node_genes[i]))
        else
            # copies extra genes from the fittest parent
            push!(child.node_genes, deepcopy(ng1))
        end
    end

    child.node_order = parent1.node_order[:]
end

#----------  Recurrent Mutation  ----------

function mutate_add_node!(ch::Chromosome, g::Global,::Recurrent)
    # Choose a random connection to split
    ks = collect(keys(ch.connection_genes))
    toSpilt = rand(1:length(ks))
    conn_to_split = ch.connection_genes[ks[toSpilt]]

    ng = NodeGene(length(ch.node_genes)+1,:HIDDEN, 0., 1., g.cf.nn_activation, 1.0)
#     println("To Split: $toSpilt \n   $conn_to_split\n   $ng")
    push!(ch.node_genes, ng)
    new_conn1, new_conn2 = split(g, conn_to_split, ng.id)
    ch.connection_genes[new_conn1.key] = new_conn1
    ch.connection_genes[new_conn2.key] = new_conn2
    return (ng, conn_to_split) # the return is only used in genome_feedforward
end

function mutate_add_connection!(ch::Chromosome, g::Global, ::Recurrent)
    # Only for recurrent networks
    total_possible = (length(ch.node_genes) - ch.inputCnt)  * length(ch.node_genes)
    remaining_conns = total_possible - length(ch.connection_genes)

    # Check if new connection can be added:
    if remaining_conns > 0
        n = rand(1:remaining_conns)
        count = 1
        # Count connections
        for in_node in ch.node_genes
            for out_node in ch.node_genes[ch.inputCnt+1:end]
                if !haskey(ch.connection_genes,(in_node.id, out_node.id)) # if fDree connection
                    if count == n # Connection to create
                        weight = randn() * g.cf.weight_stdev
                        cg = ConnectionGene(g, in_node.id, out_node.id, weight, true)
                        ch.connection_genes[cg.key] = cg
#                         println(cg)
                        return
                    end
                    count += 1
                end
            end
        end
    end
end

#----------  FeedForward Mutation  ----------

function mutate_add_node!(ch::Chromosome, g::Global,::FeedForward)

    ng, split_conn = mutate_add_node!(ch::Chromosome, g::Global, Recurrent())

    if length(ch.node_order) == 0
        push!(ch.node_order, ng.id)
        return ng, split_conn
    end
    # Add node to node order list: after the presynaptic node of the split connection
    # and before the postsynaptic node of the split connection
    mini = ch.node_genes[split_conn.inId].ntype == :HIDDEN?
        findfirst(ch.node_order, split_conn.inId)+1:1

    maxi = ch.node_genes[split_conn.outId].ntype == :HIDDEN?
        findfirst(ch.node_order, split_conn.outId):length(ch.node_order)

    idx = mini <= maxi? rand(mini:maxi):mini # unnecessary?
    insert!(ch.node_order, idx, ng.id)
#     assert(length(ch.node_order) == length([n for n in ch.node_genes if n.type == :HIDDEN]))
    return ng, split_conn

end

function mutate_add_connection!(ch::Chromosome, g::Global,::FeedForward)

    # Only for feedforwad networks
    num_hidden = length(ch.node_order)
    num_output = length(ch.node_genes) - ch.inputCnt - num_hidden

    total_possible_conns = (num_hidden + num_output)*(ch.inputCnt + num_hidden) - sum([1:num_hidden+1])

    remaining_conns = total_possible_conns - length(ch.connection_genes)
    # Check if new connection can be added:
    if remaining_conns > 0
        n = rand(0:remaining_conns - 1)
        count = 0
        # Count connections
        for in_node in vcat(ch.node_genes[1:ch.inputCnt],ch.node_genes[end-num_hidden+1:end])
            for out_node in ch.node_genes[ch.inputCnt+1:end]
                if !haskey(ch.connection_genes,(in_node.id, out_node.id)) &&　is_connection_feedforward(ch, in_node, out_node)
                    # Free connection
                    if count == n # Connection to create
                        weight = randn() * g.cf.weight_stdev
                        cg = ConnectionGene(g, in_node.id, out_node.id, weight, true)
                        ch.connection_genes[cg.key] = cg
                        return
                    else
                        count += 1
                    end
                end
            end
        end
    end

end

function is_connection_feedforward(ch::Chromosome, in_node::NodeGene, out_node::NodeGene)
    return in_node.ntype == :INPUT || out_node.ntype == :OUTPUT ||
        findfirst(ch.node_order, in_node.id) < findfirst(ch.node_order, out_node.id)
end

#----------  End of Mutation  ----------

# compatibility function
function distance(self::Chromosome, other::Chromosome, cf::Config)

    # Returns the distance between this chromosome and the other.
    chromo1, chromo2 = length(self.connection_genes) > length(other.connection_genes)? (self,other):(other,self)

    weight_diff = 0
    matching = 0
    disjoint = 0
    excess = 0

    max_cg_chromo2 = maxInnov(chromo2.connection_genes)

    for (k, cg1) in chromo1.connection_genes
        if haskey(chromo2.connection_genes, k)
            # Homologous genes
            cg2 = chromo2.connection_genes[cg1.key]
            weight_diff += abs(cg1.weight - cg2.weight)
            matching += 1
        else
            if cg1.innovNumber > max_cg_chromo2.innovNumber
                excess += 1
            else
                disjoint += 1
            end
        end
    end

    disjoint += length(chromo2.connection_genes) - matching
    d = cf.excess_coeficient * excess + cf.disjoint_coeficient * disjoint

    return matching > 0? d + cf.weight_coeficient * weight_diff / matching : d
end

function size(ch::Chromosome)
    # Defines chromosome 'complexity': number of hidden nodes plus
    # number of enabled connections (bias is not considered)
    num_hidden = length(ch.node_genes) - ch.inputCnt - ch.outputCnt
    conns_enabled = sum(map(cg->ch.connection_genes[cg].enable==true? 1:0, collect(keys(ch.connection_genes))))

    return num_hidden, conns_enabled
end

#     def __cmp__(self, other)
#         """ First compare chromosomes by their fitness and then by their id.
#             Older chromosomes (lower ids) should be prefered if newer ones
#             performs the same.
#         """
#         #return cmp(self.fitness, other.fitness) or cmp(other.id, self.id)
#         return cmp(self.fitness, other.fitness)

function Base.show(io::IO, ch::Chromosome)
    s = "$(ch.node_gene_type) type\nNodes:\n"
    for ng in ch.node_genes
        s = "$(s)\t$(ng)"
    end
    s = "$(s)Order:\n\t$(ch.node_order)"

    s = "$(s)\nConnections:"
    ord = map(x->(ch.connection_genes[x].innovNumber,x), collect(keys(ch.connection_genes)))
    sort!(ord, by=x->x[1])
    for (innov, (inId ,outId)) in ord
        s = "$(s)\n\t$(ch.connection_genes[(inId ,outId)])"
    end
    @printf(io,"\n%6s", s)
    return
end

function add_hidden_nodes!(g::Global, ch::Chromosome, num_hidden::Int64, ::Recurrent)

    id = length(ch.node_genes)+1
    for i in 1:num_hidden
        node_gene = NodeGene(id, :HIDDEN, 0., 1., g.cf.nn_activation)
        push!(ch.node_genes, node_gene)
        id += 1
        # Connect all nodes to it
        for pre in ch.node_genes
            weight = randn() * g.cf.weight_stdev
            cg = ConnectionGene(g, pre.id, node_gene.id, weight)
            ch.connection_genes[cg.key] = cg
        end
        # Connect it to all nodes except input nodes
        for post in ch.node_genes[ch.inputCnt+1:end]
            weight = randn() * g.cf.weight_stdev
            cg = ConnectionGene(g, node_gene.id, post.id, weight)
            ch.connection_genes[cg.key] = cg
        end
    end
end

function create_unconnected(g::Global)

    # Creates a chromosome for an unconnected feedforward network with no hidden nodes.
    c = Chromosome(g, 0, 0, (g.cf.feedforward? FeedForward():Recurrent()), :ConnectionGene)

    id = 1
    # Create node genes
    for i = 1:c.inputCnt
        push!(c.node_genes, NodeGene(id, :INPUT, 0., 1., :none))
        id += 1
    end
#         #c._input_nodes += num_input
    for i in 1:c.outputCnt
        push!(c.node_genes, NodeGene(id, :OUTPUT, 0., 1., g.cf.nn_activation))
        id += 1
    end
    @assert id == length(c.node_genes) + 1
    return c
end

function create_minimally_connected(g::Global)

    # Creates a chromosome for a minimally connected feedforward network with no hidden nodes.
    # That is, each output node will have a single connection from a randomly chosen input node.
    ch = create_unconnected(g)
    for node_gene in ch.node_genes[end-ch.outputCnt+1:end] # each output node
        # Connect it to a random input node
        input_node = ch.node_genes[rand(1:ch.inputCnt)]
        weight = randn() * g.cf.weight_stdev
        cg = ConnectionGene(g, input_node.id, node_gene.id, weight)
        ch.connection_genes[cg.key] = cg
    end
    return ch
end

function create_fully_connected(g::Global)

    # Creates a chromosome for a fully connected feedforward network with no hidden nodes.
    ch = create_unconnected(g)
    for node_gene in ch.node_genes[end-ch.outputCnt+1:end] # each output node
        # Connect it to all input nodes
        for input_node in ch.node_genes[1:ch.inputCnt]
            # TODO: review the initial weights distribution
            # weight = random.uniform(-1, 1)*Config.random_range
            weight = randn() * g.cf.weight_stdev
            cg = ConnectionGene(g, input_node.id, node_gene.id, weight)
            ch.connection_genes[cg.key] = cg
        end
    end
    return ch
end

function add_hidden_nodes!(g::Global, ch::Chromosome, num_hidden::Int64, ::FeedForward)
    id = length(ch.node_genes)+1
    for i in 1:num_hidden
        node_gene = NodeGene(id, :HIDDEN, 0., 1., g.cf.nn_activation)
        push!(ch.node_genes, node_gene)
        push!(ch.node_order,node_gene.id)
        id += 1

        # Connect all input nodes to it
        for pre in ch.node_genes[1:ch.inputCnt]
            weight = randn() * g.cf.weight_stdev
            cg = ConnectionGene(g, pre.id, node_gene.id, weight, true)
            ch.connection_genes[cg.key] = cg
            @assert is_connection_feedforward(ch, pre, node_gene) == true
        end

        # Connect all previous hidden nodes to it
        # Comment: Makes for one wierd network! Omit hidden->hidden??
#         for pre_id in ch.node_order[1:end-1]
#             @assert pre_id != node_gene.id
#             weight = randn() * g.cf.weight_stdev
#             cg = ConnectionGene(g, pre_id, node_gene.id, weight, true)
#             ch.connection_genes[cg.key] = cg
#         end

        # Connect it to all output nodes
        for post in ch.node_genes[ch.inputCnt+1:(ch.inputCnt + ch.outputCnt)]
            @assert post.ntype == :OUTPUT
            weight = randn() * g.cf.weight_stdev
            cg = ConnectionGene(g, node_gene.id, post.id, weight, true)
            ch.connection_genes[cg.key] = cg
#             assert self.__is_connection_feedforward(node_gene, post)
        end
    end
end

# if __name__ == '__main__':
#     # Example
#     import visualize
#     # define some attributes
#     node_gene_type = genome.NodeGene         # standard neuron model
#     conn_gene_type = genome.ConnectionGene   # and connection link
#     Config.nn_activation = 'exp'             # activation function
#     Config.weight_stdev = 0.9                # weights distribution

#     Config.input_nodes = 2                   # number of inputs
#     Config.output_nodes = 1                  # number of outputs

#     # creates a chromosome for recurrent networks
#     #c1 = Chromosome.create_fully_connected()

#     # creates a chromosome for feedforward networks
#     c2 = FFChromosome.create_fully_connected()
#     # add two hidden nodes
#     c2.add_hidden_nodes(2)
#     # apply some mutations
#     #c2._mutate_add_node()
#     #c2._mutate_add_connection()

#     # check the result
#     #visualize.draw_net(c1) # for recurrent nets
#     visualize.draw_ff(c2)   # for feedforward nets
#     # print the chromosome
#     print  c2

max_stars_repo_name: Andy-P/NEAT

max_issues_repo_name: Andy-P/NEAT

max_forks_repo_name: Andy-P/NEAT


text:
const KN_var_type_map = Dict(
    KN_VARTYPE_CONTINUOUS => :Cont,
    KN_VARTYPE_INTEGER => :Int,
    KN_VARTYPE_BINARY => :Bin
)

const KN_rev_var_type_map = Dict(
    :Cont => KN_VARTYPE_CONTINUOUS,
    :Int => KN_VARTYPE_INTEGER,
    :Bin => KN_VARTYPE_BINARY
)

# This list was obtained through AWK with Knitro 11.1 with the following command:
# grep "#define" knitro.h | grep "KN_PARAM_" | awk '{ printf("\"%s\" => Int32(%s),\n",$2,$3) }'
# For versions up to 10.3, use the folowing command:
# grep "#define" knitro.h | grep "KTR_PARAM_" | awk '{ printf("\"%s\" => Int32(%s),\n",$2,$3) }'
const KN_paramName2Indx = Dict(
    "KN_PARAM_NEWPOINT" => Int32(1001),
    "KN_PARAM_HONORBNDS" => Int32(1002),
    "KN_PARAM_ALGORITHM" => Int32(1003),
    "KN_PARAM_ALG" => Int32(1003),
    "KN_PARAM_BAR_MURULE" => Int32(1004),
    "KN_PARAM_BAR_FEASIBLE" => Int32(1006),
    "KN_PARAM_GRADOPT" => Int32(1007),
    "KN_PARAM_HESSOPT" => Int32(1008),
    "KN_PARAM_BAR_INITPT" => Int32(1009),
    "KN_PARAM_ACT_LPSOLVER" => Int32(1012),
    "KN_PARAM_CG_MAXIT" => Int32(1013),
    "KN_PARAM_MAXIT" => Int32(1014),
    "KN_PARAM_OUTLEV" => Int32(1015),
    "KN_PARAM_OUTMODE" => Int32(1016),
    "KN_PARAM_SCALE" => Int32(1017),
    "KN_PARAM_SOC" => Int32(1019),
    "KN_PARAM_DELTA" => Int32(1020),
    "KN_PARAM_BAR_FEASMODETOL" => Int32(1021),
    "KN_PARAM_FEASTOL" => Int32(1022),
    "KN_PARAM_FEASTOLABS" => Int32(1023),
    "KN_PARAM_MAXTIMECPU" => Int32(1024),
    "KN_PARAM_BAR_INITMU" => Int32(1025),
    "KN_PARAM_OBJRANGE" => Int32(1026),
    "KN_PARAM_OPTTOL" => Int32(1027),
    "KN_PARAM_OPTTOLABS" => Int32(1028),
    "KN_PARAM_LINSOLVER_PIVOTTOL" => Int32(1029),
    "KN_PARAM_XTOL" => Int32(1030),
    "KN_PARAM_DEBUG" => Int32(1031),
    "KN_PARAM_MULTISTART" => Int32(1033),
    "KN_PARAM_MSENABLE" => Int32(1033),
    "KN_PARAM_MSMAXSOLVES" => Int32(1034),
    "KN_PARAM_MSMAXBNDRANGE" => Int32(1035),
    "KN_PARAM_MSMAXTIMECPU" => Int32(1036),
    "KN_PARAM_MSMAXTIMEREAL" => Int32(1037),
    "KN_PARAM_LMSIZE" => Int32(1038),
    "KN_PARAM_BAR_MAXCROSSIT" => Int32(1039),
    "KN_PARAM_MAXTIMEREAL" => Int32(1040),
    "KN_PARAM_CG_PRECOND" => Int32(1041),
    "KN_PARAM_BLASOPTION" => Int32(1042),
    "KN_PARAM_BAR_MAXREFACTOR" => Int32(1043),
    "KN_PARAM_LINESEARCH_MAXTRIALS" => Int32(1044),
    "KN_PARAM_BLASOPTIONLIB" => Int32(1045),
    "KN_PARAM_OUTAPPEND" => Int32(1046),
    "KN_PARAM_OUTDIR" => Int32(1047),
    "KN_PARAM_CPLEXLIB" => Int32(1048),
    "KN_PARAM_BAR_PENRULE" => Int32(1049),
    "KN_PARAM_BAR_PENCONS" => Int32(1050),
    "KN_PARAM_MSNUMTOSAVE" => Int32(1051),
    "KN_PARAM_MSSAVETOL" => Int32(1052),
    "KN_PARAM_PRESOLVEDEBUG" => Int32(1053),
    "KN_PARAM_MSTERMINATE" => Int32(1054),
    "KN_PARAM_MSSTARTPTRANGE" => Int32(1055),
    "KN_PARAM_INFEASTOL" => Int32(1056),
    "KN_PARAM_LINSOLVER" => Int32(1057),
    "KN_PARAM_BAR_DIRECTINTERVAL" => Int32(1058),
    "KN_PARAM_PRESOLVE" => Int32(1059),
    "KN_PARAM_PRESOLVE_TOL" => Int32(1060),
    "KN_PARAM_BAR_SWITCHRULE" => Int32(1061),
    "KN_PARAM_HESSIAN_NO_F" => Int32(1062),
    "KN_PARAM_MA_TERMINATE" => Int32(1063),
    "KN_PARAM_MA_MAXTIMECPU" => Int32(1064),
    "KN_PARAM_MA_MAXTIMEREAL" => Int32(1065),
    "KN_PARAM_MSSEED" => Int32(1066),
    "KN_PARAM_MA_OUTSUB" => Int32(1067),
    "KN_PARAM_MS_OUTSUB" => Int32(1068),
    "KN_PARAM_XPRESSLIB" => Int32(1069),
    "KN_PARAM_TUNER" => Int32(1070),
    "KN_PARAM_TUNER_OPTIONSFILE" => Int32(1071),
    "KN_PARAM_TUNER_MAXTIMECPU" => Int32(1072),
    "KN_PARAM_TUNER_MAXTIMEREAL" => Int32(1073),
    "KN_PARAM_TUNER_OUTSUB" => Int32(1074),
    "KN_PARAM_TUNER_TERMINATE" => Int32(1075),
    "KN_PARAM_LINSOLVER_OOC" => Int32(1076),
    "KN_PARAM_BAR_RELAXCONS" => Int32(1077),
    "KN_PARAM_MSDETERMINISTIC" => Int32(1078),
    "KN_PARAM_BAR_REFINEMENT" => Int32(1079),
    "KN_PARAM_DERIVCHECK" => Int32(1080),
    "KN_PARAM_DERIVCHECK_TYPE" => Int32(1081),
    "KN_PARAM_DERIVCHECK_TOL" => Int32(1082),
    "KN_PARAM_LINSOLVER_INEXACT" => Int32(1083),
    "KN_PARAM_LINSOLVER_INEXACTTOL" => Int32(1084),
    "KN_PARAM_MAXFEVALS" => Int32(1085),
    "KN_PARAM_FSTOPVAL" => Int32(1086),
    "KN_PARAM_DATACHECK" => Int32(1087),
    "KN_PARAM_DERIVCHECK_TERMINATE" => Int32(1088),
    "KN_PARAM_BAR_WATCHDOG" => Int32(1089),
    "KN_PARAM_FTOL" => Int32(1090),
    "KN_PARAM_FTOL_ITERS" => Int32(1091),
    "KN_PARAM_ACT_QPALG" => Int32(1092),
    "KN_PARAM_BAR_INITPI_MPEC" => Int32(1093),
    "KN_PARAM_XTOL_ITERS" => Int32(1094),
    "KN_PARAM_LINESEARCH" => Int32(1095),
    "KN_PARAM_OUT_CSVINFO" => Int32(1096),
    "KN_PARAM_INITPENALTY" => Int32(1097),
    "KN_PARAM_ACT_LPFEASTOL" => Int32(1098),
    "KN_PARAM_CG_STOPTOL" => Int32(1099),
    "KN_PARAM_RESTARTS" => Int32(1100),
    "KN_PARAM_RESTARTS_MAXIT" => Int32(1101),
    "KN_PARAM_BAR_SLACKBOUNDPUSH" => Int32(1102),
    "KN_PARAM_CG_PMEM" => Int32(1103),
    "KN_PARAM_BAR_SWITCHOBJ" => Int32(1104),
    "KN_PARAM_OUTNAME" => Int32(1105),
    "KN_PARAM_OUT_CSVNAME" => Int32(1106),
    "KN_PARAM_ACT_PARAMETRIC" => Int32(1107),
    "KN_PARAM_ACT_LPDUMPMPS" => Int32(1108),
    "KN_PARAM_ACT_LPALG" => Int32(1109),
    "KN_PARAM_ACT_LPPRESOLVE" => Int32(1110),
    "KN_PARAM_ACT_LPPENALTY" => Int32(1111),
    "KN_PARAM_BNDRANGE" => Int32(1112),
    "KN_PARAM_BAR_CONIC_ENABLE" => Int32(1113),
    "KN_PARAM_CONVEX" => Int32(1114),
    "KN_PARAM_OUT_HINTS" => Int32(1115),
    "KN_PARAM_EVAL_FCGA" => Int32(1116),
    "KN_PARAM_MIP_METHOD" => Int32(2001),
    "KN_PARAM_MIP_BRANCHRULE" => Int32(2002),
    "KN_PARAM_MIP_SELECTRULE" => Int32(2003),
    "KN_PARAM_MIP_INTGAPABS" => Int32(2004),
    "KN_PARAM_MIP_INTGAPREL" => Int32(2005),
    "KN_PARAM_MIP_MAXTIMECPU" => Int32(2006),
    "KN_PARAM_MIP_MAXTIMEREAL" => Int32(2007),
    "KN_PARAM_MIP_MAXSOLVES" => Int32(2008),
    "KN_PARAM_MIP_INTEGERTOL" => Int32(2009),
    "KN_PARAM_MIP_OUTLEVEL" => Int32(2010),
    "KN_PARAM_MIP_OUTINTERVAL" => Int32(2011),
    "KN_PARAM_MIP_OUTSUB" => Int32(2012),
    "KN_PARAM_MIP_DEBUG" => Int32(2013),
    "KN_PARAM_MIP_IMPLICATNS" => Int32(2014),
    "KN_PARAM_MIP_GUB_BRANCH" => Int32(2015),
    "KN_PARAM_MIP_KNAPSACK" => Int32(2016),
    "KN_PARAM_MIP_ROUNDING" => Int32(2017),
    "KN_PARAM_MIP_ROOTALG" => Int32(2018),
    "KN_PARAM_MIP_LPALG" => Int32(2019),
    "KN_PARAM_MIP_TERMINATE" => Int32(2020),
    "KN_PARAM_MIP_MAXNODES" => Int32(2021),
    "KN_PARAM_MIP_HEURISTIC" => Int32(2022),
    "KN_PARAM_MIP_HEUR_MAXIT" => Int32(2023),
    "KN_PARAM_MIP_HEUR_MAXTIMECPU" => Int32(2024),
    "KN_PARAM_MIP_HEUR_MAXTIMEREAL" => Int32(2025),
    "KN_PARAM_MIP_PSEUDOINIT" => Int32(2026),
    "KN_PARAM_MIP_STRONG_MAXIT" => Int32(2027),
    "KN_PARAM_MIP_STRONG_CANDLIM" => Int32(2028),
    "KN_PARAM_MIP_STRONG_LEVEL" => Int32(2029),
    "KN_PARAM_MIP_INTVAR_STRATEGY" => Int32(2030),
    "KN_PARAM_MIP_RELAXABLE" => Int32(2031),
    "KN_PARAM_MIP_NODEALG" => Int32(2032),
    "KN_PARAM_MIP_HEUR_TERMINATE" => Int32(2033),
    "KN_PARAM_MIP_SELECTDIR" => Int32(2034),
    "KN_PARAM_PAR_NUMTHREADS" => Int32(3001),
    "KN_PARAM_PAR_CONCURRENT_EVALS" => Int32(3002),
    "KN_PARAM_PAR_BLASNUMTHREADS" => Int32(3003),
    "KN_PARAM_PAR_LSNUMTHREADS" => Int32(3004),
    "KN_PARAM_PAR_MSNUMTHREADS" => Int32(3005),
)


text:
#=##############################################################################
# DESCRIPTION
    Methods for geometric manipulation.
# AUTHORSHIP
  * Author    : Eduardo J Alvarez
  * Email     : Edo.AlvarezR@gmail.com
  * Created   : Nov 2017
  * License   : MIT License
=###############################################################################

################################################################################
# MESHING AND DISCRETIZATION
################################################################################
"""
  `discretize(f, xlow, xhigh, n::Int64, r::Float64; central::Bool=false)`

Discretizes the continuous function `f` between the range `xlow` and `xhigh`
into `n` intervals, with `r` the ratio between first and last interval if
`central=false` or between first and central interval if `central=true`.

  **Examples**

  The following lines show the discretization of a semi-circumference perimeter
  into 100 intervals of uniform length:

  ```julia
    julia> f(x) = (x, sqrt(1-round(x,digits=8)^2), 0)  # Semi-circunference of radius 1
    julia> discretize(f, -1, 1, 100, 1.0)
  ```
"""
function discretize(f, xlow, xhigh, n::Int64, r::Float64; central::Bool=false,
                      check::Bool=true)

  # ERROR CASES
  if n == 0
    return Any[]
  elseif n < 0
    error("Invalid number of intervals (n < 0)")
  elseif r <= 0
    error("Invalid expansion ratio (r <= 0)")
  end

  out = Any[f(xlow)]

  l = xhigh - xlow        # Length of domain to be discretized
  cumlen = 0              # Cumulative length already walked

  for i in 1:n

    # Case of uniform discretization
    if r==1.0
      len = l/n

    # Linear increment discretization (see notebook entry 20170519)
    # Case of no central expansion
    elseif !central
      p = l/( (n*(n-1)/2)*(r+1)/(r-1) )
      d1 = p*(n-1)/(r-1)
      len = d1 + p*(i-1)
      # println("i=$i\tp=$p\td1=$d1\tlen=$len")

    # Case of central expansion
    else
      _central = 0.5
      # Left of the center
      if i<=floor(n*_central)
        _l = l*_central
        _n = floor(n*_central)
        _r = r
        _i = i
      # Right of the center
      else
        _l = l*(1-_central)
        _n = n-floor(n*_central)
        _r = 1/r
        _i = i-floor(n*_central)
      end
      p = _l/( (_n*(_n-1)/2)*(_r+1)/(_r-1) )
      d1 = p*(_n-1)/(_r-1)
      len = d1 + p*(_i-1)
      # println("$i\t$_r\t$cumlen\t_n=$_n\t_i=$_i")
    end

    cumlen += len
    this_x = xlow + (cumlen/l)*(xhigh-xlow)
    # println("i=$i\tx=$this_x\tlen=$len")
    this_f = f(this_x)
    push!(out, this_f)
  end

  # Verifies correct discretization
  if check && abs((cumlen-l)/l)>0.0001
    error("Critical logic error! cumlen!=l ($cumlen!=$l)")
  end

  return out
end

"""
  `multidiscretize(f, xlow, xhigh, sections)`

Discretizes the continuous function `f` between the range `xlow` and `xhigh`
into multiple sections of refinement as specified in `sections`.

  ** Arguments **
  * `f`         : Continuous function of the form `f(x)` to be discretized
                  between `xlow` and `xhigh`,
  * `xlow`      : Lower bound.
  * `xhigh`     : Upper bound.
  * `sections`  : Array `[sec1, sec2, ...]`specifying the
                  sections of discretization in the format
                  `sec = (c::Float64, n::Int64, r::Float64, central::Bool)`,
                  with `c` the normalized length of this section (the sum of all
                  c must equal one), `n` the number of intervals in this section
                  , `r` the increment ratio between first and last interval if
                  `central=false` or between first and central interval if
                  `central=true`.

  **Examples**

  The following lines show the discretization of a semi-circumference perimeter
  into 90 intervals done in three sections of discretization:

  ```julia
    julia> f(theta) = (cos(theta), sin(theta), 0)
    julia> sec = (1/3, 30, 1/8, true)
    julia> points = multidiscretize(f, 0, pi, [sec, sec, sec])
  ```
"""
function multidiscretize(f, xlow, xhigh, sections::multidisctype;
                                                              check::Bool=true)
  out = Any[f(xlow)]
  ctot = sum([sec[1] for sec in sections]) # Sum of all `c`s

  # Iterates over sections
  prev_xhigh = xlow
  for (c,n,r,central) in sections

    this_c = c/ctot # Makes sure that the sum of `c`s will add 1
    this_xlow = prev_xhigh
    this_xhigh = prev_xhigh + this_c*(xhigh-xlow)

    this_out = discretize(f, this_xlow, this_xhigh, n, r; central=central,
                                  check=check)
    out = vcat(out, this_out[2:end])

    prev_xhigh = this_xhigh
  end

  return out
end

"""
  Receives a closed 2D contour and splits it up in upper and lower surfaces as
  divided by lowest and highest x values. It returns `(upper, lower)` with
  `upper=(x,y)` the points of the upper surface, ditto for `lower`. Both `upper`
  and `lower` are given in increasing order in x. This function is useful for
  splitting up a closed contour into two sections that are injective in x, and
  that can be received by `parameterize()`.
"""
function splitcontour(x,y)

  # Flag indicating whether the contour start at the trailing or leading edge
  start_TE = x[1]==maximum(x)

  # Find the opposite end of the contour
  end_i = -1
  for (i, xi) in enumerate(x)
    if i==1
      nothing
    # Case of starting from the trailing edge
    elseif start_TE && xi > x[i-1]
      end_i = i-1
      break
    # Case of leading edge
    elseif !start_TE  && xi < x[i-1]
      end_i = i-1
      break
    end
  end

  # ERROR CASE
  if end_i==-1
    error("Logic error! End of contour not found!")
  end

  # Splits them up
  x_sec1, y_sec1 = x[1:end_i], y[1:end_i]
  x_sec2, y_sec2 = x[end_i:end], y[end_i:end]

  # Sorts them from LE to TE
  if x_sec1[1] > minimum(x); reverse!(x_sec1); reverse!(y_sec1); end;
  if x_sec2[1] > minimum(x); reverse!(x_sec2); reverse!(y_sec2); end;

  # Determines upper and lower surfaces
  if Statistics.mean(y_sec1) > Statistics.mean(y_sec2)
    upper = [x_sec1, y_sec1]
    lower = [x_sec2, y_sec2]
  else
    upper = [x_sec2, y_sec2]
    lower = [x_sec1, y_sec1]
  end

  return upper, lower
end

##### END OF MESHING ###########################################################

################################################################################
# PARAMETRIZATION
################################################################################
"""
  `parameterize(x, y, z; inj_var::Int64=1, s=0.0001, debug=false)`

Receives a contour (line) and returns a parametrization function f(s) of the
contour. The parametrization is done on the path of the splined
contour such that f(0.5) returns the point (x,y,z) where half of the entire
contour has been walked, and f(1.0) returns the last point in the contour. To
perform this parametrization, the contour must be injective in least one of the
variables (x, y, or z) (i.e., all values of the variable are unique), and that
variable must be given in increasing order.

  **Arguments**
  * `x::Array{Float64,1}`    : x-coordinates of the contour.
  * `y::Array{Float64,1}`    : y-coordinates of the contour.
  * `z::Array{Float64,1}`    : z-coordinates of the contour.

  **Optional Arguments**
  * `inj_var::Int64`         : Indicates the variable that is injective, with
                                  1=x, 2=y, and 3=z.
  * `s::Float64`             : Spline smoothness.

See Parametrization section in documentation for more details.
"""
function parameterize(x, y, z; inj_var::Int64=1, s=0.0001, debug=false, atol=0,
                                            kspl="automatic", bc="extrapolate")
  # ERROR CASES
  if size(x)!=size(y)!=size(z)  # Invalid contour
    error("Invalid contour. "*
              "$(size(x))!=$(size(y))!=$(size(z)) (size(x)!=size(y)!=size(z))")
  elseif !(inj_var in [1,2,3])  # Invalid input
    error("Invalid `inj_var`=$inj_var (!(inj_var in [1,2,3]))")
  end

  # Identifies injective and dependant variables
  inj = nothing   # Injective variable
  dep = []        # Dependant variables
  for (i,var) in enumerate([x, y, z])
    if i==inj_var
      inj = var
    else
      push!(dep, var)
    end
  end

  # Splines
  # bc="extrapolate"           # Out of boundary case
  # s=0.0001                 # Spline smoothness
  k = kspl=="automatic" ? min(size(x)[1]-1, 5) : kspl  # Spline order
  spl = []                   # Spline of each variable respect the injective

  for var in dep
    this_spl = Dierckx.Spline1D(inj, var; k=k, bc=bc, s=s)
    push!(spl, this_spl)
  end

  # Defines the path function
  dfdx1(x) = Dierckx.derivative(spl[1], x)    # Derivative of f respect x1
  dfdx2(x) = Dierckx.derivative(spl[2], x)    # Derivative of f respect x2
  fun(x) = sqrt.(1 .+ (dfdx1(x)).^2 .+ (dfdx2(x)).^2)   # Integrand
                                                  # Integral between xmin and x
  fun_s(this_x) = QuadGK.quadgk(fun, inj[1], this_x; atol=atol)[1]

  # Defines the normalized path function
  stot = fun_s(inj[end])      # Total length of the path
  norm_s(x) = fun_s(x)/stot         # Normalized path function

  # Defines the inverse normalized function
  function inv_norm_s(sstar; debug=false)
      # ERROR CASES
      if sstar<-0.001 || sstar>1.001 # Outside of domain
          error("Invalid normalized path length $(sstar) (sstar<0 || sstar>1)")
      end

      # Finds the x that matches the target (sstar)
      this_fun(x) = sstar - norm_s(x)   # Zeroed function
      bracket = [ inj[1]*(1-0.01*sign(inj[1])),# Bracket of the zero
                  inj[end]*(1+0.01*sign(inj[end]))]

      if debug
          println("sstar=$sstar\tbracket=$bracket")
          println("flow=$(this_fun(bracket[1]))")
          println("fup=$(this_fun(bracket[2]))")
      end

      # Zero
      this_x = Roots.find_zero(this_fun, bracket, Roots.Bisection())

      return this_x
  end

  # Redefine the original function in terms of the new parametrization s*
  function fstar(sstar)
      x = inv_norm_s(sstar; debug=debug)       # Finds the x that matches the path
      out = []                    # Output
      dep_i = 1                   # Next spline to evaluate
      for i in 1:3
        if i==inj_var
          push!(out, x)           # Injective variable
        else
          push!(out, spl[dep_i](x))   # Dependant variables
          dep_i += 1
        end
      end
      return out
  end

  return fstar
end

##### END OF PARAMETRIZATION ###################################################

################################################################################
# LINEAR ALGEBRA TRANSFORMATIONS
################################################################################
"""
  `transform(V::Array{Float64,1}, M::Array{Float64,2}, T::Array{Float64,1})`

Rotates and translates the vector V. Receives the i', j', k' unit vectors of an
euclidean system with origin T, and returns ``V'=M(V-T)`` (In this version, the
unit vectors have been organized as a matrix `M=[i'; j'; k']`).
"""
function transform(V::Array{T1,1}, M::Array{T2,2}, T::Array{T3,1}
                                        ) where{T1<:Real, T2<:Real, T3<:Real}
    return [
                M[1,1]*(V[1]-T[1]) + M[1,2]*(V[2]-T[2]) + M[1,3]*(V[3]-T[3]),
                M[2,1]*(V[1]-T[1]) + M[2,2]*(V[2]-T[2]) + M[2,3]*(V[3]-T[3]),
                M[3,1]*(V[1]-T[1]) + M[3,2]*(V[2]-T[2]) + M[3,3]*(V[3]-T[3])
            ]
end

function transform!(out::Array{T1, 1}, V::Array{T2,1},
                            M::Array{T3,2}, T::Array{T4,1}
                           ) where{T1<:Real, T2<:Real, T3<:Real, T4<:Real}
    out[1] = M[1,1]*(V[1]-T[1]) + M[1,2]*(V[2]-T[2]) + M[1,3]*(V[3]-T[3])
    out[2] = M[2,1]*(V[1]-T[1]) + M[2,2]*(V[2]-T[2]) + M[2,3]*(V[3]-T[3])
    out[3] = M[3,1]*(V[1]-T[1]) + M[3,2]*(V[2]-T[2]) + M[3,3]*(V[3]-T[3])
end

function transform(Vs::AbstractArray{Array{P,1},1}, M::AbstractArray{P,2}, T::AbstractArray{P,1}
                                                                ) where{P<:Real}
  out = Array{Float64,1}[]
  for V in Vs
    push!(out, transform(V, M, T))
  end
  return out
end

"""
  `countertransform(Vp::Array{Float64,1}, invM::Array{Float64,2},
T::Array{Float64,1})`

Rotates and translates back a vector V' that had been rotated and translated
into the system (i', j', k') with origin T, and returns the original
``V=M^{-1}V' + T``. To ease repetitive computation, instead of giving the unit
vectors, give the inverse of their matrix.
"""
function countertransform(Vp::Array{T1,1}, invM::Array{T2,2}, T::Array{T3,1}
                                        ) where{T1<:Real, T2<:Real, T3<:Real}
    return [
                invM[1,1]*Vp[1] + invM[1,2]*Vp[2] + invM[1,3]*Vp[3] + T[1],
                invM[2,1]*Vp[1] + invM[2,2]*Vp[2] + invM[2,3]*Vp[3] + T[2],
                invM[3,1]*Vp[1] + invM[3,2]*Vp[2] + invM[3,3]*Vp[3] + T[3]
            ]
end

function countertransform!(out::Array{T1, 1}, Vp::Array{T2,1},
                            invM::Array{T3,2}, T::Array{T4,1}
                           ) where{T1<:Real, T2<:Real, T3<:Real, T4<:Real}
   out[1] = invM[1,1]*Vp[1] + invM[1,2]*Vp[2] + invM[1,3]*Vp[3] + T[1]
   out[2] = invM[2,1]*Vp[1] + invM[2,2]*Vp[2] + invM[2,3]*Vp[3] + T[2]
   out[3] = invM[3,1]*Vp[1] + invM[3,2]*Vp[2] + invM[3,3]*Vp[3] + T[3]
end

function countertransform(Vps::AbstractArray{Array{P,1},1}, invM::AbstractArray{P,2},
                                                  T::AbstractArray{P,1}) where{P<:Real}
  out = Array{Float64,1}[]
  for Vp in Vps
    push!(out, countertransform(Vp, invM, T))
  end
  return out
end

"""
  `check_coord_sys(M::Array{Float64,2}; raise_error::Bool=true)`

Checks that the unit vectors given as the matrix `M=[i;j;k]` define a coordinate
system
"""
function check_coord_sys(M::Array{T,2}; raise_error::Bool=true) where{T<:Real}
  # Checks normalization
  for i in 1:size(M)[1]
    if abs(norm(M[i,:])-1) > 0.00000001
      println(M)
      if raise_error
        error("Not unitary axis: $(M[i,:])")
      else
        return false
      end
    end
  end

  # Checks ortogonality
  for i in size(M)[1]
    xi = M[i, :]
    xip1 = M[(i%size(M)[1])+1, :]
    proj = abs(dot(xi, xip1))
    if proj>0.00000001
      if raise_error
        error("Non-ortogonal system $M")
      else
        return false
      end
    end
  end
  return true
end

function check_coord_sys(M::Array{Array{T,1},1}; raise_error::Bool=true
                                                                ) where{T<:Real}
  dims = 3
  newM = Array{Float64}(undef, dims, dims)
  for i in 1:dims
    newM[i, :] = M[i]
  end
  return check_coord_sys(newM; raise_error=raise_error)
end

"""
  `axis_rotation(r::Array{Float64, 1}, angle_deg::Float64)`

Returns the transformation matrix of rotation around an arbitrary axis of unit
vector `r`
"""
function axis_rotation(r::Array{T, 1}, angle_deg::Real) where{T<:Real}
  ux, uy, uz = r
  C = cos(angle_deg*pi/180)
  S = sin(angle_deg*pi/180)
  t = 1 - C
  M = [t*ux^2+C t*ux*uy-S*uz t*ux*uz+S*uy;
        t*ux*uy+S*uz t*uy^2+C t*uy*uz-S*ux;
        t*ux*uz-S*uy t*uy*uz+S*ux t*uz^2+C]
  return M
end

"""
  `rotation_matrix(yaw::Real, pitch::Real, roll::Real)`

Receives yaw, pitch, and roll angles (in degrees) and returns the rotation
matrix corresponding to this rotation.
(see http://planning.cs.uiuc.edu/node102.html)

NOTE: Naming follows aircraft convention, with
* yaw:    rotation about z-axis.
* pitch:  rotation about y-axis.
* roll:   rotation about x-axis.
"""
function rotation_matrix(yaw::Real, pitch::Real, roll::Real)
  a, b, g = yaw*pi/180, pitch*pi/180, roll*pi/180
  Rz = [cos(a) -sin(a) 0; sin(a) cos(a) 0; 0 0 1]
  Ry = [cos(b) 0 sin(b); 0 1 0; -sin(b) 0 cos(b)]
  Rx = [1 0 0; 0 cos(g) -sin(g); 0 sin(g) cos(g)]
  return Rz*Ry*Rx
end

"""
  `rotation_matrix2(roll::Real, pitch::Real, yaw::Real)`

Receives yaw, pitch, and roll angles (in degrees) and returns the rotation
matrix corresponding to this rotation.
(see http://planning.cs.uiuc.edu/node102.html)

NOTE: Naming follows aircraft convention, with
* roll:   rotation about x-axis.
* pitch:  rotation about y-axis.
* yaw:    rotation about z-axis.

**Examples**
```jldoctest
julia> M = gt.rotation_matrix(0, 0, 0)
3×3 Array{Float64,2}:
  1.0  0.0  0.0
  0.0  1.0  0.0
  0.0  0.0  1.0

julia> M = gt.rotation_matrix(90, 0, 0)
3×3 Array{Float64,2}:
  1.0  0.0   0.0
  0.0  0.0  -1.0
  0.0  1.0   0.0

julia> X = [0, 1, 0];
julia> Xp = M*X
3-element Array{Float64,1}:
  0.0
  0.0
  1.0

julia> M = gt.rotation_matrix(0, 90, 0)
3×3 Array{Float64,2}:
  0.0  0.0  1.0
  0.0  1.0  0.0
 -1.0  0.0  0.0

julia> X = [1, 0, 0];
julia> Xp = M*X
3-element Array{Float64,1}:
  0.0
  0.0
 -1.0

julia> M = gt.rotation_matrix(0, 0, 90)
3×3 Array{Float64,2}:
  0.0  -1.0  0.0
  1.0   0.0  0.0
  0.0   0.0  1.0

julia> X = [1, 0, 0];
julia> Xp = M*X
3-element Array{Float64,1}:
  0.0
  1.0
  0.0

julia> M = gt.rotation_matrix(0, 45, 0)
3×3 Array{Float64,2}:
  0.707  0.0  0.707
  0.0    1.0  0.0
 -0.707  0.0  0.707

julia> X = [1, 0, 0];
julia> Xp = M*X
3-element Array{Float64,1}:
  0.707
  0.0
 -0.707

julia> M = gt.rotation_matrix(0, 45, 45)
3×3 Array{Float64,2}:
  0.5    -0.707  0.5
  0.5     0.707  0.5
 -0.707   0.0    0.707

julia> X = [1, 0, 0];
julia> Xp = M*X
3-element Array{Float64,1}:
  0.5
  0.5
 -0.707
```
"""
function rotation_matrix2(roll::Real, pitch::Real, yaw::Real)
  return rotation_matrix(yaw, pitch, roll)
end
##### END OF ALGEBRA ###########################################################

################################################################################
# ORTHOGONAL SPACE TRANSFORMATIONS
################################################################################

# ORTHOGONAL SYSTEMS (See https://en.wikipedia.org/wiki/Orthogonal_coordinates)
function cylindrical3D(X)
    r, theta, z = X
    return [ r*cos(theta), r*sin(theta), z ]
end
function cylindrical2D(X)
    r, theta = X
    return [ r*cos(theta),  r*sin(theta)]
end
function spherical3D(X)
    r, theta, phi = X
    return [  r*sin(theta)*cos(phi),
              r*sin(theta)*sin(phi),
              r*cos(theta) ]
end
function parabolic3D(X)
    u,v,z = X
    return [ 1/2*(u^2-v^2), u*v, z ]
end
function paraboloidal3D(X)
    u,v,phi = X
    return [u*v*cos(phi), u*v*sin(phi), 1/2*(u^2-v^2)]
end
function elliptic3D(X; a=1)
    u,v,z = X
    return [a*cosh(u)*cos(v), a*sinh(u)*sin(v), z]
end
function prolate3D(X; a=1)
    xi,etha,phi = X
    return [a*sinh(xi)*sin(etha)*cos(phi),
            a*sinh(xi)*sin(etha)*sin(phi),
            a*cosh(xi)*cos(etha)]
end
function oblate3D(X; a=1)
    xi,etha,phi = X
    return [a*cosh(xi)*cos(etha)*cos(phi),
            a*cosh(xi)*cos(etha)*sin(phi),
            a*sinh(xi)*sin(etha)]
end
function bipolar3D(X; a=1)
    u,v,z = X
    return [a*sinh(v)/(cosh(v)-cos(u)),
            a*sin(u)/(cosh(v)-cos(u)),
            z]
end
function toroidal3D(X; a=1)
    u,v,phi = X
    return [a*sinh(v)*cos(phi)/( cosh(v)-cos(u) ),
            a*sinh(v)*sin(phi)/( cosh(v)-cos(u) ),
            a*sin(u)/( cosh(v)-cos(u) )]
end
function conical3D(X; b=2, c=1)
    r,mu,nu = X
    return [r*mu*nu/(b*c),
            r/b*sqrt( (mu^2-b^2)*(nu^2-b^2) / (b^2-c^2) ),
            r/c*sqrt( (mu^2-c^2)*(nu^2-c^2) / (c^2-b^2) )]
end
##### END OF SPACE TRANSFORMATIONS #############################################


text:
module Day4
export day4

function parseLine(line)
    if line == "" return nothing end
    Dict(map(p->Pair(split(p, ':'; limit=2)...), split(line, ' ')))
end

function readInput(inputfile="inputs/input04.txt")
    parsed = [parseLine(l) for l in eachline(inputfile)]

    data = Any[]
    current = Dict{String, String}()
    for p in parsed
        if p === nothing && length(current) > 0
            push!(data, current)
            current = Dict{String, String}()
            continue
        end
        current = merge(current, p)
    end
    if length(current) > 0
        push!(data, current)
    end
    data
end

function day4(inputfile="inputs/input04.txt")
    values = readInput(inputfile)
    validPart1 = filter(day4Part1, values)
    (length(validPart1), count(day4Part2, validPart1))
end

function day4Part1(fields)
    (length(fields) == 8) || (length(fields) == 7 && !haskey(fields, "cid"))
end

function day4Part2(values)
    validateBYR(values["byr"]) &&
    validateIYR(values["iyr"]) &&
    validateEYR(values["eyr"]) &&
    validateHGT(values["hgt"]) &&
    validateHCL(values["hcl"]) &&
    validateECL(values["ecl"]) &&
    validatePID(values["pid"])
end

inRange(v, lower, upper) = lower <= parse(Int, v) <= upper

validateBYR(v) = occursin(r"^\d{4}$", v) && inRange(v, 1920, 2002)
validateIYR(v) = occursin(r"^\d{4}$", v) && inRange(v, 2010, 2020)
validateEYR(v) = occursin(r"^\d{4}$", v) && inRange(v, 2020, 2030)
function validateHGT(v)
    m = match(r"^(\d+)(in|cm)$", v)
    if m === nothing return false end
    num, unit = m.captures
    if unit == "cm"
        inRange(num, 150, 196)
    elseif unit == "in"
        inRange(num, 59, 76)
    else
        false
    end
end
validateHCL(v) = occursin(r"^#[0-9a-f]{6}$", v)
validateECL(v) = v in ("amb", "blu", "brn", "gry", "grn", "hzl", "oth")
validatePID(v) = occursin(r"^[0-9]{9}$", v)

end #module