text:
__precompile__()
"""
Module to apply coordinate transformations.
"""
module Transformations

using LazySets, ..Utils, MathematicalSystems

export transform

"""
    transform(S; [method])

Interface function that calls the respective transformation function.

### Input

- `S`      -- discrete or continuous system
- `method` -- (optional, default: `'schur'`) transformation method name; valid
              otions are:

    * `'schur'`

### Output

A tuple containing:

- transformed discrete or continuous system
- inverse transformation matrix for reverting the transformation

### Notes

The functions that are called in the background should return a the transformed
system components `A`, `X0`, and `U`, and also an inverse transformation matrix `M`.
"""
function transform(S::InitialValueProblem; method::String="schur")

    if method == "schur"
        return schur_transform(S)
    else
        error("the transformation method $method is undefined")
    end
end

"""
    schur_transform(S)

Applies a Schur transformation to a discrete or continuous system.

### Input

- `S` -- discrete or continuous system

### Output

A tuple containing:

- transformed discrete or continuous system
- inverse transformation matrix for reverting the transformation

### Algorithm

We use Julia's default `schurfact` function to compute a
[Schur decomposition](https://en.wikipedia.org/wiki/Schur_decomposition)
of the coefficients matrix ``A``.
"""
function schur_transform(S::InitialValueProblem)

    A_new, T_new = schur(full(S.A)) # full (dense) matrix is required

    # recall that for Schur matrices, inv(T) == T'
    Z_inverse = T_new.'

    # apply transformation to the initial states
    X0_new = Z_inverse * S.x0

    # apply transformation to the inputs
    U_new = Z_inverse * S.s.U

    # obtain the transformation matrix for reverting the transformation again
    T_inverse = F[:vectors]

    if S.s isa DiscreteSystem
        return (DiscreteSystem(A_new, X0_new, U_new), T_inverse)
    elseif S.s isa ContinuousSystem
        return (ContinuousSystem(A_new, X0_new, U_new), T_inverse)
    end
end

end # module


text:
### A Pluto.jl notebook ###
# v0.12.20

using Markdown
using InteractiveUtils

# ╔═╡ b9fc511c-6007-11eb-0c6b-1f6871a40710
using Pkg, DrWatson

# ╔═╡ 20d377b2-6008-11eb-364a-617b6934ecb2
begin
	@quickactivate "StatsModelComparisons"
	using StatsModelComparisons
	using StanSample, StatsFuns, StatsPlots
	using DataFrames, CSV, RDatasets
end

# ╔═╡ af6b0b20-6008-11eb-2fa1-2f61145ab7db
md"
!!! note

	This script assumes that RDatasets.jl is loaded. RDatasets.jl is not included in the dependencies of PSIS.jl (as it would require R to be installed)."

# ╔═╡ 20d43cda-6008-11eb-09f0-53489a26110d
df = RDatasets.dataset("datasets", "cars");

# ╔═╡ 20e22570-6008-11eb-1565-f13541b41861
cars_stan = "
data {
    int N;
    vector[N] speed;
    vector[N] dist;
}
parameters {
    real a;
    real b;
    real sigma;
}
transformed parameters{
    vector[N] mu;
    mu = a + b * speed;
}
model {
    a ~ normal(0, 100);
    b ~ normal(0, 10);
    sigma ~ exponential(1);
    dist ~ normal(mu, sigma)    ;
}
generated quantities {
    vector[N] log_lik;
    for (i in 1:N)
        log_lik[i] = normal_lpdf(dist[i] | mu[i], sigma);
}
";

# ╔═╡ 20e2ebb6-6008-11eb-34f5-61f1ec3d024c
begin
	cars_stan_model = SampleModel("cars.model", cars_stan)
	data = (N = size(df, 1), speed = df.Speed, dist = df.Dist)
	rc = stan_sample(cars_stan_model; data)

	if success(rc)
		cars_df = read_samples(cars_stan_model; output_format=:dataframe)
		isdefined(Main, :StatisticalRethinking) &&
            PRECIS(cars_df[:, [:a, :b, :sigma]])
	end
end

# ╔═╡ 5fc59200-6008-11eb-3e06-1d0bcdf11d7d
if success(rc)
	nt_cars = read_samples(cars_stan_model);
	log_lik = nt_cars.log_lik'
end;

# ╔═╡ 20ed768a-6008-11eb-13f4-458ca1a29592
begin
    ns, n = size(log_lik)
	lppd = reshape(logsumexp(log_lik .- log(ns); dims=1), n)
	pwaic = [var(log_lik[:, i]) for i in 1:n]
	-2(sum(lppd) - sum(pwaic))
end

# ╔═╡ efe960f0-600a-11eb-1df4-5be83899715a
begin
	waic_vec = -2(lppd - pwaic)
	sqrt(n*var(waic_vec))
end

# ╔═╡ 20f5a31e-6008-11eb-2ce1-075893273872
waic(log_lik)

# ╔═╡ 6b67e38e-6009-11eb-3d9f-afc517b7d9fb
begin
	loo, loos, pk = psisloo(log_lik)
	loo
end

# ╔═╡ 82232dfe-6009-11eb-3b32-dbc487c6e4e7
pk_qualify(pk)

# ╔═╡ 8ddf22b0-6009-11eb-08da-5198ba046628
pk_plot(pk)

# ╔═╡ Cell order:
# ╟─af6b0b20-6008-11eb-2fa1-2f61145ab7db
# ╠═b9fc511c-6007-11eb-0c6b-1f6871a40710
# ╠═20d377b2-6008-11eb-364a-617b6934ecb2
# ╠═20d43cda-6008-11eb-09f0-53489a26110d
# ╠═20e22570-6008-11eb-1565-f13541b41861
# ╠═20e2ebb6-6008-11eb-34f5-61f1ec3d024c
# ╠═5fc59200-6008-11eb-3e06-1d0bcdf11d7d
# ╠═20ed768a-6008-11eb-13f4-458ca1a29592
# ╠═efe960f0-600a-11eb-1df4-5be83899715a
# ╠═20f5a31e-6008-11eb-2ce1-075893273872
# ╠═6b67e38e-6009-11eb-3d9f-afc517b7d9fb
# ╠═82232dfe-6009-11eb-3b32-dbc487c6e4e7
# ╠═8ddf22b0-6009-11eb-08da-5198ba046628

ext: jl


text:
using JuMP
using CPLEX
include("parser.jl")

function dual_relax(MyFileName::String)
  n, s, t, S, d1, d2, p, ph, d, D = read_instance(MyFileName)

  # Create the model
  m = Model(CPLEX.Optimizer)

  ## Variables
  # x relaxe, y relaxe
  @variable(m, x[1:n,1:n] >= 0) # certains arcs n'existent pas, ils ne sont pas utilises dans les contraintes
  @variable(m, y[1:n] >= 0)
  @variable(m, t1 >= 0)
  @variable(m, t2 >= 0)
  @variable(m, z[1:n,1:n] >= 0, Int) # meme remarque
  @variable(m, zprim[1:n] >= 0, Int)

  ## Constraints
  @constraint(m,[i in 1:n], y[i] <= 1) # relaxation
  @constraint(m,[i in 1:n, j in 1:n], x[i,j] <= 1) # relaxation
  @constraint(m,sum(x[s,j] for j in 1:n if d[s,j] != 0) == 1)
  @constraint(m,sum(x[i,t] for i in 1:n if d[i,t] != 0) == 1)
  @constraint(m,[v in 1:n; v != s && v != t], sum(x[i,v] for i in 1:n if d[i,v] != 0) == sum(x[v,j] for j in 1:n if d[v,j] != 0))
  @constraint(m,[v in 1:n; v != t], y[v] == sum(x[v,j] for j in 1:n if d[v,j] != 0))
  @constraint(m,[i in 1:n; i != t], y[t] == sum(x[i,t] for i in 1:n if d[i,t] != 0))
  @constraint(m, t2*d2 + sum(p[i]*y[i] + 2*zprim[i] for i in 1:n) <= S)
  @constraint(m, [i in 1:n, j in 1:n ; d[i,j] != 0], t1 + z[i,j] >= d[i,j]*x[i,j])
  @constraint(m, [i in 1:n], t2 + zprim[i] >= ph[i]*y[i])

  ## Objective
  @objective(m, Min, d1*t1 + sum(d[i,j]*x[i,j] + D[i,j]*z[i,j] for i in 1:n, j in 1:n if d[i,j] != 0))

  #resolution
  optimize!(m)

  return JuMP.objective_value.(m) # enlever . et JuMP
end

function Trouve_min(P,Q,d)
  min1 = -1
  min2 = -1
  dist = -1
  for p in P
    for q in Q
      if d[p,q] > 0 && (dist == -1 || d[p,q] < dist)
        min1 = p
        min2 = q
        dist = d[p,q]
      end
    end
  end
  return min1,min2,dist
end

function Dijkstra(n,s,t,d)
  M = Array{Int64,2}(zeros(n,3))
  # 1:Distance à s
  # 2:Prédecesseur
  # 3:Nombre de prédecesseurs
  P = [s]
  Q = [i for i in 1:n if i!=s]
  for i in 1:n
    if i != s
      M[i,1] = -1
      M[i,2] = -1
    end
  end
  for q in Q
    if d[s,q] > 0
      M[q,1] = d[q,s]
      M[q,2] = s
      M[q,3] = 1
    end
  end
  while size(Q)[1] > 0
    p,q,dist = Trouve_min(P,Q,d)
    push!(P,q)
    j = 1
    while Q[j] != q
      j = j+1
    end
    deleteat!(Q,j)
    for r in Q
      if (d[q,r] > 0) && (M[q,1] + d[q,r] < M[r,1] || M[r,1] == -1)
        M[r,1] = M[q,1] + d[q,r]
        M[r,2] = q
        M[r,3] = M[q,3]+1
      end
    end
  end
  return M
end

function sortdistance(n,d)
  Sorted_d = Array{Int64}(zeros(0))
  Sorted_i = Array{Int64}(zeros(0))
  Sorted_j = Array{Int64}(zeros(0))
  for i in 1:n
    for j in 1:n
      k = 1
      dist = d[i,j]
      done = false
      while k <= size(Sorted_d)[1] && done == false
        if dist < Sorted_d[k]
          k += 1
        else
          done = true
        end
      end
      insert!(Sorted_d,k,dist)
      insert!(Sorted_i,k,i)
      insert!(Sorted_j,k,j)
    end
  end
  return Sorted_i,Sorted_j
end

function evaluatedist(L,Sorted_i,Sorted_j,d,D,d1)
  distance = 0
  for l in 1:(size(L)[1]-1)
    distance += d[L[l],L[l+1]]
  end
  remaining = d1
  k = size(L)[1] - 1
  acc = 1
  while k > 0 && remaining > 0
    i,j = Sorted_i[acc],Sorted_j[acc]
    for l in 1:(size(L)[1]-1)
      if L[l] == i && L[l+1] == j
        k -= 1
        if D[i,j] < remaining
          distance += d[i,j]*D[i,j]
          remaining -= D[i,j]
        else
          distance += d[i,j] * remaining
          remaining = 0
        end
      end
    end
    acc = acc + 1
  end
  return distance
end

function sortweight(n,ph)
  Sorted = []
  for i in 1:n
    k = 1
    done = false
    while k <= size(Sorted)[1] && done == false
      if ph[i] < ph[Sorted[k]]
        k += 1
      else
        done =true
      end
    end
    insert!(Sorted,k,i)
  end
  return Sorted
end

function evaluateweight(L,Sorted,ph,d2)
  k = size(L)[1]
  acc = 1
  remaining = d2
  weight = 0
  while k > 0 && remaining > 0
    i = Sorted[acc]
    if i in L
      k -= 1
      if 2 < remaining
        weight += 2*ph[i]
        remaining -= 2
      else
        weight += remaining * ph[i]
        remaining = 0
      end
    end
    acc += 1
  end
  return weight
end

function TotalValue(L,Value,Weight,Relativeweight,Order,Limit)
  total = 0
  weight = Array{Int64}(zeros(0))
  relativeweight = Array{Int64}(zeros(0))
  for l in L
    total += value[l]
    push!(weight,Weight[l])
    push!(relativeweight,Relativeweight[l])
  end
  total += evaluate(L,weight,relativeweight,Order,Limit)
  return total
end

function coeff(i,j)
  return 1
end

function heuristic(MyFileName::String)
  n, s, t, S, d1, d2, p, ph, d, D = read_instance(MyFileName)
  global_tol = 5
  instant_tol = 2
  timelimit = 60
  relax_dual_value = dual_relax(MyFileName)
#function heuristic(instant_tol,global_tol,n, s, t, S, d1, d2, p, ph, d, D,timelimit)
  Ones = Array{Int64,2}(zeros(n,n))
  for i in 1:n
    for j in 1:n
    Ones[i,j] = 1
    end
  end
  H1 = Dijkstra(n,s,t,Ones)
  H2 = Dijkstra(n,t,s,Ones)
  Length1 = Array{Int64}(zeros(n))
  for i in [1:n]
    Length1[i] = H1[i,3] + H2[i,3]
  end

  New_d = Array{Int64,2}(zeros(n,n))
  for j in 1:n
    for i in 1:n
      if d[i,j] > 0
        length = (Length1[i]+Length1[j])/2
        extra = min(floor(d1/length),coeff(i,j)*D[i,j]*d[i,j])
        New_d[i,j] = d[i,j] + extra
      end
    end
  end
  H3 = Dijkstra(n,s,t,New_d)
  H4 = Dijkstra(n,t,s,New_d)
  Length2 = Array{Int64}(zeros(n))
  for i in [1:n]
    Length2[i] = H3[i,3] + H4[i,3]
  end

  Sortedweight = sortweight(n,p)
  println(Sortedweight)
  I,J = sortdistance(n,d)
  pmin = p[Sortedweight[n]]
  pmax = p[Sortedweight[1]]
  P_lim = [S - H2[i,3]*pmin for i in 1:n]
  Distance = n * d[I[1],J[1]]
  Path = [s]
  Weight = 0

  DecisionTree = [[s,p[s],global_tol,false]]
  DecisionMatrix = [[] for i in 1:n]
  Current = s
  t0 = time()
  Started = false
  iter = 0
  while (time() - t0 < timelimit) && (Started == false || size(DecisionTree)[1] > 1)
    iter+=1
    k = size(DecisionTree)[1]
    u,tol = DecisionTree[k][1],DecisionTree[k][3]
    for v in 1:n
      if New_d[u,v] > 0 && ((H4[v,3] - H4[u,3]) <= (-1 + min(tol,instant_tol)))
        i = 1
        done = false
        while i < size(DecisionMatrix[u])[1] && done == false
          if (H4[v,3],New_d[u,v]) > (H4[DecisionMatrix[u][i],3],New_d[u,DecisionMatrix[u][i]])
            i += 1
          else
            done = true
          end
        end
        insert!(DecisionMatrix[u],i,v)
      end
    end
    println(DecisionMatrix[s])
    forward = false
    while forward == false && (k > 1 || Started == false)
      Started = true
      while size(DecisionMatrix[u])[1] == 0 && k > 1
        pop!(DecisionTree)
        k -= 1
        if k > 0
          u,tol = DecisionTree[k][1],DecisionTree[k][3]
          if size(DecisionMatrix[u])[1] > 0
            popfirst!(DecisionMatrix[u])
          end
        else
          Current = s
        end
      end
      if size(DecisionMatrix[u])[1] > 0
        v = DecisionMatrix[u][1]
        w = DecisionTree[k][2] + p[v]
        checkedweight = DecisionTree[k][4]
        if checkedweight == false && S - w < 2*d2*pmax
          checkedweight = true
          L = [DecisionTree[i][1] for i in 1:k]
          push!(L,v)
          w += evaluateweight(L,Sortedweight,ph,d2)
        end
        path = [DecisionTree[i][1] for i in 1:size(DecisionTree)[1]]
        if w > P_lim[v] || v in path
         popfirst!(DecisionMatrix[u])
        elseif v == t
          push!(path,v)
          weight = evaluateweight(path,Sortedweight,ph,d2)
          for i in path
            weight += p[i]
          end
          distance = evaluatedist(path,I,J,d,D,d1)
          if weight <= S && distance < Distance
            Path = path
            Weight = weight
            Distance = distance
          end
          popfirst!(DecisionMatrix[u])
        else
          push!(DecisionTree,[v,w,tol-(H4[v,3]-H4[u,3]+1),checkedweight])
          forward = true
          Current = v
        end
      end
    end
  end
  println(Path," ",Weight," ",Distance," ",(Distance/relax_dual_value)-1," ",iter)
  return Path,Weight,Distance,(Distance/relax_dual_value)-1
end

max_stars_repo_name: boumpteryx/ECMA

max_issues_repo_name: boumpteryx/ECMA

max_forks_repo_name: boumpteryx/ECMA


text:
function handle_callbacks!(integrator)
    discrete_callbacks = integrator.opts.callback.discrete_callbacks
    continuous_callbacks = integrator.opts.callback.continuous_callbacks
    atleast_one_callback = false

    continuous_modified = false
    discrete_modified = false
    saved_in_cb = false
    if !(typeof(continuous_callbacks) <: Tuple{})
        time, upcrossing, event_occured, event_idx, idx, counter =
            DiffEqBase.find_first_continuous_callback(integrator, continuous_callbacks...)
        if event_occured
            integrator.event_last_time = idx
            integrator.vector_event_last_time = event_idx
            continuous_modified, saved_in_cb = DiffEqBase.apply_callback!(
                integrator,
                continuous_callbacks[idx],
                time,
                upcrossing,
                event_idx,
            )
        else
            integrator.event_last_time = 0
            integrator.vector_event_last_time = 1
        end
    end
    if !(typeof(discrete_callbacks) <: Tuple{})
        discrete_modified, saved_in_cb =
            DiffEqBase.apply_discrete_callback!(integrator, discrete_callbacks...)
    end

    integrator.u_modified = continuous_modified || discrete_modified
    if integrator.u_modified
        handle_callback_modifiers!(integrator)
    end

    if !saved_in_cb
        savevalues!(integrator)
    end

    integrator.u_modified = false
end

function DiffEqBase.savevalues!(
    integrator::AbstractSundialsIntegrator,
    force_save = false,
)::Tuple{Bool, Bool}
    saved, savedexactly = false, false
    !integrator.opts.save_on && return saved, savedexactly
    uType = eltype(integrator.sol.u)
    # The call to first is an overload of Base.first implemented in DataStructures
    while !isempty(integrator.opts.saveat) &&
        integrator.tdir * first(integrator.opts.saveat) < integrator.tdir * integrator.t
        saved = true
        curt = pop!(integrator.opts.saveat)

        tmp = integrator(curt)
        save_value!(integrator.sol.u, tmp, uType, integrator.sizeu, integrator.opts.save_idxs, Val{false})
        push!(integrator.sol.t, curt)
        if integrator.opts.dense
            tmp = integrator(curt, Val{1})
            save_value!(integrator.sol.interp.du, tmp, uType, integrator.sizeu, integrator.opts.save_idxs, Val{false})
        end
    end

    if force_save || integrator.opts.save_everystep || (integrator.opts.save_everystep && (isempty(integrator.sol.t) ||
                                         (integrator.t !== integrator.sol.t[end])))
        saved = true
        save_value!(integrator.sol.u, integrator.u, uType, integrator.sizeu, integrator.opts.save_idxs)
        push!(integrator.sol.t, integrator.t)
        if integrator.opts.dense
            tmp = integrator(integrator.t, Val{1})
            save_value!(integrator.sol.interp.du, tmp, uType, integrator.sizeu, integrator.opts.save_idxs)
        end
    end
    savedexactly = !isempty(integrator.sol.t) && last(integrator.sol.t) == integrator.t
    return saved, savedexactly
end

function save_value!(
    save_array,
    val,
    ::Type{T},
    sizeu, save_idxs,
    make_copy::Type{Val{bool}} = Val{true},
) where {T <: Number, bool}
    push!(save_array, first(val))
end
function save_value!(
    save_array,
    val,
    ::Type{T},
    sizeu, save_idxs,
    make_copy::Type{Val{bool}} = Val{true},
) where {T <: Vector, bool}
    save = if save_idxs !== nothing
        val[save_idxs]
    else
        bool ? copy(val) : val
    end
    push!(save_array, save)
end
function save_value!(
    save_array,
    val,
    ::Type{T},
    sizeu, save_idxs,
    make_copy::Type{Val{bool}} = Val{true},
) where {T <: AbstractArray, bool}
    save = if save_idxs !== nothing
        reshape(val, sizeu)[save_idxs]
    else
        x = bool ? copy(val) : val
        reshape(x, sizeu)
    end
    push!(save_array, save)
end

function handle_callback_modifiers!(integrator::CVODEIntegrator)
    CVodeReInit(integrator.mem, integrator.t, integrator.u)
end

function handle_callback_modifiers!(integrator::ARKODEIntegrator)
    ARKStepReInit(integrator.mem, integrator.t, integrator.u)
end

function handle_callback_modifiers!(integrator::IDAIntegrator)
    IDAReInit(integrator.mem, integrator.t, integrator.u, integrator.du)
    DiffEqBase.initialize_dae!(integrator)
end

function DiffEqBase.add_tstop!(integrator::AbstractSundialsIntegrator, t)
    integrator.tdir * (t - integrator.t) < 0 &&
        error("Tried to add a tstop that is behind the current time. This is strictly forbidden")
    push!(integrator.opts.tstops, t)
end

function DiffEqBase.add_saveat!(integrator::AbstractSundialsIntegrator, t)
    integrator.tdir * (t - integrator.t) < 0 &&
        error("Tried to add a saveat that is behind the current time. This is strictly forbidden")
    push!(integrator.opts.saveat, t)
end

DiffEqBase.get_tmp_cache(integrator::AbstractSundialsIntegrator) = (integrator.tmp,)

@inline function DiffEqBase.u_modified!(integrator::AbstractSundialsIntegrator, bool::Bool)
    integrator.u_modified = bool
end

function DiffEqBase.terminate!(
    integrator::AbstractSundialsIntegrator,
    retcode = :Terminated,
)
    integrator.sol = DiffEqBase.solution_new_retcode(integrator.sol, retcode)
    integrator.opts.tstops.valtree = typeof(integrator.opts.tstops.valtree)()
end

@inline function DiffEqBase.get_du(integrator::CVODEIntegrator)
    integrator(integrator.t, Val{1})
end

@inline function DiffEqBase.get_du!(out, integrator::CVODEIntegrator)
    integrator(out, integrator.t, Val{1})
end

@inline function DiffEqBase.get_du(integrator::IDAIntegrator)
    reshape(integrator.du, integrator.sizedu)
end

@inline function DiffEqBase.get_du!(out, integrator::IDAIntegrator)
    out .= reshape(integrator.du, integrator.sizedu)
end

function DiffEqBase.change_t_via_interpolation!(integrator::AbstractSundialsIntegrator, t)
    integrator.t = t
    integrator(integrator.u, integrator.t)
    return nothing
end

@inline function Base.getproperty(integrator::AbstractSundialsIntegrator, sym::Symbol)
    if sym == :dt
        return integrator.t - integrator.tprev
    else
        return getfield(integrator, sym)
    end
end

DiffEqBase.reeval_internals_due_to_modification!(integrator::AbstractSundialsIntegrator) =
    nothing
DiffEqBase.reeval_internals_due_to_modification!(integrator::IDAIntegrator) =
    handle_callback_modifiers!(integrator::IDAIntegrator)

# Required for callbacks
DiffEqBase.set_proposed_dt!(i::AbstractSundialsIntegrator, dt) = nothing

DiffEqBase.initialize_dae!(integrator::AbstractSundialsIntegrator) = nothing
function DiffEqBase.initialize_dae!(integrator::IDAIntegrator)
    integrator.f(integrator.tmp, integrator.du, integrator.u, integrator.p, integrator.t)
    if any(abs.(integrator.tmp) .>= integrator.opts.reltol)
        if integrator.sol.prob.differential_vars === nothing && !integrator.alg.init_all
            error("Must supply differential_vars argument to DAEProblem constructor to use IDA initial value solver.")
        end
        if integrator.alg.init_all
            init_type = IDA_Y_INIT
        else
            init_type = IDA_YA_YDP_INIT
            integrator.flag =
                IDASetId(integrator.mem, integrator.sol.prob.differential_vars)
        end
        integrator.flag = IDACalcIC(integrator.mem, init_type, integrator.dt)
    end
end


text:
include(joinpath(@__DIR__,"..","src","open_file.jl"))

using ArgParse, LibGit2

s = ArgParseSettings()
@add_arg_table! s begin
    "--diff"
        help = "print to STDOUT the output of `diff`"
        action = :store_true
    "--two"
        help = "display both files"
        action = :store_true
    "--bw"
        help = "generate, save and display a B&W difference image for PNG and SVG files.  requires Rsvg, Cairo, and Images"
        action = :store_true
    "filter"
        help = "a regular expression describing the filenames to compare"
        default=""
end

args = parse_args(s)

if args["bw"]
    using Rsvg, Cairo, Images

    function svg2img(filename)
        r = Rsvg.handle_new_from_file(filename)
        d = Rsvg.handle_get_dimensions(r)
        cs = Cairo.CairoImageSurface(d.width,d.height,Cairo.FORMAT_ARGB32)
        c = Cairo.CairoContext(cs)
        Rsvg.handle_render_cairo(c,r)
        fout = tempname()
        Cairo.write_to_png(cs,fout)
        png = load(fout)
        rm(fout)
        png
    end
end

# delete diffedoutput/
repo = LibGit2.GitRepo(dirname(@__DIR__))
options = LibGit2.StatusOptions(flags=LibGit2.Consts.STATUS_OPT_INCLUDE_IGNORED |
                                      LibGit2.Consts.STATUS_OPT_RECURSE_IGNORED_DIRS)
status = LibGit2.GitStatus(repo, status_opts=options)
for i in 1:length(status)
    entry = status[i]
    entry.index_to_workdir == C_NULL && continue
    index_to_workdir = unsafe_load(entry.index_to_workdir)
    if index_to_workdir.status == Int(LibGit2.Consts.DELTA_IGNORED)
        filepath = unsafe_string(index_to_workdir.new_file.path)
        startswith(filepath,joinpath("test/diffed-output")) || continue
        rm(filepath)
    end
end

function display_two(master,devel)
    open_file(master)
    open_file(devel)
end

# Compare with cached output
masterout = joinpath((@__DIR__), "master-output")
develout  = joinpath((@__DIR__), "devel-output")
diffedout = joinpath((@__DIR__), "diffed-output")
ndifferentfiles = 0
const creator_producer = r"(Creator|Producer)"
filter_mkdir_git(x) = !mapreduce(y->x==y,|,[".mkdir","git.log","git.status"])
filter_regex(x) = occursin(Regex(args["filter"]), x)
master_files = filter(x->filter_mkdir_git(x) && filter_regex(x), readdir(masterout))
devel_files = filter(x->filter_mkdir_git(x) && filter_regex(x), readdir(develout))
cached_notin_genned = setdiff(master_files, devel_files)
isempty(cached_notin_genned) ||
      @warn string("files in master-output/ but not in devel-output/: ", join(cached_notin_genned,", "))
genned_notin_cached = setdiff(devel_files, master_files)
isempty(genned_notin_cached) ||
      @warn string("files in devel-output/ but not in master-output/: ", join(genned_notin_cached,", "))
for file in intersect(master_files,devel_files)
    print("Comparing ", file, " ... ")
    cached = open(readlines, joinpath(masterout, file))
    genned = open(readlines, joinpath(develout, file))
    same = length(cached) == length(genned)
    if same
        lsame = Bool[cached[i] == genned[i] for i = 1:length(cached)]
        if !all(lsame)
            for idx in findall(lsame.==false)
                # Don't worry about lines that are due to
                # Creator/Producer (e.g., Cairo versions)
                if findfirst(creator_producer, cached[idx]) !== nothing
                    lsame[idx] = true
                end
            end
        end
        same = same & all(lsame)
    end
    if same
        println("same!")
    else
        global ndifferentfiles
        ndifferentfiles +=1
        println("different :(")
        if args["diff"]
            diffcmd = `diff $(joinpath(masterout, file)) $(joinpath(develout, file))`
            run(ignorestatus(diffcmd))
        end
        args["two"] && display_two(joinpath(masterout,file), joinpath(develout,file))
        if args["bw"] && (endswith(file,".svg") || endswith(file,".png"))
            wait_for_user = false
            if endswith(file,".svg")
                gimg = svg2img(joinpath(develout,file));
                cimg = svg2img(joinpath(masterout,file));
            elseif endswith(file,".png")
                gimg = load(joinpath(develout,file));
                cimg = load(joinpath(masterout,file));
            end
            if size(gimg)==size(cimg)
                dimg = convert(Matrix{Gray}, gimg.==cimg)
                if any(dimg.==0)
                    fout = joinpath(diffedout,file*".png")
                    Images.save(fout, dimg)
                    wait_for_user = true
                    open_file("$fout")
                else
                    println("files are different but PNGs are the same")
                end
            else
                wait_for_user = true
                println("PNGs are different sizes :(")
            end
        end
        args["diff"] || args["two"] || (args["bw"] &&
                (endswith(file,".svg") || endswith(file,".png")) && wait_for_user) || continue
        println("Enter 'two' to display both files, nothing to continue, or press CTRL-C to quit")
        while true
          resp = readline()
          resp=="" && break
          resp=="two" && display_two(joinpath(masterout,file), joinpath(develout,file))
        end
    end
end

result = string("# different images = ",ndifferentfiles)
if ndifferentfiles==0
  @info result
else
  @warn result
end

repo = options = status = 0;  GC.gc()   # see https://github.com/JuliaLang/julia/issues/28306


text:
using MPStates, Test, LinearAlgebra

@testset "Operations with Mpo and Mps" begin
# Define a Hubbard Hamiltonian that we will use for testing.
L = 6
# Physical dimension.
d = 2
t = 1.5
U = 1.5
J = t .* Symmetric(diagm(1 => ones(L-1)))
V = U .* diagm(1 => ones(L-1))

H = Mpo(L, d)
# Add fermionic operators.
add_ops!(H, "c+", "c", J, ferm_op="Z")
# Add interaction terms.
add_ops!(H, "n", "n", V)

# Define "W", "GHZ", and "product" states.
W = Mps(L, "W")
GHZ = Mps(L, "GHZ")
product = Mps(L, "product")

# Make a spin Hamiltonian to test.
Hs = Mpo(L, 3)
Jzz = diagm(1 => ones(L-1))
Jpm = 0.3*Symmetric(diagm(1 => ones(L-1)))
add_ops!(Hs, "Sz", "Sz", Jzz)
add_ops!(Hs, "S+", "S-", Jpm)

rtest3 = MPStates.testMps("rtest3")

@testset "expectation value of Mps" begin
    @test expected(H, product) ≈ 0. atol=1e-15
    @test expected(H, GHZ) ≈ U/2*(L-1)
    @test expected(H, W) ≈ 2t*(L-1)/L
    @test expected(Hs, rtest3) ≈ 2/3 - 1/36 - 0.3*8/3
end

@testset "expectation value of squared Mps" begin
    @test MPStates.norm_of_apply(H, product) ≈ 0. atol=1e-15
    @test MPStates.norm_of_apply(H, GHZ) ≈ U^2/2*(L-1)^2
    @test MPStates.norm_of_apply(H, W) ≈ t^2*(4L-6)/L
    @test MPStates.norm_of_apply(Hs, rtest3) ≈ norm(apply!(Hs, rtest3))
end

@testset "variance of Mpo" begin
    @test m_variance(H, product) ≈ 0. atol=1e-15
    @test m_variance(H, GHZ) ≈ U^2/4*(L-1)^2
    @test m_variance(H, W) ≈ t^2*(4L-6)/L - (2t*(L-1)/L)^2
end

@testset "apply Mpo to Mps" begin
    L = 6
    J = zeros(L, L)
    J[1, 2] = 1.
    J[4, 5] = 1.

    rtest1 = MPStates.testMps("rtest1")
    ctest2 = MPStates.testMps("ctest2")

    Op = Mpo(L, 2)
    add_ops!(Op, "c+", "c", J, ferm_op="Z")

    apply!(Op, rtest1)
    MPStates.make_left_canonical!(rtest1, false)
    @test norm(rtest1) ≈ 4/9
    @test expected(rtest1, "n", 1) ≈ 4/9
    @test expected(rtest1, "n", 2) ≈ 0. atol=1e-15
    @test expected(rtest1, "n", 3) ≈ 0. atol=1e-15
    @test expected(rtest1, "n", 4) ≈ 4/9
    @test expected(rtest1, "n", 5) ≈ 4/9*0.64
    @test expected(rtest1, "n", 6) ≈ 4/9

    cOp = Mpo(ComplexF64, L, 2)
    add_ops!(cOp, "c+", "c", convert.(ComplexF64, J), ferm_op="Z")

    apply!(cOp, ctest2)
    MPStates.make_left_canonical!(ctest2, false)
    @test norm(ctest2) ≈ 4/9
    @test expected(ctest2, "n", 1) ≈ 4/9
    @test expected(ctest2, "n", 2) ≈ 0. atol=1e-15
    @test expected(ctest2, "n", 3) ≈ 4/9
    @test expected(ctest2, "n", 4) ≈ 4/9
    @test expected(ctest2, "n", 5) ≈ 4/9
    @test expected(ctest2, "n", 6) ≈ 4/9*0.64
end
end # @testset "Operations with Mpo and Mps"


text:
#
# Initialize simulation data
#

function init(;psf=nothing,
               dcd=nothing,
               log=nothing,
               pdb=nothing,
               parfiles=nothing,
               vmd="vmd")

  #
  # Read PSF file data
  #

  if psf == nothing
    error(" At least a PSF must be provided with psf=filename.psf ")
  end
  atoms, bonds, angles, dihedrals, impropers = readpsf(psf)
  natoms = length(atoms)
  nbonds = length(bonds)
  nangles = length(angles)
  ndihedrals = length(dihedrals)
  nimpropers = length(impropers)

  # Reading coordinates from PDB file, if provided

  if pdb != nothing
    getpdbcoords!(pdb,atoms)
  else
    pdb = "none"
  end

  # Reading parameter files, if provided

  if parfiles != nothing
    readprm!(parfiles,atoms,bonds,angles,dihedrals,impropers)
  else
    parfiles = ["none"]
  end

  #
  # Reads DCD file header, returns nframes (correctly, if set) and ntotat
  #

  if dcd != nothing
    FortranDCD = FortranFile(dcd)
    nframes, read_natoms, dcdaxis = dcdheader(FortranDCD)
    if read_natoms != natoms
      error(" Number of atoms in DCD file is different from that of PSF file.")
    end
    if nframes < 1
      println(" WARNING: The DCD header does not inform the number of frames. Will read the DCD to fetch it. ")
      nframes = getnframes(FortranDCD, dcdaxis)
    end
  else
    dcd = "none"
    nframes = 0
    dcdaxis = false
    FortranDCD = FortranFile(psf)
    println(" Warning: no DCD file provided, so many functions won't work. ")
  end

  #
  # Read data from log file, if provided
  #

  if log != nothing

    logfile = strip(log)
    println(" Reading data from LOG file: ", logfile,"\n" )

    timestep = 0
    dcdfreq = 0
    dcdfirststep = 0
    nsteps = 0
    for line in eachline(logfile)
      if line[1:min(length(line),7)] == "ENERGY:"
        nsteps = nsteps + 1
      end
      try
        if line[1:14] == "Info: TIMESTEP"
          data = split(line)
          timestep = parse(Int64,data[3])
        end
      catch msg
      end
      try
        if line[1:19] == "Info: DCD FREQUENCY"
          data = split(line)
          dcdfreq = parse(Int64,data[4])
        end
      catch msg
      end
      try
        if line[1:19] == "Info: DCD FIRST STEP"
          data = split(line)
          dcdfirststep = parse(Int64,data[5])
        end
      catch msg
      end
    end
    if timestep == 0
      error("Could not find TIMESTEP in log file.")
    end

    time = Vector{Float64}(undef,nsteps)
    ts = Vector{Int64}(undef,nsteps)
    bond = Vector{Float64}(undef,nsteps)
    angle = Vector{Float64}(undef,nsteps)
    dihed = Vector{Float64}(undef,nsteps)
    imprp = Vector{Float64}(undef,nsteps)
    elect = Vector{Float64}(undef,nsteps)
    vdw = Vector{Float64}(undef,nsteps)
    boundary = Vector{Float64}(undef,nsteps)
    misc = Vector{Float64}(undef,nsteps)
    kinetic = Vector{Float64}(undef,nsteps)
    total = Vector{Float64}(undef,nsteps)
    temperature = Vector{Float64}(undef,nsteps)
    potential = Vector{Float64}(undef,nsteps)
    total3 = Vector{Float64}(undef,nsteps)
    tempavg = Vector{Float64}(undef,nsteps)
    pressure = Vector{Float64}(undef,nsteps)
    gpressure = Vector{Float64}(undef,nsteps)
    volume = Vector{Float64}(undef,nsteps)
    pressavg = Vector{Float64}(undef,nsteps)
    gpressavg = Vector{Float64}(undef,nsteps)

    nsteps = 0
    for line in eachline(logfile)
      if line[1:min(length(line),7)] == "ENERGY:"
        nsteps = nsteps + 1
        data = split(line)
        ts[nsteps] = parse(Int64,data[2])
        time[nsteps] = ts[nsteps] / ( 1.e6 * timestep ) # in ns
        bond[nsteps] = parse(Float64,data[3])
        angle[nsteps] = parse(Float64,data[4])
        dihed[nsteps] = parse(Float64,data[5])
        imprp[nsteps] = parse(Float64,data[6])
        elect[nsteps] = parse(Float64,data[7])
        vdw[nsteps] = parse(Float64,data[8])
        boundary[nsteps] = parse(Float64,data[9])
        misc[nsteps] = parse(Float64,data[10])
        kinetic[nsteps] = parse(Float64,data[11])
        total[nsteps] = parse(Float64,data[12])
        temperature[nsteps] = parse(Float64,data[13])
        potential[nsteps] = parse(Float64,data[14])
        total3[nsteps] = parse(Float64,data[15])
        tempavg[nsteps] = parse(Float64,data[16])
        pressure[nsteps] = parse(Float64,data[17])
        gpressure[nsteps] = parse(Float64,data[18])
        volume[nsteps] = parse(Float64,data[19])
        pressavg[nsteps] = parse(Float64,data[20])
        gpressavg[nsteps] = parse(Float64,data[21])
      end
    end
    logdata = LogData(timestep,dcdfreq,dcdfirststep,time,ts,bond,angle,dihed,imprp,
                      elect,vdw,boundary,misc,kinetic,total,temperature,potential,total3,tempavg,
                      pressure,gpressure,volume,pressavg,gpressavg)
  else
    logfile = "none"
    logdata = LogData()
  end

  simulation = Simulation(psf,pdb,dcd,
                          natoms,atoms,
                          nbonds,bonds,
                          nangles,angles,
                          ndihedrals,dihedrals,
                          nimpropers,impropers,
                          nframes,dcdaxis,
                          vmd,
                          FortranDCD,
                          logfile,logdata,parfiles)

  return simulation
end


text:
using ColorTypes

"""
utilities for dealing data structs like FullScrollableDat or SingleSliceDat
"""
module StructsManag
using  Logging, Setfield,   ..ForDisplayStructs,   ..DataStructs, Rocket
export getThreeDims,addToforUndoVector,cartTwoToThree,getHeightToWidthRatio,threeToTwoDimm,modSlice!, threeToTwoDimm,modifySliceFull!,getSlicesNumber,getMainVerticies

```@doc
given two dim dat it sets points in given coordinates in given slice to given value
coords - coordinates in a plane of chosen slice to modify
value - value to set for given points
return reference to modified slice
```
function modSlice!(data::TwoDimRawDat{T}
                ,coords::Vector{CartesianIndex{2}}
                ,value::T ) where {T}
   data.dat[coords].=value
   data
end#modSlice

```@doc
gives access to the slice of intrest - way of slicing is defined at the begining
typ - type of data
slice - slice we want to access
sliceDim - on the basis of what dimension we are slicing
return 2 dimensional array  wrapper -TwoDimRawDat  object representing slice of given 3 dimensional array
!! important returned TwoDimRawDat holds view to the original 3 dimensional data
```
function threeToTwoDimm(typ::Type{T}
                ,sliceInner::Int
                ,sliceDim::Int
                ,threedimDat::ThreeDimRawDat{T})::TwoDimRawDat{T} where {T}

maxSlice =size(threedimDat.dat)[sliceDim]
slice= sliceInner
if(sliceInner>maxSlice) slice=maxSlice   end
               return TwoDimRawDat{T}(typ,threedimDat.name,selectdim(threedimDat.dat, sliceDim, slice)   )
end#ThreeToTwoDimm

modifySliceFull!Str= """
modifies given slice in given coordinates of given data - queried by name
data - full data we work on and modify
coords - coordinates in a plane of chosen slice to modify (so list of x and y coords)
value - value to set for given points
return reference to modified slice
"""
@doc modifySliceFull!Str
function modifySliceFull!(data::FullScrollableDat
                        ,slice::Int
                        ,coords::Vector{CartesianIndex{2}}
                        ,name::String
                        ,value)

     threeDimDat=data.nameIndexes[name] |>
     (ind)-> data.dataToScroll[ind]
     if(typeof(value)!=threeDimDat.type )  throw(DomainError(value, "supplied value should be of compatible type - $(threeDimDat.type )"))  end #if

     return threeToTwoDimm(threeDimDat.type,slice,data.dimensionToScroll,threeDimDat  ) |>
     (twoDimDat)-> modSlice!(twoDimDat,coords,value)
end#modifySliceFull!

```@doc
Return number of slices present in on slice data - takes into account slices dimensions
```
function getSlicesNumber(data::FullScrollableDat)::Int32
return Int32(size(data.dataToScroll[1].dat)[data.dimensionToScroll])
end#getSlicesNumber

```@doc
Based on DataToScrollDims it will enrich passed CalcDimsStruct texture width, height and  heightToWithRatio
based on data passed from DataToScrollDims
```
function getHeightToWidthRatio( calcDim::CalcDimsStruct ,dataToScrollDims::DataToScrollDims)::CalcDimsStruct
  toSelect= filter(it-> it!=dataToScrollDims.dimensionToScroll , [1,2,3] )# will be used to get texture width and height

  return setproperties(calcDim, (imageTextureWidth= dataToScrollDims.imageSize[toSelect[1]]
                              ,imageTextureHeight =dataToScrollDims.imageSize[toSelect[2]]
                              ,heightToWithRatio =  dataToScrollDims.voxelSize[toSelect[1]]/dataToScrollDims.voxelSize[toSelect[2]]
                              ,textTextureZeros= calcDim.textTextureZeros
  ))
end#getHeightToWidthRatio

```@doc
Based on DataToScrollDims ,2 dim cartesian coordinate and  slice number it gives 3 dimensional coordinate of mouse position
```
function cartTwoToThree(dataToScrollDims::DataToScrollDims
                           ,sliceNumber::Int
                           ,cartIn::CartesianIndex{2})::CartesianIndex{3}

 toSelect= filter(it-> it!=dataToScrollDims.dimensionToScroll , [1,2,3] )# will be used to get texture width and height
resArr= [1,1,1]

resArr[toSelect[1]] = cartIn[1]
resArr[toSelect[2]] = cartIn[2]
resArr[dataToScrollDims.dimensionToScroll]=Int64(sliceNumber)

  return CartesianIndex(resArr[1],resArr[2],resArr[3]  )
end#cartTwoToThree

```@doc
Given function and actor it passes the function to forUndoVector -
   in case the length of the vector is too big the last element woill be removed
```
function addToforUndoVector(actor::SyncActor{Any, ActorWithOpenGlObjects}
                           ,fun)

   push!(actor.actor.forUndoVector,fun)

    if(length(actor.actor.forUndoVector) >actor.actor.maxLengthOfForUndoVector )
      popfirst!(actor.actor.forUndoVector)
    end

end#addToforUndoVector

```@doc
utility function to create series of ThreeDimRawDat from list of tuples where
first entry is String and second entry is 3 dimensional array with data
```
function getThreeDims(list )

  return map(tupl->ThreeDimRawDat{typeof(tupl[2][1])}(typeof(tupl[2][1]),tupl[1],tupl[2]) ,list)

end#getThreeDims

# parameter_type(x::TextureSpec) = parameter_type(typeof(x))

```@doc
calculates proper dimensions form main quad display on the basis of data stored in CalcDimsStruct
some of the values calculated will be needed for futher derivations for example those that will  calculate mouse positions
reurn CalcDimsStruct enriched by new data
    ```
function getMainVerticies(calcDimStruct::CalcDimsStruct)::CalcDimsStruct
    #corrections that will be added on both sides (in case of height correction top and bottom in case of width correction left and right)
    # to achieve required ratio
    widthCorr=0.0
    heightCorr=0.0

    correCtedWindowQuadHeight= calcDimStruct.avWindHeightForMain
    correCtedWindowQuadWidth = calcDimStruct.avWindWidtForMain
    isWidthToBeCorrected = false
    isHeightToBeCorrected = false

    #first one need to check weather current height to width ratio is as we want it  and if not  is it to high or to low
    if(calcDimStruct.avMainImRatio>calcDimStruct.heightToWithRatio )
        #if we have to big height to width ratio we need to reduce size of acual quad from top and bottom
        # we know that we would not need to change width  hence we will use the width to calculate the quad height
        correCtedWindowQuadHeight= calcDimStruct.heightToWithRatio * calcDimStruct.avWindWidtForMain
        isHeightToBeCorrected= true
    end# if to heigh

    if(calcDimStruct.avMainImRatio<calcDimStruct.heightToWithRatio )
        #if we have to low height to width ratio we need to reduce size of acual quad from left and right
        # we know that we would not need to change height  hence we will use height to calculate the quad height
        correCtedWindowQuadWidth = calcDimStruct.avWindHeightForMain/calcDimStruct.heightToWithRatio
        isWidthToBeCorrected= true
    end# if to wide

    # now we still need ratio of the resulting quad window size after corrections relative to  total window size
    quadToTotalHeightRatio= correCtedWindowQuadHeight/calcDimStruct.windowHeight
    quadToTotalWidthRatio= correCtedWindowQuadWidth/calcDimStruct.windowWidth

    # original ratios of available space of main image to total window dimensions
    avQuadToTotalHeightRatio= calcDimStruct.avWindHeightForMain/calcDimStruct.windowHeight
    avQuadToTotalWidthRatio= calcDimStruct.avWindWidtForMain/calcDimStruct.windowWidth
    # if those would be equal to corrected ones we would just start from the  bottom left corner and create quad from there - yet now we need to calculate corrections based on the diffrence of quantities just above and corrected ones

    if(isHeightToBeCorrected)
        heightCorr=abs(quadToTotalHeightRatio-avQuadToTotalHeightRatio)
    end # if isHeightToBeCorrected

    if(isWidthToBeCorrected)
        widthCorr=abs(quadToTotalWidthRatio-avQuadToTotalWidthRatio)
    end # if isWidthToBeCorrected

    correctedWidthForTextAccounting = (-1+ calcDimStruct.fractionOfMainIm*2)
    #as in OpenGl we start from -1 and end at 1 those ratios needs to be doubled in order to translate them in the OPEN Gl coordinate system yet we will achieve this doubling by just adding the corrections from both sides
    #hence we do not need  to multiply by 2 becose we get from -1 to 1 so total is 2
      res =  Float32.([
      # positions                  // colors           // texture coords
      correctedWidthForTextAccounting-widthCorr,  1.0-heightCorr, 0.0,   1.0, 0.0, 0.0,   1.0, 1.0,   # top right
      correctedWidthForTextAccounting-widthCorr, -1.0+heightCorr, 0.0,   0.0, 1.0, 0.0,   1.0, 0.0,   # bottom right
      -1.0+widthCorr, -1.0+heightCorr, 0.0,                0.0, 0.0, 1.0,   0.0, 0.0,   # bottom left
      -1.0+widthCorr,  1.0-heightCorr, 0.0,                1.0, 1.0, 0.0,   0.0, 1.0    # top left
    ])

    windowWidthCorr=Int32(round( (widthCorr/2)*calcDimStruct.windowWidth))
    windowHeightCorr= Int32(round((heightCorr/2)*calcDimStruct.windowHeight))

    return setproperties(calcDimStruct, (correCtedWindowQuadHeight= Int32(round(correCtedWindowQuadHeight))
                                        ,correCtedWindowQuadWidth= Int32(round(correCtedWindowQuadWidth))
                                        ,quadToTotalHeightRatio=quadToTotalHeightRatio
                                        ,quadToTotalWidthRatio=quadToTotalWidthRatio
                                        ,widthCorr=widthCorr
                                        ,heightCorr=heightCorr
                                        ,mainImageQuadVert=res
                                        ,mainQuadVertSize = sizeof(res)
                                        ,windowWidthCorr=windowWidthCorr
                                        ,windowHeightCorr=windowHeightCorr
                                        ))

    end #getMainVerticies

end#StructsManag


text:
using AcronymGenerator
using Test

@testset "AcronymGenerator.jl" begin
    # Write your tests here.
end


text:
using QXTools
using Test
using TestSetExtensions

@testset ExtendedTestSet "All QXTools.jl tests" begin
    @includetests ARGS
end


text:
using StochasticSearch, FactCheck

facts("[unit_value]") do
    context("EnumParameter") do
        p = EnumParameter([StringParameter("value", "i0"),
                           FloatParameter(1.3, 56.3, 3.3, "i1"),
                           IntegerParameter(1, 40, 3, "i3")],
                           "test")
        c = p.value
        u = unit_value(p)
        v = unit_value!(p, u)
        @fact (unit_value!(p, u) == c)     --> true
        @fact (1 <= v <= length(p.values)) --> true
    end
    context("IntegerParameter") do
        p = IntegerParameter(1, 40, 23, "i0")
        c = p.value
        u = unit_value(p)
        v = unit_value!(p, u)
        @fact (unit_value!(p, u) == c)     --> true
        @fact (p.min <= v <= p.max)        --> true
    end
    context("FloatParameter") do
        p = FloatParameter(1.0, 40.0, 23.0, "i0")
        c = p.value
        u = unit_value(p)
        v = unit_value!(p, u)
        # Inexact conversion.
        @fact (p.min <= v <= p.max)        --> true
        @fact (unit_value!(p, u) == c)     --> true
    end
end


text:
using HelloJuliaHubHTTPApp

@info("Starting HelloJuliaHub...")

HelloJuliaHubHTTPApp.serve()

@info("Finished HelloJuliaHub")

max_stars_repo_name: tanmaykm/HelloJuliaHubHTTPApp

max_issues_repo_name: tanmaykm/HelloJuliaHubHTTPApp

max_forks_repo_name: tanmaykm/HelloJuliaHubHTTPApp


text:
module RPConstraints

include("./../../constraints/geometric/constants.jl")
using JuMP
using .Constants

export rp_constraints

function rp_constraints(m)

    #AND TOP
    @variable(m, robot_on_top[i=1:length(boxes)-2, j=(i+1):length(boxes)-1], Bin)
    @constraint(m, c_robot_on_top[i=1:length(boxes)-2, j=(i+1):length(boxes)-1],
        m[:ifz_on_top][i, j] + m[:rp_ifx_or][i, j] + m[:rp_ify_or][i, j] >=
        3 * m[:robot_on_top][i, j]
    )

    #AND RIGHT
    @variable(m, robot_right[i=1:length(boxes)-2, j=(i+1):length(boxes)-1], Bin)
    @constraint(m, c_robot_right[i=1:length(boxes)-2, j=(i+1):length(boxes)-1],
        m[:ifx_right][i, j] + m[:rp_ify_or][i, j] + m[:rp_ifz_or][i, j] >=
        3 * m[:robot_right][i, j]
    )

    #AND behind
    @variable(m, robot_behind[i=1:length(boxes)-2, j=(i+1):length(boxes)-1], Bin)
    @constraint(m, c_robot_behind[i=1:length(boxes)-2, j=(i+1):length(boxes)-1],
        m[:ify_behind][i, j] + m[:rp_ifx_or][i, j] + m[:rp_ifz_or][i, j] >=
        3 * m[:robot_behind][i, j]
    )

    @expression(m, robot_positions[i=1:length(boxes)-2, j=(i+1):length(boxes)-1],
        m[:robot_on_top][i, j] + m[:robot_right][i, j] + m[:robot_behind][i, j]
    )

    for k in 2:length(boxes)-1
        positions = []
        for j in k:k
            for i in 1:j
                if i!=j
                    push!(positions, [i, j])
                end
            end
        end
        expressions = []
        for l in 1:length(positions)
            push!(expressions, m[:robot_positions][value(positions[l][1]), value(positions[l][2])])
        end
        @show(sum(expressions[n] for n in 1:length(expressions)))
        @constraint(m,
            sum(expressions[n] for n in 1:length(expressions)) == 1
        )
    end

    # The first box is placed at (0, 0, 0) for this pattern to work.
    @constraint(m, robot_box_one_x, m[:flbx_front_left_bot_x][1] == 0)
    @constraint(m, robot_box_one_y, m[:flby_front_left_bot_y][1] == 0)
    @constraint(m, robot_box_one_z, m[:flbz_front_left_bot_z][1] == 0)

    return m
end

end


text:
#========================================================================================#
#	Laboratory 8
#
# Genetic algorithms and refactoring.
#
# Author: Niall Palfreyman, 16/03/2022
#========================================================================================#
[
	Activity(
		"""
		In this laboratory we will use our Casinos module to implement a genetic algorithm
		(GA). I have already written most of the code for you, but you will see that there
		is still some work left for you to do, so let's get started ...

		Our GA will locate the global minima of a function that we shall call the objective
		of the GA. Objective functions can be very complicated, so we want to be able to view
		them graphically to understand them. I have written a module Objectives to define and
		view objective functions. Run the simple demo1() method of the Objectives module now,
		and tell me to the nearest integer value the x-coordinate of the global minimum of
		the 7-th function in our test-suite of Objective functions:

			include("src/Development/SimpleGAs/Objectives.jl")
			using .Objectives
			Objectives.demo1()
		""",
		"The low values in this contour map are colourd deep blue - look to the top-left! :)",
		x -> x==9
	),
	Activity(
		"""
		Objectives.demo2() shows you how to insert your own function into an Objective. How
		many hills does the Valleys function possess within the displayed domain?
		""",
		"The hills are coloured green-yellow",
		x -> x==3
	),
	Activity(
		"""
		Go to the source code of the Objectives module and look at the Objective type definition.
		the first field of this type is the function itself - whether you take this from the test
		suite or insert your own function. What is the name of the second field?
		""",
		"",
		x -> occursin("domain",lowercase(x))
	),
	Activity(
		"""
		Create several different Objective functions and study the structure of their domains.
		What do you think is the dimension of test suite function 1?
		""",
		"The dimension is the number of number-pairs [lwb,upb] in the domain",
		x -> x==1
	),
	Activity(
		"""
		Use the evaluate() method to find the value of Objective(7) at various locations. Note
		that evaluate() evaluates the objective function at several different points at once,
		so it expects a vector of points. So to check the value for the point [1,2], you will
		need to enter:

			evaluate( obj, [[1,2]])

		What do you now think is the y-coordinate of the global minimum of Objective(7)?
		""",
		"You should find that your earlier answer from graphical inspection is correct",
		x -> x==9
	),
	Activity(
		"""
		Create an instance of the 16-th test suite function. Find out its dimensionality:
		""",
		"You can either inspect the domain or you can use the dim() method",
		x -> x==128
	),
	Activity(
		"""
		Clearly, that is a very large number of dimensions! This makes it very difficult to find
		the global minimum of the function. Indeed, this function is the mepi() (maximally
		epistatic) function created by Robert Watson in 2007 to be maximally difficult for GAs
		to solve. To see how difficult it is to find this function's minimum, you can explore the
		mepi function in any number of dimensions (here 3):

			Objectives.mepi([1,0,1])						# Returns the value 8
			Objectives.mepi.([rand(3) for _ in 1:9])		# Returns nine different mepi values

		Play around with the mepi function to discover its minimum value in 4 dimensions:
		""",
		"Try confining your search to vectors with components either 0 or 1",
		x -> x==4
	),
	Activity(
		"""
		To find out why it is so difficult to find mepi's global minimum, try depicting the
		mepi function in all 128 dimensions. Clearly, you cannot produce a 128-dimensional
		graphic, but you can plot a cross-section of mepi's values in 128 dimensions using the
		command

			depict(Objective(16),[])

		Modify the code of the demo1() method to display this cross-section. Suppose we only had
		8 dimensions. In this case the two minima occur at the vector arguments [0,0,0,0,0,0,0,0]
		and [1,1,1,1,1,1,1,1]. The next lowest local minimum occurs at a vector maximally distant
		from both of these two vectors. What is this vector?
		""",
		"How many 0/1 switches are necessary to convert one of these vectors into the other?",
		x -> x==[0,0,0,0,1,1,1,1] || x==[1,1,1,1,0,0,0,0]
	),
	Activity(
		"""
		OK, now let's look at the module Decoders. Our GA will work with genes consisting of
		0/1 values, but we want to optimise objective functions whose arguments consist of
		real numbers. We therefore need a Decoder to convert our GA chromosomes into input
		numbers for an Objective.

		The method Decoders.demo() demonstrates how to use a Decoder to decode the two
		bit-vectors [1,0,1,0,1,0,1,0,0,0,0,1,0,1,1] and [0,1,0,1,0,1,0,1,1,0,1,0,1,1,0]. How
		many bits does each of these two bit-vectors contain?
		""",
		"Count the bits in the vectors! :)",
		x -> x==15
	),
	Activity(
		"""
		Each bit-vector encodes three floating-point values (for example, the three arguments
		of an objective function). How many bits are used to encode each floating-point number?
		""",
		"Length of the bit-vector divided by the number of float values",
		x -> x==5
	),
	Activity(
		"""
		If you look at where the decoder is constructed in the demo() method, you can see that
		it actually specifies a domain with three dimensions and then the number of bits (5)
		that we want to use for encoding each float value. If you then go up to the definition
		of the Decoder constructor, you will see how it converts these two values into the fields
		inside a Decoder.

		Notice in particular how it uses a functional expression (x->x[1]) to pull out all the
		lower-bound values from the domain and store them in the field lwb. What is the lower
		bound of the third dimension of the domain in our demonstration example?
		""",
		"",
		x -> x==-7
	),
	Activity(
		"""
		Now study the code inside the Decoder constructor - it creates a matrix that will be
		stored in the field coeffs of the Decoder. At the Julia prompt, enter the lines that
		calculate the vector qlevels (quantum levels), using the value 5 instead of nbits.

		Notice how qlevels first contains the sequence [1/2,1/4,1/8,1/16,...], and then these
		numbers are divided by their sum, so that finally the sum of all values in qulevels is
		equal to 1. Is qlevels a row vector or a column vector?
		""",
		"",
		x -> occursin("col",lowercase(x))
	),
	Activity(
		"""
		The vector span contains the span-width of each domain dimension, which we will need
		to decode the bits. Investigate whether span is a row vector or a column vector:
		""",
		"",
		x -> occursin("row",lowercase(x))
	),
	Activity(
		"""
		The Decoder constructor calculates the coeffs matrix by multiplying qlevels * span.
		Notice that qlevels has size (5,1), while span has size (1,3). These numbers are
		important: 5 is the accuracy of our bit-coding of numbers, and 3 is the number of
		different dimensions we need to decode. What is the size of coeffs?
		""",
		"",
		x -> x==(5,3)
	),
	Activity(
		"""
		Now comes your work. The result of decoding our two bit-vectors should be the two
		floating-point vectors [1.67742, 3.51613, -2.74194] and [1.32258, 4.41935,  1.51613].
		But if you run demo(), you will see that it calculates a silly result. The reason for
		this is that I have changed the two important lines of code in the method decode(). it
		is your job now to write the necessary code. You will find hints and information in
		the comments. When you are getting the right answers, come back here and reply
		"Yay!!!" ...
		""",
		"",
		x -> x=="Yay!!!"
	),
	Activity(
		"""
		Great - well done! Now it's time to look at the GA itself - the source code is in the
		module SimpleGAs. Once again, the method demo() shows you how to use the SimpleGA, but
		again it won't give you correct answers yet because I've changed the code in a couple
		of places. So let's just start exploring the code.

		As you can see, in demo(), I simply create an Objective and attach it to a SimpleGA,
		then let it run. The main work of finding the minimum of the Objective is done by the
		method evolve!(), so let's go there next.

		By the way, how many iterations of the GA does demo() perform?
		""",
		"It is the second argument of evolve!()",
		x -> x==5000
	),
	Activity(
		"""
		evolve!() is basically just a loop over lots of generations, and in the loop you can
		see that each generation of the SimpleGA (sga) first measures the current fitness of
		sga's population, then uses this information to (a) sort out any elitists that will
		survive unchanged into the next generation, (b) perform selection (i.e. fitness-
		dependent recombination/crossover) and (c) mutate the children in the next generation.

		If you have already successfully implemented the Casinos module, you should now be
		getting a compile-error inside the method SimpleGAs.mutate!(). Go there now and look at
		the first line of code. Here, we call Casinos.draw() to create a random Boolean matrix
		that will tell us which genes to mutate. However, we haven't yet implemented a draw()
		method with 4 arguments. How many arguments does your current draw() method have?
		""",
		"Look at the signature of the method Casinos.draw()",
		x -> x==3
	),
	Activity(
		"""
		First let's look at how to create a Bernoulli matrix of Boolean values with a particular
		(Bernoulli) probability mu of being true. Enter this line at the Julia prompt:

			rand(5,5) < 0.2

		This will give you a random Boolean matrix. Find the number of true elements like this:

			sum(ans)

		Now enter sum(rand(5,5) < 0.2) many times - what is the average result that you get?
		""",
		"You'll need to do it about 10 times to estimate the average value",
		x -> x==5
	),
	Activity(
		"""
		OK, so we can calculate a Bernoulli matrix from a matrix of random floats between 0
		and 1. But we already have a draw() method that produceds random float matrices, so
		you now need to create a new Casinos.draw() method with 4 arguments that calls the
		3-argument method, compares it with the Bernoulli probability mu, and returns the
		resulting Boolean matrix.

		Do this now, then recompile and run SimpleGAs.demo(). Then give me a "Hi 5!" when
		it runs without errors.
		""",
		"",
		x -> occursin('5',x) || occursin("five",lowercase(x))
	),
	Activity(
		"""
		Notice how good this feels: You have brought a broken program to the stage where it
		runs without errors! Of course, the answer you get is not yet correct, but we will
		fix that now ...

		Go to the method SimpleGAs.recombine!(). First notice the code that selects parents
		for recombination according to their fitness. THIS CODE IS TOTALLY COOL! It is not
		easy to understand, but you will need it if you ever want to implement roulette-wheel
		selection for yourself!

		Further down is a loop that actually generates the individual progeny (children)
		sally and billy from their parents mummy[m] and daddy[m]. You will see there are five
		lines of code that you need to write yourself. Do this now according to the explanations
		in the code.

		You can tell when you have succeeded by looking at the results of the GA: The winning
		individual should be close to [9,9], which was your estimate for the location of the
		global minimum of Objective(7). What is the value (to the nearest integer) of the
		objective function at this point?
		""",
		"",
		x -> x==-16
	),
	Activity(
		"""
		Well done! You have now completed the set work for this laboratory, and you have
		created an up-and-running genetic algorithm capable of optimising multi-variable
		functions. Now try out SimpleGA with various other objective functions from our
		test suite. In particular, find out how close you can get to the global minimum of
		mepi() - or Objective(16) - by tuning sga using the following methods:

			mu!(), temperature!(), elite!() and px!()
		""",
		"",
		x -> true
	),
]


text:
using ApproxFun, SingularIntegralEquations, Test
import ApproxFunBase: choosedomainspace, promotedomainspace, ConstantSpace, interlace,
                        testraggedbelowoperator, testblockbandedoperator, blocklengths

@testset "Ideal Fluid Flow" begin
    k = 50
    Γ = Segment(0.,1+0.5im)
    z = Fun(Γ)
    α = exp(-π*k/50im)

    @test choosedomainspace(Hilbert(),z) == JacobiWeight(-0.5,-0.5,ChebyshevDirichlet{1,1}(Γ))

    Ai=[1 Hilbert()]

    @test choosedomainspace(Ai.ops[1],space(z)) isa ConstantSpace

    S = choosedomainspace(Ai,space(z))
    @test component(S,1) isa ConstantSpace
    @test component(S,2) isa JacobiWeight{ChebyshevDirichlet{1,1,Segment{Complex{Float64}},Float64},
                                                            Segment{Complex{Float64}},Float64,Float64}

    AiS = promotedomainspace(Ai,S)
    @test domainspace(AiS) == S

    S = JacobiWeight(0.5,0.5,Ultraspherical(1,Γ))
    @time c,ui = [1 Hilbert(S)]\[imag(α*z)]

    u =(x,y) -> α*(x+im*y)+2cauchy(ui,x+im*y)

    @test u(0.1,0.2) ≈ 0.039532462109794025-0.3188804984805561im # empirical

    k = 227;
    Γ = 0.5+exp(im*Segment(0.1,-42))
    z = Fun(Γ)
    α = exp(-k/50im)

    Ai = [1 PseudoHilbert()]

    @test isa(choosedomainspace(Ai.ops[1],space(z)),ConstantSpace)
    @test isa(component(choosedomainspace(Ai,space(z)),1),ConstantSpace)
    @test isa(component(choosedomainspace(Ai,space(z)),2),
                JacobiWeight{ChebyshevDirichlet{1,1,Arc{Float64,Float64,Complex{Float64}},Float64},
                                                                        Arc{Float64,Float64,Complex{Float64}},Float64})

    S = choosedomainspace(Ai,space(z))
    AiS = promotedomainspace(Ai,S)
    @test domainspace(AiS) == S

    k = 227;
    Γ = 0.5+exp(im*Segment(0.1,-42))
    z = Fun(Γ)
    α = exp(-k/50im)
    S = JacobiWeight(0.5,0.5,Γ)
    @time c,ui = [1 PseudoHilbert(S)]\[imag(α*z)]

    u = (x,y) -> α*(x+im*y)+2pseudocauchy(ui,x+im*y)

    @test u(0.1,0.2) ≈ 0.6063720775017964 - 0.6382733554119975im # empirical

    Γ = Circle()
    z = Fun(Fourier(Γ))

    Ai = [0 DefiniteLineIntegral();
          1 real(Hilbert())]

    S = choosedomainspace(Ai,space([Fun(0.);z]))

    @test isa(component(S,1),ConstantSpace)
    @test component(S,2) == Fourier(Γ)

    @test domainspace(promotedomainspace(Ai,S)) == S

    A=promotedomainspace(Ai,S)

    testraggedbelowoperator(A)

    k=239;
    α=exp(-k/45im)

    c,ui=[0 DefiniteLineIntegral();
          1 real(Hilbert())]\[Fun(0.);imag(α*z)]

    u =(x,y)->α*(x+im*y)+2cauchy(ui,x+im*y)
    @test u(2.,1.1) ≈ 2.426592437403252-0.8340542386599383im

    ## Curve
    Γ=Curve(Fun(x->exp(0.8im)*(x+x^2-1+im*(x-4x^3+x^4)/6)))
        z=Fun(Γ)
        α=im
        S=JacobiWeight(0.5,0.5,Γ)
        @time c,ui=[1 real(Hilbert(S))]\[imag(α*z)]

    u=(x,y)->α*(x+im*y)+2cauchy(ui,x+im*y)

    @test u(0.1,0.2) ≈ (-1.1657816742288978-0.21306668168680534im)

    ## 2 intervals
    Γ=Segment(-1.,-0.5) ∪ Segment(-0.3,1.)
    z=Fun(Γ)

    S=PiecewiseSpace(map(d->JacobiWeight(0.5,0.5,Ultraspherical(1,d)),components(Γ)))

    k=114;
        α=exp(k/50*im)

    Ai = [ones(component(Γ,1))+zeros(component(Γ,2)) zeros(component(Γ,1))+ones(component(Γ,2)) Hilbert(S)]
    testblockbandedoperator(Ai)

    QR = qrfact(Ai)
          resizedata!(QR, :, Block(10)) # check for bug in resizedata!

    @time a,b,ui = [ones(component(Γ,1))+zeros(component(Γ,2)) zeros(component(Γ,1))+ones(component(Γ,2)) Hilbert(S)] \ [imag(α*z)]

    u=(x,y)->α*(x+im*y)+2cauchy(ui,x+im*y)

    @test u(0.1,0.2) ≈ (-0.5762722129639637 + 0.04348367760228282im) # empirical
end


text:
using FractalTools
using Plots
using Cubature

# Definitions
ti, ts, tf = 0, 0.1, 1.
t = ti : ts : tf
f(t) = sin(2π * t)
g(t) = cos(2π * t)
# f(t) = t * (1 - t)

# Construct interpolation
interp = interpolate(t, f.(t), d=0.01)

# Construct hidden interpolation
hiddeninterp = hiddenfif(t, f.(t), g.(t))

# Quadratute integration
hquad(tk) = hquadrature(f, ti, tk, abstol=1e-9)[1]

# Fractal integration
hfrac(tk) = integrate(ti:ts:tk, f.(ti:ts:tk), d=0.01)

# Fractal hidden integration
hhiddenfrac(tk) = integrate(ti:ts:tk, f.(ti:ts:tk), g.(ti:ts:tk))

# Plots
td = 0 : 0.01 * ts : tf
plot(td, f.(td), label="f")
plot!(td, hquad.(td), label="intfquad")
plot!(td, hfrac.(td), label="intffrac")
plot!(td, hhiddenfrac.(td), label="intfhiddenfrac")
plot!(td, getindex.(hiddeninterp.(td), 1), label="fhiddenfrac")
scatter!(t, f.(t), marker=(:circle, 2), label="data")


text:
#**************************************************************************************
# test_Output_Interface.jl
# ======================== part of the GeoEfficiency.jl package.
#
#
#
#**************************************************************************************

using GeoEfficiency
using Base.Test


text:
"""
    Interface(uϵ::Symbol,uσ::Symbol,from_ϵσ::Function,to_ϵσ::Function)

The `Interface` struct contain information to map stress and strain to other related variables such as force and displacement in experimental systems such as atomic force microscopes or tweezers.

"""
struct Interface
		  uϵ::Symbol
		  uσ::Symbol
		  from_ϵσ::Function
		  to_ϵσ::Function
end

"""
    AFM(R::Real)

Create and return an `Interface` struct using the Hertz contact model with indenter radius `R`.
"""
function AFM(R::Real)
  Interface(:d, :f, (ϵ,σ)->(d = (R.*(3. /4.).^(2. /3)) .* ϵ.^(2. ./3), f = R^2 .* σ ), (d,f)->(ϵ = 4 ./(3 .*R.^1.5) .* d.^1.5, σ = f./R^2) )
end

"""
    Tweezers(R::Real, slip::Real)

Create and return an `Interface` struct for a spherical probe of radius `R` embeded in a material, with probe/material viscosity ratio `slip`.
"""
function Tweezers(R::Real, slip::Real=Inf)
  # slip can go between 0 and Inf
  # no slip is at Inf, frictionless slip is at 0
  R >= 0 || error("sphere must have a positive radius")
  slip >= 0 || error("slip less than 0 is not physical")

  if slip == Inf
    hadamard = 3. / 2.
  else
    hadamard = (((3 * slip) + 2) / (2 * (slip + 1)))
  end

  stress_conversion = hadamard * 4 * π * R^2

  Interface(:d, :f, (ϵ,σ)->(d = R .* ϵ, f = stress_conversion .* σ ), (d,f)->(ϵ = d ./ R, σ = f ./ stress_conversion) )
end

function Base.getindex(d::RheoTimeData, i::Interface)
	return( NamedTuple{(i.uϵ, i.uσ),Tuple{Vector{RheoFloat},Vector{RheoFloat}}}( i.from_ϵσ(d.ϵ,d.σ) )   )
end

function RheoTimeData(interface::Interface ; t::Vector{T} = RheoFloat[], comment="Created from generic constructor", savelog = true, log = savelog ? RheoLogItem(comment) : nothing, kwargs...)  where {T<:Real}

   if interface.uϵ in keys(kwargs)
	  uϵ = convert(Vector{RheoFloat}, kwargs[interface.uϵ])
   else
	  uϵ = RheoFloat[]
   end

   if interface.uσ in keys(kwargs)
	  uσ = convert(Vector{RheoFloat}, kwargs[interface.uσ])
   else
	  uσ = RheoFloat[]
   end

  typecheck = check_time_data_consistency(t,uϵ,uσ)
  ϵ,σ = interface.to_ϵσ(uϵ, uσ)

  RheoTimeData(convert(Vector{RheoFloat},σ), convert(Vector{RheoFloat},ϵ), convert(Vector{RheoFloat},t),
  log === nothing ? nothing : [ RheoLogItem(log.action,merge(log.info, (type=typecheck, interface = interface)))]     )

end

"""
    importcsv(filepath::String, interface::Interface; delimiter = ',', header = false, comment = "Imported from csv file", savelog = true, kwargs...)

Import function for raw data to be transformed to stress and strain using an `Interface`.
Column numbers in the csv containging the force/displacement data need to be indicated as keyword arguments using the corresponding symbols in the `interface`..

# Examples
```julia-repl
julia> importcsv("myfile.csv", AFM(2e-6), t=1, d=2, f=3)
```
"""
function importcsv(filepath::String, interface::Interface; delimiter = ',', header = false, comment = "Imported from csv file", savelog = true, kwargs...)
   cols=symbol_to_unicode(kwargs.data)

    # legacy keyword arguments in place until v 1.0
   uϵ_col_sym = Symbol(string(interface.uϵ) * "_col")
   uσ_col_sym = Symbol(string(interface.uσ) * "_col")

   if uϵ_col_sym in keys(cols)
	  uϵ_col = cols[uϵ_col_sym]
   elseif interface.uϵ in keys(cols)
	  uϵ_col = cols[interface.uϵ]
   else
	  uϵ_col = nothing
   end

   if uσ_col_sym in keys(cols)
	  uσ_col = cols[uσ_col_sym]
    elseif interface.uσ in keys(cols)
      uσ_col = cols[interface.uσ]
    else
	  uσ_col = nothing
   end

   data = importcsv(filepath, t = cols.t, σ = uσ_col, ϵ = uϵ_col, delimiter = delimiter, header = header, comment = comment, savelog = false)

   # need to make sure that missing columns are properly transfered.
   ϵ,σ = interface.to_ϵσ(data.ϵ, data.σ)

   log = if savelog
			 info = (comment=comment, folder=pwd(), stats=(t_min=data.t[1],t_max=data.t[end], n_sample=size(data.t[:])))
			 kwds = NamedTuple{(:t, interface.uϵ, interface.uσ, :header),Tuple{Any,Any,Any,Bool}}( (cols.t, uϵ_col, uσ_col, header) )
			 RheoLogItem( (type=:source, funct=:importcsv, params=(filepath=filepath, interface=interface), keywords=kwds), info )
		 else
			 RheoLogItem(type = nothing, info = "nothing")
		 end

   return RheoTimeData(σ, ϵ, data.t, [log])

end


text:
# Types

type Mixture{C}  # C is the component type
	K::Int
	components::Vector{C}
	pi::Vector{Float64}
end

function Mixture{C}(components::Vector{C}, pi::Vector{Float64})
	K = length(components)
	@check_argdims K == length(pi)
	Mixture{C}(K, components, pi)
end

ncomponents(m::Mixture) = m.K


text:
resample(uv::AbstractVector{<:Number}) = uv

# No constraints
function resample(uv::DT) where {DT <: AbstractVector{<:AbstractUncertainValue}}
    resample.(uv)
end

function resample(uv::DT, n::Int) where {DT <: AbstractVector{<:AbstractUncertainValue}}
    [resample.(uv) for i = 1:n]
end

# Vectors of sampling constraints mapped to unique sampling constraints
function resample(uv::DT, constraint::Vector{<:SamplingConstraint}) where {
        DT <: AbstractVector{<:AbstractUncertainValue}}
    [resample(uv[i], constraint[i]) for i in 1:length(uv)]
end

function resample(uv::DT, constraint::Vector{<:SamplingConstraint}, n::Int) where {
        DT <: AbstractVector{<:AbstractUncertainValue}}
    [[resample(uv[i], constraint[i]) for i in 1:length(uv)] for k = 1:n]
end

# Vectors of sampling constraints mapped to unique sampling constraints
function resample(uv::DT, constraint::SamplingConstraint) where {
        DT <: AbstractVector{<:AbstractUncertainValue}}
    [resample(uv[i], constraint) for i in 1:length(uv)]
end

function resample(uv::DT, constraint::SamplingConstraint, n::Int) where {
        DT <: AbstractVector{<:AbstractUncertainValue}}
    [[resample(uv[i], constraint) for i in 1:length(uv)] for k = 1:n]
end

function resample(uv::DT, constraint::NoConstraint) where {
    DT <: AbstractVector{<:AbstractUncertainValue}}

    return uv
end

function resample(uv::DT, constraint::NoConstraint, n::Int) where {
    DT <: AbstractVector{<:AbstractUncertainValue}}

    return [uv for i = 1:n]
end


text:
# ----------------------------------------------------------------------------------- #
# Copyright (c) 2016 Varnerlab
# Robert Frederick School of Chemical and Biomolecular Engineering
# Cornell University, Ithaca NY 14850
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
# THE SOFTWARE.
# ----------------------------------------------------------------------------------- #
#
# ----------------------------------------------------------------------------------- #
# Balances: Evaluates model equations given time, state and the data_dictionary.
# Type: GRN-JULIA
# Version: 1.0
#
# Input arguments:
# t  - current time
# x  - state array
# data_dictionary  - Data dictionary instance (holds model parameters)
#
# Return arguments:
# dxdt - derivative array at current time step
# ----------------------------------------------------------------------------------- #
function Balances(t,x,data_dictionary)

  # correct for negatives -
  idx_small = find(x.<0)
  x[idx_small] = 0.0

  # Get model matricies and other required data from the data_dictionary -
  stoichiometric_matrix = data_dictionary["stoichiometric_matrix"]
	dilution_matrix = data_dictionary["dilution_matrix"]
	degradation_matrix = data_dictionary["degradation_matrix"]
  mugmax = data_dictionary["maximum_specific_growth_rate"]

  # Calculate the kinetics -
  transcription_rate_array = calculate_transcription_rates(t,x,data_dictionary)
  translation_rate_array = calculate_translation_rates(t,x,data_dictionary)
  mRNA_degradation_rate_array = calculate_mRNA_degradation_rates(t,x,data_dictionary)
  protein_degradation_rate_array = calculate_protein_degradation_rates(t,x,data_dictionary)
  background_transcription_rate_array = calculate_background_transcription_rates(t,x,transcription_rate_array,data_dictionary)
  input_array = calculate_input_array(t,x,data_dictionary)

  # Call my control function -
  control_array = Control(t,x,data_dictionary)

  # Compute the modified rate -
  transcription_rate_array = transcription_rate_array.*control_array;

  # rate array -
  txtl_rate_array = [transcription_rate_array ; translation_rate_array]
  degradation_rate_array = [mRNA_degradation_rate_array ; protein_degradation_rate_array]

  # Evaluate the balance equations -
  dxdt_array = stoichiometric_matrix*txtl_rate_array+degradation_matrix*degradation_rate_array+mugmax*dilution_matrix*x+background_transcription_rate_array+input_array

  # correct for trigger and cRAF dilution -
  dxdt_array[53] = 0.0;
  dxdt_array[54] = 0.0;

  # return -
  return dxdt_array
end


text:
#=
Levis, C., Johnson, J. T., & Teixeira, F. L. (2010). Radiowave propagation: physics and applications. John Wiley & Sons.
=#
@doc """
 Flat Earth Model
 freq: signal frequency
 txH: base station's height
 rxH: mobile's height
 """
type FlatEarthModel<:Model
  freq::Float64
  txH::Float64
  rxH::Float64
 FlatEarthModel() = new(800,50,1.5)
end

"Calculate Path loss for distance dist (km) with Flat Earth Model"
function pathloss_op(model::FlatEarthModel, dist::AbstractFloat)
  hm::Float64 = model.rxH
  hb::Float64 = model.txH
  d::Float64 = dist
  f::Float64 = model.freq

   L1 = -20*log10(hb)
   L2 = -20*log10(hm)
   Lo = 120 +10*4*log10(d)
   L = Lo + L1 + L2
   return L
end # pathloss_op end