text:
"""
For an initial guess from the index of a Fold bifurcation point located in ContResult.bifpoint, returns a point which will be refined using `newtonFold`.
"""
function FoldPoint(br::ContResult, index::Int64)
	@assert br.foldpoint[index].type == :fold "The provided index does not refer to a Fold point"
	bifpoint = br.foldpoint[index]
	return BorderedArray(_copy(bifpoint.x), bifpoint.param)
end

####################################################################################################
# Method using Minimally Augmented formulation
@with_kw struct FoldProblemMinimallyAugmented{TF, TJ, TJa, Td2f, Tl <: Lens, vectype, S <: AbstractLinearSolver, Sa <: AbstractLinearSolver, Sbd <: AbstractBorderedLinearSolver}
	F::TF								# Function F(x, p) = 0
	J::TJ								# Jacobian of F wrt x
	Jadjoint::TJa						# Adjoint of the Jacobian of F
	d2F::Td2f = nothing					# Hessian of F
	lens::Tl							# parameter axis for the Fold point
	a::vectype							# close to null vector of J^T
	b::vectype							# close to null vector of J
	linsolver::S						# linear solver
	linsolverAdjoint::Sa = linsolver	# linear solver for the jacobian adjoint
	linbdsolver::Sbd					# bordered linear solver
end

hasHessian(pb::FoldProblemMinimallyAugmented{TF, TJ, TJa, Td2f, Tl, vectype, S, Sa, Sbd}) where {TF, TJ, TJa, Td2f, Tp, Tl, vectype, S, Sa, Sbd} = Td2f != Nothing

hasAdjoint(pb::FoldProblemMinimallyAugmented{TF, TJ, TJa, Td2f, Tl, vectype, S, Sa, Sbd}) where {TF, TJ, TJa, Td2f, Tl, vectype, S, Sa, Sbd} = TJa != Nothing

FoldProblemMinimallyAugmented(F, J, Ja, d2F, lens::Lens, a, b, linsolve) = FoldProblemMinimallyAugmented(F, J, Ja, d2F, lens, a, b, linsolve, linsolve, BorderingBLS(linsolve))

FoldProblemMinimallyAugmented(F, J, Ja, lens::Lens, a, b, linsolve) = FoldProblemMinimallyAugmented(F, J, Ja, nothing, lens, a, b, linsolve)

function (fp::FoldProblemMinimallyAugmented)(x::vectype, p::T, par) where {vectype, T}
	# These are the equations of the minimally augmented (MA) formulation of the Fold bifurcation point
	# input:
	# - x guess for the point at which the jacobian is singular
	# - p guess for the parameter at which the jacobian is singular
	# The jacobian of the MA problem is solved with a bordering method
	a = fp.a
	b = fp.b

	# we solve Jv + a σ1 = 0 with <b, v> = n
	# the solution is v = -σ1 J\a with σ1 = -n/<b, J^{-1}a>
	n = T(1)
	v, _, _ = fp.linsolver(fp.J(x, set(par, fp.lens, p)), a)
	bv = dot(b, v)
	bv == T(0) && @error "Error when using Minimally Augmented formulation for Fold bifurcation. The dot product should non-zero."

	σ1 = -n / bv
	return fp.F(x, set(par, fp.lens, p)), σ1
end

# this function is for the functional
function (foldpb::FoldProblemMinimallyAugmented)(x::BorderedArray, params)
	res = foldpb(x.u, x.p, params)
	return BorderedArray(res[1], res[2])
end

# Struct to invert the jacobian of the fold MA problem.
struct FoldLinearSolverMinAug <: AbstractLinearSolver; end

function foldMALinearSolver(x, p::T, pbMA::FoldProblemMinimallyAugmented, par,
							rhsu, rhsp;
							debug_::Bool = false) where T
	# The jacobian should be passed as a tuple as Jac_fold_MA(u0, pb::FoldProblemMinimallyAugmented) = (return (u0, pb, d2F::Bool))
	# The Jacobian J of the vector field is expressed at (x, p)
	# We solve here Jfold⋅res = rhs := [rhsu, rhsp]
	# The Jacobian expression of the Fold problem is Jfold = [J dpF ; σx σp] where σx := ∂_xσ
	# We recall the expression of σx = -< w, d2F(x,p)[v, x2]> where (w, σ2) is solution of J'w + b σ2 = 0 with <a, w> = n
	############### Extraction of function names #################
	# N = length(du) - 1

	F = pbMA.F
	J = pbMA.J

	d2F = pbMA.d2F
	a = pbMA.a
	b = pbMA.b

	# parameter axis
	lens = pbMA.lens

	# we define the following jacobian. It is used at least 3 times below. This avoid doing 3 times the possibly costly building of J(x, p)
	J_at_xp = J(x, set(par, lens, p))

	# we do the following in order to avoid computing J_at_xp twice in case pbMA.Jadjoint is not provided
	if hasAdjoint(pbMA)
		JAd_at_xp = pbMA.Jadjoint(x, set(par, lens, p))
	else
		JAd_at_xp = transpose(J_at_xp)
	end

	n = T(1)
	# we solve Jv + a σ1 = 0 with <b, v> = n
	# the solution is v = -σ1 J\a with σ1 = -n/<b, J\a>
	v = pbMA.linsolver(J_at_xp, a)[1]
	σ1 = -n / dot(b, v)
	rmul!(v, -σ1)

	# we solve J'w + b σ2 = 0 with <a, w> = n
	# the solution is w = -σ2 J'\b with σ2 = -n/<a, J'\b>
	w = pbMA.linsolverAdjoint(JAd_at_xp, b)[1]
	σ2 = -n / dot(a, w)
	rmul!(w, -σ2)

	δ = T(1e-8)
	ϵ1, ϵ2, ϵ3 = T(δ), T(δ), T(δ)
	################### computation of σx σp ####################
	dpF = minus(F(x, set(par, lens, p + ϵ1)), F(x, set(par, lens, p - ϵ1))); rmul!(dpF, T(1) / T(2ϵ1))
	dJvdp = minus(apply(J(x, set(par, lens, p + ϵ3)), v), apply(J(x, set(par, lens, p - ϵ3)), v)); rmul!(dJvdp, T(1) / T(2ϵ3))
	σp = -dot(w, dJvdp) / n

	if hasHessian(pbMA) == false
		# We invert the jacobian of the Fold problem when the Hessian of x -> F(x, p) is not known analytically. We thus rely on finite differences which can be slow for large dimensions
		prod(size(x)) > 1e4 && @warn "You might want to pass the Hessian, you are using finite differences with $(prod(size(x))) unknowns to compute a gradient"
		# we first need to compute the value of ∂_xσ written σx
		σx = zero(x)
		e  = zero(x)

		# this is the part which does not work if x is not AbstractArray. We use CartesianIndices to support AbstractArray as a type for the solution we are looking for
		for ii in CartesianIndices(x)
			e[ii] = T(1)
			# d2Fve := d2F(x,p)[v,e]
			d2Fve = (apply(J(x + ϵ2 * e, set(par, lens, p)), v) - apply(J(x - ϵ2 * e, set(par, lens, p)), v)) / T(2ϵ2)
			σx[ii] = -dot(w, d2Fve) / n
			e[ii] = T(0)
		end

		########## Resolution of the bordered linear system ########
		dX, dsig, flag, it = pbMA.linbdsolver(J_at_xp, dpF, σx, σp, rhsu, rhsp)

	else
		# We invert the jacobian of the Fold problem when the Hessian of x -> F(x, p) is known analytically. Much faster than the previous case

		# we solve it here instead of calling linearBorderedSolver because this removes the need to pass the linear form associated to σx
		x1, x2, _, it = pbMA.linsolver(J_at_xp, rhsu, dpF)

		d2Fv = d2F(x, set(par, lens, p), x1, v)
		σx1 = -dot(w, d2Fv ) / n

		copyto!(d2Fv, d2F(x, set(par, lens, p), x2, v))
		σx2 = -dot(w, d2Fv ) / n

		dsig = (rhsp - σx1) / (σp - σx2)

		# dX = x1 .- dsig .* x2
		dX = copyto!(similar(x1), x1); axpy!(-dsig, x2, dX)
	end

	if debug_
		return dX, dsig, true, sum(it), 0., [J(x, set(par, lens, p)) dpF ; σx' σp]
	else
		return dX, dsig, true, sum(it), 0.
	end
end

function (foldl::FoldLinearSolverMinAug)(Jfold, du::BorderedArray{vectype, T}, debug_::Bool = false) where {vectype, T}
	out =  foldMALinearSolver((Jfold.x).u,
				 (Jfold.x).p,
				 Jfold.fldpb,
				 Jfold.param,
				 du.u, du.p;
				 debug_ = debug_)

	if debug_ == false
		return BorderedArray(out[1], out[2]), out[3], out[4], out[5]
	else
		return BorderedArray(out[1], out[2]), out[3], out[4], out[5], out[6]
	end
end

################################################################################################### Newton / Continuation functions
"""
	newtonFold(F, J, foldpointguess, par, lens::Lens, eigenvec, options::NewtonPar; Jt = nothing, d2F = nothing, normN = norm, kwargs...)

This function turns an initial guess for a Fold point into a solution to the Fold problem based on a Minimally Augmented formulation. The arguments are as follows
- `F   = (x, p) -> F(x, p)` where `p` is a set of parameters.
- `dF  = (x, p) -> d_xF(x, p)` associated jacobian
- `foldpointguess` initial guess (x_0, p_0) for the Fold point. It should be a `BorderedArray` as returned by the function `FoldPoint`
- `par` parameters used for the vector field
- `lens` parameter axis used to locate the Fold point.
- `eigenvec` guess for the 0 eigenvector
- `options::NewtonPar` options for the Newton-Krylov algorithm, see [`NewtonPar`](@ref).

# Optional arguments:
- `Jt = (x, p) -> transpose(d_xF(x, p))` jacobian adjoint, it should be implemented in an efficient manner. For matrix-free methods, `transpose` is not readily available and the user must provide a dedicated method. In the case of sparse based jacobian, `Jt` should not be passed as it is computed internally more efficiently, i.e. it avoid recomputing the jacobian as it would be if you pass `Jt = (x, p) -> transpose(dF(x, p))`
- `d2F = (x, p, v1, v2) ->  d2F(x, p, v1, v2)` a bilinear operator representing the hessian of `F`. It has to provide an expression for `d2F(x,p)[v1,v2]`.
- `normN = norm`
- `kwargs` keywords arguments to be passed to the regular Newton-Krylov solver

# Simplified call
Simplified call to refine an initial guess for a Fold point. More precisely, the call is as follows

	newtonFold(F, J, br::ContResult, ind_fold::Int64, par, lens::Lens; Jt = nothing, d2F = nothing, options = br.contparams.newtonOptions, kwargs...)

where the optional argument `Jt` is the jacobian transpose and the Hessian is `d2F`. The parameters / options are as usual except that you have to pass the branch `br` from the result of a call to `continuation` with detection of bifurcations enabled and `index` is the index of bifurcation point in `br` you want to refine. You can pass newton parameters different from the ones stored in `br` by using the argument `options`.

!!! tip "Jacobian tranpose"
    The adjoint of the jacobian `J` is computed internally when `Jt = nothing` by using `tranpose(J)` which works fine when `J` is an `AbstractArray`. In this case, do not pass the jacobian adjoint like `Jt = (x, p) -> transpose(d_xF(x, p))` otherwise the jacobian will be computed twice!

!!! warning "Hessian"
    The hessian of `F`, when `d2F` is not passed, is computed with Finite differences. This can be slow for many variables, e.g. ~1e6
"""
function newtonFold(F, J, foldpointguess, par, lens::Lens, eigenvec, options::NewtonPar; Jt = nothing, d2F = nothing, normN = norm, kwargs...) where {T, vectype}
	foldproblem = FoldProblemMinimallyAugmented(
		F, J, Jt, d2F, lens,
		_copy(eigenvec), #copy(eigenvec),
		_copy(eigenvec), #copy(eigenvec),
		options.linsolver)

	# Jacobian for the Fold problem
	Jac_fold_MA = (x, param) -> (x = x, param = param, fldpb = foldproblem)

	opt_fold = @set options.linsolver = FoldLinearSolverMinAug()

	# solve the Fold equations
	return newton(foldproblem, Jac_fold_MA, foldpointguess, par, opt_fold; normN = normN, kwargs...)
end

function newtonFold(F, J, br::ContResult, ind_fold::Int64, par, lens::Lens; Jt = nothing, d2F = nothing, options = br.contparams.newtonOptions, kwargs...)
	foldpointguess = FoldPoint(br, ind_fold)
	bifpt = br.foldpoint[ind_fold]
	eigenvec = bifpt.tau

	# solve the Fold equations
	return newtonFold(F, J, foldpointguess, par, lens, eigenvec, options; Jt = Jt, d2F = d2F, kwargs...)
end

"""
$(SIGNATURES)

Codim 2 continuation of Fold points. This function turns an initial guess for a Fold point into a curve of Fold points based on a Minimally Augmented formulation. The arguments are as follows
- `F = (x, p) ->	F(x, p)` where `p` is a set of parameters
- `J = (x, p) -> d_xF(x, p)` associated jacobian
- `foldpointguess` initial guess (x_0, p1_0) for the Fold point. It should be a `BorderedArray` as returned by the function `FoldPoint`
- `par` set of parameters
- `lens1` parameter axis for parameter 1
- `lens2` parameter axis for parameter 2
- `eigenvec` guess for the 0 eigenvector at p1_0
- `options_cont` arguments to be passed to the regular [`continuation`](@ref)

# Optional arguments:

- `Jt = (x, p) -> transpose(d_xF(x, p))` associated jacobian transpose
- `d2F = p -> ((x, p, v1, v2) -> d2F(x, p, v1, v2))` this is the hessian of `F` computed at `(x, p)` and evaluated at `(v1, v2)`.
- `kwargs` keywords arguments to be passed to the regular [`continuation`](@ref)

# Simplified call
The call is as follows

	continuationFold(F, J, br::ContResult, ind_fold::Int64, par, lens1::Lens, lens2::Lens, options_cont::ContinuationPar ; Jt = nothing, d2F = nothing, kwargs...)

where the parameters are as above except that you have to pass the branch `br` from the result of a call to `continuation` with detection of bifurcations enabled and `index` is the index of Fold point in `br` you want to continue.

!!! tip "Jacobian tranpose"
    The adjoint of the jacobian `J` is computed internally when `Jt = nothing` by using `tranpose(J)` which works fine when `J` is an `AbstractArray`. In this case, do not pass the jacobian adjoint like `Jt = (x, p) -> transpose(d_xF(x, p))` otherwise the jacobian would be computed twice!

!!! warning "Hessian"
    The hessian of `F`, when `d2F` is not passed, is computed with Finite differences. This can be slow for many variables, e.g. ~1e6
"""
function continuationFold(F, J, foldpointguess::BorderedArray{vectype, T}, par, lens1::Lens, lens2::Lens, eigenvec, options_cont::ContinuationPar ; Jt = nothing, d2F = nothing, kwargs...) where {T,vectype}

	# options for the Newton Solver inheritated from the ones the user provided
	options_newton = options_cont.newtonOptions

	foldPb = FoldProblemMinimallyAugmented(
			F, J, Jt, d2F,
			lens1,
			_copy(eigenvec), #copy(eigenvec),
			_copy(eigenvec), #copy(eigenvec),
			options_newton.linsolver)

	# Jacobian for the Fold problem
	Jac_fold_MA = (x, param) -> (x = x, param = param, fldpb = foldPb)

	opt_fold_cont = @set options_cont.newtonOptions.linsolver = FoldLinearSolverMinAug()

	# solve the Fold equations
	branch, u, tau = continuation(
		foldPb, Jac_fold_MA,
		foldpointguess, par, lens2,
		opt_fold_cont,
		printSolution = (u, p) -> u.p,
		plotSolution = (x, p; kwargs...) -> (xlabel!("p2", subplot=1); ylabel!("p1", subplot=1)  ); kwargs...)
	return setproperties(branch; type = :FoldCodim2, functional = foldPb), u, tau
end

function continuationFold(F, J, br::ContResult, ind_fold::Int64, par, lens1::Lens, lens2::Lens, options_cont::ContinuationPar ; Jt = nothing, d2F = nothing, kwargs...)
	foldpointguess = FoldPoint(br, ind_fold)
	bifpt = br.foldpoint[ind_fold]
	eigenvec = bifpt.tau
	return continuationFold(F, J, foldpointguess, par, lens1, lens2, eigenvec, options_cont ; Jt = Jt, d2F = d2F, kwargs...)
end


text:
include("./single_bh_superradiance_param.jl")

s = [ p2 ]

for x in s
    solve_wave_equation_2D(x)
end

ext: jl


text:
module SBMs

abstract SBM

export SBM
export loglikelihood, computeparams
export considermoves, enactmove
export getindices, getgroups
export setgroups

"""

    loglikelihood(sbm::SBM, g=nothing, ps=nothing)

  Computes the loglikelihood for a model at the cluster state `g`, or given the
model parameters `ps`.

  If a `ps` argument is passed, the `g` variable is ignored.
If neither a `g` nor a `ps` is passed, g is taken from the internal SBM state,
and the parameters are computed from that state.
"""
function loglikelihood(sbm::SBM, g=nothing, ps=nothing)
  error("loglikelihood not implemented for type $(typeof(sbm))")
end

"""

    computeparams(sbm::SBM, g=nothing)

  Computes the model parameters for a model at the cluster state `g`.

  If no `g` variable is passed, it's taken from the internal model state.
"""
function computeparams(sbm::SBM, g=nothing, ps=nothing)
  error("computeparams not implemented for type $(typeof(sbm))")
end

"""

    updateparams!(sbm::SBM, ps, old_g, g=nothing)

  Sometimes, it can be quicker to update parameters given the group changes
instead or recomputing them. This performs an update step, or defaults to
recomputation.

  If no `g` variable is passed, it's taken from the internal model state.
"""
function updateparams!(sbm::SBM, ps, old_g, g=nothing)
  computeparams(sbm,g)
end

"""

    considermoves(sbm::sbm, i, g=nothing, ps=nothing; forcemove=true)

  Computes the loglikelihood for all possible moves of node `i` starting
from cluster state `g` or parameters `ps`.

  If a `ps` argument is passed, the `g` variable is ignored.
If neither a `g` nor a `ps` is passed, g is taken from the internal SBM state,
and the parameters are computed from that state.

  If `forcemove`, then only the cluster assignments not equal to the current
setting are considered.
"""
function considermoves(sbm::SBM, i, g=nothing, ps=nothing; forcemove=true)
  error("considermove not implemented for type $(typeof(sbm))")
end

"""

    bestmove(sbm::SBM, i, g=nothing, ps=nothing; forcemove=true)

Selects the move which maximizes the loglikelihood for a given model.

  If `forcemove`, then only the cluster assignments not equal to the current
setting are considered, and the move selected could be one which least decreases
the ll.
"""
function bestmove(sbm::SBM, i, g=nothing, ps=nothing; forcemove=true)
  moves, logliks = considermoves(sbm, i, g, ps; forcemove=forcemove)
  maxll, maxi = findmax(logliks)

  moves[maxi], maxll
end

"""

    enactmove(sbm::SBM, i, new_group)

Changes the internal cluster state of `sbm` by a single move.

See `setgroups` for batch editing
"""
enactmove(sbm::SBM, i, new_group) = sbm.g[i] = new_group

getindices(sbm::SBM) = 1:size(sbm.G,1)
getgroups(sbm::SBM) = copy(sbm.g)
getgroups(sbm::SBM, i) = sbm.g[i]

"""
    setgroups(sbm::SBM, g)

  Sets the cluster state to `g`
"""
function setgroups(sbm::SBM, g)

  @assert typeof(g) == typeof(sbm.g)
  @assert length(g) == length(sbm.g)

  sbm.g = copy(g)

end

#==========================================
General SBM Utility Functions
==========================================#

"""

    fill_G!(G::SparseMatrixCSC, g, dists, dir)

  Fills the graph matrix with edges depending upon the group
distributions within dists. Sampling from the distribution at
index `[i,j]` should instantiate an edge/multiedge between
a node in group `i` and a node in group `j`.
"""
function fill_G!(G::SparseMatrixCSC, g, dists, dir)

  @assert size(G,1) == size(G,2) "G must be a square matrix"
  num_groups = length(unique(g))
  @assert size(dists) == (num_groups,num_groups)

  n = size(G,1)

  for i in 1:n, j in 1:n

    gi,gj = g[i], g[j]

    dist = dists[gi,gj]

    v = rand(dist)
    if v == 0 continue end

    G[i,j] = v
    if !dir  G[j,i] = v  end
  end

end

"""

    fill_degcorr_G!(G::SparseMatrixCSC, g, dists, dir)

  Fills the graph matrix with edges depending upon the group
distributions within dists. Sampling from the distribution at
index `[i,j]` should instantiate an edge/multiedge between
a node in group `i` and a node in group `j`.
"""
function fill_degcorr_G!(G::SparseMatrixCSC, g, dists, dir)
  0#stub
end

function randomize_g(g, group_types)

  rand_g = zeros(Int,length(g))

  group_to_type = Dict();
  for (i,t) in enumerate(group_types)
    for group in t
      group_to_type[group] = i
    end
  end

  for i in eachindex(g)
    rand_g[i] = rand(group_types[group_to_type[g[i]]])
  end

  rand_g
end

randomize_g!(sbm::SBM) = setgroups(sbm, randomize_g(sbm.g, sbm.t))

"""
Finds the valid cluster pairs to evaluate for the likelihood fn
"""
function assemble_pairs(gtypes, directed)

  @assert length(gtypes) in [1,2] "only full or bipartite graphs supported"

  pairs = Tuple{Int,Int}[] #init
  if length(gtypes) == 2 && directed #directed bipartite graph
    pairs12 = [(r,s) for r in gtypes[1], s in gtypes[2]][:]
    pairs21 = [(r,s) for r in gtypes[2], s in gtypes[1]][:]
    pairs = vcat(pairs12,pairs21)

  elseif length(gtypes) == 2 #undirected bipartite
    pairs = [(r,s) for r in gtypes[1], s in gtypes[2]][:]

  elseif directed #directed full graph
    pairs = [(r,s) for r in gtypes[1], s in gtypes[1]][:]

  else #undirected full graph
    gs = gtypes[1]; ng = length(gs)
    for i in 1:ng, j in i+1:ng
      push!(pairs,(gs[i],gs[j]))
    end
  end

  pairs
end

end #module SBMs


text:
@testset "RainYday.jl" begin
    data_surv = CSO.dataset("cso_raw")
    data_trudeau = CSO.dataset("Montreal-Trudeau_15min_precipitations")
    data = innerjoin(data_surv, data_trudeau, on=:Date)

    notInc = 0
    if data[1, :d24h] == 1
        notInc = 1
    end
    for i in 2:nrow(data)
        if data[i-1, :d24h] .> 0 && year(data[i,:].Date) .!= year(data[i-1,:].Date)
            notInc = notInc + 1
        end
    end

    transformedDF = addColRainYday!(data, "d24h", "RainYday", 5)
    nb1 = sum(transformedDF.RainYday)

    #New column contains only 0's and 1'
    @test unique(transformedDF.RainYday) in [0, 1, [0,1]]

    #The number of 1 equals the number of positive precipitations but doesn't include the first day of the year
    @test count(transformedDF[:, :d24h] .> 0.0) == (nb1 + notInc)

    #Ancient DF has been changed
    @test nrow(data) == nrow(transformedDF)
end

max_stars_repo_name: AlexLsn/CSO

max_issues_repo_name: AlexLsn/CSO

max_forks_repo_name: AlexLsn/CSO


text:
struct Sphere{T} <: AbstractRotationalSymmetricSurface{T}
    curvature::T
end

function Sphere(curvature, aperture, n, t, id = nothing)
    OpticalSurface(
        Sphere(curvature),
        aperture,
        n,
        t,
        id,
    )
end

Plano(aperture, n, t, id = nothing) = Sphere(0.0, aperture, n, t, id)

function sag(radius, s::Sphere)
    radius2 = radius^2
    return (s.curvature * radius2 / (1 + sqrt(1 - (s.curvature^2 * radius2))))
end

# Raytrace für sphärische Flächen
function transfer_to_intersection!(ray, t, s::Sphere)
    cv = s.curvature
    e = t * ray.cz - (ray.x * ray.cx + ray.y * ray.cy + ray.z * ray.cz)
    M1z = ray.z + e * ray.cz - t
    M1squared = ray.x^2 + ray.y^2 + ray.z^2 - e^2 + t^2 - 2 * t * ray.z
    E1Arg = ray.cz^2 - cv * (cv * M1squared - 2 * M1z)
    if E1Arg < 0
        return ray_miss_error()
    end
    E1 = sqrt(E1Arg)
    L = e + (cv * M1squared - 2 * M1z) / (ray.cz + E1)
    z1 = ray.z + L * ray.cz - t
    y1 = ray.y + L * ray.cy
    x1 = ray.x + L * ray.cx

    ray.x = x1
    ray.y = y1
    ray.z = z1
    ray._e1 = E1
    return ray
end

function refract!(ray, m, s::Sphere, m1, λ)
    E1 = ray._e1
    if m == m1
        X1 = ray.cx
        Y1 = ray.cy
        Z1 = ray.cz
    else
        n0 = refractive_index(m, λ)
        n1 = refractive_index(m1, λ)
        EprimeArg = 1 - (n0 / n1)^2 * (1 - E1^2)
        if EprimeArg < 0
            return total_internal_reflection_error()
        end
        Eprime = sqrt(EprimeArg)
        g1 = Eprime - n0 / n1 * E1
        Z1 = n0 / n1 * ray.cz - g1 * s.curvature * ray.z + g1
        Y1 = n0 / n1 * ray.cy - g1 * s.curvature * ray.y
        X1 = n0 / n1 * ray.cx - g1 * s.curvature * ray.x
    end

    ray.cx = X1
    ray.cy = Y1
    ray.cz = Z1
    return ray
end


text:
using ImageFiltering

@testset "SSIM" begin
    @info "test: SSIM"

    iqi_a = SSIM()
    iqi_b = SSIM(KernelFactors.gaussian(1.5, 11))
    iqi_c = SSIM(KernelFactors.gaussian(1.5, 11), (1.0, 1.0, 1.0))
    iqi_d = SSIM(KernelFactors.gaussian(1.5, 11), (1.0, 1.0, 1.0); crop=false)
    @test iqi_a == iqi_b == iqi_c == iqi_d
    @test SSIM(; crop=true) == SSIM(KernelFactors.gaussian(1.5, 11), (1.0, 1.0, 1.0); crop=true)

    iqi = SSIM()
    sz_img_3 = (3, 3, 3)

    # numerical test
    img1 = testimage("cameraman")
    img2 = testimage("lena_gray_512")
    @test assess_ssim(img1, img2; crop=false) ≈ 0.3595 atol=1e-4 # MATLAB built-in ssim result
    @test assess_ssim(img1, img2; crop=true) ≈ 0.3526 atol=1e-4 # MATLAB original ssim result & tensorflow implementation

    # this calls the general implementation
    iqi_δ1 = SSIM(KernelFactors.gaussian(1.5, 11), (1.0+1e-5, 1.0, 1.0))
    @test assess(iqi_δ1, img1, img2) ≈ assess(SSIM(), img1, img2) atol = 1e-4
    # this calls the general implementation with max.(s, 0)
    iqi_δ2 = SSIM(KernelFactors.gaussian(1.5, 11), (1.0, 1.0, 1.0-1e-5))
    @test assess(iqi_δ2, img1, img2) ≈ assess(SSIM(), img1, img2) atol = 1e-2

    # non-standard powers
    iqi_γ = SSIM(KernelFactors.gaussian(1.5, 11), (0.5, 0.5, 0.5))
    @test iqi_γ(img1, img2) ≈ 0.5261 atol=1e-4
    # non standard parameters, result may differ from other implementations

    # images should have the same size, but their axes can differs
    @test assess_ssim(img1, OffsetArray(img2, -2, -2)) == assess_ssim(img1, img2)
    @test_throws ArgumentError assess_ssim(img1, restrict(img2))

    # Gray image
    type_list = generate_test_types([Bool, Float32, N0f8], [Gray])
    A = [1.0 1.0 1.0; 1.0 1.0 1.0; 0.0 0.0 0.0]
    B = [1.0 1.0 1.0; 0.0 0.0 0.0; 1.0 1.0 1.0]
    for T in type_list
        test_ndarray(iqi, sz_img_3, T)

        a = A .|> T
        b = B .|> T

        @test assess_ssim(a, b) == assess_ssim(a, b; crop = false) == assess(SSIM(), a, b) == SSIM()(a, b)
        @test assess_ssim(a, b; crop = true) == assess(SSIM(crop = true), a, b) == SSIM(crop = true)(a, b)
        @test assess_ssim(a, a) ≈ 1.0

        # FIXME: the result of Bool type is not strictly equal to others
        eltype(T) <: Bool && continue
        test_numeric(iqi, a, b, T; filename="references/SSIM_2d_Gray")
        test_numeric(iqi, channelview(a), channelview(b), T; filename="references/SSIM_2d_Gray")
    end
    test_cross_type(iqi, A, B, type_list)

    # RGB image
    img1 = testimage("mandril_color")
    img2 = testimage("lena_color_512")
    # this differs from MATLAB built-in result as our implementation don't slide the window in the channel dimension
    @test assess_ssim(img1, img2; crop=false) ≈ 0.11069226443828077 atol=1e-4 # MATLAB built-in: 0.0664
    @test assess_ssim(img1, img2; crop=true) ≈ 0.10967952100784095 atol=1e-4 # the original implementation: 0.1047

    # gray images are promoted to RGB images before calculation
    img3 = testimage("lena_gray_512")
    @test assess_ssim(img2, img3) == assess_ssim(img2, RGB.(img3))
    @test assess_ssim(img2, img3) ≠ assess_ssim(Gray.(img2), img3)

    type_list = generate_test_types([Float32, N0f8], [RGB])
    A = [RGB(0.0, 0.0, 0.0) RGB(0.0, 1.0, 0.0) RGB(0.0, 1.0, 1.0)
        RGB(0.0, 0.0, 1.0) RGB(1.0, 0.0, 0.0) RGB(1.0, 1.0, 0.0)
        RGB(1.0, 1.0, 1.0) RGB(1.0, 0.0, 1.0) RGB(0.0, 0.0, 0.0)]
    B = [RGB(0.0, 0.0, 0.0) RGB(0.0, 0.0, 1.0) RGB(1.0, 1.0, 1.0)
        RGB(0.0, 1.0, 0.0) RGB(1.0, 0.0, 0.0) RGB(1.0, 0.0, 1.0)
        RGB(0.0, 1.0, 1.0) RGB(1.0, 1.0, 0.0) RGB(0.0, 0.0, 0.0)]
    for T in type_list
        test_ndarray(iqi, sz_img_3, T)

        a = A .|> T
        b = B .|> T

        @test assess_ssim(a, b; crop = false) == assess_ssim(a, b) == assess(SSIM(), a, b) == SSIM()(a, b)
        @test assess_ssim(a, b; crop = true) == assess(SSIM(crop = true), a, b) == SSIM(crop = true)(a, b)
        @test assess_ssim(a, b) ≠ assess_ssim(channelview(a), channelview(b))
        @test assess_ssim(a, a) ≈ 1.0

        # RGB is treated as 3d gray image
        test_numeric(iqi, a, b, T; filename="references/SSIM_2d_RGB")
    end
    type_list = generate_test_types([Float32, N0f8], [RGB, BGR])
    test_cross_type(iqi, A, B, type_list)

    # general Color3 images
    type_list = generate_test_types([Float32], [Lab, HSV])
    A = [RGB(0.0, 0.0, 0.0) RGB(0.0, 1.0, 0.0) RGB(0.0, 1.0, 1.0)
        RGB(0.0, 0.0, 1.0) RGB(1.0, 0.0, 0.0) RGB(1.0, 1.0, 0.0)
        RGB(1.0, 1.0, 1.0) RGB(1.0, 0.0, 1.0) RGB(0.0, 0.0, 0.0)]
    B = [RGB(0.0, 0.0, 0.0) RGB(0.0, 0.0, 1.0) RGB(1.0, 1.0, 1.0)
        RGB(0.0, 1.0, 0.0) RGB(1.0, 0.0, 0.0) RGB(1.0, 0.0, 1.0)
        RGB(0.0, 1.0, 1.0) RGB(1.0, 1.0, 0.0) RGB(0.0, 0.0, 0.0)]
    for T in type_list
        a = A .|> T
        b = B .|> T

        @test_nowarn assess_ssim(A, b), assess_ssim(a, B)

        @test assess_ssim(a, b; crop = false) == assess_ssim(a, b) == assess(SSIM(), a, b) == SSIM()(a, b)
        @test assess_ssim(a, b; crop = true) == assess(SSIM(crop = true), a, b) == SSIM(crop = true)(a, b)
        @test assess_ssim(a, a) ≈ 1.0

        # conversion to RGB first differs from no conversion
        @test assess_ssim(a, b) ≠ assess_ssim(channelview(a), channelview(b))
    end
    @test assess_ssim(A, B) ≈ assess_ssim(Lab.(A), B) atol=1e-4
end


text:
using MyQuadtree
#using Base.Test

## Quadtree used
quad = Quadtree(0,0,160,160)

## Shapes definition
shape2      = Shape([1 1; 1 79; 79 79; 79 1])
shape3333   = Shape([1 159; 1 151; 9 151; 9 159])
shape33333  = Shape([2 158; 2 157; 8 157; 8 158])
shape31422  = Shape([61 101; 61 104; 64 101; 64 104])
shape3143   = Shape([61 111; 61 119; 69 111; 69 119])
shape3      = Shape([10 90; 70 90; 70 150; 10 150])
shape1      = Shape([150 10; 90 10; 90 70; 150 70])
shape111    = Shape([141 10; 141 15; 150 15; 150 10])
shape12     = Shape([81 1; 119 1; 119 5; 81 5])
shape_big   = Shape([1 1; 1 159; 159 159; 159 1])
shape_      = Shape([71 79; 79 79; 79 71; 71 71])

## Point used to find if there are colisions
point_test      = Shape([62 103])
line_test       = Shape([23 85; 93 115])
triangle_test   = Shape([40 140; 70 140; 60 70])

## Insert the shapes onto the quadtree
#insert!(quad, shape2)     ## quad.nw.objects          [1 1;    1 79;   79 79;  79 1]
#insert!(quad, shape31422) ##                          [61 101; 61 104; 64 101; 64 104]
#insert!(quad, shape3333)  ## quad.sw.sw.sw.sw.objects [1 159;  1 151;  9 159;  9 151]
#insert!(quad, shape3143)  ## quad.sw.ne.se.sw.objects [61 111; 61 119; 69 111; 69 119]
#insertOn!(quad, shape3)     ## quad.sw.objects          [10 90;  70 90;  10 150; 70 150]
#insert!(quad, shape1)
#insert!(quad, shape111)
#insert!(quad, shape12)
#insert!(quad, triangle_test)
#insert!(quad, shape_big)

## Insert random objects inside the tree and retrieve them
const MAX_SIZE = 5
k = 5
number_of_objects = 10^k
shapes_test = Array(Shape, number_of_objects)

@time for i = 1:number_of_objects
    size2 = rand(1:MAX_SIZE)

    x1 = rand(1:(159-size2))
    x2 = x1+size2
    y1 = rand(1:(159-size2))
    y2 = y1+size2
    a = [x1 y1;
         x1 y2;
         x2 y2;
         x2 y1
        ]
    s = Shape(a)
    shapes_test[i] = s
end

@time for i = 1:number_of_objects
    insertOn!(quad, shapes_test[i])
end

@time o = retrieve!(quad, triangle_test)

return


text:
using RoME
using Test
using TensorCast

##

@testset "test tree based autoinit on small bimodal point2 trilateration..." begin

##

N=150
fg = initfg()
getSolverParams(fg).graphinit=false
getSolverParams(fg).treeinit=false
getSolverParams(fg).useMsgLikelihoods = true
getSolverParams(fg).gibbsIters = 5
# getSolverParams(fg).inflation = 10.0

addVariable!(fg, :x0, Point2, N=N)
addFactor!(fg, [:x0], PriorPoint2(MvNormal([100.0;0], [1.0 0; 0 1])), graphinit=false)

addVariable!(fg, :x1, Point2, N=N)
addFactor!(fg, [:x1], PriorPoint2(MvNormal([0.0;100.0], [1.0 0; 0 1])), graphinit=false)

addVariable!(fg, :l1, Point2, N=N)
addFactor!(fg, [:x0;:l1], Point2Point2Range(Normal(100.0, 1.0)), graphinit=false)
addFactor!(fg, [:x1;:l1], Point2Point2Range(Normal(100.0, 1.0)), graphinit=false)

## what would

eo = [:x1; :x0; :l1]

## Mock clique solve

# tree = buildTreeReset!(fg, eo)
# # from clique 2
# sfg2 = buildCliqSubgraph(fg, tree, :x1)

# prior first

# pts, = predictbelief(sfg2,:x1, :)
# initManual!(sfg2, :x1, pts)
# pts, = predictbelief(sfg2, :l1, :)
# initManual!(sfg2, :l1, pts)

##

# plotKDE(getBelief(sfg2,:l1))

##

i = 3
while 0 < i
##
  tree = solveTree!(fg, eliminationOrder=eo);
##
  try
    #
    @cast l1_vals[j,i] := getVal(fg, :l1)[i][j]

    @test 5 < sum( 80 .< l1_vals[1,:] .< 120 )
    @test 5 < sum( -20 .< l1_vals[2,:] .< 20 )

    @test 5 < sum( -20 .< l1_vals[1,:] .< 20 )
    @test 5 < sum( 80 .< l1_vals[2,:] .< 120 )

    @cast x0_vals[j,i] := getVal(fg, :x0)[i][j]

    @test 80 < sum( 80 .< x0_vals[1,:] .< 120 )
    @test 80 < sum( -20 .< x0_vals[2,:] .< 20 )

    @cast x1_vals[j,i] := getVal(fg, :x1)[i][j]

    @test 80 < sum( -20 .< x1_vals[1,:] .< 20 )
    @test 80 < sum( 80 .< x1_vals[2,:] .< 120 )

    break
    #
  catch
    i -= 1
    i < 1 ? error("3 attempts failed") : nothing
  end
end

## testset
end

##

# using RoMEPlotting
# Gadfly.set_default_plot_size(25cm,20cm)

##

# plotKDE(fg, [:x0;:x1]) # |> SVG("/tmp/test.svg") || @async run(`eog /tmp/test.svg`)
# plotKDE(fg, :l1) # |> SVG("/tmp/test.svg") || @async run(`eog /tmp/test.svg`)

##

# IIF._getCCW(fg, :x1l1f1).inflation = 100.0
# pts = approxConv(fg, :x1l1f1, :l1)
# initManual!(fg, :l1, pts)

# plotKDE(manikde!(pts, Point2))

##


text:
# Example of training with learning rates computed from Jacobians
# Author: Gabrio Rizzuti, grizzuti3@gatech.edu
# Date: October 2020

using InvertibleNetworks, LinearAlgebra, Statistics, Flux, PyPlot, Images, TestImages
using Random; Random.seed!(123)
flag_gpu = true
# flag_gpu = false
flag_gpu && (using CUDA)

# Load image
Y = Float32.(testimage("mandril_gray")); Y = (Y.-mean(Y))./sqrt(var(Y))
nx, ny = size(Y)
Y = reshape(Y, nx, ny, 1, 1)
for i = 1:2
    global Y = wavelet_squeeze(Y)
end
nx, ny, n_ch = size(Y)[1:3]
flag_gpu && (Y = Y |> gpu)

# Initialize HINT layer
n_hidden = 2*n_ch
batchsize = 1
N = CouplingLayerHINT(n_ch, n_hidden; logdet=false, permute="full")∘
    CouplingLayerHINT(n_ch, n_hidden; logdet=false, permute="full")∘
    CouplingLayerHINT(n_ch, n_hidden; logdet=false, permute="full")
flag_gpu && (N = N |> gpu)

# Fixed input
X = randn(Float32, nx, ny, n_ch, batchsize)
flag_gpu && (X = X |> gpu)

# Loss function
loss(Y_) = 0.5f0*norm(Y-Y_)^2f0
∇loss(Y_) = Y_-Y

# Training
lr = 0.5f0
maxiter = 1000
fval = zeros(Float32, maxiter)
for i = 1:maxiter

    # Evaluate network
    Y_ = N.forward(X)

    # Evaluate objective
    fval[i] = loss(Y_)
    (mod(i, 10) == 0 || i == 1) && (print("Iteration: ", i, "; err_rel = ", sqrt(2f0*fval[i]/norm(Y)^2f0), "\n"))

    # Compute gradient
    ΔY = ∇loss(Y_)
    ΔX, Δθ, _ = N.adjointJacobian(ΔY, Y_)

    # Computing quasi-optimal step-length
    JΔθ, _ = N.jacobian(0f0.*X, Δθ, X)
    α = dot(JΔθ, ΔY)/dot(JΔθ, JΔθ)

    # Update parameters
    θ = get_params(N)
    set_params!(N, θ-lr*α.*Δθ)

end

Y_ = N.forward(X)
for i = 1:2
    global Y_ = wavelet_unsqueeze(Y_)
    global Y = wavelet_unsqueeze(Y)
end
flag_gpu && (Y = Y |> cpu; Y_ = Y_ |> cpu)


text:
#Test EasyData's read/write capabilities using sample/ plots
#-------------------------------------------------------------------------------
#=
Look for EasyData._write() & EasyData._read() near end of file.
=#

module CMDimData_SampleUsage

using CMDimData
using CMDimData.EasyPlot
using CMDimData.EasyData
CMDimData.@includepkg EasyPlotInspect

#==Constants
===============================================================================#
demolist = EasyPlot.demofilelist()

#==Helper functions
===============================================================================#
printsep(label, sep="-") = println("\n", label, "\n", repeat(sep, 80))
printheader(label) = printsep(label, "=")

#==Main Code
===============================================================================#
pdisp = EasyPlot.GUIDisplay(:InspectDR)
for filepath in demolist
	corename = splitext(basename(filepath))[1]
	filename = "sample_$corename.hdf5"
	printsep("Executing $filepath...")
	pcoll = evalfile(filepath)
	display(pdisp, pcoll)

	#EasyData portion:
	@info("Writing $filename...")
	EasyData.writeplot(filename, pcoll)
	@info("Reading back $filename...")
	plot2 = EasyData.readplot(filename)
	plot2.title = "Compare Results: " * pcoll.title
	display(pdisp, plot2)
end

end
:SampleCode_Executed


text:
"""
    struct ElementFEAInfo{dim, T}
        Kes::AbstractVector{<:AbstractMatrix{T}}
        fes::AbstractVector{<:AbstractVector{T}}
        fixedload::AbstractVector{T}
        cellvolumes::AbstractVector{T}
        cellvalues::CellValues{dim, T}
        facevalues::FaceValues{<:Any, T}
        metadata::Metadata
        black::AbstractVector
        white::AbstractVector
        varind::AbstractVector{Int}
        cells
    end

An instance of the `ElementFEAInfo` type stores element information such as:
- `Kes`: the element stiffness matrices,
- `fes`: the element load vectors,
- `cellvolumes`: the element volumes,
- `cellvalues` and `facevalues`: two `Ferrite` types that facilitate cell and face iteration and queries.
- `metadata`: that stores degree of freedom (dof) to node mapping, dof to cell mapping, etc.
- `black`: a `BitVector` such that `black[i]` is 1 iff element `i` must be part of any feasible design.
- `white`: a `BitVector` such that `white[i]` is 1 iff element `i` must never be part of any feasible design.
- `varind`: a vector such that `varind[i]` gives the decision variable index of element `i`.
- `cells`: the cell connectivities.
"""
@params struct ElementFEAInfo{dim,T}
    Kes::AbstractVector{<:AbstractMatrix{T}}
    fes::AbstractVector{<:AbstractVector{T}}
    fixedload::AbstractVector{T}
    cellvolumes::AbstractVector{T}
    cellvalues::CellValues{dim,T,<:Any}
    facevalues::FaceValues{<:Any,T,<:Any}
    metadata::Metadata
    black::AbstractVector
    white::AbstractVector
    varind::AbstractVector{Int}
    cells::Any
end

function Base.show(io::Base.IO, ::MIME"text/plain", efeainfo::ElementFEAInfo)
    return print(
        io,
        "ElementFEAInfo: Kes |$(length(efeainfo.Kes))|, fes |$(length(efeainfo.fes))|, fixedload |$(length(efeainfo.fixedload))|, cells |$(length(efeainfo.cells))|",
    )
end

"""
    ElementFEAInfo(sp, quad_order=2, ::Type{Val{mat_type}}=Val{:Static}) where {mat_type}

Constructs an instance of `ElementFEAInfo` from a stiffness problem `sp` using a Gaussian quadrature order of `quad_order`. The element matrix and vector types will be:
1. `SMatrix` and `SVector` if `mat_type` is `:SMatrix` or `:Static`, the default,
2. `MMatrix` and `MVector` if `mat_type` is `:MMatrix`, or
3. `Matrix` and `Vector` otherwise.

The static matrices and vectors are more performant and GPU-compatible therefore they are used by default.
"""
function ElementFEAInfo(
    sp, quad_order=2, ::Type{Val{mat_type}}=Val{:Static}
) where {mat_type}
    Kes, weights, dloads, cellvalues, facevalues = make_Kes_and_fes(
        sp, quad_order, Val{mat_type}
    )
    element_Kes = convert(
        Vector{<:ElementMatrix},
        Kes;
        bc_dofs=sp.ch.prescribed_dofs,
        dof_cells=sp.metadata.dof_cells,
    )
    fixedload = Vector(make_cload(sp))
    assemble_f!(fixedload, sp, dloads)
    cellvolumes = get_cell_volumes(sp, cellvalues)
    cells = sp.ch.dh.grid.cells
    return ElementFEAInfo(
        element_Kes,
        weights,
        fixedload,
        cellvolumes,
        cellvalues,
        facevalues,
        sp.metadata,
        sp.black,
        sp.white,
        sp.varind,
        cells,
    )
end

"""
    struct GlobalFEAInfo{T, TK<:AbstractMatrix{T}, Tf<:AbstractVector{T}, Tchol}
        K::TK
        f::Tf
        cholK::Tchol
    end

An instance of `GlobalFEAInfo` hosts the global stiffness matrix `K`, the load vector `f` and the cholesky decomposition of the `K`, `cholK`.
"""
@params mutable struct GlobalFEAInfo{T}
    K::AbstractMatrix{T}
    f::AbstractVector{T}
    cholK::Any
    qrK::Any
end
function Base.show(::IO, ::MIME{Symbol("text/plain")}, ::GlobalFEAInfo)
    return println("TopOpt global FEA information")
end

"""
    GlobalFEAInfo(::Type{T}=Float64) where {T}

Constructs an empty instance of `GlobalFEAInfo` where the field `K` is an empty sparse matrix of element type `T` and the field `f` is an empty dense vector of element type `T`.
"""
GlobalFEAInfo(::Type{T}=Float64) where {T} = GlobalFEAInfo{T}()
function GlobalFEAInfo{T}() where {T}
    return GlobalFEAInfo(sparse(zeros(T, 0, 0)), zeros(T, 0), cholesky(one(T)), qr(one(T)))
end

"""
    GlobalFEAInfo(sp::StiffnessTopOptProblem)

Constructs an instance of `GlobalFEAInfo` where the field `K` is a sparse matrix with the correct size and sparsity pattern for the problem instance `sp`. The field `f` is a dense vector of the appropriate size. The values in `K` and `f` are meaningless though and require calling the function `assemble!` to update.
"""
function GlobalFEAInfo(sp::StiffnessTopOptProblem)
    K = initialize_K(sp)
    f = initialize_f(sp)
    return GlobalFEAInfo(K, f)
end
function GlobalFEAInfo(
    K::Union{AbstractSparseMatrix,Symmetric{<:Any,<:AbstractSparseMatrix}}, f
)
    chol = cholesky(spdiagm(0 => ones(size(K, 1))))
    qrfact = qr(spdiagm(0 => ones(size(K, 1))))
    return GlobalFEAInfo{eltype(K),typeof(K),typeof(f),typeof(chol),typeof(qrfact)}(
        K, f, chol, qrfact
    )
end

"""
    GlobalFEAInfo(K, f)

Constructs an instance of `GlobalFEAInfo` with global stiffness matrix `K` and load vector `f`.
"""
function GlobalFEAInfo(K, f)
    chol = cholesky(Matrix{eltype(K)}(I, size(K)...))
    qrfact = qr(Matrix{eltype(K)}(I, size(K)...))
    return GlobalFEAInfo(K, f, chol, qrfact)
end

"""
    get_cell_volumes(sp::StiffnessTopOptProblem{dim, T}, cellvalues)

Calculates an approximation of the element volumes by approximating the volume integral of 1 over each element using Gaussian quadrature. `cellvalues` is a `Ferrite` struct that facilitates the computation of the integral. To initialize `cellvalues` for an element with index `cell`, `Ferrite.reinit!(cellvalues, cell)` can be called. Calling `Ferrite.getdetJdV(cellvalues, q_point)` then computes the value of the determinant of the Jacobian of the geometric basis functions at the point `q_point` in the reference element. The sum of such values for all integration points is the volume approximation.
"""
function get_cell_volumes(sp::StiffnessTopOptProblem{dim,T}, cellvalues) where {dim,T}
    dh = sp.ch.dh
    cellvolumes = zeros(T, getncells(dh.grid))
    for (i, cell) in enumerate(CellIterator(dh))
        reinit!(cellvalues, cell)
        cellvolumes[i] = sum(
            Ferrite.getdetJdV(cellvalues, q_point) for
            q_point in 1:Ferrite.getnquadpoints(cellvalues)
        )
    end
    return cellvolumes
end


text:
################################################################################
#------------------------------ Init Functions ---------------------------------
################################################################################
#Init the Net
function init_θ_combined(net_s::Array{Int,1})
    x      = net_s[1]
    hidden = net_s[2:end-1]
    ysize  = net_s[end]

    w = []
    for y in [hidden...]
        push!(w, xavier(y, x))
        push!(w, zeros(y, 1))
        x = y
    end
    push!(w, xavier(ysize, x)) # Prob actions
    push!(w, zeros(ysize, 1))
    push!(w, xavier(1, x))     # Value function
    push!(w, zeros(1, 1))

    return w
end

#Init the Net
function init_θ_single(net_s::Array{Int,1})
    x      = net_s[1]
    hidden = net_s[2:end-1]
    ysize  = net_s[end]

    w = []
    for y in [hidden...]
        push!(w, xavier(y, x))
        push!(w, xavier(y, 1))
        x = y
    end
    push!(w, xavier(ysize, x)) # Prob actions
    push!(w, xavier(ysize, 1))

    return w
end

################################################################################
#------------------------------ Output Functions -------------------------------
################################################################################
#Evaluate the Network (Actions + Value Function)
@everywhere function mlp_nlayer(w::Any,x::AbstractVector,params::Params)
    σ       = params.σ
    u_range = params.u_range

    for i = 1:2:length(w)-4
        x = relu.(w[i] * x .+ w[i+1])
    end
    mu_act   = vec((1 .- σ) .* u_range.*tanh.(w[end-3] * x .+ w[end-2]))
    value    = (w[end-1] * x .+ w[end])[1]
    return mu_act, value
end

#Evaluate the Policy-Network
@everywhere function mlp_nlayer_policy(w::Any,x::AbstractVector,params::Params)
    σ       = params.σ
    u_range = params.u_range

    for i = 1:2:length(w)-2
        x = relu.(w[i] * x .+ w[i+1])
    end
    mu_act   = (1 .- σ) .* u_range.*tanh.(w[end-1] * x .+ w[end])
    return vec(mu_act)
end

#Evaluate the Policy-Network
@everywhere function mlp_nlayer_policy_skeleton(w::Any,x::AbstractVector,params::Params)
    σ       = params.σ
    u_range = params.u_range

    for i = 1:2:length(w)-2
        x = relu.(w[i] * x .+ w[i+1])
    end
    mu_act   = 0.5 .* (1 .- 1 .* σ) .* u_range.*tanh.(w[end-1] * x .+ w[end]) .+ 0.5 * u_range
    return vec(mu_act)
end

#Evaluate the Value-Network
@everywhere function mlp_nlayer_value(w::Any,x::AbstractVector,params::Params)
    for i = 1:2:length(w)-2
        x = relu.(w[i] * x .+ w[i+1])
    end
    value    = (w[end-1] * x .+ w[end])[1]
    return value
end

#Get the log(pi(a|s))
@everywhere function logpdf(a_t,out_t,params::Params)
    fac = Float64(-log(sqrt(2pi)))
    r   = (a_t - out_t) ./ params.σ
    return -r .* r .* Float64(0.5) - (log.(params.σ)) .+ fac
end

################################################################################
#---------------------------- Gradient Functions -------------------------------
################################################################################

#Idea 1: Do everything in one go -> But bad performance!
@everywhere function lossfunc_do_all_combined(θ::Any,params::Params)
    #generate from struct:
    l_traj = params.l_traj
    log_pi = Any[]
    r_t    = Float64[]
    a_t    = Float64[]
    val_t  = Any[]
    s_t1   = copy(params.s_init)
    totalr = 0

    #Sum over one Trajectory and then build the gradient:
    for i_t = 1: l_traj
        out_t, loc_val  = mlp_nlayer(θ,s_t1,params) #evaluate
        loc_a_t         = sample_action(out_t, params.σ)[1]
        s_t, loc_r_t    = environment_call(s_t1,loc_a_t,params)

        loc_log_pi = logpdf(loc_a_t,out_t[1],params)

        #Extend the running storage
        push!(val_t ,loc_val)
        push!(r_t   ,loc_r_t[1])
        push!(a_t   ,loc_a_t)
        push!(log_pi,loc_log_pi)
        totalr += loc_r_t

    end
    println(totalr)

    #Get the Accumulated Reward:
    r_t = discount(r_t,params.γ_disc)
    A_t = r_t - val_t
    A_t = reshape(A_t, 1, length(A_t))

    #Return the final Sum over one Trajectory
    sum(A_t .* hcat(log_pi)) + L2Reg(A_t)
end

@everywhere ĝ_traj = grad(lossfunc_do_all_combined,1) #Gradient construction

#Better: Do it after the episode was built:
function loss_a3c_combined(θ::Any,params::Params,traj::Trajectory)
    #Take params from struct:
    l_traj = traj.len
    γ      = params.γ_disc
    λ      = params.λ_critic

    #Init
    Err_actor  = 0
    Err_critic = 0
    s_t = traj.s_t[1,:]
    val_old    = 0;
    A_t        = 0;

    #Hope: improved peroformance by going backwards...
    for i_t = l_traj-1:-1:1
        #Take from Struct
        s_t = traj.s_t[i_t,:]
        a_t = traj.a_t[i_t]
        r_t = traj.r_t[i_t]

        out_t, val_new = mlp_nlayer(θ,s_t,params)          #Evaluate the Network
        d_val          = r_t + γ*val_old - val_new         #Get local TD-error
        val_old        = val_new                           #Update old error
        A_t            = d_val + (γ*λ)*A_t                 #GAE update -> Advantage
        Err_actor     += sum(logpdf(a_t,out_t,params)*A_t) #Actor correction
        Err_critic    += 0.5*(A_t * A_t)                   #Critic correction
    end

    #Return the final Error_Sum
    return Err_actor + (Err_critic / l_traj)
end

a3cgrad_combined = grad(loss_a3c_combined,1) #Gradient construction

#Decompose into separate Value und Policy Function:
#-------------------------------------------------------------------------------

#Better: Do it after the episode was built:
@everywhere function loss_policy(θa::Any,params::Params,traj::Trajectory)
    #Take params from struct:
    l_traj = traj.len

    #Init
    Err_actor  = 0
    s_t = traj.s_t[1,:]

    #Hope: improved peroformance by going backwards...
    for i_t = 1:l_traj
        #Take from Struct
        s_t = traj.s_t[i_t,:]
        a_t = traj.a_t[i_t,:]
        r_t = traj.r_t[i_t]
        A_t = traj.A_t[i_t]

        #Only for Segway
        s_mod = s_t
        #=---------------------------------
        s_mod        = zeros(8)
        s_mod[1:2]   = s_t[1:2]
        s_mod[3]     = sin(s_t[3])
        s_mod[4]     = cos(s_t[3])
        s_mod[5:end] = s_t[4:end]
        #-------------------------------=#

        out_t = mlp_nlayer_policy_norm(θa,s_mod,params)  #.+ 0.5 #Evaluate the Network
        Err_actor += sum(logpdf(a_t,out_t,params) * A_t)  #Actor correction
    end

    #Return the final Error_Sum
    return Err_actor
end

@everywhere grad_policy = grad(loss_policy,1) #Gradient construction

#Better: Do it after the episode was built:
@everywhere function loss_value(θv::Any,params::Params,traj::Trajectory)
    #Take params from struct:
    l_traj = traj.len
    γ      = params.γ_disc
    λ      = params.λ_critic

    #Init
    Err_critic = 0
    s_t = traj.s_t[1,:]
    V_t_true   = 0;
    e_t        = 0;

    #Hope: improved peroformance by going backwards...
    for i_t = l_traj-1:-1:1
        #Take from Struct
        s_t = traj.s_t[i_t,:]
        r_t = traj.r_t[i_t]
        val_t2 = traj.val_t[i_t+1]

        #Only for Segway
        s_mod = s_t
        #=---------------------------------
        s_mod        = zeros(8)
        s_mod[1:2]   = s_t[1:2]
        s_mod[3]     = sin(s_t[3])
        s_mod[4]     = cos(s_t[3])
        s_mod[5:end] = s_t[4:end]
        #-------------------------------=#

        V_t_est     = mlp_nlayer_value_norm(θv,s_mod,params)[1] #Estimate the Value
        V_t_true    = traj.R_t[i_t]                      #Get the true Value (TD1)
        #V_t_true    = r_t + γ*val_t2                     #Get the Bellman Error (TD0)
        e_t         = V_t_est - V_t_true                 #Calculate the error
        Err_critic += (e_t * e_t)                        #Collect the error
    end

    #Return the final Error_Sum
    return Err_critic / l_traj
end

@everywhere grad_value = grad(loss_value,1) #Gradient construction

#Better: Do it after the episode was built:
@everywhere function loss_value_bs(θv::Any,params::Params,traj::Trajectory)
    #Take params from struct:
    l_traj = traj.len
    γ      = params.γ_disc
    λ      = params.λ_critic

    #Init
    Err_critic = 0
    s_t = traj.s_t[1,:]
    V_t_true   = 0;
    e_t        = 0;

    #Hope: improved peroformance by going backwards...
    for i_t = l_traj-1:-1:1
        #Take from Struct
        s_t   = traj.s_t[i_t,:]
        r_t   = traj.r_t[i_t]
        val_t = traj.val_t[i_t]

        V_t_est     = mlp_nlayer_value(θv,s_t,params)[1] #Estimate the Value
        e_t         = V_t_est - r_t - γ * traj.val_t[i_t+1]
        Err_critic += (e_t * e_t)                        #Collect the error

        #Update eligibility
        #V_t_true    = (1-λ) * val_t + λ * (r_t + γ * V_t_true) #Get the true Value - bootstrapped
    end

    #Return the final Error_Sum
    return Err_critic / l_traj
end

@everywhere grad_value_bs = grad(loss_value_bs,1) #Gradient construction

function supervised_training!(all_traj::Traj_Data,θv::Array{Any,1}, params::Params,opt_para_v)
    test_num = Int(floor(0.7 * all_traj.n_traj))
    test_err_old = Inf
    test_err = 1e10
    ĝ_v = copy(θv)
    iter = 0
    println()
    println("----start-training-------")
    while test_err_old > test_err
        test_err_old = test_err
        test_err = 0
        iter += 1

        #Do a Training Batch:
        for i_train = 1: test_num
            ĝ_v = grad_value(θv,params,all_traj.traj[i_train])
            Knet.update!(θv, ĝ_v, opt_para_v)
        end

        #Evaluate the Error
        for i_test = test_num+1:all_traj.n_traj
            test_err += loss_value(θv,params,all_traj.traj[i_test])
        end

        if (((test_err_old/ test_err)^2) < 1.05) && iter > 3
            break
        end

        println(test_err)
    end
    println("------end-training-------")
end

################################################################################
#------------------------------ Util Functions ---------------------------------
################################################################################

#Filter NaN
@everywhere function filterNet_NaN!(θ::Any)
    for i_layer = 1: length(θ)
        for i_element = 1: length(θ[i_layer])
             if isnan(θ[i_layer][i_element])
                 θ[i_layer][i_element] = 0
             end
         end
     end
end

#Build the Multiplikation of the Neural Network
@everywhere function mult_net(θ1::Any,θ2::Any)
    mult_θ = similar(θ1)
    for i_layer = 1:length(θ1)
        mult_θ[i_layer] = θ1[i_layer] .* θ2[i_layer]
    end
    return mult_θ
end

#Build the Divison of two Nets
@everywhere function division_net(θ1::Any,θ2::Any)
    div_θ = similar(θ1)
    for i_layer = 1:length(θ1)
        div_θ[i_layer] = θ1[i_layer] ./ θ2[i_layer]
    end
    return div_θ
end

#create net to vector
@everywhere function net_to_vec(θ::Any)
    net_vec =  Array{Float64}(undef,0)
    for i_layer = 1:length(θ)
        append!(net_vec,vec(θ[i_layer]))
    end
    return net_vec
end

#create vector to net
@everywhere function vec_to_net(net_vec::Array{Float64,1}, θ::Any)
    vec_net = copy(θ) * 0
    global_count = 1
    for i_layer = 1:length(θ)
        for i_loc = 1:length(vec(θ[i_layer]))
            vec_net[i_layer][i_loc] = net_vec[global_count]
            global_count += 1
        end
    end
    return vec_net
end

#"""
#Return the Probability of a_t given the output of the network and action a_t
#"""
@everywhere function get_π_prob(net_out::Any,σ::Array{Float64,1},a_t::Array{Float64,1})
    return (1 ./ (σ.*sqrt(2*pi))) .* exp.(((net_out .- a_t) .* (net_out .- a_t)) ./ (-2 .* σ .* σ))
end

max_stars_repo_name: NikonPic/REINFORCEMENT-JULIA

max_issues_repo_name: NikonPic/REINFORCEMENT-JULIA

max_forks_repo_name: NikonPic/REINFORCEMENT-JULIA


text:
function hull_reformulation{HR<:HRepresentation}(model::Model, regions::Vector{HR}, x::AbstractVector{Variable})
    n = length(x)
    R = length(regions)
    @boundscheck begin
        for region in regions
            # ensure that dimension of region matches number of variables
            @assert fulldim(region) == n

            # ensure that there are no rays, so that setting the bounds of the inequalities to zero will
            # result in only the origin being inside the polyhedron
            vr = vrep(polyhedron(region))
            @assert nrays(vr) == 0
        end
    end

    Hrs = @variable(model, [r = 1 : R], Bin) # indicator variables, 4.2c
    @constraint(model, sum(Hrs) == 1) # 4.2b

    xrs = [@variable(model, [k = 1 : n]) for r = 1 : R] # auxiliary variables
    @constraint(model, x .== sum(xrs)) # 4.4a

    for (region, Hr, xr) in zip(regions, Hrs, xrs)
        for ineq in ineqs(region)
            @constraint(model, ineq.a ⋅ xr <= Hr * ineq.β) # 4.4b
        end
        for eq in eqs(region)
            @constraint(model, eq.a ⋅ xr == Hr * eq.β) # 4.4b
        end
    end
    Hrs
end


text:
function MMPolya(alpha0::AbstractArray,delta::AbstractArray,alpha::AbstractArray,Y::AbstractArray,TS::AbstractArray,MY::AbstractArray,MTS::AbstractArray,mu::AbstractArray,L::Int64,I::AbstractArray,J::Int64,K::Int64,MMNum::Int64)
    Newdelta=similar(delta)
    Newalpha0=similar(alpha0)
    for t2 in 1:MMNum
        DS=zeros(J,K,L)
        DS2=zeros(J,K,L)
        DS1=zeros(K,L)
    @inbounds  for l in 1:L
                    @inbounds for k in 1:K
                        salpha=sum(view(alpha,:,k,l))

                        n2=maximum(view(mu[l],k,:)[TS[l].>=1])
                        @inbounds for c1 in 0:MTS[l]-1
                            DS1[k,l]+=sum(exp.(view(mu[l],k,:)[TS[l].>=c1+1].-n2))/(salpha+c1)
                        end
                        DS1[k,l]=n2+log(DS1[k,l])

                        @inbounds   for j in 1:J
                            if length(view(mu[l],k,:)[Y[l][j,:].>=1])==0
                                DS2[j,k,l]=-10000
                                DS[j,k,l]=0
                            else
                                n3=maximum(view(mu[l],k,:)[Y[l][j,:].>=1])
                                @inbounds   for c2 in 0:MY[j,l]-1
                                    DS2[j,k,l]+=sum(exp.(view(mu[l],k,:)[Y[l][j,:].>=c2+1].-n3))/(view(alpha,j,k,l).+c2)
                                end
                                DS2[j,k,l]=n3+log(DS2[j,k,l])
                                DS[j,k,l]=exp(DS2[j,k,l]-DS1[k,l])
                             end

                        end
                    end
                end
        Newdelta=delta.*DS
        Newalpha0=alpha0.*(sum(exp.(DS2),dims=3)./sum(exp.(DS1),dims=2)')
        delta=copy(Newdelta)
        alpha0=copy(Newalpha0)
        alpha0[alpha0.<1e-100].=1e-100
        for l in 1:L
            alpha[:,:,l]=alpha0+delta[:,:,l]
        end
    end
    alpha0, delta, alpha
end


text:
#####################################################################
# Base functions.
#####################################################################
sphere(x) = sum(abs2, x)

"""
Schwefel's ellipsoid.
"""
function ellipsoid(x)
    res = 0.0
    cumsum = 0.0
    for xx in x
        cumsum += xx
        res += abs2(cumsum)
    end
    return res
end

function elliptic(x)
    D = length(x)
    condition = 1e+6
    coefficients = condition .^ range(0, stop=1, length=D)
    return sum(coefficients .* abs2.(x))
end

#= FIXME this should be more efficient implementation, but switching to it would reset the benchmarks
function elliptic(x)
    condition = 1e+6
    powers = range(0, stop=1, length=length(x))
    @inbounds return sum(i -> (condition^powers[i])*abs2(x[i]), eachindex(x))
end
=#

function rastrigin(x)
    D = length(x)
    10 * D + sum(abs2, x) - 10 * sum(xx -> cos(2π * xx), x)
end

function ackley(x)
    D = length(x)
    try
        return 20 - 20*exp(-0.2.*sqrt(sum(abs2, x)/D)) - exp(sum(xx -> cos(2π * xx), x)/D) + MathConstants.e
    catch ex
        # Sometimes we have gotten a DomainError from the cos function so we protect this call
        println(ex)
        println("For input x = ", x)
        # Return a very large fitness value to indicate that this is NOT the solution we want.
        # TODO: Fix better way to handle this!
        return 9.99e100
    end
end

function schwefel1_2(x)
    D = length(x)
    partsums = zeros(D)
    partsum = 0
    for i in 1:D
        partsum += x[i]
        partsums[i] = partsum
    end
    return sum(abs2, partsums)
end

function rosenbrock(x)
    n = length(x)
    return sum(100*(view(x, 2:n) .- view(x, 1:(n-1)).^2).^2 .+ (view(x, 1:(n-1)) .- 1).^2)
end

#= FIXME this should be more efficient implementation, but switching to it would reset the benchmarks
function rosenbrock(x)
    @inbounds return sum(i -> 100.0*abs2(x[i+1] - x[i]^2) + abs2(x[i] - 1.0),
                         1:length(x)-1)
end
=#

step(x) = sum(xx -> ceil(xx + 0.5), x)

function griewank(x)
    n = length(x)
    1 + (1/4000)*sum(abs2, x) - prod(cos.(x ./ sqrt.(1:n)))
end

schwefel2_22(x) = sum(abs, x) + prod(abs, x)

schwefel2_21(x) = maximum(abs, x)

# I'm unsure about this one since it does not return the expected minima at
# [1.0, 1.0].
function schwefel2_26(x)
    D = length(x)
    418.98288727243369 * D - sum(xx -> xx * sin(sqrt(abs(xx))), x)
end

cigar(x) = x[1]^2 + 1e6 * sum(abs2, view(x, 2:length(x)))

cigtab(x) = x[1]^2 + 1e8 * x[end]^2 + 1e4 * sum(abs2, view(x, 2:length(x)-1))

"""
Generic function to define `ShekelN` problems.
"""
function shekel(x, a, c)
    @assert length(x) == size(a, 2)
    @assert length(c) == size(a, 1)
    sum = 0.0
    @inbounds for i in eachindex(c)
        den = 0.0
        for j in 1:size(a, 2)
            den += abs2(x[j] - a[i,j])
        end
        sum = sum - 1 / (den + c[i])
    end
    return sum
end

# constants for `Shekel10`
Shekel10_A = [4 4 4 4; 1 1 1 1; 8 8 8 8; 6 6 6 6; 3 7 3 7; 2 9 2 9; 5 5 3 3; 8 1 8 1; 6 2 6 2; 7 3.6 7 3.6]
Shekel10_C = [0.1, 0.2, 0.2, 0.4, 0.4, 0.6, 0.3, 0.7, 0.5, 0.5]

"""
`Shekel10` is a 4D, multi-minima, non-separable test problem. Our implementation
is based on the C code in:
http://www.math.ntu.edu.tw/~wwang/cola_lab/test_problems/multiple_opt/multiopt_prob/Shekel10/Shekel10.c
"""
shekel10(x) = shekel(x, Shekel10_A, Shekel10_C)

# constants for `Shekel7`
Shekel7_A = [4 4 4 4; 1 1 1 1; 8 8 8 8; 6 6 6 6; 3 7 3 7; 2 9 2 9; 5 5 3 3]
Shekel7_C = [0.1, 0.2, 0.2, 0.4, 0.4, 0.6, 0.3]

"""
`Shekel7` is a 4D, multi-minima, non-separable test problem. Our implementation
is based on the C code in:
http://www.math.ntu.edu.tw/~wwang/cola_lab/test_problems/multiple_opt/multiopt_prob/Shekel7/Shekel7.c
"""
shekel7(x) = shekel(x, Shekel7_A, Shekel7_C)

# constants for `Shekel5`
Shekel5_A = [4 4 4 4; 1 1 1 1; 8 8 8 8; 6 6 6 6; 3 7 3 7]
Shekel5_C = [0.1, 0.2, 0.2, 0.4, 0.4]

"""
`Shekel5` is a 4D, multi-minima, non-separable test problem. Our implementation
is based on the C code in:
http://www.math.ntu.edu.tw/~wwang/cola_lab/test_problems/multiple_opt/multiopt_prob/Shekel5/Shekel5.c
"""
shekel5(x) = shekel(x, Shekel5_A, Shekel5_C)

"""
Generic function for `Hartman N` problem family.
"""
function hartman(x, alpha, A, P)
    @assert size(A) == size(P)
    @assert length(x) == size(P, 2)
    @assert length(alpha) == size(P, 1)
    res = 0.0
    @inbounds for i in 1:length(alpha)
        isum = 0.0
        for j in 1:size(A, 2)
            isum += A[i,j] * (x[j] - P[i,j])^2
        end
        res -= alpha[i] * exp(-isum)
    end
    res
end

# constants for `Hartman6`
Hartman6_alpha = [1.0 1.2 3.0 3.2]
Hartman6_A = [10 3 17 3.50 1.7 8; 0.05 10 17 0.1 8 14; 3 3.5 1.7 10 17 8; 17 8 0.05 10 0.1 14]
Hartman6_P = 1e-4 * [1312 1696 5569 124 8283 5886; 2329 4135 8307 3736 1004 9991; 2348 1451 3522 2883 3047 6650; 4047 8828 8732 5743 1091 381]

"""
`Hartman 6D` is a multi-minima, non-separable test problem. Our implementation is based on:
http://www.sfu.ca/~ssurjano/hart6.html
"""
hartman6(x) = hartman(x, Hartman6_alpha, Hartman6_A, Hartman6_P)

# constants for `Hartman3`
Hartman3_alpha = [1.0 1.2 3.0 3.2]
Hartman3_A = [3.0 10 30; 0.1 10 35; 3.0 10 30; 0.1 10 36]
Hartman3_P = 1e-4 * [3689 1170 2673; 4699 4387 7470; 1091 8732 5547; 381 5743 8828]

"""
`Hartman 3D` is a multi-minima, non-separable test problem. Our implementation is based on:
http://www.sfu.ca/~ssurjano/hart3.html
However, we get a different global minima than the one stated on that page.
The minima should be -3.86278 but we get a different one so use that in the problem spec:
"""
hartman3(x) = hartman(x, Hartman3_alpha, Hartman3_A, Hartman3_P)

#####################################################################
# S2 functions in addition to the base functions above. As stated
# in table II of the JADE paper: http://150.214.190.154/EAMHCO/pdf/JADE.pdf
#####################################################################

quartic(x) = sum((i,x) -> i*x^4, enumerate(x))

noisy_quartic(x) = quartic(x) + rand()

s2_step(x) = sum(xx -> abs2(ceil(xx + 0.5)), x)

#####################################################################
# Misc other interesting optimization functions and families.
#####################################################################

# Schwefel 2.13 is a hard one...
#function f=schwefel_213(x)
#global initial_flag
#persistent a b A alpha
#[ps,D]=size(x);
#if initial_flag==0
#    initial_flag=1;
#    load schwefel_213_data
#    if length(alpha)>=D
#        alpha=alpha(1:D);a=a(1:D,1:D);b=b(1:D,1:D);
#    else
#        alpha=-3+6*rand(1,D);
#        a=round(-100+200.*rand(D,D));
#        b=round(-100+200.*rand(D,D));
#    end
#    alpha=repmat(alpha,D,1);
#    A=sum(a.*sin(alpha)+b.*cos(alpha),2);
#end
#for i=1:ps
#    xx=repmat(x(i,:),D,1);
#    B=sum(a.*sin(xx)+b.*cos(xx),2);
#    f(i,1)=sum((A-B).^2,1);
#end

"""
Generator for the family of deceptive functions from the
Cuccu2011 paper on novelty-based restarts. We have vectorized it to allow
more than 1D versions. The Cuccu2011 paper uses the following values for
```
(l, w) = [(5, 0),  (15, 0),  (30, 0),
          (5, 2),  (15, 2),  (30, 2),
          (5, 10), (15, 10), (30, 10)]
```
and notes that `(15, 2)` and `(30, 2)` are the most difficult instances.
"""
function deceptive_cuccu2011(l, w)
    (x) -> begin
        sumabsx = sum(abs, x)
        if sumabsx <= 1
            return sum(abs2, x)
        elseif sumabsx >= l+1
            return sum(xx -> abs2(abs(xx) - l), x)
        else
            return (1 - 0.5 * sum(xx -> abs2(sin( (π * w * (abs(xx) - 1)) / l )), x))
        end
    end
end

# Deceptive/hardest instances:
deceptive_cuccu2011_15_2 = deceptive_cuccu2011(15, 2)
deceptive_cuccu2011_30_2 = deceptive_cuccu2011(30, 2)

"""
From section 3, page 7, of Tsallis1996 paper:
    Tsallis and Stariolo, "Generalized simulated annealing", Physica A, 1996.
available from http://www.if.ufrgs.br/~stariolo/publications/TsSt96_PhysA233_395_1996.pdf
the original paper used this as a 4-dimensional problem but here it is generalized.
"""
energy_tsallis1996(x::AbstractArray) = sum(xx -> abs2(xx^2 - 8.0), x) + 5.0*sum(x) + 57.3276


text:
# write your own tests here
module PkgTests

using Base.Test

import BEAST

include("test_fourier.jl")

include("test_basis.jl")
#include("test_directproduct.jl")
include("test_raviartthomas.jl")
include("test_rt.jl")
include("test_dvg.jl")
include("test_bcspace.jl")
include("test_trace.jl")
include("test_timebasis.jl")

include("test_gram.jl")
include("test_vector_gram.jl")

include("test_wiltonints.jl")
include("test_nitsche.jl")
include("test_nitschehh3d.jl")
include("test_farfield.jl")
include("test_rtports.jl")
include("test_tdop_scaling.jl")
include("test_tdrhs_scaling.jl")

include("test_assemblerow.jl")

try
    Pkg.installed("BogaertInts10")
    info("`BogaertInts10` detected. Including relevant tests.")
    include("test_bogaertints.jl")
catch
    info("`Could not load BogaertInts10`. Related tests skipped.")
end

end


text:
__precompile__()

module Websocket
using HTTP, Base64, Sockets, MbedTLS
import Sockets: listen

export WebsocketServer, WebsocketClient, WebsocketConnection, WebsocketError, serve, send, ping, listen, emit, throwWSerror, logWSerror

include("opt/vars.jl")
include("opt/utils.jl")
include("lib/WebsocketConnection.jl")
include("lib/WebsocketClient.jl")
include("lib/WebsocketServer.jl")

"""
    logWSerror(err::Union{Exception, WebsocketError})
A convenience method to lift errors into a scope and log them.
# Example
```julia
using Websocket

server = WebsocketServer()
ended = Condition()

listen(server, :client) do client end
listen(server, :connectError) do err::WebsocketError
    notify(ended, err)
end

@async serve(server, 8080, "notahost")

reason = wait(ended)
reason isa Exception && logWSerror(reason)
```
"""
function logWSerror(err::WebsocketError)
    err.log()
end
function logWSerror(err::Exception)
    @error err
end
"""
    throwWSerror(err::Union{Exception, WebsocketError})
A convenience method to lift fatal errors into a scope and throw them.
# Example
```julia
using Websocket

server = WebsocketServer()
ended = Condition()

listen(server, :client) do client end
listen(server, :connectError) do err::WebsocketError
    notify(ended, err)
end

@async serve(server, 8080, "notahost")

reason = wait(ended)
reason isa Exception && throwWSerror(reason)
```
"""
function throwWSerror(err::WebsocketError)
    err.log()
    throw(err)
end
function throwWSerror(err::Exception)
    throw(err)
end

end

max_stars_repo_name: citkane/Websocket

max_issues_repo_name: citkane/Websocket

max_forks_repo_name: citkane/Websocket


text:
immutable Binary <: DiscreteUnivariateDistribution
  a::Int # Failure
  b::Int # Success
  p::Float64

  function Binary(a::Real, b::Real, p::Real)
    # @check_args(Binary, a < b && zero(p) <= p <= one(p)) # Possibly a bug in Distributions, this line is not working
    new(a, b, p)
  end

  function Binary(a::Real, b::Real)
    # @check_args(Binary, a < b) # Possibly a bug in Distributions, this line is not working
    new(a, b, 0.5)
  end

  Binary() = new(0, 1, 0.5)
end

@distr_support Binary d.a d.b

### Parameters

succprob(d::Binary) = d.p
failprob(d::Binary) = 1.0 - d.p

params(d::Binary) = (d.a, d.b, d.p)

### Properties

mean(d::Binary) = d.a*failprob(d)+d.b*succprob(d)

### Evaluation

pdf(d::Binary, x::Int) = x == d.a ? failprob(d) : x == d.b ? succprob(d) : 0.0

pdf(d::Binary) = Float64[failprob(d), succprob(d)]

logpdf(d::Binary, x::Int) = x == d.a ? log(failprob(d)) : x == d.b ? log(succprob(d)) : -Inf

### Sampling

rand(d::Binary) = rand() <= succprob(d) ? d.a : d.b


text:
f = open("cities.txt", "r")

#a = [split(s, ":")[2] for s in filter(x -> length(x) > 1, a)]
#a = [Int(strip(s)) for s in a]
sm = 0
for s in eachline(f)
  st = match(r"(\d+)", s)
  if st != nothing
    sm += round(Int, parse(st.match))
  end
end
println(sm)

close(f)


text:
function variable_analysis_ex(ex, syms::SymbolTable)
    use!(x) = usevar!(syms, x)
    allocate!(x) = allocatevar!(syms, x)
    deallocate!(x) = deallocatevar!(syms, x)
    @smatch ex begin
        :($x[$key] ← $val) || :($x[$key] → $val) => (use!(x); use!(key); use!(val))
        :($x ← $val) => allocate!(x)
        :($x → $val) => deallocate!(x)
        :($x ↔ $y) => swapvars!(syms, x, y)
        :($a += $f($(b...))) || :($a -= $f($(b...))) ||
        :($a *= $f($(b...))) || :($a /= $f($(b...))) ||
        :($a .+= $f($(b...))) || :($a .-= $f($(b...))) ||
        :($a .*= $f($(b...))) || :($a ./= $f($(b...))) ||
        :($a .+= $f.($(b...))) || :($a .-= $f.($(b...))) ||
        :($a .*= $f.($(b...))) || :($a ./= $f.($(b...))) => begin
            ex.args[1] = render_arg(a)
            b .= render_arg.(b)
            use!(a)
            use!.(b)
            check_args(Any[a, b...])
        end
        :($a ⊻= $f($(b...))) || :($a .⊻= $f($(b...)))  || :($a .⊻= $f.($(b...))) => begin
            ex.args[1] = render_arg(a)
            b .= render_arg.(b)
            use!(a)
            use!(b)
        end
        :($a ⊻= $b || $c) || :($a ⊻= $b && $c) => begin
            ex.args[1] = render_arg(a)
            ex.args[2].args .= render_arg.(ex.args[2].args)
            use!(a)
            use!(b)
        end
        Expr(:if, _...) => variable_analysis_if(ex, syms)
        :(while $condition; $(body...); end) => begin
            julia_usevar!(syms, condition)
            localsyms = SymbolTable(Symbol[], copy(syms.deallocated), Symbol[])
            variable_analysis_ex.(body, Ref(localsyms))
            checksyms(localsyms)
        end
        :(begin $(body...) end) => begin
            variable_analysis_ex.(body, Ref(syms))
        end
        # TODO: allow ommit step.
        :(for $i=$range; $(body...); end) => begin
            julia_usevar!(syms, range)
            localsyms = SymbolTable(Symbol[], copy(syms.deallocated), Symbol[])
            variable_analysis_ex.(body, Ref(localsyms))
            checksyms(localsyms)
            ex
        end
        :(@safe $line $subex) => julia_usevar!(syms, subex)
        :(@cuda $line $(args...)) => variable_analysis_ex(args[end], syms)
        :(@launchkernel $line $(args...)) => variable_analysis_ex(args[end], syms)
        :(@inbounds $line $subex) => variable_analysis_ex(subex, syms)
        :(@simd $line $subex) => variable_analysis_ex(subex, syms)
        :(@threads $line $subex) => variable_analysis_ex(subex, syms)
        :(@avx $line $subex) => variable_analysis_ex(subex, syms)
        :(@invcheckoff $line $subex) => variable_analysis_ex(subex, syms)
        # 1. precompile to expand macros
        # 2. get dual expression
        # 3. precompile to analyze vaiables
        :($f($(args...))) => begin
            args .= render_arg.(args)
            check_args(args)
            use!.(args)
        end
        :($f.($(args...))) => begin
            args .= render_arg.(args)
            check_args(args)
            use!.(args)
        end
        :(nothing) => nothing
        ::LineNumberNode => nothing
        ::Nothing => nothing
        _ => error("unsupported statement: $ex")
    end
end

function variable_analysis_if(ex, exsyms)
    syms = copy(exsyms)
    julia_usevar!(exsyms, ex.args[1])
    variable_analysis_ex.(ex.args[2].args, Ref(exsyms))
    checksyms(exsyms, syms)
    if length(ex.args) == 3
        if ex.args[3].head == :elseif
            variable_analysis_if(ex.args[3], exsyms)
        elseif ex.args[3].head == :block
            syms = copy(exsyms)
            variable_analysis_ex.(ex.args[3].args, Ref(exsyms))
            checksyms(exsyms, syms)
        else
            error("unknown statement following `if` $ex.")
        end
    end
end

usevar!(syms::SymbolTable, arg) = @smatch arg begin
    ::Number || ::String => nothing
    ::Symbol => _isconst(arg) || operate!(syms, arg)
    :(@skip! $line $x) => julia_usevar!(syms, x)
    :($x.$k) => usevar!(syms, x)
    :($a |> subarray($(ranges...))) => (usevar!(syms, a); julia_usevar!.(Ref(syms), ranges))
    :($x |> tget($f)) || :($x |> $f) || :($x .|> $f) || :($x::$f) => (usevar!(syms, x); julia_usevar!(syms, f))
    :($x') || :(-$x) => usevar!(syms, x)
    :($t{$(p...)}($(args...))) => begin
        usevar!(syms, t)
        usevar!.(Ref(syms), p)
        usevar!.(Ref(syms), args)
    end
    :($a[$(x...)]) => begin
        usevar!(syms, a)
        usevar!.(Ref(syms), x)
    end
    :(($(args...),)) => usevar!.(Ref(syms), args)
    _ => julia_usevar!(syms, arg)
end

julia_usevar!(syms::SymbolTable, ex) = @smatch ex begin
    ::Symbol => _isconst(ex) || operate!(syms, ex)
    :($a:$b:$c) => julia_usevar!.(Ref(syms), [a, b, c])
    :($a:$c) => julia_usevar!.(Ref(syms), [a, c])
    :($a && $b) || :($a || $b) || :($a[$b]) => julia_usevar!.(Ref(syms), [a, b])
    :($a.$b) => julia_usevar!(syms, a)
    :(($(v...),)) || :(begin $(v...) end) => julia_usevar!.(Ref(syms), v)
    :($f($(v...))) || :($f[$(v...)]) => begin
        julia_usevar!(syms, f)
        julia_usevar!.(Ref(syms), v)
    end
    :($args...) => julia_usevar!(syms, args)
    Expr(:parameters, targets...) => julia_usevar!.(Ref(syms), targets)
    Expr(:kw, tar, val) => julia_usevar!(syms, val)
    ::LineNumberNode => nothing
    _ => nothing
end

# push a new variable to variable set `x`, for allocating `target`
allocatevar!(st::SymbolTable, target) = @smatch target begin
    ::Symbol => allocate!(st, target)
    :(($(tar...),)) => begin
        for t in tar
            allocatevar!(st, t)
        end
    end
    :($tar = $y) => allocatevar!(st, y)
    :($tar...) => allocatevar!(st, tar)
    :($tar::$tp) => allocatevar!(st, tar)
    Expr(:parameters, targets...) => begin
        for tar in targets
            allocatevar!(st, tar)
        end
    end
    Expr(:kw, tar, val) => begin
        allocatevar!(st, tar)
    end
    _ => _isconst(target) || error("unknown variable expression $(target)")
end

# pop a variable from variable set `x`, for deallocating `target`
deallocatevar!(st::SymbolTable, target) = @smatch target begin
    ::Symbol => deallocate!(st, target)
    :(($(tar...),)) => begin
        for t in tar
            deallocatevar!(st, t)
        end
    end
    _ => error("unknow variable expression $(target)")
end

function swapvars!(st::SymbolTable, x, y)
    e1 = isemptyvar(x)
    e2 = isemptyvar(y)
    # check assersion
    for (e, v) in ((e1, x), (e2, y))
        e && dosymbol(v) do sv
            if sv ∈ st.existing || sv ∈ st.unclassified
                throw(InvertibilityError("can not assert variable to empty: $v"))
            end
        end
    end
    if e1 && e2
    elseif e1 && !e2
        dosymbol(sx -> allocate!(st, sx), x)
        dosymbol(sy -> deallocate!(st, sy), y)
        usevar!(st, x)
    elseif !e1 && e2
        dosymbol(sx -> deallocate!(st, sx), x)
        dosymbol(sy -> allocate!(st, sy), y)
        usevar!(st, y)
    else  # both are nonempty
        sx = dosymbol(identity, x)
        sy = dosymbol(identity, y)
        if sx === nothing || sy === nothing  # e.g. x.y ↔ k.c
            usevar!(st, x)
            usevar!(st, y)
        elseif sx isa Symbol && sy isa Symbol  # e.g. x ↔ y
            swapsyms!(st, sx, sy)
        elseif sx isa Vector && sy isa Vector  # e.g. (x, y) ↔ (a, b)
            @assert length(sx) == length(sy)
            swapsyms!.(Ref(st), sx, sy)
        elseif sx isa Vector && sy isa Symbol  # e.g. (x, y) ↔ args
            swapsyms_asymetric!(st, sx, sy)
        elseif sx isa Symbol && sy isa Vector  # e.g. args ↔ (x, y)
            swapsyms_asymetric!(st, sy, sx)
        end
    end
end
isemptyvar(ex) = @smatch ex begin
    :($x[end+1]) => true
    :($x::∅) => true
    _ => false
end
dosymbol(f, ex) = @smatch ex begin
    x::Symbol => f(x)
    :(@fields $line $sym) => dosymbol(f, sym)
    :($x::$T) => dosymbol(f, x)
    :(($(args...),)) => dosymbol.(Ref(f), args)
    _ => nothing
end

_isconst(x) = @smatch x begin
    ::Symbol => x ∈ Symbol[:im, :π, :Float64, :Float32, :Int, :Int64, :Int32, :Bool, :UInt8, :String, :Char, :ComplexF64, :ComplexF32, :(:), :end, :nothing]
    ::QuoteNode || ::Bool || ::Char || ::Number || ::String => true
    :($f($(args...))) => all(_isconst, args)
    :(@const $line $ex) => true
    _ => false
end

# avoid share read/write
function check_args(args)
    args_kernel = []
    for i=1:length(args)
        out = memkernel(args[i])
        if out isa Vector
            for o in out
                if o !== nothing
                    push!(args_kernel, o)
                end
            end
        elseif out !== nothing
            push!(args_kernel, out)
        end
    end
    # error on shared read or shared write.
    for i=1:length(args_kernel)
        for j in i+1:length(args_kernel)
            if args_kernel[i] == args_kernel[j]
                throw(InvertibilityError("$i-th argument and $j-th argument shares the same memory $(args_kernel[i]), shared read and shared write are not allowed!"))
            end
        end
    end
end

# Returns the memory `identifier`, it is used to avoid shared read/write.
memkernel(ex) = @smatch ex begin
    ::Symbol => ex
    :(@const $line $x) => memkernel(x)
    :($a |> subarray($(inds...))) || :($a[$(inds...)]) => :($(memkernel(a))[$(inds...)])
    :($x.$y) => :($(memkernel(x)).$y)
    :($a |> tget($x)) => :($(memkernel(a))[$x])
    :($x |> $f) || :($x .|> $f) || :($x') || :(-$x) || :($x...) => memkernel(x)
    :($t{$(p...)}($(args...))) || :(($(args...),)) => memkernel.(args)
    _ => nothing   # Julia scope, including `@skip!`, `f(x)` et. al.
end

# Modify the argument, e.g. `x.[1,3:5]` is rendered as `x |> subarray(1,3:5)`.
render_arg(ex) = @smatch ex begin
    ::Symbol => ex
    :(@skip! $line $x) => ex
    :(@const $line $x) => Expr(:macrocall, Symbol("@const"), line, render_arg(x))
    :($a.[$(inds...)]) => :($(render_arg(a)) |> subarray($(inds...)))
    :($a |> subarray($(inds...))) => :($(render_arg(a)) |> subarray($(inds...)))
    :($a[$(inds...)]) => :($(render_arg(a))[$(inds...)])
    :($x.$y) => :($(render_arg(x)).$y)
    :($a |> tget($x)) => :($(render_arg(a)) |> tget($x))
    :($x |> $f) => :($(render_arg(x)) |> $f)
    :($x .|> $f) => :($(render_arg(x)) .|> $f)
    :($x') => :($(render_arg(x))')
    :(-$x) => :(-$(render_arg(x)))
    :($ag...) => :($(render_arg(ag))...)
    :($t{$(p...)}($(args...))) => :($t{($p...)}($(render_arg.(args)...)))
    :(($(args...),)) => :(($(render_arg.(args)...),))
    _ => ex   # Julia scope, including `@skip!`, `f(x)` et. al.
end


text:
using TomoForward
using Images
using Plots
using XfromProjections
using StaticArrays
using PyCall
using Logging
using Random
using Printf

#include("../tv_primaldual_flow.jl")
#include("../optical_flow.jl")
include("../../phantoms/sperm_phantom.jl")

H,W = 609, 609
detmin, detmax = -38.0, 38.0
function radon_operator(height, width, detcount, angles)
    proj_geom = ProjGeom(0.5, detcount, angles)
    A = fp_op_parallel2d_line(proj_geom, height, width, detmin,detmax, detmin,detmax)
    return A
end
grid = collect(detmin:0.1:detmax)
r(s) = 1.0

images, tracks = get_sperm_phantom(301,r,grid)
frames = images[:,:,collect(1:31)]
@info size(frames)

w_tv = 0.1
w_flow  = 0.05
c=1.1

detcount = Int(floor(H*1.4))

cwd = @__DIR__
path = normpath(joinpath(@__DIR__, "result"))
!isdir(path) && mkdir(path)
cd(path)

#4 angles evenly spaced between 0 and pi
#angles = collect(range(0.0,π,length=5))[1:end-1]#angles = rand(0.0:0.001:π, 4)#collect(range(π/2,3*π/2,length=ang_nr+1))[1:end-1]
#@info angles
nangles = 3
a =rand(0.0:0.001:π, nangles*size(frames)[3]) #range(0.0,4*pi+pi/4, length=nangles*size(frames)[3])#
angles = map(t -> a[(t-1)*nangles+1:(t-1)*nangles+nangles], 1:size(frames)[3])
@info angles
#match size of input image (generating data)
As = map(t -> radon_operator(size(frames[:,:,1],1),size(frames[:,:,1],2),detcount, angles[t]),1:size(frames)[3])
bsex = zeros(size(As[1])[1],size(frames)[3])

map(t -> bsex[:,t] = As[t]*vec(frames[:,:,t]), 1:size(frames)[3])
@info "adding gaussion noise at level 0.01"
rho = 0.01
e = randn(size(bsex));
e = rho*norm(bsex)*e/norm(e);
bs = bsex + e;

u0s = zeros(H,W,size(frames)[3])
mask = ones(size(u0s)...)

As  = map(t -> radon_operator(H,W, detcount, angles[t]),1:size(u0s)[3])

@info "Reconstruction using tv regularization frame by frame"
niter=450
us_tv = zeros(H,W,size(frames)[3])
for t = 1:size(frames)[3]
    A = As[t]
    p = bs[:,t]
    u0 = u0s[:,:,t]
    us_tv[:,:,t] .= recon2d_tv_primaldual!(us_tv[:,:,t], A, p, niter, w_tv, c)
end
p1 = heatmap(grid, grid, Gray.(us_tv[:,:,1]), yflip=false, aspect_ratio=:equal, framestyle=:none, legend=false)
plot!(tracks[271][:,1], tracks[271][:,2], aspect_ratio=:equal, linewidth=5)
p2 = heatmap(grid, grid, Gray.(us_tv[:,:,2]), yflip=false, aspect_ratio=:equal, framestyle=:none, legend=false)
plot!(tracks[272][:,1], tracks[272][:,2], aspect_ratio=:equal, linewidth=5)
p3 = heatmap(grid, grid, Gray.(us_tv[:,:,3]), yflip=false, aspect_ratio=:equal, framestyle=:none, legend=false)
plot!(tracks[273][:,1], tracks[273][:,2], aspect_ratio=:equal, linewidth=5)
p4 = heatmap(grid, grid, Gray.(us_tv[:,:,4]), yflip=false, aspect_ratio=:equal, framestyle=:none, legend=false)
plot!(tracks[274][:,1], tracks[274][:,2], aspect_ratio=:equal, linewidth=5)
l = @layout [a b c d]
plot(p1, p2, p3, p4, layout = l, size=(2000,600), linewidth=5)
savefig(@sprintf "result_all_tv_2_%d" length(angles))

@info "Reconstructing using joint motion estimation and reconstruction"
niter1=100
niter2=500
u0s = deepcopy(us_tv)
us_flow = recon2d_tv_primaldual_flow(As, bs, u0s, niter1, niter2, w_tv, w_flow, c, mask, 0.012)

p1 = heatmap(grid, grid, Gray.(us_flow[:,:,1]), yflip=false, aspect_ratio=:equal, framestyle=:none, legend=false)
plot!(tracks[271][:,1], tracks[271][:,2], aspect_ratio=:equal, linewidth=5)
p2 = heatmap(grid, grid, Gray.(us_flow[:,:,2]), yflip=false, aspect_ratio=:equal, framestyle=:none, legend=false)
plot!(tracks[272][:,1], tracks[272][:,2], aspect_ratio=:equal, linewidth=5)
p3 = heatmap(grid, grid, Gray.(us_flow[:,:,3]), yflip=false, aspect_ratio=:equal, framestyle=:none, legend=false)
plot!(tracks[273][:,1], tracks[273][:,2], aspect_ratio=:equal, linewidth=5)
p4 = heatmap(grid, grid, Gray.(us_flow[:,:,4]), yflip=false, aspect_ratio=:equal, framestyle=:none, legend=false)
plot!(tracks[274][:,1], tracks[274][:,2], aspect_ratio=:equal, linewidth=5)
l = @layout [a b c d]
plot(p1, p2, p3, p4, layout = l, size=(2000,600), linewidth=5)
savefig(@sprintf "result_all_flow_2_%d" length(angles))

@info "Preparing results in human readable format"
anim = @animate for t=1:size(frames)[3]
    l = @layout [a b c]
    p1 = heatmap(grid, grid, Gray.(frames[:,:,t]),yflip=false, aspect_ratio=:equal, framestyle=:none, legend=false, title="Ground Truth")
    p2 = heatmap(grid, grid, Gray.(us_flow[:,:,t]),yflip=false, aspect_ratio=:equal, framestyle=:none, legend=false, title="Flow")
    plot!(tracks[t+270][:,1], tracks[t+270][:,2], aspect_ratio=:equal, linewidth=1)
    p3 = heatmap(grid, grid, Gray.(us_tv[:,:,t]),yflip=false, aspect_ratio=:equal, framestyle=:none, legend=false, title="TV")
    plot!(tracks[t+270][:,1], tracks[t+270][:,2], aspect_ratio=:equal, linewidth=1)
    plot(p1, p2, p3, layout = l)
end

cwd = @__DIR__
path = normpath(joinpath(@__DIR__, "result"))
!isdir(path) && mkdir(path)
cd(path)
gif(anim, "reconstruction_flow_tv_2.gif", fps = 1)
cd(cwd)

max_stars_repo_name: JuliaTomo/DynamicTomographyWithFlow

max_issues_repo_name: JuliaTomo/DynamicTomographyWithFlow

max_forks_repo_name: JuliaTomo/DynamicTomographyWithFlow


text:
using DataDeps

register(DataDep(
  "horse-colic",
  "http://archive.ics.uci.edu/ml/datasets/Horse+Colic",
  "http://archive.ics.uci.edu/ml/machine-learning-databases/horse-colic/horse-colic.data",
  "270db486c78796143ffa84b4808118b3d4b96d1bd0f70bf22b3bb0bfca21e547",
  post_fetch_method=(path -> begin
    UCIData.process_dataset(path,
        target_index=24,
        feature_indices=1:22,
        categoric_indices=[7:9; 11:12; 15; 17:18; 21],
        delim=' ',
    )
  end),
))

max_stars_repo_name: JackDunnNZ/uci-data

max_issues_repo_name: JackDunnNZ/uci-data

max_forks_repo_name: JackDunnNZ/uci-data


text:
function CND(X)

  y = abs(X)

  if y > 37
    cum = 0
  else
    Exponential = exp(-y^2/2)

    if y < 7.07106781186547
      SumA = 0.0352624965998911 * y + 0.700383064443688
      SumA = SumA * y + 6.37396220353165
      SumA = SumA * y + 33.912866078383
      SumA = SumA * y + 112.079291497871
      SumA = SumA * y + 221.213596169931
      SumA = SumA * y + 220.206867912376
      SumB = 0.0883883476483184 * y + 1.75566716318264
      SumB = SumB * y + 16.064177579207
      SumB = SumB * y + 86.7807322029461
      SumB = SumB * y + 296.564248779674
      SumB = SumB * y + 637.333633378831
      SumB = SumB * y + 793.826512519948
      SumB = SumB * y + 440.413735824752
      cum = Exponential * SumA / SumB
    else
      SumA = y + 0.65
      SumA = y + 4 / SumA
      SumA = y + 3 / SumA
      SumA = y + 2 / SumA
      SumA = y + 1 / SumA
      cum = Exponential / (SumA * 2.506628274631)

      end
  end

  if X > 0
    cum = 1 - cum
  else
          cum
  end

end


text:
module JSONL

import JSON3,
    Mmap,
    Tables

export readfile,
    writefile

include("helpers.jl")

"""
    readfile(file::AbstractString; kwargs...) => Vector{JSON3.Object}

Read (parts of) a JSONLines file.

* `file`: Path to JSONLines file
* Keyword Arguments:
    * `structtype = nothing`: StructType passed to JSON3.read for each row of the file
    * `nrows = nothing`: Number of rows to load
    * `skip = nothing`: Number of rows to skip before loading
    * `usemmap::Bool = (nrows !== nothing || skip !=nothing)`: Memory map file (required for nrows and skip)
"""
function readfile(file; structtype = nothing, nrows = nothing, skip = nothing, usemmap::Bool = (nrows !== nothing || skip !== nothing))
    ff = getfile(file, nrows, skip, usemmap)
    length(ff) == 0 && return JSON3.Object[]
    if isnothing(structtype)
        rows = JSON3.read.(lstrip.(String.(ff)))
    else
        rows = JSON3.read.(lstrip.(String.(ff)), structtype)
    end
    return rows
end

"""
    writefile(file, data, mode = "w")

Write `data` to `file` in the JSONLines format.

* `file`: Path to target JSONLines file
* `data`: `Tables.jl` compatible data source
* `mode = w`: Mode to open the file in [see I/O and Network](https://docs.julialang.org/en/v1/base/io-network/)
"""
function writefile(file, data, mode = "w")
	if !Tables.istable(data)
		throw(ArgumentError("data needs to be compatible with the Tables interface"))
	end
	fi = open(file, mode)
	for row in Tables.rowtable(data)
		writerow(fi, row)
	end
	close(fi)
end

end # Module


text:
function clear( a::Array{T, N} )::Array{T, N} where {T, N}
    return Array{T, N}(undef, size(a) .- size(a))
end

include("RedBlack/RedBlack.jl")
include("GeometricalWeights/GeometricalWeights.jl")
include("PotentialSimulationSetups/PotentialSimulationSetups.jl")

include("ChargeDensity.jl")
include("DielectricDistribution.jl")
include("ElectricPotential.jl")
include("WeightingPotential.jl")

include("SimulationAlgorithms/Cylindrical.jl")
include("SimulationAlgorithms/Cartesian3D.jl")
include("ConvergenceAndRefinement.jl")


text:
"""
    μ = lowerboundSSV(M,blkstruct)

Compute a lower bound on the (largest) structured singular value (SSV, also μ, also mu) of a complex matrix for a provided structured.

"""
function lowerboundSSV(M::Matrix{S},B::BlockStructure{T}) where {S<:Number,T<:Int}
    if B.real_structured = 0 && B.complex_structured == 0
        μ = opnorm(M)   # the largest singular value
    elseif B.complex_structured_blocks == 0 && B.complex_unstructured_blocks == 0
        μ = 0.0
    else
        μ = 0.0
    end
    return μ
end

max_stars_repo_name: hurak/StructuredSingularValues

max_issues_repo_name: hurak/StructuredSingularValues

max_forks_repo_name: hurak/StructuredSingularValues


text:
using Oceananigans.Architectures: CPU
using Oceananigans.Fields: AbstractField, ReducedField

import Oceananigans.Fields: set!

const CubedSphereCPUField = CubedSphereField{X, Y, Z, <:CPU} where {X, Y, Z}
const CubedSphereGPUField = CubedSphereField{X, Y, Z, <:GPU} where {X, Y, Z}

const CubedSphereCPUReducedField = CubedSphereReducedField{X, Y, Z, <:CPU} where {X, Y, Z}
const CubedSphereGPUReducedField = CubedSphereReducedField{X, Y, Z, <:GPU} where {X, Y, Z}

const CubedSphereCPUFields = Union{CubedSphereCPUField, CubedSphereCPUReducedField}
const CubedSphereGPUFields = Union{CubedSphereGPUField, CubedSphereGPUReducedField}

# We need to define the function once for CPU fields then again for GPU fields to avoid the method
# ambiguity with Fields.set!.

function set!(u::CubedSphereCPUFields , v::CubedSphereCPUFields)
    for (u_face, v_face) in zip(faces(u), faces(v))
        @. u_face.data.parent = v_face.data.parent
    end
    return nothing
end

function set!(u::CubedSphereGPUField , v::CubedSphereGPUField)
    for (u_face, v_face) in zip(faces(u), faces(v))
        @. u_face.data.parent = v_face.data.parent
    end
    return nothing
end

set!(field::CubedSphereCPUFields, f::Function) = [set_face_field!(field_face, f) for field_face in faces(field)]
set!(field::CubedSphereCPUFields, f::Number) = [set_face_field!(field_face, f) for field_face in faces(field)]
set!(field::CubedSphereGPUFields, f::Function) = [set_face_field!(field_face, f) for field_face in faces(field)]
set!(field::CubedSphereGPUFields, f::Number) = [set_face_field!(field_face, f) for field_face in faces(field)]

set_face_field!(field, a) = set!(field, a)


text:
?convert;
# ?parse

days = 365
convert(Float64, days)

# Esto regresa el código ascii (un entero) asociado con el caracter '1'
convert(Int64, '1')

# Esto regresa un entero encapsulado entre comillas
parse(Int64, '1')

"hola"^1000

m, n = 1, 1
"$m + $n = $(m + n)"

a_ray = [[0], [0]]
push!(a_ray[1], 1)
push!(a_ray[2], 1)
a_ray

miagenda = Dict("Jenny" => "867-5309", "Ghostbusters" => "555-2368")

miagenda["Emergency"] = 911
#=

Julia infiere que "miagenda" toma ambos llaves
y valores del tipo "String". Podemos ver que miagenda
es un Dict{String,String} con 2 entradas. Esto significa que
Julia no va a aceptar enteros como valores en miagenda.

=#

agenda_flexible = Dict("Jenny" => "867-5309", "Ghostbusters" => 5552368)

flexible_phonebook["Emergency"] = 911

#=

Julia infiere que miagenda_flexible toma valores del tipo
Any. A diferencia de miagenda, miagenda_flexible es un
Dict{String,Any} con 2 entradas.

Para evitar esto, podemos inicializar miagenda a un
diccionario vacío y agregamos entradas después. O podemos
decirle a Julia explícitamente que queremos un diccionario
que acepte objectos del tipo Any como valores.

¡Ve los ejemplos!

=#

miagenda = Dict()

miagenda = Dict{String, Any}("Jenny" => "867-5309", "Ghostbusters" => "555-2368")

squares = Dict()
iterable = range(1, 2, 50)
#  Otra opción
# iterable = 1:2:99
for key in iterable
    squares[key] = "$(key^2)"
    squares[key + 1] = (key + 1)^2
end

@show squares[10]
@show squares[11]

A = fill(0, (10, 10))
for i in 1:10
    A[i] = i
end
A
# ¡Julia usa order columnar! Los primeros 10 elementos van a poblar la primera columna.

N = 16
if (N % 3 == 0) & (N % 5 == 0)
    println("FizzBuzz")
else
    if (N % 3 == 0)
        println("Fizz")
    else
        if (N % 5 == 0)
            println("Buzz")
        else
            println(N)
        end
    end
end

((N % 3 == 0) & (N % 5 == 0)) ? println("FizzBuzz") : ((N % 3 == 0) ? println("Fizz") : ((N % 5 == 0) ? println("Buzz") : println(N)))

fill(0, (3, 3)) .+ 1

f(x) = x^2
v = [1, 2, 3]
map(f, v)

caesar(input_string, shift) = map(x -> x + shift, input_string)

#Pkg.add("Primes")
using Primes
maximum(primes(1000000))

using Plots
pyplot()
x = 1:10
y = x .^ 2
plot(x, y)

import Base: +
+(x::String, y::Int) = map(x -> x + y, x)

"Gr#qrw#phggoh#lq#wkh#diidluv#ri#gudjrqv#iru#|rx#duh#fuxqfk|#dqg#wdvwh#jrrg#zlwk#nhwfkxs1" + -3

circshift(A, (0, 4))

v = [1, 2, 3]
v * v'

v' * v

A =
[
 140   97   74  168  131
  97  106   89  131   36
  74   89  152  144   71
 168  131  144   54  142
 131   36   71  142   36
]

eigdec = eigfact(A)
eigdec[:values]

Diagonal(eigdec[:values])

F = hessfact(A)

isapprox(A, F[:Q] * F[:H] * F[:Q]')

max_stars_repo_name: miguelraz/JuliaTutorials

max_issues_repo_name: miguelraz/JuliaTutorials

max_forks_repo_name: miguelraz/JuliaTutorials


text:
function TricubicInterpolatingFunction(arg0::Vector{jdouble}, arg1::Vector{jdouble}, arg10::Vector{Vector{jdouble}}, arg2::Vector{jdouble}, arg3::Vector{Vector{jdouble}}, arg4::Vector{Vector{jdouble}}, arg5::Vector{Vector{jdouble}}, arg6::Vector{Vector{jdouble}}, arg7::Vector{Vector{jdouble}}, arg8::Vector{Vector{jdouble}}, arg9::Vector{Vector{jdouble}})
    return TricubicInterpolatingFunction((Vector{jdouble}, Vector{jdouble}, Vector{Vector{jdouble}}, Vector{jdouble}, Vector{Vector{jdouble}}, Vector{Vector{jdouble}}, Vector{Vector{jdouble}}, Vector{Vector{jdouble}}, Vector{Vector{jdouble}}, Vector{Vector{jdouble}}, Vector{Vector{jdouble}}), arg0, arg1, arg10, arg2, arg3, arg4, arg5, arg6, arg7, arg8, arg9)
end

function is_valid_point(obj::TricubicInterpolatingFunction, arg0::jdouble, arg1::jdouble, arg2::jdouble)
    return jcall(obj, "isValidPoint", jboolean, (jdouble, jdouble, jdouble), arg0, arg1, arg2)
end

function value(obj::TricubicInterpolatingFunction, arg0::jdouble, arg1::jdouble, arg2::jdouble)
    return jcall(obj, "value", jdouble, (jdouble, jdouble, jdouble), arg0, arg1, arg2)
end


text:
#univariate nadaraya-watson estimate

function localconstant(xdata::RealVector, ydata::RealVector; xeval::RealVector=xdata, kernel::Function=gaussiankernel, h::Real=bwlocalconstant(xdata,ydata,kernel))

    n=length(xdata)
    length(ydata) == n || error("length(ydata) != length(xdata)")
    w=ones(n)
    pre = zeros(length(xeval))
    for i in 1:length(xeval)
        kernel(xeval[i], xdata, h, w, n)
        pre[i] = wsum(w, ydata)/sum(w)
    end
    pre
end
localconstant(xdata::RealVector, ydata::RealVector, xeval::Real; kernel::Function = gaussiankernel, h::Real = bwlocalconstant(xdata,ydata,kernel)) = localconstant(xdata, ydata, xeval = [xeval;], kernel=kernel, h = h)

function wsumsqdiff(w::RealVector, xdata::RealVector, xeval::Real, n::Int)
    res = 0.0
    for i in 1:n
        @inbounds res += w[i]*(xdata[i]-xeval).^2
    end
    res
end
function wsumyxdiff(w::RealVector, xdata::RealVector, xeval::Real, ydata::RealVector, n::Int)
    res = 0.0
    for i in 1:n
        @inbounds res += w[i]*ydata[i]*(xeval-xdata[i])
    end
    res
end

##univariate local linear
function locallinear(xdata::RealVector, ydata::RealVector; xeval::RealVector=xdata, kernel::Function=gaussiankernel, h::Real=bwlocallinear(xdata, ydata, kernel))
    n=length(xdata)
    length(ydata) == n || error("length of ydata not the same with xdata")
    w = ones(n)
    pre = zeros(length(xeval))
    for i in 1:length(xeval)
        kernel(xeval[i], xdata, h, w, n)
        s0 = sum(w)
        s1 = s0*xeval[i] - wsum(w, xdata)
        s2 = wsumsqdiff(w, xdata, xeval[i], n)
        sy0 = wsum(w, ydata)
        sy1 = wsumyxdiff(w, xdata, xeval[i], ydata, n)
        pre[i] = (s2 * sy0 - s1 * sy1) /(s2 * s0 - s1 * s1)
    end
    pre
end
locallinear(xdata::RealVector, ydata::RealVector, xeval::Real; kernel::Function = gaussiankernel, h::Real = bwlocallinear(xdata,ydata,kernel)) = locallinear(xdata, ydata, xeval=[xeval;], kernel=kernel, h=h)

function boundit(xdata::RealVector, xeval::RealVector, kernel::Function, lb::Real, ub::Real)
    if (lb == -Inf) && (ub == Inf)
        return (xdata, xeval, kernel)
    elseif (lb > -Inf) && (ub < Inf)
        (all(lb .<= xeval .<= ub) && all(lb .<= xdata .<= ub)) || error("Your data are not in [lb,ub]")
        xeval = (xeval .- lb)./(ub - lb)
        xdata = (xdata .- lb)./(ub - lb)
        if kernel != betakernel
            warn("Kernel is set to betakernel")
            kernel = betakernel
        end
    elseif (lb > -Inf) && (ub == Inf)
        (all(xeval .>= lb) && all(xdata .>= lb)) || error("lb should be less than your data")
        xeval = xeval .- lb
        xdata = xdata .- lb
        if kernel != gammakernel
            warn("Kernel is set to gammakernel")
            kernel = gammakernel
        end
    elseif (lb == -Inf) && (ub < Inf)
        (all(xeval .<= ub) && all(xdata .<= ub)) || error("ub should be greater than your data")
        xeval = ub .- xeval
        xdata = ub .- xdata
        if kernel != gammakernel
             warn("Kernel is set to gamma kernel")
             kernel = gammakernel
        end
    end
    (xdata, xeval, kernel)
end

function npr(xdata::RealVector, ydata::RealVector; xeval::RealVector=xdata,
        reg::Function=locallinear, lb::Real=-Inf, ub::Real=Inf, kernel::Function=gaussiankernel, h::Real=-Inf)

    xdata, xeval, kernel = boundit(xdata, xeval, kernel, lb, ub)
    if h <= 0
        h = bwreg(xdata, ydata, reg, kernel)
    end
    reg(xdata, ydata, xeval=xeval, kernel=kernel, h=h)
end

# #multi-variate nadaraya-watson regression or local linear
function localconstant(xdata::RealMatrix, ydata::RealVector; kernel::Array{Function, 1}=Function[gaussiankernel for i in 1:size(xdata)[2]], xeval::RealMatrix=xdata,  h::RealVector=bwlocalconstant(xdata, ydata, kernel))

    m, p = size(xeval)
    n, p1 = size(xdata)
    if p1 != p || length(h) !=p
        error("xeval, xdata and h should have same dimension!")
    end
    if any(h.<=0.0)
        error("Bandwidth should be positive")
    end
    pre=zeros(m)

    for i=1:m
        w = ones(n)
        wtmp = ones(n)
        for j in 1:p
            kernel[j](xeval[i, j], xdata[:, j], h[j], wtmp, n)
            for k in 1:n
                w[k] = w[k] * wtmp[k]
            end
        end
        divide!(w, sum(w))
        pre[i] = wsum(w, ydata)
    end
    pre
end

lp0=localconstant
lp1=locallinear