text:
using jInv.Mesh
using Multigrid
using SparseArrays
using LinearAlgebra
println("************************************************* GMG RAP for Elasticity 2D ******************************************************");

domain = [0.0, 1.0, 0.0, 1.0];
n = [128,128];
Mr = getRegularMesh(domain,n)
mu = 1.0*ones(prod(Mr.n));
lambda = 1.0*mu;
Ar = GetLinearElasticityOperator(Mr,mu,lambda);
Ar = Ar + 1e-3*opnorm(Ar,1)*sparse(1.0I,size(Ar,2),size(Ar,2));

levels      = 5;
numCores 	= 2;
maxIter     = 5;
relativeTol = 1e-10;
relaxType   = "SPAI";
relaxParam  = 0.75;
relaxPre 	= 2;
relaxPost   = 2;
cycleType   ='V';
coarseSolveType = "Julia";
transferOperatorType = "SystemsFacesLinear";

MG = getMGparam(Float64,Int64,levels,numCores,maxIter,relativeTol,relaxType,relaxParam,relaxPre,relaxPost,cycleType,coarseSolveType,0.5,0.0,transferOperatorType);

N = size(Ar,2);
b = Ar*rand(N,2); b = b/norm(b);
x = zeros(N,2);
MGsetup(Ar,Mr,MG,size(b,2),true);
println("****************************** Stand-alone GMG RAP: ******************************")
solveMG(MG,b,x,true);
@test norm(Ar*x - b) < 0.05;
println("****************************** GMG + CG: ******************************")
x[:].=0.0;
solveCG_MG(Ar,MG,b,x,true)
@test norm(Ar*x - b) < 0.01;

println("************************************************* GMG RAP for Elasticity 3D ******************************************************");

domain = [0.0, 1.0, 0.0, 1.0,0.0,1.0];
n = [16,16,12];
Mr = getRegularMesh(domain,n)
mu = 2.0*ones(prod(Mr.n));
lambda = mu;
Ar = GetLinearElasticityOperator(Mr,mu,lambda);
Ar = Ar + 1e-2*opnorm(Ar,1)*sparse(1.0I,size(Ar,2),size(Ar,2));

N = size(Ar,2);
b = Ar*rand(N,2); b = b/norm(b);
x = zeros(N,2);
clear!(MG)
MGsetup(Ar,Mr,MG,size(b,2),true);
println("****************************** Stand-alone GMG RAP: ******************************")
solveMG(MG,b,x,true);
@test norm(Ar*x - b) < 0.05;
println("****************************** GMG + CG: ******************************")
x[:].=0.0;
solveCG_MG(Ar,MG,b,x,true)
@test norm(Ar*x - b) < 0.01;


text:
@testset "init" begin
    _initialize_download_cache = DiagnosisClassification._initialize_download_cache
    @test _initialize_download_cache() isa Nothing
end

ext: jl


text:
module ARMSampling

import Base: length, copy, insert!
export TrapezoidalProposal, sample, sample!, overrelaxation!

#From StatsBase.jl/src/sampling.jl
function _sample(w)
    t = rand() * sum(w)
    n = length(w)
    i = 1
    cw = w[1]
    while cw < t && i < n
        i += 1
        @inbounds cw += w[i]
    end
    return i
end

slope(x0, x1, y0, y1) = (y1-y0)/(x1-x0)
intercept(x0, x1, y0, y1) = (x1*y0 - x0*y1)/(x1-x0)
slope_intercept(x0, x1, y0, y1) = (y1-y0)/(x1-x0), (x1*y0 - x0*y1)/(x1-x0)
straight_line_fun(x, x0, x1, y0, y1) = slope(x0, x1, y0, y1)*x + intercept(x0, x1, y0, y1)

function integrate_straight_line_exp(x0, x1, y0, y1)
    a, b = slope_intercept(x0, x1, y0, y1)
    if abs(a) < eps(typeof(a))
        return exp(b)*(x1-x0)
    end

    return (exp(a*x1+b)-exp(a*x0+b))/a
end

function trapezoidalweight(x::AbstractVector{T}, y::AbstractVector{T}) where T
    n = length(x) - 2
    w = T[integrate_straight_line_exp(x[i], x[i+1], y[i], y[i+1]) for i in 1:(n+1)]
    return w
end

struct TrapezoidalProposal{T,V<:AbstractVector{T}}
    x::V
    y::V
    w::V

    function TrapezoidalProposal{T,V}(x::V, y::V, w::V) where {T,V<:AbstractVector{T}}
        @assert issorted(x, lt = (x,y) -> x <= y)
        @assert length(x) >= 3
        @assert length(y) == length(x)
        @assert length(w) == length(x) - 1
        new{T,V}(x, y, w)
    end
end

function TrapezoidalProposal(x::AbstractVector{T}, y::AbstractVector{T}) where T
    w = trapezoidalweight(x, y)
    return TrapezoidalProposal{T,typeof(x)}(copy(x), copy(y), copy(w))
end

function TrapezoidalProposal(x::AbstractVector{T}, logf::Function) where T
    y = T[logf(xx) for xx in x]
    return TrapezoidalProposal(x, y)
end

length(proposal::TrapezoidalProposal) = length(proposal.x) - 2

function copy(proposal::TrapezoidalProposal{T, V}) where {T,V<:AbstractVector{T}}
    newx = copy(proposal.x)
    newy = copy(proposal.y)
    neww = copy(proposal.w)
    TrapezoidalProposal{T,V}(newx, newy, neww)
end

function insert!(proposal::TrapezoidalProposal{T}, x0::T, y0::T) where T
    n = length(proposal)
    x = proposal.x
    y = proposal.y
    w = proposal.w

    index = searchsorted(x, x0)
    left = index.stop
    right = index.start

    if left == right
        return proposal
    elseif left == 0
        w0 = integrate_straight_line_exp(x0, x[1], y0, y[1])
        insert!(w, 1, w0)
    elseif right == n + 3
        w0 = integrate_straight_line_exp(x[n+2], x0, y[n+2], y0)
        insert!(w, n+2, w0)
    else
        w_left  = integrate_straight_line_exp(x[left], x0, y[left], y0)
        w_right = integrate_straight_line_exp(x0, x[right], y0, y[right])
        splice!(w, left, [w_left, w_right])
    end
    insert!(x, right, x0)
    insert!(y, right, y0)
    return proposal
end

insert!(proposal::TrapezoidalProposal{T}, x0::T, logf::Function) where T = insert!(proposal, x0, logf(x0))

function (proposal::TrapezoidalProposal)(x0)
    if x0 < proposal.x[1] || x0 > proposal.x[end]
        error("x out of range of proposal: x = $x0")
    end
    x = proposal.x
    y = proposal.y
    index = searchsorted(x, x0)
    left = index.stop
    right = index.start
    if left == right
        return y[left]
    end

    return straight_line_fun(x0, x[left], x[right], y[left], y[right])
end

function sample(proposal::TrapezoidalProposal)
    x = proposal.x
    y = proposal.y
    w = proposal.w

    n = length(proposal)
    k = _sample(w)

    a, b = slope_intercept(x[k], x[k+1], y[k], y[k+1])
    if abs(a) < eps(typeof(a))
        return rand()*w[k]*exp(-b) + x[k]
    end

    #out = log(rand()*w[k]*a*exp(-b)+exp(a*x[k]))/a
    out = (log(rand()*w[k]*a + exp(y[k]))-b)/a
    if !isfinite(out)
        println()
        error("Got infinity! k = $k, a = $a, b = $b, w = $(w[k])")
    end
    return out
end

function sample!(proposal::TrapezoidalProposal{T},
                 logf::Function,
                 x0::T;
                 max_n = 2*length(proposal)) where T
    logfx0 = logf(x0)
    logfx = zero(T)
    x = zero(T)
    while true
        x = sample(proposal)
        logfx = logf(x)
        if rand() > exp(logfx - proposal(x))
            if length(proposal) < max_n
                insert!(proposal, x, logfx)
            end
        else
            break
        end
    end

    if rand() > min(1, exp(logfx - logfx0)*min(exp(logfx0), exp(proposal(x0)))/min(exp(logfx),exp(proposal(x))) )
    #if rand() > min(1, exp(logfx - logfx0 + proposal(x0) - proposal(x)))

        if length(proposal) < max_n
            insert!(proposal, x0, logfx0)
        end

        return x0
    else
        return x
    end
end

function sample!(out::Vector{T},
                 nsample::Integer,
                 proposal::TrapezoidalProposal{T},
                 logf::Function;
                 kwargs...) where T
    nout = length(out)
    append!(out, zeros(T, nsample))
    for i in 1:nsample
        out[i+nout] = sample!(proposal, logf, out[i+nout-1]; kwargs...)
    end
    return out
end

function rank(out)
    x0 = out[1]
    sorted_out = sort(out)
    index = searchsorted(sorted_out, x0)

    if index.start == index.stop
        r = index.start - 1
    else
        r = (index.start - 1) + _sample(ones(length(index))) - 1
    end
    return r, sorted_out
end

function overrelaxation!(proposal::TrapezoidalProposal{T},
                         logf::Function,
                         x0::T,
                         K::Integer;
                         max_n = 2*length(proposal),
                         max_try = 10*length(proposal)) where T
    x = x0
    max_try_counter = 0
    while max_try_counter < max_try && length(proposal) < max_n
        x = sample!(proposal, logf, x; max_n = max_n)
        max_try_counter += 1
    end

    out = T[x0]
    sample!(out, 2*K, proposal, logf; max_n = max_n)
    out = out[1:2:end]
    r, sorted_out = rank(out)
    return sorted_out[K-r+1], sorted_out
end

end # module


text:
# definitions of concrete subtypes of AbstractParamType and AbstractSolutionData

mutable struct ParamType{Tsol, Tres, Tdim} <: AbstractParamType{Tdim}

  f::BufferedIO     # TODO: needed?

  t::Float64
  time::Timings

  function ParamType{Tsol, Tres, Tdim}(mesh, sbp, opts) where {Tsol, Tres, Tdim}
    myrank = mesh.myrank
    if DB_LEVEL >= 1
      f = BufferedIO("log_$myrank.dat", "a")
    else
      f = BufferedIO(DevNull)  # create a dummy IOStream
    end

    t = 0.0

    time = Timings()

    return new(f, t, time)
  end

end   # end of ParamType type def

ParamType2{Tsol, Tres} =  ParamType{Tsol, Tres, 2}
ParamType3{Tsol, Tres} =  ParamType{Tsol, Tres, 3}
const ParamTypes = Union{ParamType2, ParamType3}

@doc """
### SimpleODEMod.SimpleODEData_

  This type is an implementation of the abstract SimpleODEData.  It is
  parameterized by the residual type Tres and the mesh type Tmsh
  because it stores some arrays of those types.  Tres is the 'maximum' type of
  Tsol and Tmsh, where Tsol is the type of the conservative variables.

"""->
mutable struct SimpleODEData_{Tsol, Tres, Tdim, Tmsh} <: SimpleODEData{Tsol, Tres, Tdim}

  # params::ParamType{Tdim}
  params::ParamType{Tsol, Tres, Tdim}

  comm::MPI.Comm
  commsize::Int
  myrank::Int

  t::Float64
  res_type::DataType  # type of res
  q::Array{Tsol, 3}
  q_face::Array{Tsol, 4}  # store solution values interpolated to faces
  res::Array{Tres, 3}      # result of computation
  res_vec::Array{Tres, 1}  # result of computation in vector form
  res_edge::Array{Tres, 4} # edge based residual storage
  q_vec::Array{Tres,1}     # initial condition in vector form
  q_bndry::Array{Tsol, 3}  # store solution variables interpolated to
                          # the boundaries with boundary conditions
  shared_data::Array{SharedFaceData{Tsol}, 1}
  M::Array{Float64, 1}       # mass matrix
  Minv::Array{Float64, 1}    # inverse mass matrix
  Minv3D::Array{Float64, 3}    # inverse mass matrix for application to res, not res_vec
  multiplyA0inv::Function       # multiply an array by inv(A0), where A0
                                # is the coefficient matrix of the time
                                # derivative
  majorIterationCallback::Function # called before every major (Newton/RK) itr

  function SimpleODEData_{Tsol, Tres, Tdim, Tmsh}(eqn::SimpleODEData_) where {Tsol, Tres, Tdim, Tmsh}

    return new()

  end

  function SimpleODEData_{Tsol, Tres, Tdim, Tmsh}(mesh::AbstractMesh, sbp::AbstractSBP, opts) where {Tsol, Tres, Tdim, Tmsh}
    println("\nConstruction SimpleODEData object")
    println("  Tsol = ", Tsol)
    println("  Tres = ", Tres)
    println("  Tdim = ", Tdim)
    println("  Tmsh = ", Tmsh)
    numfacenodes = mesh.numNodesPerFace

    @assert mesh.commsize == 1

    eqn = new()  # incomplete initialization
    eqn.comm = mesh.comm
    eqn.commsize = mesh.commsize
    eqn.myrank = mesh.myrank

    eqn.params = ParamType{Tsol, Tres, Tdim}(mesh, sbp, opts)
    eqn.t = 0.0
    eqn.res_type = Tres
    eqn.majorIterationCallback = majorIterationCallback
    eqn.M = calcMassMatrix(mesh, sbp, eqn)
    eqn.Minv = calcMassMatrixInverse(mesh, sbp, eqn)
    eqn.Minv3D = calcMassMatrixInverse3D(mesh, sbp, eqn)
    eqn.q = zeros(Tsol, 1, sbp.numnodes, mesh.numEl)
    eqn.res = zeros(Tsol, 1, sbp.numnodes, mesh.numEl)
    eqn.res_edge = Array{Tres}(0, 0, 0, 0)
    if mesh.isDG
      eqn.q_vec = reshape(eqn.q, mesh.numDof)
      eqn.res_vec = reshape(eqn.res, mesh.numDof)
    else
      eqn.q_vec = zeros(Tres, mesh.numDof)
      eqn.res_vec = zeros(Tres, mesh.numDof)
    end
    eqn.multiplyA0inv = matVecA0inv

    if mesh.isDG
      eqn.q_face = zeros(Tsol, 1, 2, numfacenodes, mesh.numInterfaces)
      eqn.q_bndry = zeros(Tsol, 1, numfacenodes, mesh.numBoundaryFaces)

      for i=1:mesh.npeers
        eqn.flux_sharedface[i] = zeros(Tres, 1, numfacenodes,
                                       mesh.peer_face_counts[i])
      end
    else
      eqn.q_face = Array{Tres}(0, 0, 0, 0)
      eqn.q_bndry = Array{Tsol}(0, 0, 0)
    end

    # send and receive buffers
    # TODO: update this for parallel
    eqn.shared_data = Array{SharedFaceData{Tsol}}(0)

    return eqn
  end # ends the constructer SimpleODEData_

end # End type SimpleODEData_

import ODLCommonTools.getAllTypeParams

@doc """
### SimpleODEMod.getAllTypeParameters

Gets the type parameters for mesh and equation objects.

**Input**

* `mesh` : Object of abstract meshing type.
* `eqn`  : Euler Equation object.
* `opts` : Options dictionary

**Output**

* `tuple` : Tuple of type parameters. Ordering is same as that of the concrete eqn object within this physics module.

"""->
function getAllTypeParams(mesh::AbstractMesh{Tmsh}, eqn::SimpleODEData_{Tsol, Tres, Tdim, Tmsh}, opts) where {Tmsh, Tsol, Tres, Tdim}

  tuple = (Tsol, Tres, Tdim, Tmsh)

  return tuple
end


text:
#=
grid_get_subgrid:
- Julia version: 1.1.0
- Author: Jirik
- Date: 2019-03-29
=#

function data_sub_from_block_sub(block_size, margin_size, bsub)
    xst = (block_size[1] * (bsub[1] - 1)) + 1 - margin_size
    xsp = (block_size[1] * (bsub[1] + 0)) + 0 + margin_size
    yst = (block_size[2] * (bsub[2] - 1)) + 1 - margin_size
    ysp = (block_size[2] * (bsub[2] + 0)) + 0 + margin_size
    zst = (block_size[3] * (bsub[3] - 1)) + 1 - margin_size
    zsp = (block_size[3] * (bsub[3] + 0)) + 0 + margin_size
    return xst, xsp, yst, ysp, zst, zsp
end

function number_of_blocks_per_axis(seg3d_size, block_size)
    println("block_size: ", block_size)
    if typeof(block_size) == Tuple{Int64}
        dim = nfields(seg3d_size)
    elseif  typeof(block_size) == Array{Int64, 1}
        dim = length(block_size)
    else
        warn("Unknown type of block_size")
    end
    println( "dim:", dim, " seg size:", seg3d_size, " block size: ", block_size)

    blocks_number = Array{Int}(undef, dim)
    for k in 1:dim
        # print(k)
        delenec = seg3d_size[k]
        delitel = block_size[k]
        number_for_this_axis = cld(delenec, delitel)
        # print(" ", delenec, ", ", delitel, ", ", number_for_this_axis)

        blocks_number[k] = number_for_this_axis
    end
    return prod(blocks_number), blocks_number
end

function limit_by!(stop, data_size)
    for i=1:length(stop)
        if data_size[i] < stop[i]
            stop[i] = data_size[i]
        end
    end
    return stop
end

    segmentation = LarSurf.generate_slope([11,12,13])
    data_size = LarSurf.size_as_array(size(segmentation))
    data3d = segmentation
    block_size = [2,3,4]
    block_i = 91
    margin_size = 0

    block_number, blocks_number_axis = number_of_blocks_per_axis(data_size, block_size)
    a = Array{Int}(
        undef,
        blocks_number_axis[1],
        blocks_number_axis[2],
        blocks_number_axis[3]
    )
    bsub = CartesianIndices(a)[block_i]
#     faces = LarSurf.cube_in_block_surface([1,2,3], [1,1,1], [1,2,2])
    bsub_arr = [bsub[1], bsub[2], bsub[3]]

    first = (bsub_arr .== [1, 1, 1])
    last = (bsub_arr .== blocks_number_axis)

    # if margin_size == 0

    xst, xsp, yst, ysp, zst, zsp = LarSurf.data_sub_from_block_sub(
        block_size, margin_size, bsub
    )
    stop = [xsp, ysp, zsp]
    print("stop ", stop)
    limit_by!(stop, data_size)
#     new_stop = zeros(eltype(stop), size(stop)...)
    print("stop ", stop)

#     if data_size[1] < xsp
#         xsp = data_size[1]
#     end
#     if data_size[1] < xsp
#         xsp = data_size[1]
#     end


text:
using Pkg
Pkg.activate(".")
Pkg.instantiate()

using ArgParse
using CGT
using CurrencyAmounts

function parse_commandline()
	s = ArgParseSettings(description="Computes the capital gain tax (in Ireland)")

	@add_arg_table s begin
		"--verbose", "-v"
			help = "Show transactions"
			action = :store_true
		"file"
			help = "Etrade Excel file"
			required = true
	end

	return parse_args(s)
end

function main()
	parsed_args = parse_commandline()

	sales = CGT.load(parsed_args["file"])
	sales = convert(CGT.Sales{Currency{:EUR}}, sales)
	CGT.compute(sales; details=parsed_args["verbose"])
end

main()


text:
module SimpleSudoku

export solve_sudoku

# find zeros go through each one and find the ones that are not in the same row, column or box
function get_candidates(grid, row, col)
    grid_id_row = ceil(Int, row / 3)
    grid_id_col = ceil(Int, col / 3)

    candidates = setdiff(0:9, mapreduce(vcat, (
        grid[1:row-1, col],
        grid[row+1:9, col],
        grid[row, 1:col-1],
        grid[row, col+1:9], +
        # the 3x3 grid within which the cell is contained
        grid[(grid_id_row-1)*3+1:grid_id_row*3, (grid_id_col-1)*3+1:grid_id_col*3])) do items
        vec(items)
    end)

    candidates
end

get_candidates(grid, cart) = get_candidates(grid, cart[1], cart[2])

"""
    solve_sudoku(grid)

Args:

    grid a 9x9 integer array where 0 indicates unfilled poistions

Returns:
    Tuple{grid::Array{2}, is_solution_valid::Bool}

    where `grid` is a solved grid if the second return value is true
    If the second return value is false then it means there is
        no solution starting with the input grid
"""
function solve_sudoku(grid)
    @assert size(grid) == (9, 9)
    @assert all(in.(grid, Ref(0:9)))
    # these are the unfilled positions
    zero_pos = findall(==(0), grid)

    if length(zero_pos) == 0
        return grid, true
    end

    cand = get_candidates.(Ref(grid), zero_pos)

    # any of them has no possible candidates; then it's an invalid solution
    if any(length.(cand) .== 0)
        return grid, false
    end

    # find the cell with the lowest number of candidates
    idx = argmin(length.(cand))

    # try the candidates one by one
    zp = zero_pos[idx]
    for c in cand[idx]
        cgrid = copy(grid)
        cgrid[zp] = c
        new_grid, valid = solve_sudoku(cgrid)
        if valid
            return new_grid, true
        end
    end

    return grid, false
end

end


text:
module PseudoBoolean

include("hierarchy.jl")
include("variables.jl")
include("assignment.jl")
include("binarymodels/binarymodels.jl")

end # module


text:
# Ref: https://github.com/stan-dev/stan/blob/develop/src/stan/mcmc/hmc/hamiltonians/diag_e_metric.hpp

"""
    gen_grad_func(vi::VarInfo, sampler::Sampler, model)

Generate a function that takes a vector of reals `θ` and compute the logpdf and
gradient at `θ` for the model specified by `(vi, sampler, model)`.
"""
function gen_grad_func(vi::VarInfo, sampler::Sampler, model)
    return θ::AbstractVector{<:Real}->gradient(θ, vi, model, sampler)
end

"""
    gen_lj_func(vi::VarInfo, sampler::Sampler, model)

Generate a function that takes `θ` and returns logpdf at `θ` for the model specified by
`(vi, sampler, model)`.
"""
function gen_lj_func(vi::VarInfo, sampler::Sampler, model)
    return function(θ::AbstractVector{<:Real})
        vi[sampler] = θ
        return runmodel!(model, vi, sampler).logp
    end
end

"""
  gen_rev_func(vi::VarInfo, sampler::Sampler)

Generate a function on `(θ, logp)` that sets the variables referenced by `sampler` to `θ`
and the current `vi.logp` to `logp`.
"""
function gen_rev_func(vi::VarInfo, sampler::Sampler)
    return function(θ::AbstractVector{<:Real}, logp::Real)
        vi[sampler] = θ
        setlogp!(vi, logp)
    end
end

"""
    gen_log_func(sampler::Sampler)

Generate a function that takes no argument and performs logging for the number of leapfrog
steps used in `sampler`.
"""
function gen_log_func(sampler::Sampler)
    return function()
        sampler.info[:lf_num] += 1
    end
end

# TODO: improve typing for all generator functions

function gen_momentum_sampler(vi::VarInfo, spl::Sampler)
    d = length(vi[spl])
    return function()
        return randn(d)
    end
end

function gen_H_func()
    return function(θ::AbstractVector{<:Real},
                    p::AbstractVector{<:Real},
                    logp::Real)
        H = sum(abs2, p) / 2 - logp
        return isnan(H) ? Inf : H
    end
end

function gen_momentum_sampler(vi::VarInfo, spl::Sampler, ::UnitPreConditioner)
    d = length(vi[spl])
    return function()
        return randn(d)
    end
end

function gen_H_func(::UnitPreConditioner)
    return function(θ::AbstractVector{<:Real},
                    p::AbstractVector{<:Real},
                    logp::Real)
        H = sum(abs2, p) / 2 - logp
        return isnan(H) ? Inf : H
    end
end

function gen_momentum_sampler(vi::VarInfo, spl::Sampler, pc::DiagPreConditioner)
    d = length(vi[spl])
    std = pc.std
    return function()
        return randn(d) ./ std
    end
end

function gen_H_func(pc::DiagPreConditioner)
    std = pc.std
    return function(θ::AbstractVector{<:Real},
                    p::AbstractVector{<:Real},
                    logp::Real)
        H = sum(abs2, p .* std) / 2 - logp
        return isnan(H) ? Inf : H
    end
end

# NOTE: related Hamiltonian change: https://github.com/stan-dev/stan/blob/develop/src/stan/mcmc/hmc/hamiltonians/dense_e_metric.hpp
function gen_momentum_sampler(vi::VarInfo, spl::Sampler, pc::DensePreConditioner)
    d = length(vi[spl])
    A = Symmetric(pc.covar)
    C = LinearAlgebra.cholesky(A)
    return function()
        return C.U \ randn(d)
    end
end

function gen_H_func(pc::DensePreConditioner)
    A = pc.covar
    return function(θ::AbstractVector{<:Real},
                    p::AbstractVector{<:Real},
                    logp::Real)
        H = p' * A * p / 2 - logp
        return isnan(H) ? Inf : H
    end
end

###

function runmodel!(model, vi::VarInfo, spl::Union{Nothing,Sampler})
    setlogp!(vi, zero(Real))
    if spl != nothing && :eval_num ∈ keys(spl.info)
        spl.info[:eval_num] += 1
    end
    model(vi, spl)
    return vi
end

function leapfrog(θ::AbstractVector{<:Real},
                  p::AbstractVector{<:Real},
                  τ::Int,
                  ϵ::Real,
                  model,
                  vi::VarInfo,
                  sampler::Sampler,
                  )
    lp_grad = gen_grad_func(vi, sampler, model)
    rev = gen_rev_func(vi, sampler)
    logger = gen_log_func(sampler)
    return _leapfrog(θ, p, τ, ϵ, lp_grad; rev_func=rev, log_func=logger)
end

function _leapfrog(θ::AbstractVector{<:Real},
                   p::AbstractVector{<:Real},
                   τ::Int,
                   ϵ::Real,
                   lp_grad_func::Function;
                   rev_func=nothing,
                   log_func=nothing,
                   )
    _, grad = lp_grad_func(θ)
    verifygrad(grad) || (return θ, p, 0)

    p, θ, τ_valid = deepcopy(p), deepcopy(θ), 0

    p .-= ϵ .* grad ./ 2
    for t in 1:τ

        log_func != nothing && log_func()

        θ .+= ϵ .* p
        logp, grad = lp_grad_func(θ)

        # If gradients explode, tidy up and return.
        if ~verifygrad(grad)
            θ .-= ϵ .* p
            rev_func != nothing && rev_func(θ, logp)
            break
        end

        p .-= ϵ .* grad
        τ_valid += 1
    end

    # Undo half a step in the momenta.
    p .+= ϵ .* grad ./ 2

    return θ, p, τ_valid
end

function _hmc_step(θ::AbstractVector{<:Real},
                   lj::Real,
                   lj_func::Function,
                   grad_func::Function,
                   H_func::Function,
                   τ::Int,
                   ϵ::Real,
                   momentum_sampler::Function;
                   rev_func=nothing,
                   log_func=nothing,
                   )
    @debug "sampling momentums..."
    p = momentum_sampler()

    @debug "recording old values..."
    H = H_func(θ, p, lj)

    @debug "leapfrog for $τ steps with step size $ϵ"
    θ_new, p_new, τ_valid = _leapfrog(θ, p, τ, ϵ, grad_func; rev_func=rev_func, log_func=log_func)

    @debug "computing new H..."
    lj_new = lj_func(θ_new)
    H_new = (τ_valid == 0) ? Inf : H_func(θ_new, p_new, lj_new)

    @debug "deciding wether to accept and computing accept rate α..."
    is_accept, logα = mh_accept(H, H_new)

    if is_accept
        θ = θ_new
        lj = lj_new
    end

    return θ, lj, is_accept, τ_valid, exp(logα)
end

function _hmc_step(θ::AbstractVector{<:Real},
                   lj::Real,
                   lj_func::Function,
                   grad_func::Function,
                   H_func::Function,
                   ϵ::Real,
                   λ::Real,
                   momentum_sampler::Function;
                   rev_func=nothing,
                   log_func=nothing,
                   )
    τ = max(1, round(Int, λ / ϵ))
    return _hmc_step(θ, lj, lj_func, grad_func, H_func, τ, ϵ, momentum_sampler;
                     rev_func=rev_func, log_func=log_func)
end

# TODO: remove used Turing-wrapper functions

# Ref: https://github.com/stan-dev/stan/blob/develop/src/stan/mcmc/hmc/base_hmc.hpp
function find_good_eps(model, spl::Sampler{T}, vi::VarInfo) where T
    logpdf_func_float = gen_lj_func(vi, spl, model)
    momentum_sampler = gen_momentum_sampler(vi, spl)
    H_func = gen_H_func()

    @info "[Turing] looking for good initial eps..."
    ϵ = 0.1

    p = momentum_sampler()

    θ = vi[spl]
    H0 = H_func(θ, p, logpdf_func_float(θ))

    θ_prime, p_prime, τ = leapfrog(θ, p, 1, ϵ, model, vi, spl)
    h = τ == 0 ? Inf : H_func(θ_prime, p_prime, logpdf_func_float(θ_prime))

    delta_H = H0 - h
    direction = delta_H > log(0.8) ? 1 : -1

    iter_num = 1

    # Heuristically find optimal ϵ
    while (iter_num <= 12)

        p = momentum_sampler()
        H0 = H_func(vi[spl], p, logpdf_func_float(vi[spl]))

        θ_prime, p_prime, τ = leapfrog(θ, p, 1, ϵ, model, vi, spl)
        h = τ == 0 ? Inf : H_func(θ_prime, p_prime, logpdf_func_float(θ_prime))
        @debug "direction = $direction, h = $h"

        delta_H = H0 - h

        if ((direction == 1) && !(delta_H > log(0.8)))
            break
        elseif ((direction == -1) && !(delta_H < log(0.8)))
            break
        else
            ϵ = direction == 1 ? 2.0 * ϵ : 0.5 * ϵ
        end

        iter_num += 1
    end

    while h == Inf  # revert if the last change is too big
        ϵ = ϵ / 2               # safe is more important than large
        θ_prime, p_prime, τ = leapfrog(θ, p, 1, ϵ, model, vi, spl)
        h = τ == 0 ? Inf : H_func(θ_prime, p_prime, logpdf_func_float(θ_prime))
    end
    @info "\r[$T] found initial ϵ: $ϵ"

    return ϵ
end

"""
    mh_accept(H::Real, H_new::Real)
    mh_accept(H::Real, H_new::Real, log_proposal_ratio::Real)

Peform MH accept criteria. Returns a boolean for whether or not accept and the
acceptance ratio in log space.
"""
mh_accept(H::Real, H_new::Real) = log(rand()) + H_new < min(H_new, H), min(0, -(H_new - H))
function mh_accept(H::Real, H_new::Real, log_proposal_ratio::Real)
    return log(rand()) + H_new < H + log_proposal_ratio, min(0, -(H_new - H))
end


text:
fs2()

set_curpath("index.md")

@testset "Evalcode" begin
    h = raw"""
        Simple code:
        ```julia:./code/exca1
        a = 5
        print(a^2)
        ```
        then:
        \output{./code/exca1}
        done.
        """ |> seval

    h2 = raw"""
        Simple code:
        ```julia:excb1
        a = 5
        print(a^2)
        ```
        then:
        \output{excb1}
        done.
        """ |> seval

    @test h == h2

    spatha = joinpath(F.PATHS[:site], "assets", "index", "code", "exca1.jl")
    spathb = joinpath(F.PATHS[:site], "assets", "index", "code", "excb1.jl")
    @test isfile(spatha)
    @test isfile(spathb)
    @test isapproxstr(read(spatha, String), """
        $(F.MESSAGE_FILE_GEN_FMD)
        a = 5\nprint(a^2)""")

    opath = joinpath(F.PATHS[:site], "assets", "index", "code", "output", "exca1.out")
    @test isfile(opath)
    @test read(opath, String) == "25"

    @test isapproxstr(h, raw"""

Simple code:

a = 5
            print(a^2)

            then:

25

            done.
""")
end

@testset "Eval (errs)" begin
    s = raw"""
        Simple code:
        ```python:./scripts/testpy
        a = 5
        print(a**2)
        ```
        done.
        """
    h = ""
    @test_logs (:warn, "Evaluation of non-Julia code blocks is not yet supported.") (h = s |> seval)

    @test isapproxstr(h, raw"""

Simple code:

a = 5
            print(a**2)

            done.
""")
end

@testset "Eval (rinput)" begin
    h = raw"""
        Simple code:
        ```julia:/assets/scripts/test2
        a = 5
        print(a^2)
        ```
        then:
        \output{/assets/scripts/test2}
        done.
        """ |> seval

    spath = joinpath(F.PATHS[:site], "assets", "scripts", "test2.jl")
    @test isfile(spath)
    @test occursin("a = 5\nprint(a^2)", read(spath, String))

    opath = joinpath(F.PATHS[:site], "assets", "scripts", "output", "test2.out")
    @test isfile(opath)
    @test read(opath, String) == "25"

    @test isapproxstr(h, """

Simple code:

a = 5
            print(a^2)

            then:

25

            done.
""")

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

    h = raw"""
        @def fd_rpath = "pages/pg1.md"
        Simple code:
        ```julia:./code/abc2
        a = 5
        print(a^2)
        ```
        then:
        \output{./code/abc2}
        done.
        """ |> seval

    spath = joinpath(F.PATHS[:site], "assets", "pages", "pg1", "code", "abc2.jl")
    @test isfile(spath)
    @test occursin("a = 5\nprint(a^2)", read(spath, String))

    opath = joinpath(F.PATHS[:site], "assets", "pages", "pg1", "code", "output" ,"abc2.out")
    @test isfile(opath)
    @test read(opath, String) == "25"

    @test isapproxstr(h, """

Simple code:

a = 5
            print(a^2)

            then:

25
  done.
""")
end

@testset "Eval (module)" begin
    h = raw"""
        Simple code:
        ```julia:scripts/test1
        using LinearAlgebra
        a = [5, 2, 3, 4]
        print(dot(a, a))
        ```
        then:
        \output{scripts/test1}
        done.
        """ |> seval
    # dot(a, a) == 54
    @test occursin("""then:
54
 done.""", h)
end

@testset "Eval (img)" begin
    h = raw"""
        @def reeval=true
        Simple code:
        ```julia:tv2
        #hideall
        write(joinpath(@OUTPUT, "tv2.png"), "blah")
        ```
        then:
        \input{plot}{tv2}
        done.
        """ |> seval
    @test occursin("then:  done.", h)
end

@testset "Eval (throw)" begin
    s = raw"""
        @def reeval = true
        Simple code:
        ```julia:scripts/test1
        sqrt(-1)
        ```
        then:
        \output{scripts/test1}
        done.
        """
    global h
    h = ""
    @test_logs (:warn, "There was an error of type DomainError running the code.") (global h; h = s |> seval)
    # errors silently
    if VERSION >= v"1.2"
        @test occursin("""
Simple code:
sqrt(-1)
 then:
DomainError with -1.0:
    sqrt will only return a complex result if called with a complex argument. Try sqrt(Complex(x)).

 done.

    """, h)
    end
end

@testset "Eval (nojl)" begin
    h = raw"""
        Simple code:
        ```python:scripts/test1
        sqrt(-1)
        ```
        done.
        """

    @test (@test_logs (:warn, "Evaluation of non-Julia code blocks is not yet supported.") h |> seval) == "
Simple code:
sqrt(-1)
 done.
\n"
end

# temporary fix for 186: make error appear and also use `abspath` in internal include
@testset "eval #186" begin
    h = raw"""
        Simple code:
        ```julia:scripts/test186
        fn = "tempf.jl"
        write(fn, "a = 1+1")
        println("Is this a file? $(isfile(fn))")
        include(abspath(fn))
        println("Now: $a")
        rm(fn)
        ```
        done.
        \output{scripts/test186}
        """ |> seval
    @test isapproxstr(h, raw"""

Simple code:
fn = "tempf.jl"
            write(fn, "a = 1+1")
            println("Is this a file? $(isfile(fn))")
            include(abspath(fn))
            println("Now: $a")
            rm(fn)
 done.
Is this a file? true
            Now: 2

            """)
end

@testset "show" begin
    h = raw"""
        @def hascode = true
        @def reeval = true
        ```julia:ex
        a = 5
        a *= 2
        ```
        \show{ex}
        """ |> fd2html_td
    @test isapproxstr(h, """

a = 5
        a *= 2

10

        """)

    # Show with stdout
    h = raw"""
        @def hascode = true
        @def reeval = true
        ```julia:ex
        a = 5
        println("hello")
        a *= 2
        ```
        \show{ex}
        """ |> fd2html_td
    @test isapproxstr(h, """

a = 5
        println("hello")
        a *= 2

hello
        10

        """)
end


text:
"""
    sampks(X; k, metric = "eucl")
Kennard-Stone sampling.
* `X` : X-data.
* `k` : Nb. observations to sample (=> output `train`).
* `metric` : Metric used for the distance computation.
    Possible values: "eucl", "mahal".

The function divides the data X in two sets, `train` vs `test`,
using the Kennard-Stone (KS) algorithm (Kennard & Stone, 1969).
The two sets have different underlying probability distributions:
`train` has higher dispersion than `test`.

## References

Kennard, R.W., Stone, L.A., 1969. Computer aided design of experiments.
Technometrics, 11(1), 137-148.

"""
function sampks(X; k, metric = "eucl")
    if metric == "eucl"
        D = euclsq(X, X)
    elseif metric == "mahal"
        D = mahsq(X, X)
    end
    zn = 1:size(D, 1)
    # Initial 2 selections (train)
    s = findall(D .== maximum(D))
    s = [s[1][1] ; s[1][2]]
    # Candidates
    cand = zn[setdiff(1:end, s)]
    @inbounds for i = 1:(k - 2)
        u = vec(minimum(D[s, cand], dims = 1))
        zs = cand[findall(u .== maximum(u))[1]]
        s = [s ; zs]
        cand = zn[setdiff(1:end, s)]
    end
    (train = s, test = cand)
end

"""
    sampdp(X; k, metric = "eucl")
DUPLEX sampling.
* `X` : X-data.
* `k` : Nb. pairs of observations to sample. Must be <= n/2.
* `metric` : Metric used for the distance computation.
    Possible values: "eucl", "mahal".

The function divides the data `X` in two sets of equal size, `train` vs `test`,
using the DUPLEX algorithm (Snee, 1977 p.421).
The objective is to divide the data into two
sets which cover approximately the same region and
have similar statistical properties.

The user may add (a posteriori) the eventual remaining observations (`remain`)
to `train`.

## References

Kennard, R.W., Stone, L.A., 1969. Computer aided design of experiments.
Technometrics, 11(1), 137-148.

Snee, R.D., 1977. Validation of Regression Models: Methods and Examples.
Technometrics 19, 415-428. https://doi.org/10.1080/00401706.1977.10489581

"""
function sampdp(X; k, metric = "eucl")
    if(metric == "eucl")
        D = euclsq(X, X)
    elseif(metric == "mahal")
        D = mahsq(X, X)
    end
    n = size(D, 1)
    zn = 1:n
    # Initial selection of 2 pairs
    s = findall(D .== maximum(D))
    s1 = [s[1][1] ; s[1][2]]
    zD = copy(D)
    zD[s1, :] .= -Inf ;
    zD[:, s1] .= -Inf ;
    s = findall(D .== maximum(zD))
    s2 = [s[1][1] ; s[1][2]]
    # Candidates
    cand = zn[setdiff(1:end, [s1 ; s2])]
    @inbounds for i = 1:(k - 2)
        u = vec(minimum(D[s1, cand], dims = 1))
        zs = cand[findall(u .== maximum(u))[1]]
        s1 = [s1 ; zs]
        cand = zn[setdiff(1:end, s1)]
        u = vec(minimum(D[s2, cand], dims = 1))
        zs = cand[findall(u .== maximum(u))[1]]
        s2 = [s2 ; zs]
        cand = zn[setdiff(1:end, [s1 ; s2])]
    end
    (train = s1, test = s2, remain = cand)
end

"""
    sampsys(y; k)
Systematic sampling.
* `y` : Quantitative variable.
* `k` : Nb. observations to sample (=> output `train`). Must be >= 2.

Systematic sampling (regular grid) over the quantitative variable y.

The minimum and maximum of y are always sampled.
"""
function sampsys(y; k)
    y = vec(y)
    n = length(y)
    nint = k - 1            # nb. intervals
    alpha = (n - 1) / nint  # step
    z = collect(1:alpha:n)
    z = Int64.(round.(z))
    z = unique(z)
    id = sortperm(y)
    zn = collect(1:n)
    s = zn[id[z]]
    test = zn[setdiff(1:end, s)]
    (train = s, test = test)
end

"""
    sampclas(x; k)
Stratified sampling.
* `x` : Classes of the observations.
* `y` : Quantitative variable for systematic sampling.
* `k` : Nb. observations to sample in each class (=> output `train`).

The length of `k` must be either 1 (equal number of training observations to select per class)
or the number of classes in `x`.

If `y = nothing` (default), the sampling is random, else it is systematic (grid over `y`).

## References

Naes, T., 1987. The design of calibration in near infra-red reflectance analysis by clustering.
Journal of Chemometrics 1, 121-134.
"""
function sampclas(x, y = nothing; k)
    x = vec(x)
    n = length(x)
    res = tab(x)
    lev = res.keys
    nlev = length(lev)
    ni = collect(values(res))
    isequal(length(k), 1) ? k = fill(k[1], nlev) : nothing
    s = list(nlev, Vector{Int64})
    @inbounds for i in 1:nlev
        k[i] = min(k[i], ni[i])
        zs = findall(x .== lev[i])
        if isnothing(y)
            s[i] = sample(zs, k[i]; replace = false)
        else
            y = vec(y)
            u = sampsys(y[zs]; k = k[i]).train
            s[i] = zs[u]
        end
    end
    s = reduce(vcat, s)
    zn = collect(1:n)
    test = zn[setdiff(1:end, s)]
    (train = s, test = test, lev = lev, ni = ni, k = k)
end


text:
push!(LOAD_PATH, "../src/")
using Documenter , AbstractSDRs

makedocs(sitename="AbstractSDRs.jl",
		 format = Documenter.HTML(),
		 pages    = Any[
						"Introduction to AbstractSDRs"   => "index.md",
						"Function list"         => "base.md",
						"Examples"              => Any[
														 "Examples/example_setup.md"
														 "Examples/example_parameters.md"
														 "Examples/example_benchmark.md"
														 "Examples/example_mimo.md"
														 ],
						],
		 );

# makedocs(modules = [AbstractSDRs],sitename="AbstractSDRs Documentation", format = Documenter.HTML(prettyurls = false))

deploydocs(
    repo = "github.com/JuliaTelecom/AbstractSDRs.jl",
)


text:
module KungIChi

using Configurations
using Expronicon
using MLStyle
using JSON

export ValidationError,
    @schema,
    @validate, @check,
    multiple_of,
    Constraint,
    Positive,
    Negative,
    MultipleOf,
    SecretString,
    LessThan,
    LessEqual,
    GreaterThan,
    GreaterEqual,
    Interval,
    interval,
    to_json,
    from_json,
    Maybe,
    # Configurations
    @option,
    to_toml,
    no_default,
    Maybe,
    # reflection
    field_aliases,
    field_default,
    field_alias,
    type_alias,
    # traits
    is_option,
    # parse
    from_dict,
    from_kwargs,
    from_toml,
    from_toml_if_exists,
    # serialize
    to_dict

include("pirates.jl")
include("types.jl")
include("check.jl")
include("schema.jl")
include("json.jl")
include("openapi.jl")

end


text:
# loss = zeros(size(g.W))
# if fiber.loss == true
#     for i in eachindex(gases2)
#         if gases2[i].type == :O3
#             path = joinpath(pwd(), "ozone")
#             data = readdlm("$path\\full_absorption_data.txt")
#             data = data[:,1]
#             freq = readdlm("$path\\freqdata_ozone.txt")
#             freq = freq[:,1]

#             λ = 2pi*3e8./freq*1e9
#             λ_bound = 2.0 .< λ .< 400.0
#             data = data[λ_bound]
#             freq = freq[λ_bound]
#             data *= gases2[i].p*gases2[i].pp*N0
#             spl = Spline1D(freq, data)
#             loss .= spl.(g.W)
#         end
#     end
# end

# λ = 2pi*3e8./g.W*1e9
# λ_bound = 2.0 .< λ .< 400.0

# open(dir*"\\loss.txt", "w") do io
#     writedlm(io,  loss)
# end

max_stars_repo_name: SabbahMohammed/Dev

max_issues_repo_name: SabbahMohammed/Dev

max_forks_repo_name: SabbahMohammed/Dev


text:
""" Very coarse solution for gm to check for feasibility and an initial starting point. """
function surveysolve(gm::GlobalModel)
    bbcs = [bbl for bbl in gm.bbls if bbl isa BlackBoxClassifier]
    if !isempty(bbcs)
        for bbc in bbcs
            learn_constraint!(bbc)
            update_tree_constraints!(gm, bbc)
        end
    end
    bbrs = [bbl for bbl in gm.bbls if bbl isa BlackBoxRegressor]
    if !isempty(bbrs)
        for bbr in bbrs
            learn_constraint!(bbr, regression_sparsity = 0, max_depth = 6, minbucket = 2*length(bbr.vars))
            update_tree_constraints!(gm, bbr)
        end
    end
    optimize!(gm)
    return
end

"""
    function add_relaxation_variables!(gm::GlobalModel, bbl::Union{BlackBoxLearner, LinkedLearner})
    function add_relaxation_variables!(gm::GlobalModel, bbls::Array)

Populates relax_var attributes of all substructs.
"""
function add_relaxation_variables!(gm::GlobalModel, bbl::Union{BlackBoxLearner, LinkedLearner})
    if bbl.relax_var isa Real
        bbl.relax_var = @variable(gm.model)
        set_lower_bound(bbl.relax_var, 0)
    end
    if !isempty(bbl.mi_constraints)
        clear_tree_constraints!(gm, bbl)
        @info "Clearing MI constraints from $(bbl.name) due to relaxation."
    end
    if bbl isa BlackBoxLearner
        for ll in bbl.lls
            add_relaxation_variables!(gm, ll)
        end
    end
end

function add_relaxation_variables!(gm::GlobalModel, bbls::Array)
    for bbl in bbls
        add_relaxation_variables!(gm, bbl)
    end
end

add_relaxation_variables!(gm::GlobalModel) = add_relaxation_variables!(gm, gm.bbls)

"""
    $(TYPEDSIGNATURES)

Clears slack variables in GlobalModel.
"""
function clear_relaxation_variables!(gm::GlobalModel, bbl::Union{BlackBoxLearner, LinkedLearner})
    if bbl.relax_var isa JuMP.VariableRef
        is_valid(gm.model, bbl.relax_var) ||
            throw(OCTHaGOnException("$(bbl.relax_var) is not a valid variable of  $(gm.model)."))
        delete_lower_bound(bbl.relax_var)
        delete(gm.model, bbl.relax_var)
        bbl.relax_var = 0
    end
    if bbl isa BlackBoxLearner
        for ll in bbl.lls
            clear_relaxation_variables!(gm, ll)
        end
    end
end

function clear_relaxation_variables!(gm::GlobalModel, bbls::Array)
    for bbl in bbls
        clear_relaxation_variables!(gm, bbl)
    end
end

clear_relaxation_variables!(gm::GlobalModel) = clear_relaxation_variables!(gm, gm.bbls)

"""
    $(TYPEDSIGNATURES)

Includes relaxation variables in as a linear penalty in the objective function.
"""
function relax_objective!(gm::GlobalModel, M::Real = 1e8)
    no = 0
    for bbl in gm.bbls
        no += M*bbl.relax_var
        for ll in bbl.lls
            no += M*ll.relax_var
        end
    end
    @objective(gm.model, Min, gm.objective + no)
    return
end

"""
Removes relaxation variables from objective.
"""
function tight_objective!(gm::GlobalModel)
    @objective(gm.model, Min, gm.objective)
end

"""
    $(TYPEDSIGNATURES)

Performs gradient descent on the last optimal solution in gm.solution_history.
In case of infeasibility, first projects the feasible point using the local
constraint gradients.

Optional kwargs:
max_iterations: maximum number of gradient steps or projections.
step_penalty: magnitude of penalty on step size during projections.
step_size: Size of 0-1 normalized Euclidian ball we can step.
decay_rate: Exponential coefficient of step size reduction.
"""
function descend!(gm::GlobalModel; kwargs...)
    # Initialization
    clear_tree_constraints!(gm)
    bbls = gm.bbls
    vars = gm.vars
    bounds = get_bounds(vars)

    # Update descent algorithm parameters
    if !isempty(kwargs)
        for item in kwargs
            set_param(gm, item.first, item.second)
        end
    end

    # Checking for a nonlinear objective
    obj_bbl = nothing
    if gm.objective isa VariableRef
        obj_bbl = filter(x -> x.dependent_var == gm.objective,
                            [bbl for bbl in gm.bbls if bbl isa BlackBoxRegressor])
        if isempty(obj_bbl)
            obj_bbl = nothing
        elseif length(obj_bbl) == 1
            obj_bbl = obj_bbl[1]
            vars = [var for var in vars if var != obj_bbl.dependent_var]
            bbls = [bbl for bbl in gm.bbls if bbl != obj_bbl]
        elseif length(obj_bbl) > 1
            throw(OCTHaGOnException("GlobaModel $(gm.name) has more than one BlackBoxRegressor" *
                               " on the objective variable."))
        end
    end

    grad_shell = DataFrame(string.(vars) .=> [Float64[] for i=1:length(vars)])

    # Add relaxation variables
    add_relaxation_variables!(gm)

    # Final checks
    obj_gradient = copy(grad_shell)
    if isnothing(obj_bbl)
        if gm.objective isa VariableRef
            append!(obj_gradient, DataFrame(string.(gm.objective) => 1), cols = :subset)
        elseif gm.objective isa JuMP.GenericAffExpr
            append!(obj_gradient, DataFrame(Dict(string(key) => value for (key,value) in gm.objective.terms)), cols = :subset)
        else
            @warn "Type of objective $(gm.objective) is unsupported."
        end
        obj_gradient = coalesce.(obj_gradient, 0)
    end
    if isempty(obj_gradient)
        obj_gradient = DataFrame(string.(vars) .=> zeros(length(vars)))
    end

    # x0 initialization, and actual objective computation
    x0 = DataFrame(gm.solution_history[end, :])
    sol_vals = x0[:,string.(vars)]
    if !isnothing(obj_bbl)
        actual_cost = obj_bbl(x0)
        push!(gm.cost, actual_cost[1])
        x0[:, string(gm.objective)] = actual_cost
        sol_vals = x0[:,string.(vars)]
        feas_gap(gm, x0) # Checking feasibility gaps with updated objective
        append!(gm.solution_history, x0) # Pushing last solution
    end

    # Descent direction, counting and book-keeping
    d = @variable(gm.model, [1:length(vars)])
    var_diff = []
    max_diff = 0
    for key in vars
        push!(var_diff, maximum(bounds[key]) - minimum(bounds[key]))
        if !isinf(var_diff[end])
            max_diff = maximum([var_diff[end], max_diff])
        end
    end
    if any(isinf.(var_diff))
        @warn "Unbounded variables detected. " *
              "PGD may fail to converge to a local minimum."
        replace!(var_diff, Inf => max_diff)
    end
    replace!(var_diff, 0 => 1) # Sanitizing fixed values.

    ct = 0
    d_improv = 1e5
    prev_feas = false
    feas = is_feasible(gm)
    abstol = get_param(gm, :abstol)
    tighttol = get_param(gm, :tighttol)
    max_iterations = get_param(gm, :max_iterations)
    step_penalty = get_param(gm, :step_penalty)
    equality_penalty = get_param(gm, :equality_penalty)
    step_size = get_param(gm, :step_size)
    decay_rate = get_param(gm, :decay_rate)

    # WHILE LOOP
    @info("Starting projected gradient descent...")
    while (!feas || !prev_feas ||
        (prev_feas && feas && d_improv >= get_param(gm, :abstol)) ||
        ct == 0) && ct < max_iterations
        prev_obj = gm.cost[end]
        constrs = []
        ct += 1
        if step_penalty == 0
            step_penalty = maximum([1e3, 1e4*abs(prev_obj)])
        end
        if equality_penalty == 0
            equality_penalty = 100*step_penalty
        end

        # Linear objective gradient and constraints
        if !isnothing(obj_bbl)
            update_gradients(obj_bbl, [length(obj_bbl.Y)])
            obj_gradient = copy(grad_shell)
            append!(obj_gradient,
                DataFrame(obj_bbl.gradients[end,:]), cols = :subset)
            obj_gradient = coalesce.(obj_gradient, 0)
            # Update objective constraints
            push!(constrs, @constraint(gm.model, sum(Array(obj_gradient[end,:]) .* d) +
                              obj_bbl.dependent_var >= obj_bbl.Y[end]))
        end

        # Constraint evaluation
        errors = 0
        for bbl in bbls
            update_gradients(bbl, [length(bbl.Y)])
            constr_gradient = copy(grad_shell)
            append!(constr_gradient, DataFrame(bbl.gradients[end,:]), cols = :subset)
            constr_gradient = coalesce.(constr_gradient, 0)
            if bbl isa BlackBoxClassifier
                error_diff = maximum([tighttol, abs(bbl.Y[end])])
                if bbl.equality
                    push!(constrs, @constraint(gm.model, sum(Array(constr_gradient[end,:]) .* d)  + bbl.Y[end] + error_diff * bbl.relax_var >= 0))
                    push!(constrs, @constraint(gm.model, sum(Array(constr_gradient[end,:]) .* d)  + bbl.Y[end] <= error_diff * bbl.relax_var))
                    push!(constrs, @constraint(gm.model, bbl.relax_var <= 1))
                    errors += bbl.relax_var^2
                elseif !is_feasible(bbl) || (is_feasible(bbl) && bbl.feas_gap[end] <= 0)
                    push!(constrs, @constraint(gm.model, sum(Array(constr_gradient[end,:]) .* d)  + bbl.Y[end] + error_diff * bbl.relax_var >= 0))
                    push!(constrs, @constraint(gm.model, bbl.relax_var <= 1))
                    errors += bbl.relax_var^2
                else # feasible to zero tolerance
                    push!(constrs, @constraint(gm.model, sum(Array(constr_gradient[end,:]) .* d)  + bbl.Y[end] >= 0))
                end
            elseif bbl isa BlackBoxRegressor
                error_diff = maximum([tighttol, abs(bbl.optima[end] - bbl.actuals[end])])
                if bbl.equality
                    push!(constrs, @constraint(gm.model, sum(Array(constr_gradient[end,:]) .* d) + bbl.dependent_var + error_diff * bbl.relax_var >= bbl.Y[end]))
                    push!(constrs, @constraint(gm.model, sum(Array(constr_gradient[end,:]) .* d) + bbl.dependent_var <= bbl.Y[end] + error_diff * bbl.relax_var))
                    push!(constrs, @constraint(gm.model, bbl.relax_var <= 1))
                    errors += bbl.relax_var^2
                elseif !is_feasible(bbl) || (is_feasible(bbl) && bbl.feas_gap[end] <= 0)
                    push!(constrs, @constraint(gm.model, sum(Array(constr_gradient[end,:]) .* d)  + bbl.Y[end] + error_diff * bbl.relax_var >= 0))
                    push!(constrs, @constraint(gm.model, bbl.relax_var <= 1))
                    errors += bbl.relax_var^2
                else # feasible to zero tolerance
                    push!(constrs, @constraint(gm.model, sum(Array(constr_gradient[end,:]) .* d) + bbl.dependent_var >= bbl.Y[end]))
                end
            end
            if !isempty(bbl.lls)
                n_lls = length(bbl.lls)
                update_gradients(bbl, collect(length(bbl.Y) - n_lls:length(bbl.Y)-1))
                for i = 1:n_lls
                    ll = bbl.lls[i]
                    constr_gradient = copy(grad_shell)
                    append!(constr_gradient, # Changing the headers of the gradient
                        DataFrame(string.(ll.vars) .=> Array(bbl.gradients[end-n_lls-1+i,:])), cols = :subset)
                    constr_gradient = coalesce.(constr_gradient, 0)
                    Y_val = bbl.Y[end-n_lls-1+i]
                    if bbl isa BlackBoxClassifier
                        error_diff = maximum([tighttol, abs(Y_val)])
                        if bbl.equality
                            push!(constrs, @constraint(gm.model, sum(Array(constr_gradient[end,:]) .* d) + Y_val + error_diff * ll.relax_var >= 0))
                            push!(constrs, @constraint(gm.model, sum(Array(constr_gradient[end,:]) .* d) + Y_val <= error_diff * ll.relax_var))
                            push!(constrs, @constraint(gm.model, ll.relax_var <= 1))
                            errors += ll.relax_var^2
                        elseif !is_feasible(ll) || (is_feasible(ll) && ll.feas_gap[end] <= 0)
                            push!(constrs, @constraint(gm.model, sum(Array(constr_gradient[end,:]) .* d)  + Y_val + error_diff * ll.relax_var >= 0))
                            push!(constrs, @constraint(gm.model, ll.relax_var <= 1))
                            errors += ll.relax_var^2
                        else # feasible to zero tolerance
                            push!(constrs, @constraint(gm.model, sum(Array(constr_gradient[end,:]) .* d)  + Y_val >= 0))
                        end
                    elseif bbl isa BlackBoxRegressor
                        error_diff = maximum([tighttol, abs(ll.optima[end] - ll.actuals[end])])
                        if bbl.equality
                            push!(constrs, @constraint(gm.model, sum(Array(constr_gradient[end,:]) .* d) + ll.dependent_var + error_diff * ll.relax_var >= Y_val))
                            push!(constrs, @constraint(gm.model, sum(Array(constr_gradient[end,:]) .* d) + ll.dependent_var <= Y_val + error_diff * ll.relax_var))
                            push!(constrs, @constraint(gm.model, ll.relax_var <= 1))
                            errors += ll.relax_var^2
                        elseif !is_feasible(ll) || (is_feasible(ll) && ll.feas_gap[end] <= 0)
                            push!(constrs, @constraint(gm.model, sum(Array(constr_gradient[end,:]) .* d) + ll.dependent_var + error_diff * ll.relax_var >= Y_val))
                            push!(constrs, @constraint(gm.model, ll.relax_var <= 1))
                            errors += ll.relax_var^2
                        else # feasible to zero tolerance
                            push!(constrs, @constraint(gm.model, sum(Array(constr_gradient[end,:]) .* d) + ll.dependent_var >= Y_val))
                        end
                    end
                end
            end
        end

        # Initializing descent direction
        push!(constrs, @constraint(gm.model, d .== vars .- Array(sol_vals[end,:])))
        if feas # Small gradient step
            push!(constrs, @constraint(gm.model, sum((d ./ var_diff).^2) <=
                step_size/exp(decay_rate*(ct-1)/max_iterations)))
            @objective(gm.model, Min, sum(Array(obj_gradient[end,:]) .* d) + equality_penalty*errors)
        else    # Projection step
            @objective(gm.model, Min, sum(Array(obj_gradient[end,:]) .* d) + step_penalty*sum((d ./ var_diff).^2) + equality_penalty*errors)
        end

        # Optimizing the step, and finding next x0
        optimize!(gm.model)
        x0 = DataFrame(string.(vars) .=> getvalue.(vars))
        if !isnothing(obj_bbl)
            x0[!, string(gm.objective)] = obj_bbl(x0)
        end
        sol_vals = x0[:,string.(vars)]
        feas_gap(gm, x0)
        append!(gm.solution_history, x0)

        # Measuring improvements
        if !isnothing(obj_bbl)
            push!(gm.cost, gm.solution_history[end, string(gm.objective)])
        else
            push!(gm.cost,  JuMP.getvalue(gm.objective))
        end
        d_improv = gm.cost[end-1] - gm.cost[end]

        # Feasibility check
        prev_feas = feas
        feas = is_feasible(gm)

        # Saving solution dict for JuMP-style recovery.
        # TODO: do not regenerate soldict unless final solution.
        gm.soldict = Dict(key => JuMP.getvalue.(gm.model[key]) for (key, value) in gm.model.obj_dict)

        # Delete gradient constraints (TODO: perhaps add constraints to avoid cycling?)
        for con in constrs
            delete(gm.model, con)
        end
    end

    fincost = round(gm.cost[end], digits = -Int(round(log10(abstol))))
    # Termination criteria
    if ct >= max_iterations
        @info "Max iterations ($(ct)) reached."
        if prev_feas && feas
            @info "Solution is not converged to tolerance $(abstol)!"
        elseif (feas && !prev_feas) && (abs(d_improv) <= 100*abstol)
            @info "Solution is feasible and likely cycling, but the solution is close. Reduce step size and descend again. "
        elseif (!feas && prev_feas) && (abs(d_improv) <= 100*abstol)
            @info "Solution is infeasible and likely cycling, but the solution is close. Reduce step size and descend again. "
        elseif ((feas && !prev_feas) || (prev_feas && !feas))
            @info "Solution is likely cycling, with > $(abstol) changes in cost. Reduce step size and descend again."
        elseif !feas && !prev_feas
            @info "Solution is infeasible to tolerance $(tighttol)."
        end
        @info("Final cost is $(fincost).")
    else
        @info("PGD converged in $(ct) iterations!")
        @info("The optimal cost is $(fincost).")
    end

    # Reverting objective, and deleting vars
    @objective(gm.model, Min, gm.objective)
    delete(gm.model, d)
    return
end

""" Complete solution procedure for GlobalModel. """
function globalsolve!(gm::GlobalModel)
    @info "GlobalModel " * gm.name * " solution in progress..."
    for bbl in gm.bbls
        if !is_sampled(bbl)
            uniform_sample_and_eval!(bbl)
        end
    end
    set_param(gm, :ignore_accuracy, true)
    set_param(gm, :ignore_feasibility, true)
    @info "Training OptimalTreeLearners..."
    for bbl in gm.bbls
        if isempty(bbl.learners)
            learn_constraint!(bbl)
        end
        if isempty(bbl.mi_constraints)
            add_tree_constraints!(gm, bbl)
        end
    end
    @info "Solving MIP..."
    optimize!(gm)
    descend!(gm)
end

function globalsolve_and_time!(m::Union{JuMP.Model, GlobalModel})
    t1 = time()
    if m isa JuMP.Model
        optimize!(m);
    else
        globalsolve!(m)
    end
    println("Model solution time: " * string(time()-t1) * " seconds.")
end


text:
function MCTS.RandomBandit(mdp::VDPTagMDP; kwargs...)
    RandomBandit(;
        enable_action_pw=true,
        check_repeat_action=true,
        kwargs...)
end
function MCTS.DPWBandit(mdp::VDPTagMDP; kwargs...)
    DPWBandit(;
        enable_action_pw=true,
        check_repeat_action=true,
        exploration_constant=0.75,
        k_action=10.0,
        alpha_action=0.5,
        kwargs...)
end
function MCTS.CBTSBandit(mdp::VDPTagMDP; kwargs...)
    CBTSBandit(;
        enable_action_pw=true,
        check_repeat_action=true,
        exploration_constant=1.0,
        A_max=20,
        n_proposes=100,
        log_length_scale=0.0,
        log_signal_sigma=0.0,
        log_obs_noise=-1.0,
        action_dims=1,
        n_sig=2.0,  #number of standard deviations for GP-UCB
        kwargs...)
end
function MCTS.CBTSDPWBandit(mdp::VDPTagMDP; kwargs...)
    CBTSDPWBandit(;
        enable_action_pw=true,
        check_repeat_action=true,
        exploration_constant=2.0,
        k_action=10.0,
        alpha_action=0.5,
        n_proposes=500,
        log_length_scale=0.0,
        log_signal_sigma=0.0,
        log_obs_noise=-1.0,
        fit_qs=false,
        action_dims=1,
        n_sig=2.0,  #number of standard deviations for GP-UCB
        kwargs...)
end
function MCTS.ModularSolver(mdp::VDPTagMDP, b::ModularBandit, seed=0;
                            kwargs...)
    rng = MersenneTwister(seed)
    ModularSolver(;
        n_iterations=100,
        bandit=b,
        depth=10,
        k_state=2.0,
        alpha_state=0.9,
        keep_tree=true,
        check_repeat_state=true,
        rng=rng,
        init_Q=0.0,
        init_N=0,
        next_action=RandomActionGenerator(rng),
        estimate_value=RolloutEstimator(RandomSolver(rng)),
        track_best_path=true,
        kwargs...
       )
end
Base.string(::VDPTagMDP) = "VDPTag"


text:
# Author: Lukasz Bratos

# Odkomentować przy pierwszym użyciu
# import Pkg
# Pkg.add("Plots")
# Pkg.add("PyPlot")
using Plots

# Nasze dane testowe do zadania
consts = [-2, -2, -2, -1, -1, -1, -1]
x0s = [1, 2, 1.99999999999999, 1, -1, 0.75, 0.25]

# Funkcja obliczajaca wartosc wyrazenia z zadania
# n - liczba iteracji, c - pewna stała, x0 - nasz argument
function func(n, c, x0)
    arr = zeros(Float64, 0)
    xn = x0
    push!(arr, xn)
    for i in 1:n
        xn = xn ^ 2 + c
        push!(arr, xn)
    end
    return arr
end

x = 0:40
for i in 1:7
    println("Wynik dla c = ", consts[i], ", x0 = ", x0s[i], ": \n", func(40, Float64(consts[i]), Float64(x0s[i])))
    y = func(40, Float64(consts[i]), Float64(x0s[i]))
    pyplot()
    savefig(scatter(x, y, title = "$i.", xlabel = "x", ylabel = "y"), "plot$i.png")
end


text:
abstract GLM

@doc "Representation of a generalized linear model using SharedArrays" ->
type PGLM{T<:AbstractFloat,D<:UnivariateDistribution,L<:Link} <: GLM
    Xt::SharedMatrix{T}                 # transposed model matrix
    XtWX::SharedArray{T,3}
    XtWr::SharedMatrix{T}
    wt::SharedVector{T}                 # prior case weights
    y::SharedVector{T}                  # observed response vector
    β::Vector{T}                        # base value of β
    βs::Vector{T}                       # value of β + s*δβ
    δβ::Vector{T}                       # increment
    η::SharedVector{T}                  # current linear predictor vector
    μ::SharedVector{T}                  # current mean vector
    dev::SharedVector{T}
    d::D
    l::L
    fit::Bool
end

scvt(a::Array) = convert(SharedArray,a)

function PGLM{T<:AbstractFloat}(Xt::SharedMatrix{T},
                                y::SharedVector{T},
                                wt::SharedVector{T},
                                d::UnivariateDistribution,
                                l::Link)
    p,n = size(Xt)
    n == length(y) == length(wt) || throw(DimensionMismatch(""))
    pr = procs(y)
    Set(pr) == Set(procs(Xt)) == Set(procs(wt)) || error("SharedArrays must have same procs")
    β = zeros(p)
    ntot = maximum(pr)
    g = PGLM(Xt,scvt(similar(y,(p,p,ntot))),scvt(similar(y,(p,ntot))),wt,y,
             β,copy(β),copy(β),scvt(similar(y)),scvt(similar(y)),scvt(similar(y,(ntot,))),
             d,l,false)
    initμη!(g)
    updateXtW!(g)
    copy!(g.β,g.δβ)
    g
end
function PGLM{T<:AbstractFloat}(Xt::SharedMatrix{T},y::SharedVector{T},
                                d::UnivariateDistribution,l::Link)
    PGLM(Xt,y,scvt(ones(y)),d,l)
end
function PGLM{T<:AbstractFloat}(Xt::SharedMatrix{T},
                                y::SharedVector{T},
                                d::UnivariateDistribution)
    PGLM(Xt,y,d,canonical(d))
end
function PGLM{T<:AbstractFloat}(Xt::Matrix{T},y::Vector{T},d::UnivariateDistribution)
    PGLM(scvt(Xt),scvt(y),d)
end

type SGLM{T<:AbstractFloat,D<:UnivariateDistribution,L<:Link} <: GLM
    Xt::Matrix{T}                    # transposed model matrix
    XtW::Matrix{T}                   # weighted, transposed model matrix
    XtWX::Matrix{T}                  # weighted cross-product
    wt::Vector{T}                    # prior case weights
    wtres::Vector{T}                 # weighted residuals
    y::Vector{T}                     # observed response vector
    β::Vector{T}                     # base value of β
    βs::Vector{T}                    # value of β + s*δβ
    δβ::Vector{T}                    # increment
    η::Vector{T}                     # current linear predictor vector
    μ::Vector{T}                     # current mean vector
    d::D
    l::L
    blas::Bool                       # use BLAS.syrk! to evaluate XtWX
    fit::Bool
end
function SGLM{T<:AbstractFloat}(Xt::Matrix{T},
                                y::Vector{T},
                                wt::Vector{T},
                                d::UnivariateDistribution,
                                l::Link,
                                blas::Bool)
    p,n = size(Xt)
    XtW = blas ? similar(Xt) : similar(Xt,(0,0))
    n == length(y) == length(wt) || throw(DimensionMismatch(""))
    β = Array(T,(p,))
    g = SGLM(Xt,XtW,zeros(T,(p,p)),wt,similar(y),y,β,Array(T,(p,)),Array(T,(p,)),
             similar(y),similar(y),d,l,blas,false)
    initμη!(g)
    updateXtW!(g)
    copy!(g.β,g.δβ)
    g
end
function SGLM{T<:AbstractFloat}(Xt::Matrix{T},y::Vector{T},
                                d::UnivariateDistribution,l::Link,blas::Bool)
    SGLM(Xt,y,fill!(similar(y),one(T)),d,l,blas)
end
function SGLM{T<:AbstractFloat}(Xt::Matrix{T},
                                y::Vector{T},
                                d::UnivariateDistribution,
                                blas::Bool)
    SGLM(Xt,y,d,canonical(d),blas)
end
function SGLM{T<:AbstractFloat}(Xt::Matrix{T},
                                y::Vector{T},
                                d::UnivariateDistribution)
    SGLM(Xt,y,d,canonical(d),true)
end


text:
# we define our own RegexAndMatchData, similar to the definition in Base
# but we call create_match_data once in __init__
mutable struct RegexAndMatchData
    re::Regex
    match_data::Ptr{Cvoid}
    RegexAndMatchData(re::Regex) = new(re) # must create_match_data in __init__
end

function initialize!(re::RegexAndMatchData)
    re.match_data = Base.PCRE.create_match_data(re.re.regex)
    return
end

@static if isdefined(Base, :RegexAndMatchData)

"""
Execute a regular expression without the overhead of `Base.Regex`
"""
exec(re::RegexAndMatchData, bytes, offset::Int=1) =
    Base.PCRE.exec(re.re.regex, bytes, offset-1, re.re.match_options, re.match_data)

"""
`SubString` containing the bytes following the matched regular expression.
"""
nextbytes(re::RegexAndMatchData, bytes) = SubString(bytes, unsafe_load(Base.PCRE.ovec_ptr(re.match_data), 2) + 1)

"""
`SubString` containing a regular expression match group.
"""
function group(i, re::RegexAndMatchData, bytes)
    p = Base.PCRE.ovec_ptr(re.match_data)
    SubString(bytes, unsafe_load(p, 2i+1) + 1, prevind(bytes, unsafe_load(p, 2i+2) + 1))
end

function group(i, re::RegexAndMatchData, bytes, default)
    p = Base.PCRE.ovec_ptr(re.match_data)
    return unsafe_load(p, 2i+1) == Base.PCRE.UNSET ? default :
        SubString(bytes, unsafe_load(p, 2i+1) + 1, prevind(bytes, unsafe_load(p, 2i+2) + 1))
end

else # old Regex style

exec(re::RegexAndMatchData, bytes, offset::Int=1) = Base.PCRE.exec(re.re.regex, bytes, offset-1, re.re.match_options, re.re.match_data)
nextbytes(re::RegexAndMatchData, bytes) = SubString(bytes, re.re.ovec[2]+1)
group(i, re::RegexAndMatchData, bytes) = SubString(bytes, re.re.ovec[2i+1]+1, prevind(bytes, re.re.ovec[2i+2]+1))
group(i, re::RegexAndMatchData, bytes, default) = re.re.ovec[2i+1] == Base.PCRE.UNSET ? default : SubString(bytes, re.re.ovec[2i+1]+1, prevind(bytes, re.re.ovec[2i+2]+1))

end


text:
include("../src/components/RCPSSPScenario.jl")
include("../src/components/AbatementCostParameters.jl")
include("../src/components/AbatementCosts.jl")

using DataFrames
using Test

for testscen in 1:2
    valdir, scenario, use_permafrost, use_seaice = get_scenario(testscen)
    println(scenario)

    m = page_model()
    scenario = addrcpsspscenario(m, scenario)

    for gas in [:CO2, :CH4, :N2O, :Lin]
        abatementcostparameters = addabatementcostparameters(m, gas)
        abatementcosts = addabatementcosts(m, gas)

        abatementcostparameters[:yagg] = readpagedata(m,"test/validationdata/yagg_periodspan.csv")
        abatementcostparameters[:cbe_absoluteemissionreductions] = abatementcosts[:cbe_absoluteemissionreductions]

        abatementcosts[:zc_zerocostemissions] = abatementcostparameters[:zc_zerocostemissions]
        abatementcosts[:q0_absolutecutbacksatnegativecost] = abatementcostparameters[:q0_absolutecutbacksatnegativecost]
        abatementcosts[:blo] = abatementcostparameters[:blo]
        abatementcosts[:alo] = abatementcostparameters[:alo]
        abatementcosts[:bhi] = abatementcostparameters[:bhi]
        abatementcosts[:ahi] = abatementcostparameters[:ahi]
        if gas == :Lin
            abatementcosts[:er_emissionsgrowth] = scenario[:er_LGemissionsgrowth]
        else
            abatementcosts[:er_emissionsgrowth] = scenario[Symbol("er_" * String(gas) * "emissionsgrowth")]
        end
    end

    p = load_parameters(m)
    p["y_year_0"] = 2015.
    p["y_year"] = Mimi.dim_keys(m.md, :time)
    set_leftover_params!(m, p)

    run(m)

    @test !isnan(m[:AbatementCostsCO2, :tc_totalcost][10, 5])
    @test !isnan(m[:AbatementCostsCH4, :tc_totalcost][10, 5])
    @test !isnan(m[:AbatementCostsN2O, :tc_totalcost][10, 5])
    @test !isnan(m[:AbatementCostsLin, :tc_totalcost][10, 5])

    #compare output to validation data
    zc_compare_co2=readpagedata(m, "test/validationdata/zc_zerocostemissionsCO2.csv")
    zc_compare_ch4=readpagedata(m, "test/validationdata/zc_zerocostemissionsCH4.csv")
    zc_compare_n2o=readpagedata(m, "test/validationdata/zc_zerocostemissionsN2O.csv")
    zc_compare_lin=readpagedata(m, "test/validationdata/zc_zerocostemissionsLG.csv")

    @test m[:AbatementCostParametersCO2, :zc_zerocostemissions] ≈ zc_compare_co2 rtol=1e-2
    @test m[:AbatementCostParametersCH4, :zc_zerocostemissions] ≈ zc_compare_ch4 rtol=1e-3
    @test m[:AbatementCostParametersN2O, :zc_zerocostemissions] ≈ zc_compare_n2o rtol=1e-3
    @test m[:AbatementCostParametersLin, :zc_zerocostemissions] ≈ zc_compare_lin rtol=1e-3

    tc_compare_co2=readpagedata(m, "test/validationdata/$valdir/tc_totalcosts_co2.csv")
    tc_compare_ch4=readpagedata(m, "test/validationdata/$valdir/tc_totalcosts_ch4.csv")
    tc_compare_n2o=readpagedata(m, "test/validationdata/$valdir/tc_totalcosts_n2o.csv")
    tc_compare_lin=readpagedata(m, "test/validationdata/$valdir/tc_totalcosts_linear.csv")

    @test m[:AbatementCostsCO2, :tc_totalcost] ≈ tc_compare_co2 rtol=1e-2
    @test m[:AbatementCostsCH4, :tc_totalcost] ≈ tc_compare_ch4 rtol=1e-2
    @test m[:AbatementCostsN2O, :tc_totalcost] ≈ tc_compare_n2o rtol=1e-2
    @test m[:AbatementCostsLin, :tc_totalcost] ≈ tc_compare_lin rtol=1e-2
end

max_stars_repo_name: jkikstra/PAGE-VAR

max_issues_repo_name: jkikstra/PAGE-VAR

max_forks_repo_name: jkikstra/PAGE-VAR


text:
julia> @eval macro $(:do)(block, when::Symbol, condition)
           when ≠ :when && error("@do expected `when` got `$s`")
           quote
               let
                   $block
                   while $condition
                       $block
                   end
               end
           end |> esc
       end
@do (macro with 1 method)

julia> i = 0
0

julia> @do begin
           @show i
           i += 1
       end when i % 6 ≠ 0
i = 0
i = 1
i = 2
i = 3
i = 4
i = 5

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using CSV
using Test
using DecisionTrees

elt = [Int,	Float64,	Float64,	Float64,	Float64,	Float64,	Int,	String,	String,	String,	Int,	String,	String,	String,	String,	String,	String,	String,	String,	String,	String,	Int,	String,	Int,	String,	String,	Int,	String,	String,	Int,	String,	String,	Int,	String,	String,	Int,	String,	Int,	String,	String,	Int,	Int,	String,	Int,	Int,	String,	Int,	Int,	String,	Int,	String,	String,	String,	Int,	String,	String,	String,	String,	Int,	Int,	Int,	String,	String,	String,	String,	String,	Int,	Int,	Int,	Int,	Int,	Int,	Int,	Int,	Int,	Int,	Int,	Int,	Int,	Float64,	Int,	Float64,	Float64,	Float64,	Float64,	Int,	Int,	Int,	Int,	Int,	Int]

fi = "data1Medium.csv"
fi = "data1Small.csv"
@time df_tmp = CSV.read(joinpath(datadir, "GermanMotorPremiums", fi), allowmissing=:none, types=elt, categorical=false, rows_for_type_detect=10000);

selected_explanatory_vars = ["PLZ_WOHNORT","ART_DES_WOHNEIGENTUM","GEBURTSDATUM","FAMILIENSTAND","NATIONALITAET","GESCHLECHT","FINANZIERUNGSART","STADT","KENNZEICHEN"]

# for x in selected_explanatory_vars
#    @show x,eltype(df_tmp[Symbol(x)])
# end

dtmtable, sett, df_prepped = prepare_dataframe_for_dtm!(df_tmp, treat_as_categorical_variable=["PLZ_WOHNORT"], weightcol="EXPOSURE", numcol="LOSS20HALF", denomcol="PREMIUM66", independent_vars=selected_explanatory_vars);
sett.minWeight = -.2

using StatsBase
using BenchmarkTools

v = df_prepped[:PLZ_WOHNORT]
# version 1
@benchmark uq1 = StatsBase.countmap(v)

# version 2
function me(y)
    u = unique(y)
    d = Dict([(i, count(x->x == i, y)) for i in u])
    return d
end
@benchmark uq2 = me(v) # this is about 40% faster than countmap. Then again

# sorting / cutoff
function frequencyRel(x)
    d = proportionmap(x)
    srt = sortperm(collect(values(d)))
    return d, srt
end


text:
@testset "ridge-reg" begin
    rng = StableRNG(622161)
    n, p = 100, 5
    X = randn(rng, n, p)
    y = randn(rng, n)
    X1 = R.augment_X(X, true)

    λ = 0.3

    Xt = MLJBase.table(X)
    rr = RidgeRegressor(lambda=λ, penalize_intercept=true,
                        scale_penalty_with_samples = false)
    fr, = MLJBase.fit(rr, 1, Xt, y)
    ŷ = MLJBase.predict(rr, fr, Xt)

    θ = (X1'X1 + λ*I)\(X1'y)
    coefs = θ[1:end-1]
    intercept = θ[end]

    fp = MLJBase.fitted_params(rr, fr)
    @test last.(fp.coefs) ≈ coefs
    @test fp.intercept ≈ intercept
end

@testset "logistic" begin
    ((X, y, θ), (X1, y1, θ1)) = generate_binary(100, 5)

    λ = 0.7
    γ = 0.1

    Xt = MLJBase.table(X)
    yc = MLJBase.categorical(y1)

    lr = LogisticClassifier(lambda=λ, gamma=γ, scale_penalty_with_samples = false)
    fr, = MLJBase.fit(lr, 1, Xt, yc)

    ŷ = MLJBase.predict(lr, fr, Xt)
    ŷ = MLJBase.mode.(ŷ)

    mcr = MLJBase.misclassification_rate(ŷ, yc)
    @test mcr ≤ 0.2
end

@testset "multinomial" begin
    ((X, y, θ), (X1, y1, θ1)) = generate_multiclass(100,  5, 3)

    λ = 0.5
    γ = 0.2

    Xt = MLJBase.table(X)
    yc = MLJBase.categorical(y1)

    mc = MultinomialClassifier(lambda=λ, gamma=γ, fit_intercept=false,
                               scale_penalty_with_samples = false)
    fr, = MLJBase.fit(mc, 1, Xt, yc)

    mach = MLJBase.machine(mc, Xt, yc)
    MLJBase.fit!(mach)
    fp = MLJBase.fitted_params(mach)
    @test length(fp.coefs) == 5

    ŷ = MLJBase.predict(mc, fr, Xt)
    ŷ = MLJBase.mode.(ŷ)

    mcr = MLJBase.misclassification_rate(ŷ, yc)
    @test mcr ≤ 0.3
end

# see issue https://github.com/alan-turing-institute/MLJ.jl/issues/387
@testset "String-Symbol" begin
    model = LogisticClassifier(penalty="l1")
    @test model.penalty == "l1"
    gr = MLJLinearModels.glr(model, 2)
    @test gr isa GLR
    @test gr.penalty isa ScaledPenalty{L1Penalty}
end

# see issue #71
@testset "Logistic-m" begin
    X, y = MLJBase.make_blobs(centers=3)
    model = LogisticClassifier()
    mach = MLJBase.machine(model, X, y)
    MLJBase.fit!(mach)
    fp = MLJBase.fitted_params(mach)
    @test unique(fp.classes) == [1,2,3]
end


text:
using Random, Distributions, Plots, LaTeXStrings; pyplot()
Random.seed!(0)

actualAlpha, actualLambda = 2,3
gammaDist = Gamma(actualAlpha,1/actualLambda)
n = 10^2
sample = rand(gammaDist, n)

alphaGrid = 1:0.02:3
lambdaGrid = 2:0.02:5

likelihood = [prod([pdf.(Gamma(a,1/l),v) for v in sample])
                        for l in lambdaGrid, a in alphaGrid]

surface(alphaGrid, lambdaGrid, likelihood, lw=0.1,
	c=cgrad([:blue, :red]), legend=:none, camera = (135,20),
	xlabel=L"\alpha", ylabel=L"\lambda", zlabel="Likelihood")

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
#=
MIT License

Copyright (c) 2021 David-W-Wallis

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.
=#

module TwoDerrors

import Printf
import JSON
import Random
import Statistics
import JLD
import PyPlot

const C = 330

# Sensor positions
sensorA=(10,-5,0,"A")
sensorB=(40,-5,0,"B")
sensorC=(-5,45,-90,"C")
sensorsAB=(sensorA,sensorB)

# For constructing file names.
sensor2string(s) = Printf.@sprintf("%+04d%+04d%+04d",s[1],s[2],s[3])
sensors2string(s) = join(map(sensor2string,s))

# Angle from microphone array to a point on the modelled area.
theta(x,y,sensorPos) = atan.(x-sensorPos[1], y-sensorPos[2])

# Equation 8
function phi(theta, ifParams)
    f,d,C=ifParams
    2*pi*f*d.*(sin.(theta)./C)
end

# Equation 8
function phi(theta,f,d)
    return 2*pi*f*d/C * sin(deg2rad(theta))
end

function errorMapFileName(sensors,d,f)
    fmt=Printf.@sprintf("jld/%s_%02d_%05d",sensors2string(sensors),
                        convert(Int64,d*1000),convert(Int64,f))
    return "errormap_$fmt.jld"
end

coeffsFileName(froot,f,d) = return Printf.@sprintf("json/%s_%05d_%02d_%02d.json",
                                                   froot,
                                                   Int(floor(f)),
                                                   Int(d*1000),
                                                   5)

# Reads pre-calculated model coefficients from json files.
function makeCoeffs(froot)
    dict=Dict()
    for d in 1:200
        for f in [2000,4000,6000,8000]#,4000,6000,8000]
            di=Dict()
            open(coeffsFileName(froot,f,d/1000.0), "r") do f
                dicttxt = read(f,String)  # file information to string
                di=JSON.parse(dicttxt)
            end
            dict[(f,d)]=di
        end
    end
    return dict
end

BiasCoeffs=makeCoeffs("biascoeffs")
ErrorBarCoeffs=makeCoeffs("polycoeffs")

# Equation 9
# Takes radians returns degrees
function biasModel(theta,f,d)
    c=BiasCoeffs[(f,d)]
    alpha=c["alpha"]
    beta=c["beta"]
    gamma=c["gamma"]
    model = alpha/(cos(gamma*theta)^beta) - alpha
    return model*sign(theta)
end

# Equation 10
# Takes radians and returns degrees
function errorBarModel(theta,f,d)
    c=ErrorBarCoeffs[(f,d)]["coeffs"]
    return c[1]*theta^4 + c[2]*theta^3 + c[3]*theta^2 + c[4]*theta + c[1]
end

# Purturb the angle
# Takes radians returns degrees
function newAngleRad(theta,f,d)
    bias=biasModel(theta,f,d)
    err=errorBarModel(theta,f,d)
    #println("Bias: ",bias," error: ",err)
    return rad2deg(theta)-bias + err*Random.rand() - err/2
end

# Purturb the angle
# Takes degrees returns degrees
function newAngle(theta,f,d)
    return newAngleRad(deg2rad(theta),f,d)
end

# Generates and saves an error map
function errorMapN(sensors;f=2000.0,d=0.020,M=1000)
    N=length(sensors)
    GZ=51
    ifparams=(f,d,330)
    posx=Array{Float64}(undef,GZ,GZ)
    posy=Array{Float64}(undef,GZ,GZ)
    errx=Array{Float64}(undef,GZ,GZ)
    erry=Array{Float64}(undef,GZ,GZ)
    th_=zeros(N)
    th=Array{Float64}(undef,N)
    phi_=Array{Float64}(undef,N)
    slopeErr=zeros(N)
    for x = 0:GZ-1
        for y = 0:GZ-1
            for i = 1:N
                th_[i] = theta(x,y,sensors[i])
            end
            xp=zeros(M)
            yp=zeros(M)
            for test = 1:M
                th=newAngle.(th_,Int(floor(f)),Int(floor(d*1000)))
                if false
                    println("x: ",x," y: ",y)
                    println("Angle:")
                    println(map(x -> Printf.@sprintf("%.1f",x), rad2deg.(th_)))
                    println("Error:")
                    println(map(x -> Printf.@sprintf("%.1f",x), slopeErr))
                    println("New angle:")
                    println(map(x -> Printf.@sprintf("%.1f",x), rad2deg.(th)))
                    println()
                end

                A=hcat(tan.(th),ones(N))
                b=[sensors[i][1]*tan(th[i]) + sensors[i][2] for i in 1:N]
                if N > 2
                    m=GLM.lm(A,b)
                    xp[test],yp[test]=GLM.coef(m)
                else
                    xp[test],yp[test]=A\b
                end
            end
            errx[x+1,y+1]=Statistics.std(xp)
            erry[x+1,y+1]=Statistics.std(yp)
            posx[x+1,y+1]=Statistics.mean(xp)
            posy[x+1,y+1]=Statistics.mean(yp)
        end
    end
    fname=errorMapFileName(sensors,d,f)
    JLD.save(fname,"errx",errx,"erry",erry,"posx",posx,"posy",posy)
end

# Load a pre-calculated map
function loadMap(sensors,f,d)
    fname=errorMapFileName(sensors,d,f)
    return JLD.load(fname)#,"errx",errx,"erry",erry,"posx",posx,"posy",posy)
end

function plotErrorMapN(sensors,f,d;colour=false,levelsx=false,levelsy=false,PLOTSENSORS=false,contcolours="none",cmap="gist_rainbow",vmin=0,vmax=15)
    println("cmap=",cmap)
    emap=loadMap(sensors,f,d)
    errx=emap["errx"]
    erry=emap["erry"]
    posx=emap["posx"]
    posy=emap["posy"]
    fmt=Printf.@sprintf("%s_%02d_%04d",sensors2string(sensors),
                        convert(Int64,d*1000),convert(Int64,f))

    fint=convert(Int64, f)
    fig1=PyPlot.figure(1)
    fig1.clf()
    if colour
        if levelsx==false
            PyPlot.contourf(errx',cmap=cmap,vmin=vmin,vmax=vmax)
        else
            PyPlot.contourf(errx',cmap=cmap,levels=levelsx,vmin=vmin,vmax=vmax)
        end
    end

    if colour
        if contcolours == "none"
            contcolours="white"
        end
    else
        contcolours="black"
    end

    if levelsx==false
        CS = PyPlot.contour(errx',colors=contcolours)
    else
        CS = PyPlot.contour(errx',colors=contcolours,levels=levelsx)
    end

    PyPlot.clabel(CS)
    PyPlot.xlabel("Distance (meters)")
    PyPlot.ylabel("Distance (meters)")
    PyPlot.savefig("mc_error_x_$fmt.png")

    fig2=PyPlot.figure(2)
    fig2.clf()
    if colour
        if levelsy==false
            PyPlot.contourf(erry',cmap=cmap)
        else
            PyPlot.contourf(erry',cmap=cmap,levels=levelsy)
        end
    end

    if levelsy==false
        CS = PyPlot.contour(erry',colors=contcolours)
    else
        CS = PyPlot.contour(erry',colors=contcolours,levels=levelsy)
    end

    if PLOTSENSORS
        for s=sensors
            x,y=s[1],s[2]
            PyPlot.plot(x,y,"ko")
            PyPlot.annotate(s[4],xy=(x+2,y))
        end
        PyPlot.ylim([-7,50])
        PyPlot.xlim([-7,50])
    end

    PyPlot.clabel(CS)
    PyPlot.xlabel("Distance (meters)")
    PyPlot.ylabel("Distance (meters)")
    PyPlot.savefig("mc_error_y_$fmt.png")

    (errx,erry,posx,posy)
end

function paperPlot(sensors,f,d;M=600,colour=false,contcolours="black",levs=[0,1,1.5,2,3,4,5,6,8,10,15],cmap="gist_rainbow",vmin=0,vmax=15,calculate=true)
    if calculate
        errorMapN(sensors,f=f,d=d,M=M)#,levelsx=mylevs,levelsy=mylevs)#,contcolours="black")
    end
    plotErrorMapN(sensors,f,d;colour=colour,levelsx=levs,levelsy=levs,vmin=vmin,vmax=vmax,cmap=cmap,contcolours=contcolours)
    fig=PyPlot.figure(1)
    froot="$(sensors[1][4])$(sensors[2][4])$(Int(f))-$(Int(d*1000))"
    println(froot)
    PyPlot.savefig(froot*"x.png")
    fig=PyPlot.figure(2)
    PyPlot.savefig(froot*"y.png")
end

# For interactive testing
function ppRun(f,d;M=1000)
    paperPlot(TwoDerrors.sensorsAB,f,d;M=M,levs=[0,1,1.5,2,3,4,5,6,8,10,15],colour=true,cmap="jet",calculate=true)
end

# Plots for paper
function makePlots(;M=1000)
    for (f,d) in zip((2000,8000),(0.008,0.012,0.020))
        paperPlot(TwoDerrors.sensorsAB,f,d;M=M,levs=[0,1,1.5,2,3,4,5,6,8,10,15],colour=true,cmap="jet",calculate=true)
    end
    paperPlot(TwoDerrors.sensorsAB,2000,0.060;M=M,levs=[0,1,1.5,2,3,4,5,6,8,10,15],colour=true,cmap="jet",calculate=true)
end

# To generate the area coverage.
percentUnder(a,thresh) = (sum(a.<thresh) / length(a)) * 100

# To generate the area coverage.
function percentUnder(sensors,f,d,thresh)
    emap=loadMap(sensors,f,d)
    errx=emap["errx"]
    erry=emap["erry"]
    percentUnder(errx,thresh),percentUnder(erry,thresh)
end

end

max_stars_repo_name: David-W-Wallis/TwoDErrors

max_issues_repo_name: David-W-Wallis/TwoDErrors

max_forks_repo_name: David-W-Wallis/TwoDErrors


text:
using DiffEqBase: @..
using DiffEqBase
using Test
using InteractiveUtils
using MuladdMacro
function foo9(a, b, c, d, e, f, g, h, i)
    a = DiffEqBase.diffeqbc(a)
    @.. a = b + 0.1 * (0.2c + 0.3d + 0.4e + 0.5f + 0.6g + 0.6h + 0.6i)
    nothing
end
function foo26(a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t, u, v, w, x, y, z)
    @muladd @.. a = b + 0.1 * (0.2c + 0.3d + 0.4e + 0.5f + 0.6g + 0.6h + 0.6i + 0.6j + 0.6k + 0.6l + 0.6m + 0.6n + 0.6o + 0.6p + 0.6q + 0.6r + 0.6s + 0.6t + 0.6u + 0.6v + 0.6w + 0.6x + 0.6y + 0.6z)
    nothing
end
a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t, u, v, w, x, y, z = [rand(100, 100, 2) for i in 1:26];
@allocated foo9(a, b, c, d, e, f, g, h, i)
@allocated foo26(a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t, u, v, w, x, y, z)
if Base.JLOptions().check_bounds != 1
  @test @allocated(foo9(a, b, c, d, e, f, g, h, i)) == 0
  @test @allocated(foo26(a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t, u, v, w, x, y, z)) == 0
  io = IOBuffer()
  InteractiveUtils.code_llvm(io, foo9, Base.typesof(a, b, c, d, e, f, g, h, i))
  str = String(take!(io))
  @test occursin("vector.body", str)
end

fill!(a, 0.1);
@test @..(a = b + c ) == @. b + c
fill!(a, 0.1); a_ = copy(a);
@test @..(a += b + c) == @. a_ + (b + c)
fill!(a, 0.1); a_ = copy(a);
@test @..(a -= b + c) == @. a_ - (b + c)
fill!(a, 0.1); a_ = copy(a);
@test @..(a *= b + c) == @. a_ * (b + c)
fill!(a, 0.1); a_ = copy(a);
@test @..(a /= b + c) == @. a_ / (b + c)
fill!(a, 0.1); a_ = copy(a);
@test @..(a \= b + c) == @. a_ \ (b + c)
fill!(a, 0.1); a_ = copy(a);
@test @..(a ^= b + c) == @. a_ ^ (b + c)

a, b, c = [rand(Bool, 10) for i in 1:3]; a_ = copy(a)
@test @..(a &= b | c) == @. a_ & (b | c)

x = ones(10); y = similar(x)
@test_nowarn @.. y = 0.1
@test_nowarn @.. y = x


text:
# Exercise 1
using LinearAlgebra, Random
using BenchmarkTools

N = 50
A = rand(N, N)
B = rand(N, N)
@btime $A * $B

# loop by the row first
function rowloop(a, b)
    @assert size(a)[2] == size(b)[1]
    c = zeros(size(a)[1], size(b)[2])
    for j in 1:size(b)[2]
        for i in 1:size(a)[1]
            c[i, j] = a[i, :]' * b[:, j]
        end
    end
    return c
end
@btime rowloop($A, $B)

# loop by the column first
function colloop(a, b)
    @assert size(a)[2] == size(b)[1]
    c = zeros(size(a)[1], size(b)[2])
    for i in 1:size(a)[1]
        for j in 1:size(b)[2]
            c[i, j] = a[i, :]' * b[:, j]
        end
    end
    return c
end
@btime colloop($A, $B)

# loop by the dot product
function rowloopdot(a, b)
    @assert size(a)[2] == size(b)[1]
    c = zeros(size(a)[1], size(b)[2])
    for j in 1:size(b)[2]
        for i in 1:size(a)[1]
            c[i, j] = dot(a[i, :], b[:, j])
        end
    end
    return c
end
@btime rowloopdot($A, $B)

function colloopdot(a, b)
    @assert size(a)[2] == size(b)[1]
    c = zeros(size(a)[1], size(b)[2])
    for i in 1:size(a)[1]
        for j in 1:size(b)[2]
            c[i, j] = dot(a[i, :], b[:, j])
        end
    end
    return c
end
@btime colloopdot($A, $B)

# Compare colloopdot(a, b) to the built-in function
A10 = rand(10, 10)
B10 = rand(10, 10)
A1000 = rand(1000, 1000)
B1000 = rand(1000, 1000)
@btime $A10 * $B10
@btime colloopdot($A10, $B10)
@btime $A1000 * $B1000
@btime colloopdot($A1000, $B1000)

max_stars_repo_path: Problem Sets/PS4.jl

max_issues_repo_path: Problem Sets/PS4.jl

max_forks_repo_path: Problem Sets/PS4.jl


text:
using Plots
using LaTeXStrings
pgfplots()

# function save(p, fn)
#     savefig(p, fn)
#     open(fn, "r+") do io
#         contents = read(io)
#         seekstart(io)
#         preamble = "\\documentclass[tikz]{standalone}\n"*"\\usepackage{standalone}\n"*
#         "\\usepackage{tikz,pgfplots}\n"*"\\pgfplotsset{compat=newest}\n"*
#         "\\pgfplotsset{every axis legend/.append style={cells={anchor=west}}}\n"*
#         "\\usetikzlibrary{arrows}\n"*"\\tikzset{>=stealth'}\n"*"\\begin{document}\n"
#         write(io, preamble)
#         write(io, contents)
#         seekend(io)
#         write(io, "\\end{document}")
#     end
# end

P_A(n; p=0.9965) = 1/(1-p + p/n)
P_G(n; p=0.9965) = 1 - p + n*p

ps = [0.5, 0.75, 0.9, 0.95]

pa = plot([n->P_A(n, p=p) for p in ps], 1, 65535,
    label=[latexstring("p=$p") for p in ps],
    xaxis=:log2,
    xlabel="Number of Processors", ylabel="Speedup",
    legend=:topleft,
    minorgrid=true, minorgridalpha=0.1,
    framestyle=:box,
    tex_output_standalone=true
)

savefig(pa, "figures/theoretical-scaling-Amdahl.tex")

pg = plot([n->P_G(n, p=p) for p in ps], 1, 65535,
    label=[latexstring("p=$p") for p in ps],
    xaxis=:log2,
    xlabel="Number of Processors", ylabel="Speedup",
    legend=:topleft,
    minorgrid=true, minorgridalpha=0.1,
    framestyle=:box,
    tex_output_standalone=true
)

savefig(pg, "figures/theoretical-scaling-Gustafson.tex")

p1 = plot([P_A, P_G, identity], 1, 20,
    label=["strong scaling", "weak scaling", "perfect"],
    # xaxis=:log10, yaxis=:log10,
    xlabel="Number of Processors", ylabel="Speedup",
    legend=:topleft,
    minorgrid=true, minorgridalpha=0.1,
    framestyle=:box,
    tex_output_standalone=true
)

p2 = plot([P_A, P_G], 1, 65536,
    label=["strong scaling", "weak scaling"],
    xaxis=:log10, yaxis=:log10,
    xlabel="Number of Processors", ylabel="Speedup",
    legend=:topleft,
    minorgrid=true, minorgridalpha=0.1,
    framestyle=:box,
    tex_output_standalone=true
)

p = plot(p1, p2, tex_output_standalone=true)

savefig(p, "figures/theoretical-scaling-all.tex")

max_stars_repo_name: SebastianM-C/MasterThesis

max_issues_repo_name: SebastianM-C/MasterThesis

max_forks_repo_name: SebastianM-C/MasterThesis


text:
MOI.supports(::Optimizer, ::MOI.ObjectiveSense) = true
MOI.supports(::Optimizer, ::MOI.ObjectiveFunction{SVF}) = true

"""
set and get function overloads
"""
MOI.get(model::Optimizer, ::MOI.ObjectiveSense) = model.inner.sense

function MOI.set(model::Optimizer, ::MOI.ObjectiveSense, sense::MOI.OptimizationSense)
    model.inner.sense = sense
    return
end

function MOI.set(model::Optimizer, ::MOI.ObjectiveFunction, func::SVF)
    check_inbounds(model, func)
    model.inner.objective = SingleVariableObjective(CS.single_variable_objective, func.variable.value, [func.variable.value])
    return
end


text:
# Copyright 2018 IBM Corp.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

module TestRFlowGraphs
using Test

using Catlab.WiringDiagrams
using SemanticFlowGraphs

using ..IntegrationTest: db, read_semantic_graph

# Semantic flow graph
#####################

# K-means clustering on the Iris dataset using base R.
# FIXME: The encapsulated box should only have one input port.
semantic = read_semantic_graph(joinpath("r", "clustering_kmeans"))
d = WiringDiagram([], [])
file = add_box!(d, construct(pair(concepts(db, ["tabular-file", "filename"])...)))
read_file = add_box!(d, concept(db, "read-tabular-file"))
transform = add_box!(d, Box(concepts(db, ["table"]), concepts(db, ["table"])))
kmeans = add_box!(d, construct(pair(concepts(db,
  ["k-means", "clustering-model-n-clusters"])...)))
fit = add_box!(d, Hom("fit",
  otimes(concept(db, "k-means"), concept(db, "data")),
  concept(db, "k-means")))
centroids = add_box!(d, concept(db, "k-means-centroids"))
clusters = add_box!(d, concept(db, "clustering-model-clusters"))
add_wires!(d, [
  (file, 1) => (read_file, 1),
  (read_file, 1) => (transform, 1),
  (kmeans, 1) => (fit, 1),
  (transform, 1) => (fit, 2),
  (fit, 1) => (clusters, 1),
  (fit, 1) => (centroids, 1),
])
@test semantic == d

end

max_stars_repo_name: epatters/flowgraphs

max_issues_repo_name: epatters/flowgraphs

max_forks_repo_name: epatters/flowgraphs


text:
# DiagnosticConsumer
function BeginSourceFile(consumer::AbstractDiagnosticConsumer, lang::LangOptions,
                         pp::Preprocessor)
    @check_ptrs consumer lang pp
    clang_DiagnosticConsumer_BeginSourceFile(consumer, lang, pp)
    return nothing
end

function EndSourceFile(consumer::AbstractDiagnosticConsumer)
    @check_ptrs consumer
    clang_DiagnosticConsumer_EndSourceFile(consumer)
    return nothing
end

dispose(x::AbstractDiagnosticConsumer) = clang_DiagnosticConsumer_dispose(x)

# IgnoringDiagConsumer
IgnoringDiagConsumer() = IgnoringDiagConsumer(create_ignoring_diagnostic_consumer())

"""
    create_ignoring_diagnostic_consumer() -> CXDiagnosticConsumer
Return a pointer to a `clang::IgnoringDiagConsumer` object.
"""
function create_ignoring_diagnostic_consumer()
    status = Ref{CXInit_Error}(CXInit_NoError)
    consumer = clang_IgnoringDiagConsumer_create(status)
    @assert status[] == CXInit_NoError
    return consumer
end

# DiagnosticsEngine
DiagnosticsEngine() = DiagnosticsEngine(create_diagnostics_engine())

function DiagnosticsEngine(opts::DiagnosticOptions,
                           client::AbstractDiagnosticConsumer=TextDiagnosticPrinter(opts),
                           should_own_client=true)
    status = Ref{CXInit_Error}(CXInit_NoError)
    ids = create_diagnostic_ids()
    engine = clang_DiagnosticsEngine_create(ids, opts, client, should_own_client, status)
    @assert status[] == CXInit_NoError
    return DiagnosticsEngine(engine)
end

function DiagnosticsEngine(ids::DiagnosticIDs, opts::DiagnosticOptions,
                           client::AbstractDiagnosticConsumer=TextDiagnosticPrinter(opts),
                           should_own_client=true)
    status = Ref{CXInit_Error}(CXInit_NoError)
    engine = clang_DiagnosticsEngine_create(ids, opts, client, should_own_client, status)
    @assert status[] == CXInit_NoError
    return DiagnosticsEngine(engine)
end

"""
    create_diagnostics_engine() -> CXDiagnosticsEngine
Return a pointer to a `clang::DiagnosticsEngine` object.
"""
function create_diagnostics_engine()
    status = Ref{CXInit_Error}(CXInit_NoError)
    ids = create_diagnostic_ids()
    opts = DiagnosticOptions()
    client = TextDiagnosticPrinter(opts)
    should_own_client = true
    engine = clang_DiagnosticsEngine_create(ids, opts, client, should_own_client, status)
    @assert status[] == CXInit_NoError
    return engine
end

dispose(x::DiagnosticsEngine) = clang_DiagnosticsEngine_dispose(x)

function setShowColors(x::DiagnosticsEngine, should_show::Bool)
    @check_ptrs x
    return clang_DiagnosticsEngine_setShowColors(x, should_show)
end


text:
function StandaloneVertexEigenfunction(
    domain::Triangle{T},
    i::Integer;
    stride::Integer = 1,
    offset::Integer = 0,
    reversed::Bool = false,
    outside = false,
) where {T<:Union{arb,fmpq}}
    i ∈ boundaries(domain) ||
        throw(ArgumentError("attempt to get vertex $i from a $(typeof(domain))"))

    θ = domain.angles[i]

    if i == 1
        orientation = zero(θ)
    elseif i == 2
        orientation = ifelse(T == arb, domain.parent(π), 1) - θ
    elseif i == 3
        orientation = 2ifelse(T == arb, domain.parent(π), 1) - domain.angles[2] - θ
    end

    if outside
        orientation += θ
        θ = ifelse(T == arb, 2domain.parent(π), 2) - θ
    end

    return StandaloneVertexEigenfunction(
        vertex(domain, i),
        orientation,
        θ;
        stride,
        offset,
        reversed,
        domain.parent,
    )
end

function StandaloneVertexEigenfunction(
    domain::TransformedDomain,
    i::Integer;
    stride::Integer = 1,
    offset::Integer = 0,
    reversed::Bool = false,
    outside = false,
)
    u = StandaloneVertexEigenfunction(domain.original, i; stride, reversed, outside)
    # FIXME: This only works if u.θ::arb or if u.θ and
    # domain.orientation both are fmpq.
    u = StandaloneVertexEigenfunction(
        domain.map(u.vertex),
        u.orientation + domain.rotation,
        u.θ;
        stride,
        offset,
        reversed,
        domain.parent,
    )

    return u
end

function Base.show(io::IO, u::StandaloneVertexEigenfunction)
    println(io, "Standalone vertex eigenfunction")
    if !haskey(io, :compact) || !io[:compact]
        println(io, "vertex: $(u.vertex)")
        println(io, "orientation: $(u.orientation)")
        println(io, "θ: $(u.θ)")
        print(io, "number of set coefficients: $(length(u.coefficients))")
    end
end

function set_eigenfunction!(u::StandaloneVertexEigenfunction, coefficients::Vector)
    resize!(u.coefficients, length(coefficients))
    copy!(u.coefficients, u.parent.(coefficients))
end

"""
    nu(u::StandaloneVertexEigenfunction, k::Integer)

Return `k*θ` as an `arb`, the parameter used for the Bessel function.
"""
nu(u::StandaloneVertexEigenfunction{fmpq}, k::Integer = 1) = u.parent(k * inv(u.θ))
nu(u::StandaloneVertexEigenfunction{arb}, k::Integer = 1) = u.parent(k * u.parent(π) / u.θ)

"""
    coordinate_transformation(u::StandaloneVertexEigenfunction, xy::AbstractVector)

Takes a 2-element vector `xy` representing a point in the plane in
Cartesian coordinates and makes a (affine) change of coordinates so
that `u.vertex` is put at the origin and rotated `u.orientation`
clockwise.
"""
function coordinate_transformation(
    u::StandaloneVertexEigenfunction{T},
    xy::AbstractVector,
) where {T<:Union{arb,fmpq}}
    if T == fmpq
        if u.reversed
            s, c = sincospi(-u.orientation - u.θ, u.parent)
        else
            s, c = sincospi(-u.orientation, u.parent)
        end
    elseif T == arb
        s, c = sincos(-u.orientation)
    end
    M = SMatrix{2,2}(c, s, -s, c)
    res = M * (xy .- u.vertex)
    if u.reversed
        return SVector(res[1], -res[2])
    else
        return res
    end
end

function (u::StandaloneVertexEigenfunction)(
    xy::AbstractVector{T},
    λ::arb,
    ks::UnitRange{Int};
    boundary = nothing,
    notransform::Bool = false,
) where {T<:Union{arb,arb_series}}
    if !notransform
        xy = coordinate_transformation(u, xy)
    end

    r, θ = polar_from_cartesian(xy)

    # TODO: We have to choose a branch to work on. This does depend on
    # the domain but for now we only implement the one with θ on the
    # interval [0, 2π). Nemo doesn't implement mod2pi so we just do a
    # partial solution of adding 2π if it's below 0.
    if (T == arb && isnegative(θ)) || (T == arb_series && isnegative(θ[0]))
        θ += 2parent(θ)(π)
    end

    rsqrtλ = r * sqrt(λ)
    res = similar(ks, T)
    for i in eachindex(ks)
        k = 1 + (ks[i] - 1) * u.stride + u.offset
        ν = nu(u, k)

        res[i] = bessel_j(ν, rsqrtλ) * sin(ν * θ)
    end

    return res
end

# TODO: Figure out how to handle this
set_domain!(u::StandaloneVertexEigenfunction, ::AbstractDomain) = u

# TODO: Figure out how to handle this
recompute!(u::StandaloneVertexEigenfunction) = u

# Have a function that computes u from the domain?