text:
struct ODEQuadratureAdjointSensitivityFunction{rateType,uType,uType2,UF,G,JC,GC,DG,TJ,SType,CV} <: SensitivityFunction
  uf::UF
  g::G
  J::TJ
  dg_val::uType
  jac_config::JC
  g_grad_config::GC
  sensealg::QuadratureAdjoint
  f_cache::rateType
  discrete::Bool
  y::uType2
  sol::SType
  dg::DG
  colorvec::CV
end

@noinline function ODEQuadratureAdjointSensitivityFunction(g,u0,p,sensealg,discrete,sol,dg,tspan,colorvec)
  numparams = length(p)
  numindvar = length(u0)
  # if there is an analytical Jacobian provided, we are not going to do automatic `jac*vec`
  f = sol.prob.f
  isautojacvec = DiffEqBase.has_jac(f) ? false : get_jacvec(sensealg)
  J = isautojacvec ? nothing : similar(u0, numindvar, numindvar)

  if !discrete
    if dg != nothing || isautojacvec
      pg = nothing
      pg_config = nothing
    else
      pg = UGradientWrapper(g,tspan[2],p)
      pg_config = build_grad_config(sensealg,pg,u0,p)
    end
  else
    pg = nothing
    pg_config = nothing
  end

  if isautojacvec
    jac_config = nothing
    uf = nothing
  else
    uf = DiffEqDiffTools.UJacobianWrapper(f,tspan[2],p)
    jac_config = build_jac_config(sensealg,uf,u0)
  end

  y = copy(sol(tspan[1])) # TODO: Has to start at interpolation value!
  dg_val = similar(u0, numindvar) # number of funcs size
  f_cache = deepcopy(u0)

  return ODEQuadratureAdjointSensitivityFunction(uf,pg,J,dg_val,
                                       jac_config,pg_config,
                                       sensealg,f_cache,
                                       discrete,y,sol,dg,
                                       colorvec)
end

# u = λ'
function (S::ODEQuadratureAdjointSensitivityFunction)(du,u,p,t)
  @unpack y, sol, J, uf, sensealg, f_cache, jac_config, discrete, dg, dg_val, g, g_grad_config = S
  idx = length(y)
  f = sol.prob.f
  isautojacvec = DiffEqBase.has_jac(f) ? false : get_jacvec(sensealg)
  sol(y,t)
  λ     = u
  dλ    = du

  if !isautojacvec
    if DiffEqBase.has_jac(f)
      f.jac(J,y,p,t) # Calculate the Jacobian into J
    else
      uf.t = t
      jacobian!(J, uf, y, f_cache, sensealg, jac_config)
    end
    mul!(dλ',λ',J)
  else
    _dy, back = Tracker.forward(y) do u
      if DiffEqBase.isinplace(sol.prob)
        out_ = map(zero, u)
        f(out_, u, p, t)
        Tracker.collect(out_)
      else
        vec(f(u, p, t))
      end
    end
    dλ[:] = Tracker.data(back(λ)[1])
  end

  dλ .*= -one(eltype(λ))

  if !discrete
    if dg != nothing
      dg(dg_val,y,p,t)
    else
      g.t = t
      gradient!(dg_val, g, y, sensealg, g_grad_config)
    end
    dλ .+= dg_val
  end

  nothing
end

# g is either g(t,u,p) or discrete g(t,u,i)
@noinline function ODEAdjointProblem(sol,sensealg::QuadratureAdjoint,g,
                                     t=nothing,dg=nothing,
                                     callback=CallbackSet())
  f = sol.prob.f
  tspan = (sol.prob.tspan[2],sol.prob.tspan[1])
  discrete = t != nothing

  p = sol.prob.p
  p === DiffEqBase.NullParameters() && error("Your model does not have parameters, and thus it is impossible to calculate the derivative of the solution with respect to the parameters. Your model must have parameters to use parameter sensitivity calculations!")

  u0 = zero(sol.prob.u0)

  len = length(u0)
  λ = similar(u0, len)
  sense = ODEQuadratureAdjointSensitivityFunction(g,u0,
                                        p,sensealg,discrete,
                                        sol,dg,tspan,f.colorvec)

  init_cb = t !== nothing && sol.prob.tspan[2] == t[end]
  cb = generate_callbacks(sense, g, λ, t, callback, init_cb)
  z0 = vec(zero(λ))
  ODEProblem(sense,z0,tspan,p,callback=cb)
end

struct AdjointSensitivityIntegrand{pType,uType,rateType,S,AS,PF,PJC,PJT}
  sol::S
  adj_sol::AS
  p::pType
  y::uType
  λ::uType
  pf::PF
  f_cache::rateType
  pJ::PJT
  paramjac_config::PJC
  sensealg::QuadratureAdjoint
end

function AdjointSensitivityIntegrand(sol,adj_sol,sensealg)
  prob = sol.prob
  @unpack f, p, tspan, u0 = prob
  y = similar(sol.prob.u0)
  λ = similar(adj_sol.prob.u0)
  # we need to alias `y`
  pf = DiffEqDiffTools.ParamJacobianWrapper(f,tspan[1],y)
  f_cache = similar(y)
  isautojacvec = get_jacvec(sensealg)
  pJ = isautojacvec ? nothing : similar(u0,length(u0),length(p))

  if DiffEqBase.has_paramjac(f) || isautojacvec
    paramjac_config = nothing
  else
    paramjac_config = build_param_jac_config(sensealg,pf,y,p)
  end
  AdjointSensitivityIntegrand(sol,adj_sol,p,y,λ,pf,f_cache,pJ,paramjac_config,sensealg)
end

function (S::AdjointSensitivityIntegrand)(out,t)
  @unpack y, λ, pJ, pf, p, f_cache, paramjac_config, sensealg, sol, adj_sol = S
  f = sol.prob.f
  sol(y,t)
  adj_sol(λ,t)
  λ .*= -one(eltype(λ))
  isautojacvec = get_jacvec(sensealg)
  # y is aliased
  pf.t = t

  if !isautojacvec
    if DiffEqBase.has_paramjac(f)
      f.paramjac(pJ,y,p,t) # Calculate the parameter Jacobian into pJ
    else
      jacobian!(pJ, pf, p, f_cache, sensealg, paramjac_config)
    end
    mul!(out',λ',pJ)
  else
    _, back = Tracker.forward(y, p) do u, p
      if DiffEqBase.isinplace(sol.prob)
        out_ = map(zero, u)
        f(out_, u, p, t)
        Tracker.collect(out_)
      else
        vec(f(u, p, t))
      end
    end
    out[:] = vec(Tracker.data(back(λ)[2]))
  end
  out'
end

function (S::AdjointSensitivityIntegrand)(t)
  out = similar(S.p)
  S(out,t)
end

function _adjoint_sensitivities_u0(sol,sensealg::QuadratureAdjoint,args...;kwargs...)
  error("Can't get sensitivities of u0 with quadrature.")
end

function _adjoint_sensitivities(sol,sensealg::QuadratureAdjoint,alg,g,
                                t=nothing,dg=nothing;
                                abstol=1e-6,reltol=1e-3,
                                iabstol=abstol, ireltol=reltol,
                                kwargs...)
  adj_prob = ODEAdjointProblem(sol,sensealg,g,t,dg)
  adj_sol = solve(adj_prob,alg;abstol=abstol,reltol=reltol,
                               save_everystep=true,save_start=true,kwargs...)
  integrand = AdjointSensitivityIntegrand(sol,adj_sol,sensealg)

  if t === nothing
    res,err = quadgk(integrand,sol.prob.tspan[1],sol.prob.tspan[2],
                   atol=iabstol,rtol=ireltol)
  else
    res = zero(integrand.p)'
    for i in 1:length(t)-1
      res .+= quadgk(integrand,t[i],t[i+1],
                     atol=iabstol,rtol=ireltol)[1]
    end
    if t[1] != sol.prob.tspan[1]
      res .+= quadgk(integrand,sol.prob.tspan[1],t[1],
                     atol=iabstol,rtol=ireltol)[1]
    end
  end
  res
end


text:
mutable struct cairo_env
    surface::Ptr{Cvoid}
    cr::Ptr{Cvoid}
end

# init cairo layer
function initCairo(file_name::String,height::Float64,width::Float64)::cairo_env
    surface = ccall((:cairo_svg_surface_create,"libcairo"),Ptr{Cvoid},(Cstring,Float64,Float64),file_name,height,width)
    cr = ccall((:cairo_create,"libcairo"),Ptr{Cvoid},(Ptr{Cvoid},),surface)

    # set image to white
    ccall((:cairo_set_source_rgb,"libcairo"),Cvoid,(Ptr{Cvoid},Int8,Int8,Int8),cr,0,0,0)

    return cairo_env(surface,cr)
end

function destroyCairo(env::cairo_env)
    ccall((:cairo_surface_destroy,"libcairo"), Cvoid,(Ptr{Cvoid},),env.surface)
    ccall((:cairo_destroy,"libcairo"), Cvoid,(Ptr{Cvoid},),env.cr)
end

function drawTextCairo(env::cairo_env, pos_x::Float64, pos_y::Float64,font_size::Float64, font_name, text)
    ccall((:cairo_select_font_face,"libcairo"),Cvoid, (Ptr{Cvoid},Cstring, Cint,Cint), env.cr, font_name,0,0)
    ccall((:cairo_set_font_size,"libcairo"),Cvoid, (Ptr{Cvoid},Float64), env.cr, font_size)

    ccall((:cairo_move_to,"libcairo"),Cvoid,(Ptr{Cvoid},Float64,Float64),env.cr,pos_x,pos_y)

    ccall((:cairo_show_text,"libcairo"),Cvoid,(Ptr{Cvoid},Cstring),env.cr,text)
end

function drawRectCairo(env::cairo_env, pos_x::Float64, pos_y::Float64,width::Float64, height::Float64)
    ccall((:cairo_rectangle,"libcairo"),Cvoid,(Ptr{Cvoid},Float64,Float64,Float64,Float64),env.cr,pos_x,pos_y,width,height)

    # apply stroke
    ccall((:cairo_stroke,"libcairo"),Cvoid,(Ptr{Cvoid},),env.cr)
end

function setLineWidthCairo(env::cairo_env,line_width::Float64)
    #set line width first
    ccall((:cairo_set_line_width,"libcairo"),Cvoid,(Ptr{Cvoid},Float64),env.cr,line_width)
end

# pos_x and pos_y are the top of the diamond
function drawDiamondCairo(env::cairo_env,pos_x::Float64,pos_y::Float64, width::Float64,height::Float64)

    ccall((:cairo_move_to,"libcairo"),Cvoid,(Ptr{Cvoid},Float64,Float64),env.cr,pos_x,pos_y)

    ccall((:cairo_rel_line_to,"libcairo"),Cvoid,(Ptr{Cvoid},Float64,Float64),env.cr,width/2,height/2)

    ccall((:cairo_rel_line_to,"libcairo"),Cvoid,(Ptr{Cvoid},Float64,Float64),env.cr,-width/2,height/2)

    ccall((:cairo_rel_line_to,"libcairo"),Cvoid,(Ptr{Cvoid},Float64,Float64),env.cr,-width/2,-height/2)

    ccall((:cairo_close_path,"libcairo"),Cvoid,(Ptr{Cvoid},),env.cr)

    ccall((:cairo_stroke,"libcairo"),Cvoid,(Ptr{Cvoid},),env.cr)
end

function drawRelLineCairo(env::cairo_env,pos_x::Float64,pos_y::Float64,dist_x::Float64,dist_y::Float64)
    ccall((:cairo_move_to,"libcairo"),Cvoid,(Ptr{Cvoid},Float64,Float64),env.cr,pos_x,pos_y)

    ccall((:cairo_rel_line_to,"libcairo"),Cvoid,(Ptr{Cvoid},Float64,Float64),env.cr,dist_x,dist_y)

    ccall((:cairo_stroke,"libcairo"),Cvoid,(Ptr{Cvoid},),env.cr)
end

# the tip of the arrow will be at pos_x and pos_y
# width is the width of the base of the arrow
# length is the distance from tip to base
# angle is the orientation of the arrow, 0 points right, 90 points down, 180 left, and 270 up
function drawArrowCairo(env::cairo_env,pos_x::Float64,pos_y::Float64,width::Float64,height::Float64,angle::Float64)
    ccall((:cairo_move_to,"libcairo"),Cvoid,(Ptr{Cvoid},Float64,Float64),env.cr,pos_x,pos_y)

    # get the hypotenuse
    hypotenuse = ((width^2)*.25 + (height^2))^0.5

    # get the angle of the triangle
    tri_angle = atan(width*0.5/height)
    #other_angle = pi/2 - tri_angle

    # add the angle argument (convert to radians)
    arrow_angle = tri_angle + angle*pi/180
    other_angle = pi/2 - angle*pi/180

    # draw the hypotenuse on a specific angle
    ccall((:cairo_rel_line_to,"libcairo"),Cvoid,(Ptr{Cvoid},Float64,Float64),env.cr,-hypotenuse*cos(arrow_angle),hypotenuse*sin(arrow_angle))

    # draw the base line (the width line)
    ccall((:cairo_rel_line_to,"libcairo"),Cvoid,(Ptr{Cvoid},Float64,Float64),env.cr,-width*cos(other_angle),-width*sin(other_angle))

    ccall((:cairo_close_path,"libcairo"),Cvoid,(Ptr{Cvoid},),env.cr)

    # set fill rule and fill                                             CAIRO_FILL_RULE_EVEN_ODD
    ccall((:cairo_set_fill_rule,"libcairo"),Cvoid,(Ptr{Cvoid},Int32),env.cr,1)
    ccall((:cairo_fill,"libcairo"),Cvoid,(Ptr{Cvoid},),env.cr)

    # move back to where you were
    ccall((:cairo_move_to,"libcairo"),Cvoid,(Ptr{Cvoid},Float64,Float64),env.cr,pos_x,pos_y)

    #ccall((:cairo_stroke,"libcairo"),Cvoid,(Ptr{Cvoid},),env.cr)
end

# draw a line with an arrow on the end of it
function drawArrowLineCairo(env::cairo_env,
                            x::Float64,
                            y::Float64,
                            x_line_dist::Float64,
                            y_line_dist::Float64,
                            arrow_width::Float64,
                            arrow_height::Float64)

    drawRelLineCairo(env,x,y,x_line_dist,y_line_dist)

    cur_x = x+x_line_dist
    cur_y = y+y_line_dist

    if x_line_dist < 0.0
        arrow_angle = 180 -( 180*(atan(y_line_dist/x_line_dist)/pi))
    else
        arrow_angle = -( 180*(atan(y_line_dist/x_line_dist)/pi))
    end

    #put the arrow at the end of the line
    drawArrowCairo(env,cur_x,cur_y,arrow_width,arrow_height,arrow_angle)
end

ext: jl

max_stars_repo_name: aklapatch/CodeVisualized

max_issues_repo_name: aklapatch/CodeVisualized

max_forks_repo_name: aklapatch/CodeVisualized


text:
"""
This function defines how a dimension is split.
It is used everywhere dimensions are split, therefore it (and the cutoff) define the recursion pattern.

The current behavior is to try and get dimensions into primarily blocks of 8.
"""
function split_dim(N)
    N <= cutoff && return N, 0
    number_cutoff_blocks, remainder = divrem(N,cutoff)

    half_of_cutoff_blocks, odd = divrem(number_cutoff_blocks, 2)

    8(half_of_cutoff_blocks+odd), 8half_of_cutoff_blocks + remainder
end

function reasonably_square(mindim, maxdim)
    # mindim, maxdim = minmax(M, N)
    round(Int, √2 * mindim, RoundUp) >= maxdim
end

# function blocks_per_dim(N::Integer)
#     blocks = 1
#     block_size = N
#     cuts = 0
#     # bit twiddling is 50% faster
#     while cutoff < block_size
#         cuts += 1
#         blocks <<= 1 # *=2
#         block_size >>= 1 #  ÷=2
#     end
#     blocks + min(max(0, N - cutoff << cuts), blocks) #N - cutoff * 2^cuts
# end
# @generated blocks_per_dim(::Val{N}) where N = blocks_per_dim(N) # Is this needed?

# ### Implemented so that the leading blocks are larger. This is because we tend to do just slightly less with them.
# function block_position(N::Integer, i::Integer)
#     num_blocks = blocks_per_dim(N)
#     1 + num_blocks - cld(num_blocks*(N+1-i),N)
# end
# function block_position(N::Integer, i::Integer, j::Integer) #Not faster on 0.6?
#     num_blocks = blocks_per_dim(N)
#     1+num_blocks - cld(num_blocks*(N+1-i),N), 1+num_blocks - cld(num_blocks*(N+1-j),N)
# end

# function block_position(N, i)
#     position = 1
#     while N > cutoff
#         N, r = divrem(N, 2)
#         Npr = N + r
# #        @show i, N, Npr
#         if i > Npr
#             position += blocks_per_dim(Npr)#round(Int, log2(Npr/cutoff), RoundUp)
#             i -= Npr
#         else
#             N = Npr
#         end
#     end
#     position
# end

function dense_sub2ind(dims::NTuple{2}, i, j, offset = 0)
    M, N = dims
    minind, maxind = minmax(M, N)
    maxind <= cutoff && return (j-1)*M + i + offset
    if reasonably_square(minind, maxind) && minind > cutoff #divide both
        M1, M2 = split_dim(M)
        N1, N2 = split_dim(N)
        if i <= M1 && j <= N1 # upper left block
            return dense_sub2ind( (M1, N1), i, j, offset )
        elseif j <= N1 # lower left block
            return dense_sub2ind( (M2, N1), i - M1, j, offset + M1 * N1 )
        elseif i <= M1 # upper right block
            return dense_sub2ind( (M1, N2), i, j - N1, offset + M * N1 )
        else # lower right block
            return dense_sub2ind( (M2, N2), i - M1, j - N1, offset + M * N1 + M1 * N2 )
        end
    elseif M > N #divide M
        M1, M2 = split_dim(M)
        if i <= M1
            return dense_sub2ind( (M1, N), i, j, offset)
        else
            return dense_sub2ind( (M2, N), i - M1, j, offset + M1 * N)
        end
    else #divide N
        N1, N2 = split_dim(N)
        if j <= N1
            return dense_sub2ind( (M, N1), i, j, offset)
        else
            return dense_sub2ind( (M, N2), i, j - N1, offset + M * N1)
        end
    end
end

"""
Generates a recursive quote.
The function keeps track of
    M, number of rows
    N, number of columns
    linear_offset, linear offset into the data of original matrix
    row_offset, Cartesian offset of rows
    column_offset, Cartesian offset of columns
Through the recursion. When the recurion terminates, it walks the expression and replaces all instances of those symbols with their current value.
"""
function dense_recursion_quote(q, M, N, linear_offset = 0, row_offset = 0, column_offset = 0)
    minind, maxind = minmax(M, N)

    # The recursion terminates, and we substitute expressions for their value.
    if maxind <= cutoff
        return MacroTools.postwalk(x -> begin
            if isa(x, Symbol)
                if x == :M
                    return M
                elseif x == :N
                    return N
                elseif x == :linear_offset
                    return linear_offset
                elseif x == :row_offset
                    return row_offset
                elseif x == :column_offset
                    return column_offset
                else
                    return x
                end
            else
                return x
            end
        end, q)
    end

    if reasonably_square(minind, maxind) && minind > cutoff # we divide the matrix into four blocks.
        M1, M2 = split_dim(M)
        N1, N2 = split_dim(N)
        return :(if i <= $(M1 + row_offset) && j <= $(N1 + column_offset) # enter block 1
                $(dense_recursion_quote(q, M1, N1, linear_offset, row_offset, column_offset))
            elseif j <= $(N1 + column_offset) # enter block 2
                $(dense_recursion_quote(q, M2, N1, linear_offset + M1 * N1, row_offset + M1, column_offset))
            elseif i <= $(M1 + row_offset) # enter block 3
                $(dense_recursion_quote(q, M1, N2, linear_offset + M * N1, row_offset, column_offset + N1))
            else # enter block 4
                $(dense_recursion_quote(q, M2, N2, linear_offset + M * N1 + M1 * N2, row_offset + M1, column_offset + N1))
            end)
    elseif M > N # we are splitting the matrix into two blocks stacked on top of one another.
        M1, M2 = split_dim(M)
        return :(if i <= $(M1 + row_offset)
                $(dense_recursion_quote(q, M1, N, linear_offset, row_offset, column_offset))
            else
                $(dense_recursion_quote(q, M2, N, linear_offset + M1 * N, row_offset + M1, column_offset))
            end)
    else # we are splitting the matrix into two blocks, side by side.
        N1, N2 = split_dim(N)
        return :(if j <= $(N1 + column_offset)
                $(dense_recursion_quote(q, M, N1, linear_offset, row_offset, column_offset))
            else
                $(dense_recursion_quote(q, M, N2, linear_offset + M * N1, row_offset, column_offset + N1))
            end)
    end
end

function dense_sub2ind_quote(M, N)
    dense_recursion_quote(
        :( (j-(column_offset+1) ) * M + i + (linear_offset-row_offset) ),
        M, N)
end

"""
Cartesian indexing isn't recomended, but it is convenient for printing.
The approach for sub2ind here incurs branches, which would disable SIMD.
Therefore, Cartesian indexing is strongly discouraged for hot loops.
Still, it is reasonably fast -- only a handful of ns.
"""
@generated dense_sub2ind(::Val{M}, ::Val{N}, i, j) where {N,M} = dense_sub2ind_quote(M, N)
# @generated function dense_sub2ind(::Val{M}, ::Val{N}, i, j) where {N,M}
#     dense_sub2ind_quote(M, N)
# end

function triangle_sub2ind(N, i, j, offset = 0)
    N <= cutoff && return small_triangle(j) + i + offset
    N1, N2 = split_dim(N)
    if max(i,j) <= N1 # enter block 1 # triangle / symmetric block
        return triangle_sub2ind(N1, i, j, offset)
    elseif i <= N1 # enter block 3 # square block
        return dense_sub2ind((N1, N2), i, j - N1, offset + big_triangle(N1))
    elseif j <= N1 # enter block 2 #symmetric or 0 block # Not supported; check i, j before calling.
        throw("Indexing into this half not supported.")
    else # enter block 4 # triangle / symmetric block
        return triangle_sub2ind(N2, i - N1, j - N1, offset + big_triangle(N1) + N1*N2)
    end
end

"""
The triangular (or Symmetric) matrix is assumed to be square, with only upper triangular elements represented.
Lower triangular matrices are simply a view of an upper triangular matrix that pretends the storage is row-major.

This function takes in two quotes, to allow for different behavior in a triangular / diagonal block, and the offdiagonal blocks.

When the recursion terminates, the following symbols are replaced with their current value:
    M -- only used for the off diagonal blocks. See dense_recursion_quote for details.
    N -- the current size of the matrix block.
    linear_offset -- linear offset within the original data matrix.
    row_offset -- how many rows the block is offset
    column_offset -- how many columns the block has been offset
"""
function triangle_recursion_quote(triangle_q, rectangle_q, N, linear_offset = 0, row_offset = 0, column_offset = 0)
    if N <= cutoff
        return MacroTools.postwalk(x -> begin
            if isa(x, Symbol)
                if x == :N
                    return N
                elseif x == :linear_offset
                    return linear_offset
                elseif x == :row_offset
                    return row_offset
                elseif x == :column_offset
                    return column_offset
                else
                    return x
                end
            else
                return x
            end
        end, triangle_q)
    end

    N1, N2 = split_dim(N)
    :(if max( i - $row_offset, j - $column_offset) <= $N1 # enter block 1 # triangle / symmetric block
            $(triangle_recursion_quote(triangle_q, rectangle_q, N1, linear_offset, row_offset, column_offset))
        elseif i <= $(N1 + row_offset) # enter block 3 # square block
            $(dense_recursion_quote(rectangle_q, N1, N2, linear_offset + big_triangle(N1), row_offset, column_offset + N1))
        else # enter block 4 # triangle / symmetric block
            $(triangle_recursion_quote(triangle_q, rectangle_q, N2,
                linear_offset + big_triangle(N1) + N1*N2, row_offset + N1, column_offset + N1))
        end)

end

function triangle_sub2ind_quote(N)
    triangle_recursion_quote(
        :( small_triangle( j-column_offset ) + i + (linear_offset-row_offset) ),
        :( (j-(column_offset+1) ) * M + i + (linear_offset-row_offset) ),
        N
    )
end

@generated function triangle_sub2ind(::Val{N}, i, j) where N
    quote
        Base.@_inline_meta
        @boundscheck i > j && throw("Indexing into this half not supported.")
        $(triangle_sub2ind_quote(N))
    end
end

function block_dim(M,N)
    mindim, maxdim = minmax(M,N)
    # The "reasonably square" may need improvement
    # For example, a 91 x 64 breaks down into
    # 46 x 64 and 45 x 64; the 45 x 64 also fails
    # the criterion, and breaks into two 45 x 32 blocks
    # How bad is this wonky pattern?
    # Should I round in the comparison?
    if reasonably_square(mindim, maxdim) && mindim > cutoff
        M1, M2 = split_dim(M)
        N1, N2 = split_dim(N)
        out = [ (M1,N1) (M1, N2); (M2,N1) (M2,N2) ]
    elseif maxdim <= cutoff
        out = Matrix{Tuple{Int,Int}}(undef,1,1)
        out[1] = (M,N)
    elseif M > N
        out = Matrix{Tuple{Int,Int}}(undef,2,1)
        M1, M2 = split_dim(M)
        out[1] = (M1, N)
        out[2] = (M2, N)
    else
        out = Matrix{Tuple{Int,Int}}(undef,1,2)
        N1, N2 = split_dim(N)
        out[1] = (M, N1)
        out[2] = (M, N2)
    end
    out
end

# function split_dim!(n_dim::Vector{Int}, i)
#     N1, N2 = split_dim(n_dim[i])
#     n_dim[i] = N2
#     insert!(n_dim, i, N1)
# end
# function split_row(dims::Matrix{Tuple{Int}}, i)
#     N1, N2 = split_dim(n_dim[i])
#     n_dim[i] = N2
#     insert!(n_dim, i, N1)
# end
function split_col(dims::Matrix{Tuple{Int,Int}})
    num_rowblocks, num_colblocks = size(dims)
    @assert num_colblocks == 1
    vcat([block_dim(dims[i]...) for i = 1:length(dims)]...)
end
function split_row(dims::Matrix{Tuple{Int,Int}})
    num_rowblocks, num_colblocks = size(dims)
    @assert num_rowblocks == 1
    hcat([block_dim(dims[i]...) for i = 1:length(dims)]...)
end

# function block_structure(N, offset = 0)
#     bounds = Int[]
#     block_structure!(bounds, N)
#     # sort!(bounds)
# end
# function block_structure!(bounds, N, offset = 0)
#     if N > cutoff
#         N, r = divrem(N, 2)
#         Npr = N + r
#         push!(bounds, Npr + offset)
#         block_structure!(bounds, Npr, offset)
#         block_structure!(bounds, N, offset + Npr)
#         return bounds
#     else
#         return bounds
#     end
# end

# function identify_block_placement_quote(N)
#     bounds = block_structure(N)
#     num_bounds = length(bounds)
#     num_blocks = num_bounds + 1
#     q = quote
#         b_0 = 0
#         $(Symbol(:b_, num_blocks)) = $N
#     end
#     for k ∈ 1:num_bounds
#         push!(q.args, :( @inbounds $(Symbol(:b_, k)) = $(bounds[k])) )
#     end
#     push!(q.args,
#         :( Base.Cartesian.@nif $num_blocks d -> (i <= b_d)  d -> ( row_block_id = d; row_lbound = b_{d-1}; row_ubound = b_d ) d -> ( row_block_id = d; row_lbound = b_{d-1}; row_ubound = b_d ) )
#     )
#     push!(q.args,
#         :( Base.Cartesian.@nif $num_blocks d -> (j <= b_d)  d -> ( col_block_id = d; col_lbound = b_{d-1}; col_ubound = b_d ) d -> ( col_block_id = d; col_lbound = b_{d-1}; col_ubound = b_d ) )
#     )
#     push!(q.args, :(block_i = i - row_lbound))
#     push!(q.args, :(block_j = j - col_lbound))
#     q
# end
# function identify_block_placement_quote(M, N)
#     rowbounds = block_structure(M)
#     colbounds = block_structure(N)
#     num_rowbounds = length(rowbounds)
#     num_colbounds = length(colbounds)
#     num_rowblocks = num_rowbounds + 1
#     num_colblocks = num_colbounds + 1
#     q = quote
#         r_0 = 0
#         c_0 = 0
#         $(Symbol(:r_, num_rowblocks)) = $M
#         $(Symbol(:c_, num_colblocks)) = $N
#     end
#     for k ∈ 1:num_rowbounds
#         push!(q.args, :( @inbounds $(Symbol(:r_, k)) = $(rowbounds[k])) )
#     end
#     for k ∈ 1:num_colbounds
#         push!(q.args, :( @inbounds $(Symbol(:c_, k)) = $(colbounds[k])) )
#     end
#     push!(q.args,
#         :( Base.Cartesian.@nif $num_rowblocks d -> (i <= r_d)  d -> ( row_lbound = r_{d-1}; row_ubound = r_d ) d -> (row_lbound = r_{d-1}; row_ubound = r_d) )
#     )
#     push!(q.args,
#         :( Base.Cartesian.@nif $num_colblocks d -> (j <= c_d)  d -> ( col_lbound = c_{d-1}; col_ubound = c_d ) d -> (col_lbound = c_{d-1}; col_ubound = c_d) )
#     )
#     push!(q.args, :(block_i = i - row_lbound))
#     push!(q.args, :(block_j = j - col_lbound))
#     q
# end

# @generated function triangle_sub2ind(::Val{N}, i, j) where N
#     N <= cutoff && return :(small_triangle(j)+i)
#     q = identify_block_placement_quote(N)
#     push!(q.args, :(ind = big_triangle(col_lbound) + (col_ubound - col_lbound) * row_lbound))
#     push!(q.args, :(if row_block_id == col_block_id # then we're in a symmetric portion
#         ind += small_triangle(block_j) + block_i
#     else #then we are in square block
#         ind += (row_ubound-row_lbound) * (block_j-1) + block_i
#     end))
#     push!(q.args, :ind)
#     q
# end

# @generated dense_sub2ind(::Val{M}, ::Val{N}, i, j) where {M,N}
#     max(M,N) <= cutoff && return :(j*$(M-1)+i)
#     q = identify_block_placement_quote(M,N)
#     push!(q.args :( col_lbound*$N + row_lbound*(col_ubound-col_lbound) + (row_ubound - row_lbound)*(block_j-1) + block_i ))
#     q
# end

# @generated function triangle_sub2ind(::Val{N}, i, j) where N
#     N <= cutoff && return :(small_triangle(j)+i)
#     num_blocks = blocks_per_dim(N)
#     quote
#         row_block_id = $(1+num_blocks) - cld($num_blocks*($(N+1)-i),$N)
#         col_block_id = $(1+num_blocks) - cld($num_blocks*($(N+1)-j),$N)

#         row_block_lbound = cld((row_block_id-1)*$N,$num_blocks)
#         # row_block_ubound = cld( row_block_id   *$N,$num_blocks)
#         col_block_lbound = cld((col_block_id-1)*$N,$num_blocks)
#         col_block_ubound = cld( col_block_id   *$N,$num_blocks)

#         ind = big_triangle(col_block_lbound) + (col_block_ubound - col_block_lbound) * row_block_lbound

#         block_i = i - row_block_lbound
#         block_j = j - col_block_lbound

#         if row_block_id == col_block_id # then we're in a symmetric portion
#             ind += small_triangle(block_j) + block_i
#         else #then we are in square block
#             nrows_in_block = cld( row_block_id   *$N,$num_blocks)  - row_block_lbound
#             ind += nrows_in_block * (block_j-1) + block_i
#         end
#         ind
#     end
# end


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

"""
    BAT.MultimodalStudentT <: Distribution{Multivariate,Continuous}

*Experimental feature, not part of stable public API.*

The Multimodal Student-t Distribution (It's defined similarly to multimodal Cauchy
defined in [Caldwell et al.](https://arxiv.org/abs/1808.08051)).

Assumes two bimodal peaks, each in its own dimension.

Constructors:

* ```BAT.MultimodalStudentT(; μ::Real=1, σ::Float64=0.2, ν::Integer=1, n::Integer=4)```

Arguments:

* `μ::Real`: The location parameter used for the two bimodal peaks.
* `σ::Float64`: The scale parameter shared among all components.
* `ν::Int`: The degrees of freedom.
* `n::Int`: The number of dimensions.

Fields:

$(TYPEDFIELDS)

!!! note

    All fields of `MultimodalStudentT` are considered internal and subject to
    change without deprecation.
"""
struct MultimodalStudentT{M<:MixtureModel, P<:Product} <: Distribution{Multivariate,Continuous}
    bimodals::M
    σ::Float64
    n::Int
    dist::P
end

function MultimodalStudentT(;μ::Real=1., σ::Float64=0.2, ν::Int=1, n::Integer=4)
    @argcheck n > 1 "Minimum number of dimensions for MultimodalCauchy is 2"
    mixture_model = MixtureModel([LocationScale(-μ, σ, TDist(ν)), LocationScale(μ, σ, TDist(ν))])
    dist = _construct_dist(mixture_model, σ, ν, n)
    MultimodalStudentT(mixture_model, σ, n, dist)
end

Base.size(d::MultimodalStudentT) = size(d.dist)
Base.length(d::MultimodalStudentT) = length(d.dist)
Base.eltype(d::MultimodalStudentT) = eltype(d.dist)

Statistics.mean(d::MultimodalStudentT) = Distributions.mean(d.dist)
Statistics.var(d::MultimodalStudentT) = Distributions.var(d.dist)
Statistics.cov(d::MultimodalStudentT) = Distributions.cov(d.dist)

function StatsBase.params(d::MultimodalStudentT)#Check this
    (
        vcat(d.bimodals.components[1].μ, d.bimodals.components[2].μ, zeros(d.n-2)),
        [d.σ for i in 1:d.n],
        d.n
    )
end

function Distributions._logpdf(d::MultimodalStudentT, x::AbstractArray)
    convert(eltype(x), Distributions._logpdf(d.dist, x))::eltype(x)
end

function Distributions._rand!(rng::AbstractRNG, d::MultimodalStudentT, x::AbstractVector)
    Distributions._rand!(rng, d.dist, x)
end

function _construct_dist(mixture_model::MixtureModel, σ::Real, ν::Int, n::Integer)
    vector_of_dists = vcat(mixture_model, mixture_model, [LocationScale(0, σ, TDist(ν)) for i in 3:n])
    dist = product_distribution(vector_of_dists)
    return dist
end


text:
#######################################################################################################################################################################################################
#
# Changes to the function
# General
#     2021-Dec-24: migrate the function from PkgUtility to NetcdfIO
#     2022-Jan-28: remove the complicated funtion to create var and dim at the same time
#     2022-Jan-28: add global attributes to the generated file
#
#######################################################################################################################################################################################################
"""
NCDatasets.jl does not have a convenient function (1 line command) to save dataset as a file. Thus, we provide a few methods as supplements:

$(METHODLIST)

"""
function save_nc! end

"""
To save the code and effort to redefine the common attributes like latitude, longitude, and cycle index, we provide a shortcut method that handles these within the function:

    save_nc!(file::String, var_name::String, var_data::Array{T,N}, var_attribute::Dict{String,String}; compress::Int = 4, growable::Bool = false) where {T<:Union{AbstractFloat,Int,String},N}

Save the 1D, 2D, or 3D data as netcdf file, given
- `file` Path to save the dataset
- `var_name` Variable name for the data in the NC file
- `var_data` Data to save
- `var_attribute` Variable attributes for the data, such as unit and long name
- `compress` Compression level fro NetCDF, default is 4
- `growable` If true, make index growable, default is false

---
# Examples
```julia
# generate data to write into NC file
data1 = rand(12) .+ 273.15;
data2 = rand(36,18) .+ 273.15;
data3 = rand(36,18,12) .+ 273.15;

# save data as NC files (2D and 3D)
save_nc!("data1.nc", "data1", data1, Dict("description" => "Random temperature", "unit" => "K"));
save_nc!("data2.nc", "data2", data2, Dict("description" => "Random temperature", "unit" => "K"));
save_nc!("data3.nc", "data3", data3, Dict("description" => "Random temperature", "unit" => "K"));
```
"""
save_nc!(file::String, var_name::String, var_data::Array{T,N}, var_attribute::Dict{String,String}; compress::Int = 4, growable::Bool = false) where {T<:Union{AbstractFloat,Int,String},N} = (
    @assert 1 <= N <= 3 "Variable must be a 1D, 2D, or 3D dataset!";
    @assert 0 <= compress <= 9 "Compression rate must be within 0 to 9";

    # create the file
    _dset = Dataset(file, "c");

    # global title attribute
    for (_title,_notes) in ATTR_ABOUT
        _dset.attrib[_title] = _notes;
    end;

    # the case if the dimension is 1D
    if N==1
        _n_ind = (growable ? Inf : length(var_data));
        _inds  = collect(eachindex(var_data));
        add_nc_dim!(_dset, "ind", _n_ind);
        append_nc!(_dset, "ind", _inds, ATTR_CYC, ["ind"]; compress=compress);
        append_nc!(_dset, var_name, var_data, var_attribute, ["ind"]; compress=compress);

        close(_dset);

        return nothing
    end;

    # if the dimension is 2D or 3D
    _n_lon   = size(var_data, 1);
    _n_lat   = size(var_data, 2);
    _res_lon = 360 / _n_lon;
    _res_lat = 180 / _n_lat;
    _lons    = collect(_res_lon/2:_res_lon:360) .- 180;
    _lats    = collect(_res_lat/2:_res_lat:180) .- 90;
    add_nc_dim!(_dset, "lon", _n_lon);
    add_nc_dim!(_dset, "lat", _n_lat);
    append_nc!(_dset, "lon", _lons, ATTR_LON, ["lon"]; compress=compress);
    append_nc!(_dset, "lat", _lats, ATTR_LAT, ["lat"]; compress=compress);

    if N==2
        append_nc!(_dset, var_name, var_data, var_attribute, ["lon", "lat"]; compress=compress);
    elseif N==3
        _n_ind = (growable ? Inf : size(var_data,3));
        _inds  = collect(1:_n_ind);
        add_nc_dim!(_dset, "ind", _n_ind);
        append_nc!(_dset, "ind", _inds, ATTR_CYC, ["ind"]; compress=compress);
        append_nc!(_dset, var_name, var_data, var_attribute, ["lon", "lat", "ind"]; compress=compress);
    end;

    close(_dset);

    return nothing
);

"""
This method saves DataFrame as a NetCDF file to save more space (compared to a CSV file).

    save_nc!(file::String, df::DataFrame, var_names::Vector{String}, var_attributes::Vector{Dict{String,String}}; compress::Int = 4, growable::Bool = false)

Save DataFrame to NetCDF, given
- `file` Path to save the data
- `df` DataFrame to save
- `var_names` The label of data in DataFrame to save
- `var_attributes` Variable attributes for the data to save
- `compress` Compression level fro NetCDF, default is 4
- `growable` If true, make index growable, default is false

---
# Examples
```julia
df = DataFrame();
df[!,"A"] = rand(5);
df[!,"B"] = rand(5);
df[!,"C"] = rand(5);
save_nc!("dataf.nc", df, ["A","B"], [Dict("A" => "Attribute A"), Dict("B" => "Attribute B")]);
```
"""
save_nc!(file::String, df::DataFrame, var_names::Vector{String}, var_attributes::Vector{Dict{String,String}}; compress::Int = 4, growable::Bool = false) = (
    @assert 0 <= compress <= 9 "Compression rate must be within 0 to 9";
    @assert length(var_names) == length(var_attributes) "Variable name and attributes lengths must match!";

    # create the file
    _dset = Dataset(file, "c");

    # global title attribute
    for (_title,_notes) in ATTR_ABOUT
        _dset.attrib[_title] = _notes;
    end;

    # define dimension related variables
    _n_ind = (growable ? Inf : size(df)[1]);
    _inds  = collect(1:_n_ind);

    # save the variables
    add_nc_dim!(_dset, "ind", _n_ind);
    append_nc!(_dset, "ind", _inds, ATTR_CYC, ["ind"]; compress=compress);
    for _i in eachindex(var_names)
        append_nc!(_dset, var_names[_i], df[:, var_names[_i]], var_attributes[_i], ["ind"]; compress = compress);
    end;

    close(_dset);

    return nothing
);

"""
This method is a simplified version of the method above, namely when users do not want to define the attributes.

    save_nc!(file::String, df::DataFrame; compress::Int = 4, growable::Bool = false)

Save DataFrame to NetCDF, given
- `file` Path to save the data
- `df` DataFrame to save
- `notes` Global attributes (notes)
- `compress` Compression level fro NetCDF, default is 4
- `growable` If true, make index growable, default is false

---
# Examples
```julia
df = DataFrame();
df[!,"A"] = rand(5);
df[!,"B"] = rand(5);
df[!,"C"] = rand(5);
save_nc!("test.nc", df);
```
"""
save_nc!(file::String, df::DataFrame; compress::Int = 4, growable::Bool = false) = (
    _var_names = names(df);
    _var_attrs = [Dict{String,String}(_vn => _vn) for _vn in _var_names];

    save_nc!(file, df, _var_names, _var_attrs; compress=compress, growable = growable);

    return nothing
);


text:
"    SearchWindow <: GtkFrame
Search/replace panel that pops-up at the bottom of the editor.
It uses a global `GtkSourceSearchSettings` (search_settings) alongside
each `EditorTab` `GtkSourceSearchContext` (search_context).
Each tab also store the position of the current match using `GtkTextMark`'s.
"
mutable struct SearchWindow <: GtkFrame

    handle::Ptr{Gtk.GObject}
    search_entry::GtkEntry
    replace_entry::GtkEntry
    search_button::GtkButton
    replace_button::GtkButton
    replace_all_button::GtkButton
    search_settings::GtkSourceSearchSettings
    editor#not defined at this point

    function SearchWindow(editor)

        search_window = GtkFrame("") |>
            (GtkBox(:v) |>
                ((GtkBox(:h)) |>
                    (search_entry  = GtkEntry()) |>
                    (search_button = GtkButton("Search")) |>
                    (case_button = GtkToggleButton("Aa")) |>
                    (word_button = GtkToggleButton("Word"))
                ) |>
                ((GtkBox(:h)) |>
                    (replace_entry  = GtkEntry()) |>
                    (replace_button = GtkButton("Replace")) |>
                    (replace_all_button = GtkButton("Replace All"))
                )
            )

        set_gtk_property!(search_window, :height_request, 70)
        set_gtk_property!(search_entry, :hexpand, true)
        set_gtk_property!(replace_entry, :hexpand, true)

        search_settings = GtkSourceSearchSettings()
        set_gtk_property!(search_settings, :wrap_around, true)

        w = new(search_window.handle,
        search_entry, replace_entry,
        search_button, replace_button, replace_all_button,
        search_settings, editor)

        Gtk.gobject_move_ref(w, search_window)

        signal_connect(case_button_toggled_cb, case_button, "toggled", Nothing, (), false, w)
        signal_connect(word_button_toggled_cb, word_button, "toggled", Nothing, (), false, w)

        w
    end
end

import Base.open
@guarded (nothing) function open(w::SearchWindow, tab)

    visible(w, true)
    (found, it_start, it_end) = selection_bounds(tab.buffer)
    if found
        w.search_entry.text[String] = (it_start:it_end).text[String]
    end
    grab_focus(w.search_entry)
    showall(w)
    nothing
end

#first define the callbacks that are used in the constructor
function case_button_toggled_cb(widgetptr::Ptr, user_data)
#    tb = convert(GtkToggleButton, widgetptr)
    search_window = user_data
    iscase = get_gtk_property(search_window.search_settings, :case_sensitive, Bool)
    set_gtk_property!(search_window.search_settings, :case_sensitive, !iscase)
    return nothing
end
function word_button_toggled_cb(widgetptr::Ptr, user_data)
#    tb = convert(GtkToggleButton, widgetptr)
    search_window = user_data
    isword = get_gtk_property(search_window.search_settings, :at_word_boundaries, Bool)
    set_gtk_property!(search_window.search_settings, :at_word_boundaries, !isword)
    return nothing
end

import GtkSourceWidget.get_search_text
get_search_text(s::GtkSourceSearchSettings) = get_gtk_property(s, :search_text, AbstractString)

function search_entry_key_press_cb(widgetptr::Ptr, eventptr::Ptr, user_data)

    widget = convert(GtkEntry, widgetptr)
    event = convert(Gtk.GdkEvent, eventptr)
    search_window = user_data

    if event.keyval == GdkKeySyms.Escape
        set_search_text(search_window.search_settings, "")
        visible(search_window, false)
    end

    if event.keyval == GdkKeySyms.Return
        t = current_tab(search_window.editor)
        search_forward(t)
    end

    return convert(Cint, false)
end

function search_forward(t::EditorTab)

    if isnothing(t.search_mark_end)
        t.search_mark_end = create_mark(t.buffer, Gtk.GtkTextIter(t.buffer, 1))#search from the start
    end

    it = GtkTextIter(t.buffer, t.search_mark_end)
    found, its, ite = search_context_forward(t.search_context, it)

    if found
        scroll_to(t.view, its, 0, true, 0.0, 0.15)
        t.search_mark_start  = create_mark(t.buffer, its)#save the position for next search
        t.search_mark_end  = create_mark(t.buffer, ite)
    end
    return found
end

function search_entry_key_release_cb(widgetptr::Ptr, eventptr::Ptr, user_data)
    widget = convert(GtkEntry, widgetptr)
    search_window = user_data

    s = get_gtk_property(widget, :text, AbstractString)
    set_search_text(search_window.search_settings, s)

    return convert(Cint, false)
end

@guarded (INTERRUPT)  function replace_entry_key_press_cb(widgetptr::Ptr, eventptr::Ptr, user_data)
    entry = convert(GtkEntry, widgetptr)
    event = convert(Gtk.GdkEvent, eventptr)
    search_window = user_data

    if event.keyval == GdkKeySyms.Escape
        set_search_text(search_window.search_settings, "")
        visible(search_window, false)
    end

    if event.keyval == GdkKeySyms.Return
        t = current_tab(search_window.editor)
        replace_forward(t, entry, search_window)
    end

    return PROPAGATE
end

function replace_forward(t::EditorTab, entry::GtkEntry, search_window::SearchWindow)

    search_text = get_search_text(search_window.search_settings)

    do_search = false
    if t.search_mark_end != nothing && t.search_mark_start != nothing
        its = GtkTextIter(t.buffer, t.search_mark_start)
        ite = GtkTextIter(t.buffer, t.search_mark_end)

        #if it doesn't match, we search forward
        if (its:ite).text[String] != search_text
            do_search = true
        end
    else
        do_search = true
    end

    if do_search
        !search_forward(t) && return INTERRUPT
        its = GtkTextIter(t.buffer, t.search_mark_start)
        ite = GtkTextIter(t.buffer, t.search_mark_end)
    end

    s = get_gtk_property(entry, :text, AbstractString)
    search_context_replace(t.search_context, its, ite, s)
end

function replace_all(t::EditorTab, entry::GtkEntry, search_window::SearchWindow)

    search_text = get_search_text(search_window.search_settings)
    search_text == "" && return

    s = get_gtk_property(entry, :text, AbstractString)
    GtkSourceWidget.search_context_replace_all(t.search_context, s)
end

function search_button_clicked_cb(widgetptr::Ptr, user_data)
    search_window = user_data
    search_forward(current_tab(search_window.editor))
    return nothing
end

function replace_button_clicked_cb(widgetptr::Ptr, user_data)
    search_window = user_data
    t = current_tab(search_window.editor)
    replace_forward(t, search_window.replace_entry, search_window)
    return nothing
end

function init!(search_window::SearchWindow)
    signal_connect(search_entry_key_press_cb, search_window.search_entry, "key-press-event", Cint, (Ptr{Gtk.GdkEvent}, ), false, search_window)
    signal_connect(search_entry_key_release_cb, search_window.search_entry, "key-release-event", Cint, (Ptr{Gtk.GdkEvent}, ), false, search_window)
    signal_connect(replace_entry_key_press_cb, search_window.replace_entry, "key-press-event", Cint, (Ptr{Gtk.GdkEvent}, ), false, search_window)
    signal_connect(search_button_clicked_cb, search_window.search_button, "clicked", Nothing, (), false, search_window)
    signal_connect(replace_button_clicked_cb, search_window.replace_button, "clicked", Nothing, (), false, search_window)
    signal_connect(search_window_quit_cb, search_window, "delete-event", Cint, (Ptr{Gtk.GdkEvent}, ), false)
    signal_connect(replace_all_button_clicked_cb, search_window.replace_all_button, "clicked", Nothing, (), false, search_window)
end

function replace_all_button_clicked_cb(widgetptr::Ptr, user_data)
    search_window = user_data
    replace_all(current_tab(search_window.editor), search_window.replace_entry, search_window)
    return nothing
end

function search_window_quit_cb(widgetptr::Ptr, eventptr::Ptr, user_data)

    return convert(Cint, true)
end

##


text:
"""
Hadamard sense fractional derivative algorithms.
"""
abstract type HadamardDiff <: FracDiffAlg end

"""
@article{article,
author = {Yin, Chuntao and Li, Changpin and Bi, Qinsheng},
year = {2018},
month = {01},
pages = {},
title = {Approximating Hadamard Derivative and Fractional Differential Equation via the Finite Part Integral},
journal = {SSRN Electronic Journal},
doi = {10.2139/ssrn.3281675}
}
"""

"""
# Hadamard sense left rectangular approximating algorithm

    fracdiff(f, α, a, x, h, HadamardLRect())

Using finite part integral of **left rectangular formula** to compute the fractional hadamard derivative.

### Example

```julia-repl
julia> fracdiff(x->x, 0.5, 0, 1, 0.001, HadamardLRect())
0.9738507879228357
```
"""
struct HadamardLRect <: HadamardDiff end

"""
# Hadamard sense right rectangular approximating algorithm

    fracdiff(f, α, a, x, h, HadamardTrap())

Using finite part integral of **right rectangular formula** to compute the fractional hadamard derivative.

### Example

```julia-repl
julia> fracdiff(x->x, 0.5, 0, 1, 0.001, HadamardRRect())
0.9738507879228337
```
"""
struct HadamardRRect <: HadamardDiff end

"""
# Hadamard sense trapezoidal approximating algorithm

    fracdiff(f, α, a, x, h, HadamardTrap())

Using finite part integral of **trapezoidal formula** to compute the fractional hadamard derivative.

### Example

```julia-repl
julia> fracdiff(x->x, 0.5, 0, 1, 0.001, HadamardTrap())
0.9738507879228328
```
"""
struct HadamardTrap <: HadamardDiff end

################################################################
###                    Type definition done                  ###
################################################################

# FIXME: Need to verify the derivative value in the interval [x₀, x]
#=
Hadamard Left Rectangular computing algorithms
=#
function fracdiff(f::FunctionAndNumber, α, x₀, x, h::Float64, ::HadamardLRect)
    typeof(f) <: Number ? (x == 0 ? (return 0) : (return f/sqrt(pi*x))) : nothing
    x == 0 ? (return 0) : nothing
    N = round(Int, (x-x₀)/h)
    result = zero(Float64)

    @fastmath @inbounds @simd for i ∈ 0:N-1
        result += LRectCoeff(i, N, h, α, x₀)*f(x₀+i*h)
    end

    return 1/gamma(1-α)*result
end
function LRectCoeff(i, n, h, α, x₀)
    if 0 ≤ i ≤ n-2
        return ((log((x₀+n*h)/(x₀+i*h)))^(-α) - (log((x₀+n*h)/(x₀+(i+1)*h)))^(-α))
    elseif i == n-1
        return (log((x₀+n*h)/(x₀+i*h)))^(-α)
    end
end

#=
Hadamard Right Rectangular computing algorithm
=#
function fracdiff(f::FunctionAndNumber, α, x₀, x, h::Float64, ::HadamardRRect)
    typeof(f) <: Number ? (x == 0 ? 0 : f/sqrt(pi*x)) : nothing
    N = round(Int, (x-x₀)/h)
    result = zero(Float64)

    @fastmath @inbounds @simd for i ∈ 0:N-1
        result += LRectCoeff(i, N, h, α, x₀)*f(x₀+(i+1)*h)
    end

    return 1/gamma(1-α)*result
end
function RRectCoeff(i, n, h, α, x₀)
    if i == 0
        return 1/2*LRectCoeff(i, n, h, α, x₀)
    elseif 1 ≤ i ≤ n-1
        return 1/2*(LRectCoeff(i-1, n, h, α, x₀) + LRectCoeff(i, n, h, α, x₀))
    elseif i == n
        return 1/2*LRectCoeff(i-1, n, h, α, x₀)
    end
end

#=
Hadamard trapezoidal computing algorithm
=#
function fracdiff(f::FunctionAndNumber, α, x₀, x, h::Float64, ::HadamardTrap)
    typeof(f) <: Number ? (x == 0 ? 0 : f/sqrt(pi*x)) : nothing
    N = round(Int, (x-x₀)/h)

    result = zero(Float64)

    @fastmath @inbounds @simd for i ∈ 0:N
        result += RRectCoeff(i, N, h, α, x₀)*f(x₀+i*h)
    end

    return 1/gamma(1-α)*result
end


text:
## ExoplanetsSysSim/src/mission_constants.jl
## (c) 2015 Eric B. Ford

const global num_channels = 84
 const global num_modules = 42
 const global num_quarters = 17              # QUERY:  I'm favoring leaving out quarter 0, since that was engineering data.  Agree?
 const global num_cdpp_timescales = 1        # TODO SCI: Increase if incorporate CDPPs for multiple timescales, also LC/SC issue
 const global mission_data_span = 1459.789   # maximum(ExoplanetsSysSim.StellarTable.df[:dataspan])
 const global mission_duty_cycle = 0.8751    # median(ExoplanetsSysSim.StellarTable.df[:dutycycle])

 const global kepler_exp_time_internal  =  6.019802903/(24*60*60)    # https://archive.stsci.edu/kepler/manuals/archive_manual.pdf
 const global kepler_read_time_internal = 0.5189485261/(24*60*60)    # https://archive.stsci.edu/kepler/manuals/archive_manual.pdf
 const global num_exposures_per_LC = 270
 const global num_exposures_per_SC = 9
 const global LC_integration_time = kepler_exp_time_internal*num_exposures_per_LC
 const global SC_integration_time = kepler_exp_time_internal*num_exposures_per_SC
 const global LC_read_time = kepler_read_time_internal*num_exposures_per_LC
 const global SC_read_time = kepler_read_time_internal*num_exposures_per_SC
 const global LC_duration = LC_integration_time +  LC_read_time
 const global SC_duration = SC_integration_time +  SC_read_time
 const global LC_rate = 1.0/LC_duration
 const global SC_rate = 1.0/SC_duration

# Standard conversion factors on which unit system is based
const global AU_in_m_IAU2012 = 149597870700.0
 const global G_in_mks_IAU2015 = 6.67384e-11
 const global G_mass_sun_in_mks = 1.3271244e20
 const global G_mass_earth_in_mks = 3.986004e8
 const global sun_radius_in_m_IAU2015 = 6.9566e8
 const global earth_radius_eq_in_m_IAU2015 = 6.3781e6
 const global sun_mass_in_kg_IAU2010 = 1.988547e30

# Constants used by this code
const global sun_mass = 1.0
 const global earth_mass = G_mass_earth_in_mks/G_mass_sun_in_mks  # about 3.0024584e-6
 const global earth_radius = earth_radius_eq_in_m_IAU2015 / sun_radius_in_m_IAU2015 # about 0.0091705248
 const global rsol_in_au = sun_radius_in_m_IAU2015 / AU_in_m_IAU2012  # about 0.00464913034
 const global sec_in_day = 24*60*60
 const global grav_const = G_in_mks_IAU2015 * sec_in_day^2 * sun_mass_in_kg_IAU2010 / AU_in_m_IAU2012^3 # about 2.9591220363e-4 in AU^3/(day^2 Msol)


text:
export
    @newtype,
    @newimmutable,
    @stringinterface

# Given a type definition Foo <: AbstractString, or Foo{T} <: AbstractString, return the inner
# type of the new composite type..
#
# Case 1: ```Foo <: AbstractString```
#         In this case the inner type is ```AbstractString```
#
# Case 2: ```Foo{T} <: AbstractString```
#         In this case the inner type is ```T```
#
# Case 3: ```Foo{T <: AbstractString}```
#         In this case the inner type is ```T```
function inner_type_from_type_expr(typedef)
    new_typeexpr = typedef.args[1]
    if isa(new_typeexpr, Symbol)
        # Case 1:
        # First argument is a symbol, like Foo in "@newtype Foo <: AbstractsString". That means it's
        # not a generic expression. In that case we default to using the supertype (here
        # AbstractString) as inner type.
        field_type = typedef.args[3]
    else
        # First argument is an expression, like Foo{T} in "@newtype Foo{T} <: AbstractString". That
        # means it's a generic expression, so we use the T in Foo{T} or Foo{T <: AbstractString} as
        # the inner type.
        if isa(new_typeexpr.args[2], Symbol)
            # Case 2:
            field_type = new_typeexpr.args[2]
        else
            # Case 3:
            field_type = new_typeexpr.args[2].args[1]
        end
    end
    field_type
end

# Determine the type of the inner nd the supertype, given the type definition provided to the
# macro.
#
# Returns: (inner type, super type)
#
# If the supertype in the type definition is an abstract type, then that is the supertype, and we
# figure out the inner type using ```inner_type_from_type_expr()```. So, for instance, if we have
# ```@newtype Foo <: AbstractString```, then the supertype is ```AbstractString```, and by the
# rules of ```inner_type_from_type_expr()``` the inner type will also be ```AbstractString```.
# However, if the supertype is a concrete type, then we cannot subclass it, since Julia doesn't
# allow that. In this case we interpret it to mean that we want the inner type to be that cóncrete
# type, and we'll use the same supertype as the super type of the concrete type.
# For instance, ```Foo <: UTF8String``` will create a new type:
# ```
# immutable Foo <: AbstractString
#   v::UTF8String
# end
# ```
# because ```UTF8String <: AbstractString```.
function determine_types(typedef)
    local sname = typedef.args[3]
    local supertype = eval(sname)
    if isleaftype(supertype)
        return sname, super(supertype)
    else
        return inner_type_from_type_expr(typedef), sname
    end
end

macro newtype(typedef)
    local tname = typedef.args[1]
    local mytypes = determine_types(typedef)
    local inner_type = mytypes[1]
    local sname = mytypes[2]
    quote
        type $tname <: $sname
            v::$(inner_type)
        end
    end
end

macro newimmutable(typedef)
    local tname = typedef.args[1]
    local mytypes = determine_types(typedef)
    local inner_type = mytypes[1]
    local sname = mytypes[2]
    quote
        immutable $tname <: $sname
            v::$(inner_type)
        end
    end
end

macro stringinterface(typesym)
    quote
        Base.endof(x::$(esc(typesym))) = Base.endof(x.v)
        Base.next(x::$(esc(typesym)), i::Int) = Base.next(x.v, i)
    end
end


text:
module Checks
using  ...PatternStructure.Trees
import ...Helper: linesof
export EqualityCheck, TypeCheck, PredicateCheck,
       ArgsStep, QuoteStep, BlockStep,
       IterStep, SlurpStep

struct EqualityCheck{T} <: PatternCheck
  value::T
end

struct TypeCheck{T} <: PatternCheck
end

struct PredicateCheck <: PatternCheck
  predicate :: Function
end

struct ArgsStep  <: PatternStep end
struct BlockStep <: PatternStep end
struct QuoteStep <: PatternStep end
struct IterStep  <: PatternStep end
struct SlurpStep <: PatternStep end

(::ArgsStep)(ex)  = ex.args
(::BlockStep)(ex) = linesof(ex)

(::QuoteStep)(ex::Expr) = ex.args
(::QuoteStep)(ex::QuoteNode) = [ex.value]

(::IterStep)(ex)  = ex
(::SlurpStep)(ex) = error("Called the step of slurp.")

end


text:
"""julia
      rotx(θ::Real)

A matrix representing the rotation about the X-axis by an angle θ.
"""
function rotx(θ::Real)
cosθ = cos(θ);
sinθ = sin(θ);
𝐑ˣ = [1.0  0.0      0.0  ;
      0.   cosθ    -sinθ ;
      0    sinθ     cosθ ]
end

"""juliaa
      roty(θ::Real)

A matrix representing rotation about the Y-axis by an angle θ.
"""
function roty(θ::Real)
cosθ = cos(θ);
sinθ = sin(θ);
𝐑ʸ = [cosθ     0.0    sinθ ;
      0.0      1.0    0.0  ;
      -sinθ    0.0    cosθ]
end

"""julia
      rotz(θ::Real)

A matrix representing the rotation about the Z-axis by an angle θ.
"""
function rotz(θ::Real)
cosθ = cos(θ);
sinθ = sin(θ);
𝐑ᶻ = [cosθ     -sinθ    0.0 ;
      sinθ      cosθ    0.0 ;
      0.0       0.0     1.0 ]
end

"""julia
      rotz(θ::Real)

A matrix representing the rotation about the X-axis, Y-axis and Z-axis by the angles θˣ, θʸ, and θᶻ respectively.
"""
function rotxyz(θˣ::Real,θʸ::Real,θᶻ::Real)
      rotx(θˣ)*roty(θʸ)*rotz(θᶻ)
end

"""
    rotate!(camera::CameraModel, 𝐑::AbstractArray)

Rotates a camera around its optical center by a rotation matrix 𝐑.
"""
function rotate!(camera::CameraModel, 𝐑::AbstractArray)
    camera.𝐞₁ = 𝐑₂*camera.𝐞₁
    camera.𝐞₂ = 𝐑₂*camera.𝐞₂
    camera.𝐞₃ = 𝐑₂*camera.𝐞₃
    camera
end

function rodrigues2matrix(vˣ::Real,vʸ::Real,vᶻ::Real)
      𝐯 = [vˣ, vʸ, vᶻ]
      θ = norm(𝐯)
      𝐯 = θ == 0 ? 𝐯 : 𝐯/θ
      𝐈 = Matrix(1.0I, 3, 3)
      𝐖 = vec2antisym(𝐯)
      𝐑 = 𝐈 + 𝐖 * sin(θ) + 𝐖^2 * (1-cos(θ))
end


text:
"""
    L₀affine(x̄, b, bc::Tuple{BoundaryCondition, BoundaryCondition})

    Returns a vector for affine operator given an extended grid `x̄` and boundary conditions `bc`
to solve the system `L₀ * v(̄x) = b` on interior nodes.

    Returns a `length(interiornodes(x̄))`-length vector such that solving
`L₀ * v(x̄) = b` on interior nodes `x` under the boundary condition given by `bc`
is identical with `L₀bc * v(x) = b + L₀affine(x̄, bc)` The first element of `bc` is
applied to the lower bound, and second element of `bc` to the upper.

# Examples
```jldoctest; setup = :(using SimpleDifferentialOperators)
julia> x̄ = 0:5
0:5

julia> b = 1:3
1:3

julia> L₀affine(x̄, b, (Reflecting(), Reflecting()))
4-element Array{Float64,1}:
 0.0
 0.0
 0.0
 0.0

julia> L₀affine(x̄, b, (NonhomogeneousAbsorbing(1.0), Reflecting()))
4-element Array{Float64,1}:
 0.0
 0.0
 0.0
 0.0

julia> L₀affine(x̄, b, (NonhomogeneousAbsorbing(1.0), NonhomogeneousAbsorbing(2.0)))
4-element Array{Float64,1}:
 0.0
 0.0
 0.0
 0.0
```
"""
function L₀affine(x̄, b, bc::Tuple{BoundaryCondition, BoundaryCondition})
    x = interiornodes(x̄)
    M = length(x)
    return zeros(M)
end

"""
    L₁₋affine(x̄, b, bc::Tuple{BoundaryCondition, BoundaryCondition})

    Returns a vector for affine operator given an extended grid `x̄` and boundary conditions `bc`
to solve the system `L₁₋ * v(̄x) = b` on interior nodes.

    Returns a `length(interiornodes(x̄))`-length vector such that solving
`L₁₋ * v(x̄) = b` on interior nodes `x` under the boundary condition given by `bc`
is identical with `L₁₋bc * v(x) = b + L₁₋affine(x̄, bc)`. The first element of `bc` is
applied to the lower bound, and second element of `bc` to the upper.

# Examples
```jldoctest; setup = :(using SimpleDifferentialOperators)
julia> x̄ = 0:5
0:5

julia> b = 1:3
1:3

julia> L₁₋affine(x̄, b, (Reflecting(), Reflecting()))
4-element Array{Float64,1}:
 0.0
 0.0
 0.0
 0.0

julia> L₁₋affine(x̄, b, (NonhomogeneousAbsorbing(1.0), Reflecting()))
4-element Array{Float64,1}:
 1.0
 0.0
 0.0
 0.0

julia> L₁₋affine(x̄, b, (NonhomogeneousAbsorbing(1.0), NonhomogeneousAbsorbing(2.0)))
4-element Array{Float64,1}:
 1.0
 0.0
 0.0
 0.0
```
"""
function L₁₋affine(x̄, b, bc::Tuple{BoundaryCondition, BoundaryCondition})
    x = interiornodes(x̄)
    M = length(x)
    b_affine = zeros(M)

    # apply affine operator if needed
    Δ_1m = x̄[2] - x̄[1]
    b_affine[1] = (typeof(bc[1]) <: NonhomogeneousAbsorbing) ? bc[1].S / Δ_1m : b_affine[1]

    if (typeof(bc[1]) <: Absorbing && bc[1].loc > 0)
        b_affine[1:min(M, bc[1].loc)] = -b[1:min(M, bc[1].loc)]
    end

    return b_affine
end
"""
    L₁₊affine(x̄, b, bc::Tuple{BoundaryCondition, BoundaryCondition})

    Returns a vector for affine operator given an extended grid `x̄` and boundary conditions `bc`
to solve the system `L₁₊ * v(̄x) = b` on interior nodes.

    Returns a `length(interiornodes(x̄))`-length vector such that solving
`L₁₊ * v(x̄) = b` on interior nodes `x` under the boundary condition given by `bc`
is identical with `L₁₊bc * v(x) = b + L₁₊affine(x̄, bc)` The first element of `bc` is
applied to the lower bound, and second element of `bc` to the upper.

# Examples
```jldoctest; setup = :(using SimpleDifferentialOperators)
julia> x̄ = 0:5
0:5

julia> b = 1:3
1:3

julia> L₁₊affine(x̄, b, (Reflecting(), Reflecting()))
4-element Array{Float64,1}:
 0.0
 0.0
 0.0
 0.0

julia> L₁₊affine(x̄, b, (NonhomogeneousAbsorbing(1.0), Reflecting()))
4-element Array{Float64,1}:
 0.0
 0.0
 0.0
 0.0

julia> L₁₊affine(x̄, b, (NonhomogeneousAbsorbing(1.0), NonhomogeneousAbsorbing(2.0)))
4-element Array{Float64,1}:
  0.0
  0.0
  0.0
 -2.0
```
"""
function L₁₊affine(x̄, b, bc::Tuple{BoundaryCondition, BoundaryCondition})
    x = interiornodes(x̄)
    M = length(x)
    b_affine = zeros(M)

    # apply affine operator if needed
    Δ_Mp = x̄[end] - x̄[end-1]
    b_affine[end] = (typeof(bc[2]) <: NonhomogeneousAbsorbing) ? -bc[2].S / Δ_Mp : b_affine[end]

    if (typeof(bc[1]) <: Absorbing && bc[1].loc > 0)
        b_affine[1:min(M, bc[1].loc)] = -b[1:min(M, bc[1].loc)]
    end

    return b_affine
end

"""
    L₂affine(x̄, bc::Tuple{BoundaryCondition, BoundaryCondition})

    Returns a vector for affine operator given an extended grid `x̄` and boundary conditions `bc`
to solve the system `L₂ * v(̄x) = b` on interior nodes.

    Returns a `length(interiornodes(x̄))`-length vector such that solving
`L₂ * v(x̄) = b` on interior nodes `x` under the boundary condition given by `bc`
is identical with `L₂bc * v(x) = b + L₂affine(x̄, bc)` The first element of `bc` is
applied to the lower bound, and second element of `bc` to the upper.

# Examples
```jldoctest; setup = :(using SimpleDifferentialOperators)
julia> x̄ = 0:5
0:5

julia> b = 1:3
1:3

julia> L₂affine(x̄, b, (Reflecting(), Reflecting()))
4-element Array{Float64,1}:
 0.0
 0.0
 0.0
 0.0

julia> L₂affine(x̄, b, (NonhomogeneousAbsorbing(1.0), Reflecting()))
4-element Array{Float64,1}:
-1.0
 0.0
 0.0
 0.0

julia> L₂affine(x̄, b, (NonhomogeneousAbsorbing(1.0), NonhomogeneousAbsorbing(2.0)))
4-element Array{Float64,1}:
-1.0
 0.0
 0.0
-2.0
```
"""
function L₂affine(x̄, b, bc::Tuple{BoundaryCondition, BoundaryCondition})
    x = interiornodes(x̄)
    M = length(x)
    b_affine = zeros(M)

    # apply affine operator if needed
    Δ_1p = x̄[3] - x̄[2]
    Δ_1m = x̄[2] - x̄[1]
    Δ_Mp = x̄[end] - x̄[end-1]
    Δ_Mm = x̄[end-1] - x̄[end-2]
    b_affine[1] = (typeof(bc[1]) <: NonhomogeneousAbsorbing) ? -2*bc[1].S/(Δ_1m*(Δ_1m+Δ_1p)) : b_affine[1]
    b_affine[end] = (typeof(bc[2]) <: NonhomogeneousAbsorbing) ? -2*bc[2].S/(Δ_Mp*(Δ_Mm+Δ_Mp)) : b_affine[end]

    if (typeof(bc[1]) <: Absorbing && bc[1].loc > 0)
        b_affine[1:min(M, bc[1].loc)] = -b[1:min(M, bc[1].loc)]
    end

    return b_affine
end

"""
    Laffine(L, bc::Tuple{BoundaryCondition, BoundaryCondition})

    Returns a vector for affine operator given an extension operator `L` and boundary conditions `bc`
to solve the system `L * v(̄x) = b` on interior nodes.

    Returns a `size(L)[1]`-length vector such that solving
`L * v(x̄) = b` on interior nodes `x` under the boundary condition given by `bc`
is identical with `Lbc * v(x) = b + Laffine(L, bc)` The first element of `bc` is
applied to the lower bound, and second element of `bc` to the upper.

# Examples
```jldoctest; setup = :(using SimpleDifferentialOperators)
julia> x̄ = 0:5
0:5

julia> b = 1:3
1:3

julia> L = L₂(x̄)
4×6 SparseArrays.SparseMatrixCSC{Float64,Int64} with 12 stored entries:
  [1, 1]  =  1.0
  [1, 2]  =  -2.0
  [2, 2]  =  1.0
  [1, 3]  =  1.0
  [2, 3]  =  -2.0
  [3, 3]  =  1.0
  [2, 4]  =  1.0
  [3, 4]  =  -2.0
  [4, 4]  =  1.0
  [3, 5]  =  1.0
  [4, 5]  =  -2.0
  [4, 6]  =  1.0

julia> Laffine(L, b, (Reflecting(), Reflecting()))
4-element Array{Float64,1}:
 0.0
 0.0
 0.0
 0.0

julia> Laffine(L, b, (NonhomogeneousAbsorbing(1.0), Reflecting()))
4-element Array{Float64,1}:
-1.0
 0.0
 0.0
 0.0

julia> Laffine(L, b, (NonhomogeneousAbsorbing(1.0), NonhomogeneousAbsorbing(2.0)))
4-element Array{Float64,1}:
-1.0
 0.0
 0.0
-2.0
```
"""
# take M by (M+2) extended operator
function Laffine(L, b, bc::Tuple{BoundaryCondition, BoundaryCondition})
    M = size(L)[1]
    b_affine = zeros(M)

    # apply affine operator if needed
    b_affine[1] = (typeof(bc[1]) <: NonhomogeneousAbsorbing) ? -sum(L[:,1])*bc[1].S : b_affine[1]
    b_affine[end] = (typeof(bc[2]) <: NonhomogeneousAbsorbing) ? -sum(L[:,end])*bc[2].S : b_affine[end]

    if (typeof(bc[1]) <: Absorbing && bc[1].loc > 0)
        b_affine[1:min(M, bc[1].loc)] = -b[1:min(M, bc[1].loc)]
    end

    return b_affine
end


text:
"""
    demand(price, s::Student;p_neigh_parm::Int64=0)

Compute the total demand for classes based on the optimal individual bundles for a `student`.

## Arguments

Assume that N is the number of students, and C the number of classes that are available to choose from.

- `price` : a column vector of dimension C x 1, which i-th element is the price that was assigned to class i.
- `student`: instance of type `Student`

## Example

* 20 individuals
* 10 classes to choose from
* exactly 3 classes to attend
* NO CROSS PREFERENCES, i.e. Individual preferences are generated as an array of sparse diagonal matrices

```
julia> Ind_pref = []
0-element Array{Any,1}

julia> for i in 1:20
                push!(Ind_pref,
                sparse(collect(1:10),
                    collect(1:10),
                    rand(0:100, 10))
                )
       end

julia> Ind_pref
20-element Array{Any,1}:

  [1 ,  1]  =  63
  [2 ,  2]  =  51
  [3 ,  3]  =  39
  [4 ,  4]  =  85
  [5 ,  5]  =  35
  [6 ,  6]  =  60
  [7 ,  7]  =  84
  [8 ,  8]  =  96
  [9 ,  9]  =  61
  [10, 10]  =  56

  ⋮

  [1 ,  1]  =  30
  [2 ,  2]  =  9
  [3 ,  3]  =  9
  [4 ,  4]  =  58
  [5 ,  5]  =  25
  [6 ,  6]  =  82
  [7 ,  7]  =  0
  [8 ,  8]  =  75
  [9 ,  9]  =  94
  [10, 10]  =  36

julia> ind_budget = rand(150:200, 20)
20-element Array{Int64,1}:
  162
  195
  ⋮
  193

julia> price = rand(10)
10-element Array{Float32,1}:
  0.786705
  0.735436
  ⋮
  0.880273

julia> cap = fill(3, 20)
20-element Array{Int64,1}:
  3
  ⋮
  3

julia> CourseMatch.demand(price, Ind_pref, ind_budget, cap)

julia> d[:ind_demands]
20×10 Array{Int64,2}:
 0  0  0  1  1  0  1  0  0  0
 0  0  0  0  1  1  1  0  0  0
 0  1  0  0  1  1  0  0  0  0
 0  1  0  0  1  1  0  0  0  0
 0  0  0  1  0  1  0  0  0  1
 0  0  0  1  0  0  1  0  1  0
 0  1  0  0  0  0  0  1  0  1
 0  0  1  0  0  0  0  1  1  0
 1  0  0  0  0  1  0  0  1  0
 1  0  0  1  0  0  0  0  1  0
 0  0  1  0  0  0  0  1  1  0
 1  0  0  0  0  0  1  1  0  0
 0  0  1  0  1  1  0  0  0  0
 0  0  0  1  0  0  1  1  0  0
 0  0  0  0  1  1  0  1  0  0
 0  1  0  1  0  0  0  1  0  0
 1  0  0  1  0  1  0  0  0  0
 1  0  0  1  0  1  0  0  0  0
 0  1  1  0  0  0  0  0  0  1
 0  1  0  1  0  0  0  0  0  1

 julia> d[:course_demand]
1×10 Array{Int64,2}:
 5  6  4  9  6  9  5  7  5  4

```

"""
function demand!(price::Vector{Float64}, s::Student; p_neigh_parm::Int64=0)

	N = size(s.preferences)[1]  # number of choices

	# Maximization problem
	m = Model(solver=CbcSolver())

	@variable(m, x[1:N], Bin)

	# Objective: maximize utility
	@objective(m, Max, x'*s.preferences*x)

	# Constraints:

	#Should not spend more than one's budget
	@constraint(m,  dot(price, x) <= s.budget )

	#Time Constraints
	@constraint(m, x'*s.constraints.clash*x == 0 )

	#Mandatory courses
	@constraint(m, s.constraints.mandatory'*x == sum(s.constraints.mandatory))

	#Tronc commun courses
	#TC courses
	@constraint(m, s.constraints.isTC'*x .- s.req_TC .>= 0)

	#TC program constraint
	@constraint(m, s.constraints.notTCprogram'*x .== 0 )

	# TODO
	# not working from here on: need to change how those constraints are called now.
	# they are members of s.constraints now!

	#TC semester constraint
	@constraint(m, s.tc_cour_sem_const'*x .== 0 )

	#Formation commune courses
	#FC courses
	@constraint(m, s.fc_cour_const'*x .- s.fc_requirement .>= 0)

	#FC program constraint
	@constraint(m, s.fc_cour_prog_const'*x .== 0 )

	#FC semester constraint
	@constraint(m, s.fc_cour_sem_const'*x .== 0 )

	#Electives courses
	#TC courses
	@constraint(m, s.el_cour_const'*x .- s.el_requirement .>= 0)

	#TC program constraint
	@constraint(m, s.el_cour_prog_const'*x .== 0 )

	#TC semester constraint
	@constraint(m, s.el_cour_sem_const'*x .== 0 )

	#Neighboring prices constraint
	# require that course id p_neigh_parm must be chosen
	if p_neigh_parm ~= 0
		@constraint(m, x[p_neigh_parm] == 1)
	end

	# Solve problem using MIP solver
	status = solve(m)
	dem = getvalue(x)

	# store on the student
	if status == :Optimal
		s.allocation = getvalue(x)
	else
		warn("MIP problem not solved. status = $status")
	end
end

"""
	demand for all students
"""
function demand(price::Vector{Float64}, s::Array{Student})

	N = size(s[1].pref)[1]
	M = length(s)
    outdem = zeros(Int,M,N)

	for (i,stu) in enumerate(s)
		outdem[i,:] = demand(price,stu)
	end

  	return Dict(:ind_demands => out_dem,:course_demand => sum(out_dem,1),:total=>sum(out_dem) )
	# return Dict(:ind_demands => demand,:course_demand => sum(demand,1),:total=>sum(demand),:clearing_error=>err )

end

# - `price` : a column vector of dimension C x 1, which i-th element is the price that was assigned to class i.
# - `Student.pref` : an array containing S elements, and such that its n-th element is the (sparse) matrix representing student n preferences. Each of the matrices contained in that array should be a squared matrix of dimension C.
# - `Student.budget` : a column vector of dimension N x 1, which n-th element is the budget that was allocated to the n-th student.
# - `Student.time_const` : a sparse array C x C that flags with 1 if the ith element HAS a clash with jth element. 0 otherwise.
# - `Student.mand_cour_const` : a vector C x 1. the i-th element is 1 if the course must be taken, 0 otherwise
# - `Student.tc_cour_prog_const` : a vector C x 1. the i-th element is 1 if the course is NOT in the TC of the students's program. 0 otherwise.
# - `Student.tc_cour_sem_const` : a vector C x 1. the i-th element is 1 if the course is NOT in the TC of the students's semester
# - `Student.tc_cour_const` : a vector C x 1. the i-th element is 1 if the course is part of the TC for the student
# - `Student.tc_requirement` : number of TC courses that the student is required to take
# - `Student.fc_cour_prog_const` : a vector C x 1. the i-th element is 1 if the course is NOT in the FC of the students's program. 0 otherwise.
# - `Student.fc_cour_sem_const` : a vector C x 1. the i-th element is 1 if the course is NOT in the FC of the students's semester
# - `Student.fc_cour_const` : a vector C x 1. the i-th element is 1 if the course is part of the FC for the student
# - `Student.fc_requirement` : number of FC courses that the student is required to take
# - `Student.el_cour_prog_const` : a vector C x 1. the i-th element is 1 if the course is NOT in the electives of the students's program. 0 otherwise.
# - `Student.el_cour_sem_const` : a vector C x 1. the i-th element is 1 if the course is NOT in the electives of the students's semester
# - `Student.el_cour_const` : a vector C x 1. the i-th element is 1 if the course is part of the electives choice for the student
# - `Student.el_requirement` : number of electives courses that the student is required to take
# - `p_neigh_parm` : course that needs to be dropped in order to compute neigboring prices
# ## Example
# * 20 individuals
# * 15 classes to choose from
# * Random classes to attend between 4 and 6
# * Random number of mandatory classes to attend with an average of 2 per student
# * Two random time collisions per student in non?mandatory courses
# * Two random courses that are not allowed to be taken by the program of the student
# * Two random courses that are not allowed to be taken by the semester of the student
# * Two random courses are electives for each student, they are required to take 1
# * NO CROSS PREFERENCES, i.e. Individual preferences are generated as an array of sparse diagonal matrices
# ```
# num_stud = 20
# num_cour = 15
# mkt_demand = []
# for i in 1:num_stud
# 	# Individual preferences are generated as an array of sparse diagonal matrices
# 	Ind_pref = sparse(collect(1:num_cour), collect(1:num_cour), rand(-100:100, num_cour))
# 	# Individual budget
# 	ind_budget = 150+rand()
# 	# Price vector
# 	price = 100*rand(Float32, size(Ind_pref,1))
# 	# Capacity vector
# 	cap = rand(4:6)
# 	#Mandatory courses
# 	mand_cour_const = rand(Binomial(1,0.13), size(Ind_pref,1))
# 	#Flag non Mandatory courses
# 	non_mand = findn(mand_cour_const .== 0)
# 	#Schedule collisions
# 	time_const = zeros(size(Ind_pref,1), size(Ind_pref,1))
# 	time_const[non_mand[1],non_mand[2]] = time_const[non_mand[2],non_mand[1]] = 1
# 	time_const[non_mand[3],non_mand[4]] = time_const[non_mand[4],non_mand[3]] = 1
# 	#Program courses
# 	prog_cour_const = zeros(size(Ind_pref,1),1)
# 	prog_cour_const[non_mand[5]] = 1
# 	prog_cour_const[non_mand[6]] = 1
# 	#Semester courses
# 	sem_cour_const = zeros(size(Ind_pref,1),1)
# 	sem_cour_const[non_mand[7]] = 1
# 	sem_cour_const[non_mand[8]] = 1
# 	#Elective courses
# 	elec_cour_const = zeros(size(Ind_pref,1),1)
# 	elec_cour_const[non_mand[9]] = 1
# 	elec_cour_const[non_mand[10]] = 1
# 	#Number of elective courses
# 	num_elec_cour = 1
# 	#Demand computation
# 	dem = demand(price, Ind_pref, ind_budget, cap, time_const, mand_cour_const, prog_cour_const, sem_cour_const, elec_cour_const, num_elec_cour)
# 	push!(mkt_demand, dem)
# 	return mkt_demand


text:
#==============================================================================#
# XMLDict.jl
#
# Copyright Sam O'Connor 2015 - All rights reserved
#==============================================================================#

__precompile__()

module XMLDict

export parse_xml, xml_dict

using EzXML
using DataStructures
using Base.Iterators
using IterTools

#-------------------------------------------------------------------------------
# Associative type wrapper.
#-------------------------------------------------------------------------------

mutable struct XMLDictElement <: AbstractDict{Union{String,Symbol},Any}
    x
    doc
    g
end

wrap(x, doc) = XMLDictElement(x, doc, nothing)
wrap(l::Vector, doc) = [wrap(i, doc) for i in l]

Base.get(x::XMLDictElement, k, d=nothing) = XMLDict.get(x.x, x.doc, k, d)

Base.keys(x::XMLDictElement) = XMLDict.keys(x.x)
Base.length(x::XMLDictElement) = length(collect(keys(x)))

function Base.iterate(x::XMLDictElement)
    x.g = (Pair{Union{String,Symbol},Any}(n, (get(x, n))) for n in keys(x))
    iterate(x.g)
end
function Base.iterate(x::XMLDictElement, isdone::Any)
    iterate(x.g, isdone)
end
xml_dict(x, args...; options...) = xml_dict(x.x, args...; options...)

Base.show(io::IO, x::XMLDictElement) = show(io, x.x)

#-------------------------------------------------------------------------------
# XML Parsing.
#-------------------------------------------------------------------------------

# Parse "xml" string into EzXML.Document object.

function parse_xml(xml::AbstractString)
    doc = parsexml(xml)
    return wrap(doc, doc)
end

#-------------------------------------------------------------------------------
# Dynamic Associative Implementation for XMLElement
#-------------------------------------------------------------------------------

function XMLDict.keys(x::EzXML.Node)
    children = distinct(nodename(c) for c in eachelement(x))
    attribs = (Symbol(nodename(c)) for c in eachattribute(x))
    return Iterators.flatten((children, attribs))
end
XMLDict.keys(x::EzXML.Document) = keys(root(x))

# Get sub-elements that match tag.
# For leaf-nodes return element content (text).

function XMLDict.get(x::EzXML.Node, doc::EzXML.Document, tag::AbstractString, default)

    if tag == ""
        return strip(nodecontent(x))
    end

    if !haselement(x)
        return default
    end
    l = filter(el->nodename(el) == tag, elements(x))
    if isempty(l)
        return default
    end
    if countelements(l[1]) == countattributes(l[1]) == 0
        l = map(strip∘nodecontent, l)
    else
        l = wrap(l, doc)
    end
    return length(l) == 1 ? l[1] : l
end

# Get element attribute by "name".

function XMLDict.get(x::EzXML.Node, doc::EzXML.Document, name::Symbol, default)
    s = string(name)
    haskey(x, s) ? x[s] : default
end

# Wrapper for EzXML.Document.

function XMLDict.get(x::EzXML.Document, doc::EzXML.Document, tag, default)
    return XMLDict.get(root(x), doc, tag, default)
end

#-------------------------------------------------------------------------------
# Convert entire XMLDocument to OrderedDict...
#-------------------------------------------------------------------------------

# Return Dict representation of "xml" string.

function xml_dict(xml::AbstractString, dict_type::Type=OrderedDict; options...)
    doc = parse_xml(xml)
    r = xml_dict(doc, dict_type; options...)
    return r
end

# Return Dict representation of EzXML.Document.

function xml_dict(xml::EzXML.Document, dict_type::Type=OrderedDict; options...)
    r = dict_type()
    r[:version] = version(xml)
    try
        r[:encoding] = encoding(xml)
    catch
    end
    r[nodename(root(xml))] = xml_dict(root(xml), dict_type; options...)
    r
end

# Does this node have any text?

is_text(x::EzXML.Node) = istext(x) || iscdata(x)
function is_empty(x::EzXML.Node)
    c = nodecontent(x)
    return isempty(c) || all(isspace, c)
end
has_text(x::EzXML.Node) = is_text(x) && !is_empty(x)

# Return Dict representation of XMLElement.

function xml_dict(x::EzXML.Node, dict_type::Type=OrderedDict; strip_text=false)

    # Copy element attributes into dict...
    r = dict_type()
    for a in eachattribute(x)
        r[Symbol(nodename(a))] = nodecontent(a)
    end

    # Check for non-empty text nodes under this element...
    element_has_text = any(has_text, eachnode(x))

    # Check for non-contiguous repetition of sub-element tags...
    element_has_mixed_tags = false
    tags = []
    for c in eachelement(x)
        tag = nodename(c)
        if isempty(tags) || tag != tags[end]
            if tag in tags
                element_has_mixed_tags = true
                break
            end
            push!(tags, tag)
        end
    end

    # The empty-string key holds a vector of sub-elements.
    # This is necessary when grouping sub-elements would alter ordering...
    if element_has_text || element_has_mixed_tags
        r[""] = Any[]
    end

    for c in eachnode(x)

        if iselement(c)

            # Get name and sub-dict for sub-element...
            n = nodename(c)
            v = xml_dict(c, dict_type; strip_text=strip_text)

            if haskey(r, "")

                # If this is a text element, embed sub-dict in text vector...
                # "The <b>bold</b> tag" == ["The", Dict("b" => "bold"), "tag"]
                push!(r[""], dict_type(n => v))

            elseif haskey(r, n)

                # Collect sub-elements with same tag into a vector...
                # "<item>foo</item><item>bar</item>" == "item" => ["foo", "bar"]
                a = isa(r[n], Array) ? r[n] : Any[r[n]]
                push!(a, v)
                r[n] = a
            else
                r[n] = v
            end

        elseif is_text(c) && haskey(r, "")
            push!(r[""], nodecontent(c))
        end
    end

    # Collapse leaf-node vectors containing only text...
    if haskey(r, "")
        v = r[""]

        # If the vector contains a single text element, collapse the vector...
        if length(v) == 1 && isa(v[1], AbstractString)
            if strip_text
                v[1] = strip(v[1])
            end
            r[""] = v[1]

            # If "r" contains no other keys, collapse the "" key...
            if length(r) == 1
                r = r[""]
            end
        end
    end

    return r
end

#-------------------------------------------------------------------------------
# Convert Dict (produced by xml_dict) back to an XML string.
#
# dict_xml(xml_dict(xml_string)) ~= xml_string
#-------------------------------------------------------------------------------

function dict_xml(root::AbstractDict)
    string("<?xml", attr_xml(root), "?>\n", node_xml(root))
end

attrs(node::AbstractDict) = filter(pair->isa(first(pair), Symbol), node)
nodes(node::AbstractDict) = filter(pair->!isa(first(pair), Symbol), node)

function attr_xml(node::AbstractDict)
    join([" $n=\"$v\"" for (n,v) in attrs(node)])
end

attr_xml(node) = ""

node_xml(node) = join([node_xml(n,v) for (n,v) in nodes(node)])

function node_xml(name::AbstractString, value::AbstractArray)
    value_xml(name != "" ? [Dict(name=>i) for i in value] : value)
end

function node_xml(name::AbstractString, node)
    a = attr_xml(node)
    v = value_xml(node)
    name == "" ? v : string("<", name, a, v == "" ? "/>" : ">$v</$name>")
end

value_xml(value::AbstractDict) = node_xml(value)

value_xml(value::AbstractArray) = join(map(value_xml, value))

value_xml(value::AbstractString) = escape(value)

function escape(s::AbstractString)
    p = ccall((:xmlEncodeEntitiesReentrant, EzXML.libxml2), Ptr{UInt8},
              (Ptr{Cvoid}, Cstring), C_NULL, s)
    p == C_NULL ? "" : unsafe_string(p)
end

end # module XMLDict

#==============================================================================#
# End of file
#==============================================================================#


text:
using ShirleyRayTracer
using Images

#==

# launch using 40 threads
 julia -t 40 --project=.. -L RandomScene.jl  -e "main()" > random_scene.ppm

# launch force using 1 thread
 julia -t 1 --project=.. -L RandomScene.jl  -e "main()" > random_scene.ppm

# launch using thread count from environment
 julia --project=.. -L RandomScene.jl  -e "main()" > random_scene.ppm

==#

function add_random_scene!(scene::Scene)

	add!(scene, Sphere(Point3(0,-1000,0), 1000, Lambertian(0.5, 0.5, 0.5)))

	rand_material(p) = if p < 0.8
				Lambertian()
			elseif p < 0.95
				rf = randf(0.5, 1)
				Metal(rf, rf, rf, 0.5rand())
			else
				Dielectric(1.5)
			end

	for a in -11:10, b in -11:10
		center = Point3(a + 0.9rand(), 0.2, b + 0.9rand())
		if ShirleyRayTracer.magnitude(center - Point3(4, 0.2, 0)) > 0.9
			add!(scene, Sphere(center, 0.2, rand_material(rand())))
		end
	end

	add!(scene, Sphere(Point3(0, 1, 0), 1.0, Dielectric(1.5)))
	add!(scene, Sphere(Point3(-4, 1, 0), 1.0, Lambertian(0.4,0.2,0.1)))
	add!(scene, Sphere(Point3(4, 1, 0), 1.0, Metal(0.7,0.6,0.5, 0.0)))
end

function main(;filename="render.jpg", image_width=1200, aspect_ratio=16/9, samples_per_pixel=10, max_depth=50)
	image_height = round(Int, image_width / aspect_ratio)
	world = Scene(Camera(Point3(13.,2.,3.), zero(Point3), Vec3(0,1,0), 20, aspect_ratio, 0.1, 10.0))
	add_random_scene!(world)
	image = render(world, image_width, image_height, samples_per_pixel, max_depth)
	save(filename, image)
end


text:
export hwe_test, hwe

"""
    locus_chi_sq(locus::T) where T <: GenoArray
Calculate the chi square statistic and p-value for a locus
Returns a tuple with chi-square statistic, degrees of freedom, and p-value
"""
function chisq_locus(locus::T) where T <: GenoArray
    ## Get expected number of genotypes in a locus
    expected = geno_count_expected(locus)

    ## Get observed number of genotypes in a locus
    observed = geno_count_observed(locus)

    chisq_stat = expected
    @inbounds for genotype in keys(expected)
        o = get(observed, genotype, 0)
        e = get(expected, genotype, 0)
        chisq_stat[genotype] = (o - e)^2 / e
    end
    chisq_stat = values(chisq_stat) |> sum
    n_geno_exp = length(expected)
    n_alleles = length(unique_alleles(locus))
    df = n_geno_exp - n_alleles

    if df > 0
        chisq_dist = Distributions.Chisq(df)
        p_val = 1.0 - Distributions.cdf(chisq_dist, chisq_stat)
    else
        p_val = missing
    end
    return (chisq_stat, df, p_val)
end

"""
    hwe_test(data::PopData; by::String = "locus"; correction = "none")
Calculate chi-squared test of HWE for each locus and returns observed and
expected heterozygosity with chi-squared, degrees of freedom and p-values for
each locus. Use `by = "population"` to perform this separately for each population
(default: `by = "locus"`). Use `correction =` to specify a P-value
adjustment method for multiple testing.

#### example
`hwe_test(@gulfsharks, correction = "bh")` \n
`hwe_test(@gulfsharks, by = "population", correction = "bh")` \n

### `correction` methods (case insensitive)
- `"bonferroni"` : Bonferroni adjustment
- `"holm"` : Holm adjustment
- `"hochberg"` : Hochberg adjustment
- `"bh"` : Benjamini-Hochberg adjustment
- `"by"` : Benjamini-Yekutieli adjustment
- `"bl"`  : Benjamini-Liu adjustment
- `"hommel"` : Hommel adjustment
- `"sidak"` : Šidák adjustment
- `"forwardstop"` or `"fs"` : Forward-Stop adjustment
- `"bc"` : Barber-Candès adjustment
"""
@inline function hwe_test(data::PopData; by::String = "locus", correction::String = "none")
    if by == "locus"
        tmp =DataFrames.combine(
            groupby(data.loci, :locus),
            :genotype => chisq_locus => :chisq
        )
        out_table = select(tmp, :locus,
            :chisq => (i -> getindex.(i, 1)) => :chisq,
            :chisq => (i -> getindex.(i, 2)) => :df,
            :chisq => (i -> getindex.(i, 3)) => :P
        )
    else
        tmp =DataFrames.combine(
            groupby(data.loci, [:locus, :population]),
            :genotype => chisq_locus => :chisq
        )
        out_table = select(tmp, :locus, :population,
            :chisq => (i -> getindex.(i, 1)) => :chisq,
            :chisq => (i -> getindex.(i, 2)) => :df,
            :chisq => (i -> getindex.(i, 3)) => :P
        )
    end
    if correction == "none"
        return out_table
    else
        column_name = Symbol("P_"*correction)
        transform!(out_table, :P => (i -> multitest_missing(i, correction)) => column_name)
    end
end

const hwe = hwe_test


text:
@testsection "Isotropic tensors" begin
    for T ∈ (Sym, Float64), dim ∈ (2, 3)
        @testsection "type $T, dim $dim" begin
            𝟏 = tensId2(Val(dim), Val(T))
            @test opequal(𝟏, I)
            @test opequal(tr(𝟏), dim)

            if T == Sym
                α = symbols("α", real = true)
            else
                α = rand()
                while opequal(α, one(T)) || opequal(α, zero(T))
                    α = rand()
                end
            end
            @test α * 𝟏 isa TensISO{2}
            t = α * 𝟏 + (1 - α) * 𝟏
            @test opequal(t, 𝟏)
            @test t isa TensISO{2}
            @test opequal(inv(α * 𝟏), inv(α) * 𝟏)
            @test isISO(t)

            𝕀, 𝕁, 𝕂 = ISO(Val(dim), Val(T))
            @test opequal(𝕁 + 𝕂, 𝕀)
            @test opequal((𝟏 ⊗ 𝟏) / dim, 𝕁)
            if T == Sym
                α, β = symbols("α β", real = true)
            else
                α = rand()
                β = rand()
                while opequal(α + β, one(T)) || opequal(α, zero(T)) || opequal(β, zero(T))
                    α = rand()
                    β = rand()
                end
            end
            𝕋 = α * 𝕁 + β * 𝕂
            @test 𝕋 isa TensISO{4}
            @test opequal(inv(𝕋), inv(α) * 𝕁 + inv(β) * 𝕂)
            @test isISO(𝕋)

            if T == Sym && dim == 3
                E, ν = symbols("E ν", real = true)
                k = E / 3(1 - 2ν)
                μ = E / 2(1 + ν)
                λ = E * ν / ((1 + ν) * (1 - 2ν))
                ℂ = simplify(3k * 𝕁 + 2μ * 𝕂)
                @test ℂ == simplify(TensISO{dim}(3k, 2μ))
                𝕊 = simplify(inv(ℂ))
                @test simplify.(KM(𝕊)) == [
                    1/E -ν/E -ν/E 0 0 0
                    -ν/E 1/E -ν/E 0 0 0
                    -ν/E -ν/E 1/E 0 0 0
                    0 0 0 (1+ν)/E 0 0
                    0 0 0 0 (1+ν)/E 0
                    0 0 0 0 0 (1+ν)/E
                ]
                @test simplify(ℂ ⊡ 𝕊) == 𝕀

                n = 𝐞(3)
                Eᵒᵉᵈᵒ = E * (1 - ν) / ((1 + ν) * (1 - 2ν))
                Kref = simplify.([μ 0 0; 0 μ 0; 0 0 Eᵒᵉᵈᵒ])
                @test factor(n ⋅ ℂ ⋅ n) == factor(dotdot(n, ℂ, n)) == Kref
                # Hooke law
                for i ∈ 1:3, j ∈ 1:3
                    @eval $(Symbol("ε$i$j")) = symbols($"ε$i$j", real = true)
                end
                𝛆 = Tens(SymmetricTensor{2,3}((i, j) -> eval(Symbol("ε$i$j"))))
                𝛔 = ℂ ⊡ 𝛆
                @test factor(𝛔) == factor(λ * tr(𝛆) * 𝟏 + 2μ * 𝛆)
                @test factor(simplify(𝛔 ⊡ 𝛆)) == factor(simplify(λ * tr(𝛆)^2 + 2μ * 𝛆 ⊡ 𝛆))

                @test 𝕀 == 𝟏 ⊠ˢ 𝟏
                @test 3𝕁 == 𝟏 ⊗ 𝟏
                @test 𝕀 ⊙ 𝕀 == 6
                @test 𝕁 ⊙ 𝕀 == 𝕁 ⊙ 𝕁 == 1
                @test 𝕂 ⊙ 𝕀 == 𝕂 ⊙ 𝕂 == 5
                @test 𝕂 ⊙ 𝕁 == 𝕁 ⊙ 𝕂 == 0
                @test simplify(ℂ ⊙ 𝕁) == simplify(3k)
                @test simplify(ℂ ⊙ 𝕂) == simplify(10μ)

            end

        end
    end
end


text:
module TestZygote

using Zygote

include("preamble.jl")

const TEST_ALONGS = [
    (True(), False()),
    (False(), True()),
]

function makedata(alongs::NTuple{L,TypedBool}) where {L}
    psize = testdims(L)
    insize = tuple_getindex(psize, alongs)
    outsize = tuple_getindex(psize, invert(alongs))
    M = length(insize)
    N = length(outsize)
    alongs_int = tuple_getindex(ntuple(identity, L), alongs)

    flat = Array{Float64,L}(undef, psize...)
    flat[:] .= 1:length(flat)
    nested = Array{Array{Float64,M},N}(undef, outsize...)

    i = 0
    Base.mapslices(flat, dims = alongs_int) do el
        i += 1
        nested[i] = zeros(Float64, insize...)
        nested[i] .= el
        el
    end

    return (nested = nested, flat = flat)
end

showalongs(alongs) = "($(join(map(a -> a === True() ? ':' : '*', alongs), ", ")))"

@testset "slice" begin
    flat = rand(2, 3)
    let
        gs1 = Zygote.gradient(x -> prod(sum(eachcol(x))), flat)
        gs2 = Zygote.gradient(x -> prod(sum(slice(x, :, *))), flat)
        @test first(gs1) == first(gs2)
    end
    let
        gs1 = Zygote.gradient(x -> prod(sum(eachrow(x))), flat)
        gs2 = Zygote.gradient(x -> prod(sum(slice(x, *, :))), flat)
        @test first(gs1) == first(gs2)
    end
end

@testset "align/flatten" begin
    nested = [rand(2) for _=1:3]
    # Because of the splat in the below hcat, gs1 is a tuple-of-vectors, while gs2 is a
    # vector-of-vectors, so we have to convert for the comparison.
    # See also: https://github.com/FluxML/Zygote.jl/pull/501
    let # alongs = (:, *)
        gs1 = Zygote.gradient(x -> prod(sum(eachcol(hcat(x...)))), nested)
        gs2 = Zygote.gradient(x -> prod(sum(eachcol(align(x, :, *)))), nested)
        gs3 = Zygote.gradient(x -> prod(sum(eachcol(flatten(x)))), nested)
        gs1 = Tuple(Array(el) for el in first(gs1))
        gs2 = Tuple(Array(el) for el in first(gs2))
        gs3 = Tuple(Array(el) for el in first(gs3))
        @test gs1 == gs2
        @test gs1 == gs3
    end
    let # alongs = (*, :)
        gs1 = Zygote.gradient(x -> prod(sum(eachcol(hcat(x...)'))), nested)
        gs2 = Zygote.gradient(x -> prod(sum(eachcol(align(x, *, :)))), nested)
        gs1 = Tuple(Array(el) for el in first(gs1))
        gs2 = Tuple(Array(el) for el in first(gs2))
        @test gs1 == gs2
    end
    let flat = rand(2,3) # flatten(A) = A when A is not a NestedArray
       gs1 = Zygote.gradient(x -> sum(x), flat)
       gs2 = Zygote.gradient(x -> sum(flatten(x)), flat)
       @test gs1 == gs2
    end
end

end # module


text:
using OrdinaryDiffEq
import Unitful: @u_str

# w Unitful
rhsU(x,p,t) = p*x
pU = 0.5u"1/hr"
xU0 = 1.0u"g/L"
tspanU = (0.0u"hr", 5.0u"hr")
probU = ODEProblem(rhsU, xU0, tspanU, pU)
solU = solve(probU, Tsit5())
tU = range(tspanU[1], tspanU[2], length=20)
xU = solU(tU)


text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*
#*-------------------------------------------------------------------
#* Model of costs for tanks
#*--------------------------------------------------------------------
#*	Streams:
#*		* an inlet stream
#*		* an outlet stream
#*
#*	Specify:
#*		* the Inlet stream
#*		* the Outlet flow
#*		* the tank Q
#*
#*	Initial:
#*		* the tank temperature (OutletL.T)
#*		* the tank level (h)
#*		* (NoComps - 1) Outlet compositions
#*----------------------------------------------------------------------
#* Author: N�bia do Carmo Ferreira
#* $Id: tank.mso 210 2007-03-15 12:52:28Z arge $
#*--------------------------------------------------------------------
type tank_cost
	tank_cost()=begin
		new(
			VesselVolume(),
			DanaSwitcher(Dict{Symbol,Any}(
				:Valid=>["Stainless steel 316", "Stainless steel 304", "Stainless steel 347", "Nickel", "Monel", "Inconel", "Zirconium", "Titanium", "Brick_and_rubber", "Brick_and_polyester_lined steel", "Rubber", "Lead_lined steel", "Polyester" ,"Fiberglass_strengthened", "Aluminum", "Copper", "Concrete"],
				:Default=>"Stainless steel 316"
			)),
			fill(DanaReal()),
			length(Dict{Symbol,Any}(
				:Brief=>"Tank height"
			)),
			currency(Dict{Symbol,Any}(
				:Brief=>"Capital Cost"
			)),
			currency(Dict{Symbol,Any}(
				:Brief=>"Basic Cost"
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Cost Factor based on the construction material"
			)),
			volume(Dict{Symbol,Any}(
				:Brief=>"Total Volume"
			)),
			[
				:(V = Across * Height),
				:(Ce = Cb*Fm),
				:(Cb = "US\$"*exp(Cost(1,1) + Cost(1,2)*ln(V/"m^3") + Cost(1,3)*(ln(V/"m^3"))^2)),
				:(Cb = "US\$"*exp(Cost(1,1) + Cost(1,2)*ln(V/"m^3") + Cost(1,3)*(ln(V/"m^3"))^2)),
				:(Cb = "US\$"*exp(Cost(2,1) + Cost(2,2)*ln(V/"m^3") + Cost(2,3)*(ln(V/"m^3"))^2)),
				:(Cb = "US\$"*exp(Cost(2,1) + Cost(2,2)*ln(V/"m^3") + Cost(2,3)*(ln(V/"m^3"))^2)),
				:(Fm = Cost(3,1) + Cost(3,2) + Cost(3,3)),
			],
			[
				"Total Volume","Capital Cost","Basic Cost","tank bought soon","the tank manufactured in the plant","","Cost Factor based on the construction material",
			],
			[:Material,:Cost,:Height,],
			[:Ce,:Cb,:Fm,:V,]
		)
	end
	_base_1::VesselVolume
	Material::DanaSwitcher
	Cost::Array{DanaReal}
	Height::length
	Ce::currency
	Cb::currency
	Fm::positive
	V::volume
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export tank_cost
function setEquationFlow(in::tank_cost)
	addEquation(1)
	addEquation(2)
	if V < 5 * "m^3"
		addEquation(3)
		#verificar

	else
		if 5 * "m^3" > V && V < 80 * "m^3"
			addEquation(4)
		else
			if 80 * "m^3" > V && V < 45000 * "m^3"
				addEquation(5)
			else
				addEquation(6)
				#verificar

			end
		end
	end
	addEquation(7)
end
function atributes(in::tank_cost,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=true
	fields[:Icon]="icon/Tank"
	drive!(fields,_)
	return fields
end
tank_cost(_::Dict{Symbol,Any})=begin
	newModel=tank_cost()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
# Monte-Carlo Applications

# Type 1
# path like integration
# ∫₀¹ e
# todo

# Type 2
# Volume-like integration

# E.g. evaluate π

# 1. generate n points within a unit □
#    enclosed by xlim<-[0 1]
#    and         ylim<-[0 1]

n = 1000
x = rand(n)
y = rand(n)

# 2. Judge # of pts in the unit ◷
s = 0
for i in 1:n# I hate for loop , there's a simpler version below
    if (x[i]^2 + y[i]^2)<=1
        global s+=1
    end
end
s

# Area of ▧ / Area of ◍  = 4/π
π_est = 4s/n
π_est

# Vectorize? replace 2 with
#even faster?
@time size(zeros(n)[x.^2+y.^2 .<=1])
#0.000070 seconds (27 allocations: 42.813 KiB)
#Compared with

s=0
@time for i in 1:n
    if (x[i]^2 + y[i]^2)<=1
        global s+=1
    end
end
#0.000220 seconds (7.75 k allocations: 136.734 KiB)

π_est = 4s/n
π_est

plot(x,y;seriestype="scatter")
plot!(0:.01:1,sqrt.(1 .-(0:.01:1).^2))