text:
using InvertibleNetworks, LinearAlgebra, Test

# Input
nx1 = 16
nx2 = 16
nx_in = 2
nx_hidden = 4
batchsize = 2

# Observed data
nd1 = 14
nd2 = 8
nd_in = 1

# Additive SLIM layer with learned operator
CS = LearnedCouplingLayerSLIM(nx1, nx2, nx_in, nd1, nd2, nd_in, nx_hidden, batchsize; logdet=false, permute=true)

# Initializations
D = randn(Float32, nd1*nd2*nd_in, batchsize)
X = randn(Float32, nx1, nx2, nx_in, batchsize)

###################################################################################################

# Test invertibility
Y = CS.forward(X, D)
X_ = CS.inverse(Y, D)
@test isapprox(norm(X - X_)/norm(X), 0f0, atol=1e-6)

Y = CS.inverse(X, D)
X_ = CS.forward(Y, D)
@test isapprox(norm(X - X_)/norm(X), 0f0, atol=1e-6)

Y = CS.forward(X, D)
X_ = CS.backward(0f0.*Y, Y, D)[3]
@test isapprox(norm(X - X_)/norm(X), 0f0, atol=1e-6)

###################################################################################################

# Initializations
CS = LearnedCouplingLayerSLIM(nx1, nx2, nx_in, nd1, nd2, nd_in, nx_hidden, batchsize; logdet=true, permute=true)

D = randn(Float32, nd1*nd2*nd_in, batchsize) .* 10f0
D0 = randn(Float32, nd1*nd2*nd_in, batchsize) .* 10f0
dD = D - D0

X = randn(Float32, nx1, nx2, nx_in, batchsize)
X0 = randn(Float32, nx1, nx2, nx_in, batchsize)
dX = X - X0

function loss(CS, X, D)
    Y, logdet = CS.forward(X, D)
    f = -log_likelihood(Y) - logdet
    ΔY = -∇log_likelihood(Y)
    ΔX, ΔD, X_ = CS.backward(ΔY, Y, D)
    @test isapprox(norm(X - X_)/norm(X), 0f0; atol=1f-4)
    return f, ΔX, ΔD, CS.RB.W1.grad
end

# Gradient test for input
f0, gX, gD = loss(CS, X0, D0)[1:3]
h = 0.1f0
maxiter = 6
err1 = zeros(Float32, maxiter)
err2 = zeros(Float32, maxiter)

print("\nGradient test slim layer input\n")
for j=1:maxiter
    f = loss(CS, X0 + h*dX, D0 + h*dD)[1]
    err1[j] = abs(f - f0)
    err2[j] = abs(f - f0  - h*dot(gX, dX) - h*dot(gD, dD))
    print(err1[j], "; ", err2[j], "\n")
    global h = h/2f0
end

@test isapprox(err1[end] / (err1[1]/2^(maxiter-1)), 1f0; atol=1f1)
@test isapprox(err2[end] / (err2[1]/4^(maxiter-1)), 1f0; atol=1f1)

# Gradient test for weights
CS = LearnedCouplingLayerSLIM(nx1, nx2, nx_in, nd1, nd2, nd_in, nx_hidden, batchsize; logdet=true, permute=true)
CS0 = LearnedCouplingLayerSLIM(nx1, nx2, nx_in, nd1, nd2, nd_in, nx_hidden, batchsize; logdet=true, permute=true)
CS.RB.W1.data *= 10; CS0.RB.W1.data *= 10     # make weights larger
CSini = deepcopy(CS0)
dW = CS.RB.W1.data - CS0.RB.W1.data   # just test for 2 parameters

f0, gX, gD, gW = loss(CS0, X, D)
h = 0.1f0
maxiter = 5
err3 = zeros(Float32, maxiter)
err4 = zeros(Float32, maxiter)

print("\nGradient test slim layer weights\n")
for j=1:maxiter
    CS0.RB.W1.data = CSini.RB.W1.data + h*dW
    f = loss(CS0, X, D)[1]
    err3[j] = abs(f - f0)
    err4[j] = abs(f - f0 - h*dot(dW, gW))
    print(err3[j], "; ", err4[j], "\n")
    global h = h/2f0
end

@test isapprox(err3[end] / (err3[1]/2^(maxiter-1)), 1f0; atol=1f1)
@test isapprox(err4[end] / (err4[1]/4^(maxiter-1)), 1f0; atol=1f1)


text:
macro judge_var(model, variable, class, aargs, aakws)
   lag=0
   span=1000
   initial=0.0
   ub=nothing
   lb=nothing
   penalty=nothing
   state_name=:nothing

   vartype=:Con

   if length(aargs)==0
      push!(aargs,:Con)
   end

   if length(aargs)!=1
      ex = quote
         error("@"*string($class)*" macro takes at most three positional arguments")
      end
      return ex
   elseif aargs[1] ∉ [:Con,:Bin,:Int]
      ex = quote
         error("Optional third positional argument must be \'Bin\' or \'Int\'")
      end
      return ex
   end

   for (a,b) in aakws
      if a==:lag
         lag=b
      elseif a==:duration
         span=b
      elseif a==:initial
         initial=b
      elseif a==:lb
         lb=b
      elseif a==:ub
         ub=b
      elseif a==:penalty
         penalty=b
      elseif a==:state_name
         state_name=b
      else
         ex = quote
            error("Invalid keyword argument for @"*string($class)*" macro")
         end
         return ex
      end
   end

   if class==:(:state)
      if lag!=0
         @warn("'lag' keyword has been ignored")
         lag=0
      end
      if span!=1000
         @warn("'duration' keyword has been ignored")
      end
      span=1
   elseif state_name!=:nothing
      @warn("'state_name' keyword has been ignored")
      state_name=:nothing
   end

   if !(class==:state || class==:enforced) && penalty!=nothing
      @warn("'penalty' keyword has been ignored")
      penalty=nothing
   end

   tmp=nothing

   ex = quote
      if !haskey($model.ext, :expansions)
         $model.ext[:expansions] = Dict{Symbol,Any}()
         $model.ext[:options] = Dict{Symbol,Tuple}()
      end
      if length($aargs)==1
         if $aargs[1]==:Con
            tmp=@variable($model, $variable)
         elseif $aargs[1]==:Bin
            tmp=@variable($model, $variable, Bin)
         else
            tmp=@variable($model, $variable, Int)
         end
      end

      if $(Meta.quot(state_name))==:nothing
         sym=[k for (k,v) in $model.obj_dict if v===tmp][1]
         $model.ext[:expansions][sym]=tmp
         $model.ext[:options][sym]=($class,$lag,$span,$aargs[1],$lb,$ub,$initial,$penalty)
      else
         $model.ext[:expansions][$(Meta.quot(state_name))]=tmp
         $model.ext[:options][$(Meta.quot(state_name))]=($class,$lag,$span,$aargs[1],$lb,$ub,$initial,$penalty)
         $state_name=tmp
      end

   end
   return esc(ex)
end

"""
	expansion(model, variable, args...)

Defines an expansion variable `variable` within a subproblem `model`. Note that all subproblems must have the same set of expansion variables.

### Required Arguments
`model` is the JuDGE subproblem that we are adding the expansion variable to

`variable` is the name of the variable being created, this will be continuous by default; follows JuMP syntax if defining a set of variables.

### Optional Arguments
This macro has a third, unnamed, argument which can be set to Con, Bin, or Int, similar to the `@variable` macro.

`lag` is the number of nodes in the scenario between an expansion being decided, and it becoming available.

`duration` is the number of consecutive nodes in the scenario over which an expansion is available.

`lb` is the lower bound for this variable in the master problem (typically omitted).

`ub` is the upper bound for this variable in the master problem (typically omitted).

### Examples
    @expansion(model, expand[1:5], Bin) #defines an array of 5 binary variables with no lag, and unlimited lifespan
    @expansion(model, expand[1:5,1:2]>=0, lag=1) #defines a matrix of 10 continuous variables with a lag of 1, and unlimited duration
    @expansion(model, 0<=expand<=10, Int, duration=2) #defines a single integer variable with a lag of 0, and a duration of 2
"""
macro expansion(model, variable, args...)
   aargs = []
   aakws = Pair{Symbol,Any}[]
   for el in args
      if Meta.isexpr(el, :(=))
         push!(aakws, Pair(el.args...))
      else
         push!(aargs, el)
      end
   end
   ex = quote
      JuDGE.@judge_var($model, $variable, :expansion, $aargs, $aakws)
   end

   return esc(ex)
end

"""
	shutdown(model, variable, args...)

Defines an shutdown variable `variable` within a subproblem `model`. Note that all subproblems must have the same set of shutdown variables.

### Required Arguments
`model` is the JuDGE subproblem that we are adding the shutdown variable to

`variable` is the name of the variable being created, this will be continuous by default; follows JuMP syntax if defining a set of variables.

### Optional Arguments
This macro has a third, unnamed, argument which can be set to Con, Bin, or Int, similar to the `@variable` macro.

`lag` is the number of nodes in the scenario between an shutdown being decided, and it becoming unavailable.

`duration` is the number of consecutive nodes in the scenario over which the shutdown will last.

`lb` is the lower bound for this variable in the master problem (typically omitted).

`ub` is the upper bound for this variable in the master problem (typically omitted).

### Examples
    @shutdown(model, shut[1:5], Bin) #defines an array of 5 binary variables with no lag, and unlimited lifespan
    @shutdown(model, shut[1:5,1:2]>=0, lag=1) #defines a matrix of 10 continuous variables with a lag of 1, and unlimited duration
    @shutdown(model, 0<=shut<=10, Int, duration=2) #defines a single integer variable with a lag of 0, and a duration of 2
"""
macro shutdown(model, variable, args...)
   aargs = []
   aakws = Pair{Symbol,Any}[]
   for el in args
      if Meta.isexpr(el, :(=))
         push!(aakws, Pair(el.args...))
      else
         push!(aargs, el)
      end
   end
   ex = quote
      JuDGE.@judge_var($model, $variable, :shutdown, $aargs, $aakws)
   end

   return esc(ex)
end

"""
	enforced(model, variable, args...)

Defines an enforced variable `variable` within a subproblem `model`. Note that all subproblems must have the same set of enforced variables. These
variables can be used as either expansion or shutdown variables, but since the constraint in the master problem is an equality, convergence is
more difficult since there is less flexibility when solving the master problem.

### Required Arguments
`model` is the JuDGE subproblem that we are adding the expansion variable to

`variable` is the name of the variable being created, this will be continuous by default; follows JuMP syntax if defining a set of variables.

### Optional Arguments
This macro has a third, unnamed, argument which can be set to Con, Bin, or Int, similar to the `@variable` macro.

`lag` is the number of nodes in the scenario between an expansion being decided, and it becoming available.

`duration` is the number of consecutive nodes in the scenario over which an expansion is available.

`lb` is the lower bound for this variable in the master problem (typically omitted).

`ub` is the upper bound for this variable in the master problem (typically omitted).

`penalty` is a placeholder for a future feature, which may allow the violation of master/subproblem equality constraint, at a cost.

### Examples
    @expansion(model, forced[1:5], Bin) #defines an array of 5 binary variables with no lag, and unlimited lifespan
    @expansion(model, forced[1:5,1:2]>=0, lag=1) #defines a matrix of 10 continuous variables with a lag of 1, and unlimited duration
    @expansion(model, 0<=forced<=10, Int, duration=2) #defines a single integer variable with a lag of 0, and a duration of 2
"""
macro enforced(model, variable, args...)
   aargs = []
   aakws = Pair{Symbol,Any}[]
   for el in args
      if Meta.isexpr(el, :(=))
         push!(aakws, Pair(el.args...))
      else
         push!(aargs, el)
      end
   end
   ex = quote
      JuDGE.@judge_var($model, $variable, :enforced, $aargs, $aakws)
   end

   return esc(ex)
end

"""
	state(model, variable, args...)

Defines a state variable `variable` within a subproblem `model`. Note that all subproblems must have the same set of state variables. These
variables can be used to model inventory that is carried forward between the subproblems.

### Required Arguments
`model` is the JuDGE subproblem that we are adding the expansion variable to

`variable` is the name of the variable being created in the subproblem, this will be continuous by default; follows JuMP syntax if defining a set of variables.
The subproblem variable corresponds to the change in the state.

### Optional Arguments
This macro has a third, unnamed, argument which can be set to Con, Bin, or Int, similar to the `@variable` macro.

`state_name` is the name for the state variable in the master problem. If omitted, the name of the master problem variable
will match the subproblem variable (but this may cause confusion, since only the master problem variable is the state).
See the `inventory.jl` example to see how this should be implemented.

`initial` is the initial value for the master problem's state variable at the root node.

`lb` is the lower bound for the variable in the master problem (typically omitted).

`ub` is the upper bound for the variable in the master problem (typically omitted).

### Examples
    @state(sp, -50<=Δstock<=50, state_name=stock, lb=0, ub=200, initial=0) #defines a state variable called stock in the master
                                                                           #(starting at 0, and able to take values 0 to 200),
                                                                           #and Δstock in the subproblem (able to change the stock level by ±50).
"""
macro state(model, variable, args...)
   aargs = []
   aakws = Pair{Symbol,Any}[]
   for el in args
      if Meta.isexpr(el, :(=))
         push!(aakws, Pair(el.args...))
      else
         push!(aargs, el)
      end
   end
   ex = quote
      JuDGE.@judge_var($model, $variable, :state, $aargs, $aakws)
   end

   return esc(ex)
end

"""
	capitalcosts(model, expr)

Defines a linear expression specifying the capital cost of expansions and shutdowns at the current node

### Required Arguments
`model` is the JuDGE subproblem corresponding to the node in the scenario tree that we are adding specifying the costs for

`expr` is an `AffExpr` which gives the total cost of choosing expansion and shutdown variables at the current node

### Example
    @capitalcosts(model, sum(expand[i]*cost[node][i] for i in 1:5))
"""
macro capitalcosts(model, expr)
   ex = quote
      #$model.ext[:capitalcosts] = @expression($model, $expr)
      $model.ext[:capitalcosts] = Dict()
      if typeof($expr)==AffExpr
         for (term,coef) in $expr.terms
            $model.ext[:capitalcosts][term]=coef
         end
         $model.ext[:capitalcosts][:constant]=$expr.constant
      elseif typeof($expr)==VariableRef
         $model.ext[:capitalcosts][$expr]=1.0
         $model.ext[:capitalcosts][:constant]=0.0
      elseif typeof($expr)==Float64
         $model.ext[:capitalcosts][:constant]=$expr
      end
   end
   return esc(ex)
end

"""
	ongoingcosts(model, expr)

Defines a linear expression specifying the ongoing costs of expansions and shutdowns available at the current node

### Required Arguments
`model` is the JuDGE subproblem corresponding to the node in the scenario tree that we are adding specifying the costs for

`expr` is an `AffExpr` which gives the ongoing cost of expansions and shutdowns available at the current node

### Example
    @ongoingcosts(model, sum(expand[i]*ongoingcosts[node][i] for i in 1:5))
"""
macro ongoingcosts(model, expr)
   ex = quote
#      $model.ext[:ongoingcosts] = @expression($model, $expr)
      $model.ext[:ongoingcosts] = Dict()
      if typeof($expr)==AffExpr
         for (term,coef) in $expr.terms
            $model.ext[:ongoingcosts][term]=coef
         end
         $model.ext[:ongoingcosts][:constant]=$expr.constant
      elseif typeof($expr)==VariableRef
         $model.ext[:ongoingcosts][$expr]=1.0
         $model.ext[:ongoingcosts][:constant]=0.0
      elseif typeof($expr)==Float64
         $model.ext[:ongoingcosts][:constant]=$expr
      end
   end
   return esc(ex)
end

ext: jl


text:
# switch statement (with fallout)

macro switch(value, block)
  @gensym value_sym
  clauses, statements = switch(value_sym, block)
  :(for i in 1:1
      let $(esc(value_sym)) = $value
          $(clauses)
          $(esc(statements))
      end
    end)
end

@metafunction switch(value, begin :L{:case}(x); :?{xs}; :L{:case}(y); *{ys} end) begin
  @gensym label
  nextclauses, nextstatements = switch(value, quote :case($y); $(ys...) end)
  clauses    = switch_clause(value, x, label, nextclauses)
  statements = switch_statements(xs, label, nextstatements)

  clauses, statements
end

@metafunction switch(value, begin :L{:case}(x); *{xs} end) begin
  @gensym label
  clauses    = switch_clause(value, x, label)
  statements = switch_statements(xs, label)

  clauses, statements
end

switch_clause(value, x, label, nextclauses=nothing) = quote
  if $value == $x
    @goto $label
  else
    $nextclauses
  end
end

switch_statements(statements, label, nextstatements=nothing) = quote
  @label $label
  $(statements...);
  $(nextstatements)
end

# usage

r = 0
a = 2

@switch a begin
  :case(1);
    r = 1
    break

  :case(2);
    r = 2

  :case(3);
    r = 3

end

import Base.Test.@test

if a in [1,2]
  @test a != 1 || r == 1
  @test a != 2 || r == 3
end


text:
module bushing_examples
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsDeforLinear.AlgoDeforLinearModule
using FinEtools.MeshExportModule
using Statistics: mean

E = 3.0;
nu = 0.4999999; # This is the correct Poisson ratio
Ri = 0.25;
Re = 1.0;
L = Ri/2;
ang = 180/180*pi;
iuz = Ri/15
p = 0.27;
tolerance = min(L, Re-Ri, ang)/1000;
nR, nc, nt = 7, 14, 1

function bushing_h8_full()

	fens, fes  = H8block(ang, L, Re-Ri, nc, nt, nR)
	internal_fenids = selectnode(fens, box = [0 ang 0 L 0 0], inflate = tolerance);

	MR = DeforModelRed3D
	material = MatDeforElastIso(MR, E, nu)
	femm = FEMMDeforLinear(MR, IntegDomain(fes, GaussRule(3, 2)), material)

	boundaryfes = meshboundary(fes);
	ibcl = selectelem(fens, boundaryfes, box = [0 ang 0 L 0 0], inflate = tolerance);
	ebcl = selectelem(fens, boundaryfes, box = [0 ang 0 L Re-Ri Re-Ri], inflate = tolerance);

	for i in 1:count(fens)
		a = fens.xyz[i,1]; y = fens.xyz[i,2]; r = fens.xyz[i,3]+Ri;
		fens.xyz[i,:] .= (r*sin(a), y, r*cos(a));
	end

	geom = NodalField(fens.xyz)
	u = NodalField(zeros(size(fens.xyz,1), 3)) # displacement field

	x0nl = selectnode(fens, box = [0 0 -Inf Inf -Inf Inf], inflate = tolerance);
	setebc!(u, x0nl, true, 1, 0.0)

	y0nl = selectnode(fens, box = [-Inf Inf 0 0 -Inf Inf], inflate = tolerance);
	setebc!(u, y0nl, true, 2, 0.0)
	yLnl = selectnode(fens, box = [-Inf Inf L L -Inf Inf], inflate = tolerance);
	setebc!(u, yLnl, true, 2, 0.0)

	enl = connectednodes(subset(boundaryfes, ebcl))
	setebc!(u, enl, true, [1; 2; 3]; val = 0.0)
	inl = connectednodes(subset(boundaryfes, ibcl))
	setebc!(u, inl, true, [1; 2]; val = 0.0)
	setebc!(u, inl, true, [3]; val = iuz)

	applyebc!(u)
	numberdofs!(u)

	associategeometry!(femm, geom)

	K = stiffness(femm, geom, u)
	F = nzebcloadsstiffness(femm, geom, u)

	scattersysvec!(u, K\F)

	File =  "bushing_h8_full_u.vtk"
	vtkexportmesh(File, fens, fes; vectors = [("u", u.values)])
	@async run(`"paraview.exe" $File`)

    true
end # bushing_h8_full

function bushing_h8_algo_full()

	fens, fes  = H8block(ang, L, Re-Ri, nc, nt, nR)
	internal_fenids = selectnode(fens, box = [0 ang 0 L 0 0], inflate = tolerance);

	MR = DeforModelRed3D
	material = MatDeforElastIso(MR, E, nu)
	femm = FEMMDeforLinear(MR, IntegDomain(fes, GaussRule(3, 2)), material)
	regions = FDataDict[FDataDict("femm"=>femm)]

	boundaryfes = meshboundary(fes);
	ibcl = selectelem(fens, boundaryfes, box = [0 ang 0 L 0 0], inflate = tolerance);
	ebcl = selectelem(fens, boundaryfes, box = [0 ang 0 L Re-Ri Re-Ri], inflate = tolerance);

	for i in 1:count(fens)
		a = fens.xyz[i,1]; y = fens.xyz[i,2]; r = fens.xyz[i,3]+Ri;
		fens.xyz[i,:] .= (r*sin(a), y, r*cos(a));
	end

	essential_bcs = FDataDict[]
	x0nl = selectnode(fens, box = [0 0 -Inf Inf -Inf Inf], inflate = tolerance);
	push!(essential_bcs, FDataDict("displacement"=> 0.0, "component"=> 1, "node_list"=>x0nl))

	y0nl = selectnode(fens, box = [-Inf Inf 0 0 -Inf Inf], inflate = tolerance);
	push!(essential_bcs, FDataDict("displacement"=> 0.0, "component"=> 2, "node_list"=>y0nl))
	yLnl = selectnode(fens, box = [-Inf Inf L L -Inf Inf], inflate = tolerance);
	push!(essential_bcs, FDataDict("displacement"=> 0.0, "component"=> 2, "node_list"=>yLnl))

	enl = connectednodes(subset(boundaryfes, ebcl))
	push!(essential_bcs, FDataDict("displacement"=> 0.0, "component"=> 1, "node_list"=>enl))
	push!(essential_bcs, FDataDict("displacement"=> 0.0, "component"=> 2, "node_list"=>enl))
	push!(essential_bcs, FDataDict("displacement"=> 0.0, "component"=> 3, "node_list"=>enl))
	inl = connectednodes(subset(boundaryfes, ibcl))
	push!(essential_bcs, FDataDict("displacement"=> 0.0, "component"=> 1, "node_list"=>inl))
	push!(essential_bcs, FDataDict("displacement"=> 0.0, "component"=> 2, "node_list"=>inl))
	push!(essential_bcs, FDataDict("displacement"=> iuz, "component"=> 3, "node_list"=>inl))

	modeldata = FDataDict("fens"=>fens,
	"regions"=>regions,
	"essential_bcs"=>essential_bcs
	)
	modeldata = AlgoDeforLinearModule.linearstatics(modeldata)

	filebase = "bushing_h8_algo_full"
	modeldata["postprocessing"] = FDataDict("file"=>filebase * "-s", "quantity"=>:pressure, "component"=>collect(1:1), "nodevalmethod"=>:averaging, "reportat"=>:extraptrend)
	modeldata = AlgoDeforLinearModule.exportstress(modeldata)

    true
end # bushing_h8_algo_full

function bushing_h8_algo_ms()

	fens, fes  = H8block(ang, L, Re-Ri, nc, nt, nR)
	internal_fenids = selectnode(fens, box = [0 ang 0 L 0 0], inflate = tolerance);

	MR = DeforModelRed3D
	material = MatDeforElastIso(MR, E, nu)
	femm = FEMMDeforLinearMSH8(MR, IntegDomain(fes, GaussRule(3, 2)), material)
	regions = FDataDict[FDataDict("femm"=>femm)]

	boundaryfes = meshboundary(fes);
	ibcl = selectelem(fens, boundaryfes, box = [0 ang 0 L 0 0], inflate = tolerance);
	ebcl = selectelem(fens, boundaryfes, box = [0 ang 0 L Re-Ri Re-Ri], inflate = tolerance);

	for i in 1:count(fens)
		a = fens.xyz[i,1]; y = fens.xyz[i,2]; r = fens.xyz[i,3]+Ri;
		fens.xyz[i,:] .= (r*sin(a), y, r*cos(a));
	end

	essential_bcs = FDataDict[]
	x0nl = selectnode(fens, box = [0 0 -Inf Inf -Inf Inf], inflate = tolerance);
	push!(essential_bcs, FDataDict("displacement"=> 0.0, "component"=> 1, "node_list"=>x0nl))

	y0nl = selectnode(fens, box = [-Inf Inf 0 0 -Inf Inf], inflate = tolerance);
	push!(essential_bcs, FDataDict("displacement"=> 0.0, "component"=> 2, "node_list"=>y0nl))
	yLnl = selectnode(fens, box = [-Inf Inf L L -Inf Inf], inflate = tolerance);
	push!(essential_bcs, FDataDict("displacement"=> 0.0, "component"=> 2, "node_list"=>yLnl))

	enl = connectednodes(subset(boundaryfes, ebcl))
	push!(essential_bcs, FDataDict("displacement"=> 0.0, "component"=> 1, "node_list"=>enl))
	push!(essential_bcs, FDataDict("displacement"=> 0.0, "component"=> 2, "node_list"=>enl))
	push!(essential_bcs, FDataDict("displacement"=> 0.0, "component"=> 3, "node_list"=>enl))
	inl = connectednodes(subset(boundaryfes, ibcl))
	push!(essential_bcs, FDataDict("displacement"=> 0.0, "component"=> 1, "node_list"=>inl))
	push!(essential_bcs, FDataDict("displacement"=> 0.0, "component"=> 2, "node_list"=>inl))
	push!(essential_bcs, FDataDict("displacement"=> iuz, "component"=> 3, "node_list"=>inl))

	modeldata = FDataDict("fens"=>fens,
	"regions"=>regions,
	"essential_bcs"=>essential_bcs
	)
	modeldata = AlgoDeforLinearModule.linearstatics(modeldata)

	filebase = "bushing_h8_algo_ms"
	modeldata["postprocessing"] = FDataDict("file"=>filebase * "-nodal-", "quantity"=>:pressure, "component"=>collect(1:1), "nodevalmethod"=>:averaging, "reportat"=>:extraptrend)
	modeldata = AlgoDeforLinearModule.exportstress(modeldata)
	modeldata["postprocessing"] = FDataDict("file"=>filebase * "-elwise-", "quantity"=>:pressure, "component"=>collect(1:1), "nodevalmethod"=>:averaging, "reportat"=>:extraptrend)
	modeldata = AlgoDeforLinearModule.exportstresselementwise(modeldata)

    true
end # bushing_h8_algo_ms

function bushing_h8u_algo_ms()

	fens, fes  = T4block(ang, L, Re-Ri, Int.(round.((nc, nt, nR)./2).+1)...)
	fens, fes  = T4toH8(fens, fes)
	internal_fenids = selectnode(fens, box = [0 ang 0 L 0 0], inflate = tolerance);

	MR = DeforModelRed3D
	material = MatDeforElastIso(MR, E, nu)
	femm = FEMMDeforLinearMSH8(MR, IntegDomain(fes, GaussRule(3, 2)), material)
	regions = FDataDict[FDataDict("femm"=>femm)]

	boundaryfes = meshboundary(fes);
	ibcl = selectelem(fens, boundaryfes, box = [0 ang 0 L 0 0], inflate = tolerance);
	ebcl = selectelem(fens, boundaryfes, box = [0 ang 0 L Re-Ri Re-Ri], inflate = tolerance);

	for i in 1:count(fens)
		a = fens.xyz[i,1]; y = fens.xyz[i,2]; r = fens.xyz[i,3]+Ri;
		fens.xyz[i,:] .= (r*sin(a), y, r*cos(a));
	end

	essential_bcs = FDataDict[]
	x0nl = selectnode(fens, box = [0 0 -Inf Inf -Inf Inf], inflate = tolerance);
	push!(essential_bcs, FDataDict("displacement"=> 0.0, "component"=> 1, "node_list"=>x0nl))

	y0nl = selectnode(fens, box = [-Inf Inf 0 0 -Inf Inf], inflate = tolerance);
	push!(essential_bcs, FDataDict("displacement"=> 0.0, "component"=> 2, "node_list"=>y0nl))
	yLnl = selectnode(fens, box = [-Inf Inf L L -Inf Inf], inflate = tolerance);
	push!(essential_bcs, FDataDict("displacement"=> 0.0, "component"=> 2, "node_list"=>yLnl))

	enl = connectednodes(subset(boundaryfes, ebcl))
	push!(essential_bcs, FDataDict("displacement"=> 0.0, "component"=> 1, "node_list"=>enl))
	push!(essential_bcs, FDataDict("displacement"=> 0.0, "component"=> 2, "node_list"=>enl))
	push!(essential_bcs, FDataDict("displacement"=> 0.0, "component"=> 3, "node_list"=>enl))
	inl = connectednodes(subset(boundaryfes, ibcl))
	push!(essential_bcs, FDataDict("displacement"=> 0.0, "component"=> 1, "node_list"=>inl))
	push!(essential_bcs, FDataDict("displacement"=> 0.0, "component"=> 2, "node_list"=>inl))
	push!(essential_bcs, FDataDict("displacement"=> iuz, "component"=> 3, "node_list"=>inl))

	modeldata = FDataDict("fens"=>fens,
	"regions"=>regions,
	"essential_bcs"=>essential_bcs
	)
	modeldata = AlgoDeforLinearModule.linearstatics(modeldata)

	filebase = "bushing_h8u_algo_ms"
	modeldata["postprocessing"] = FDataDict("file"=>filebase * "-nodal-", "quantity"=>:pressure, "component"=>collect(1:1), "nodevalmethod"=>:averaging, "reportat"=>:extraptrend)
	modeldata = AlgoDeforLinearModule.exportstress(modeldata)
	modeldata["postprocessing"] = FDataDict("file"=>filebase * "-elwise-", "quantity"=>:pressure, "component"=>collect(1:1), "nodevalmethod"=>:averaging, "reportat"=>:extraptrend)
	modeldata = AlgoDeforLinearModule.exportstresselementwise(modeldata)

    true
end # bushing_h8u_algo_ms

function bushing_h8_export()
    nt = 1
    for nperradius in [2, 4, 8]
    	nt=nt+1;
    	fens, fes  = Q4circlen(R, nperradius)
    	fens, fes = H8extrudeQ4(fens, fes, nt, (x, k) -> [x[1],x[2],k*thickness/nt])

    	MR = DeforModelRed3D
    	material = MatDeforElastIso(MR, E, nu)
    	femm = FEMMDeforLinearMSH8(MR, IntegDomain(fes, GaussRule(3, 2)), material)

    	boundaryfes = meshboundary(fes);
    	topl = selectelem(fens,boundaryfes, box =  [-Inf Inf -Inf Inf thickness thickness], inflate =   tolerance);
    	botl = selectelem(fens,boundaryfes, box =  [-Inf Inf -Inf Inf 0.0 0.0], inflate =   tolerance);
    	x0l = selectelem(fens,boundaryfes, box =  [0.0 0.0 -Inf Inf 0.0 thickness], inflate =   tolerance);
    	y0l = selectelem(fens,boundaryfes, box =  [-Inf Inf 0.0 0.0 0.0 thickness], inflate =   tolerance);
    	cyll = setdiff(1:count(boundaryfes),topl,botl,x0l,y0l);

    	geom = NodalField(fens.xyz)
    	u = NodalField(zeros(size(fens.xyz,1), 3)) # displacement field

    	x0nl = connectednodes(subset(boundaryfes, x0l))
    	setebc!(u,x0nl,true,1,0.0)

    	y0nl = connectednodes(subset(boundaryfes, y0l))
    	setebc!(u,y0nl,true,2,0.0)

    	cylnl = connectednodes(subset(boundaryfes, cyll))
    	setebc!(u,cylnl,true,1,0.0)
    	setebc!(u,cylnl,true,2,0.0)
    	setebc!(u,cylnl,true,3,0.0)

    	applyebc!(u)
    	numberdofs!(u)

    	cnl = selectnode(fens; box = [0 0 0 0 0 thickness], inflate  =  tolerance)

    	AE = AbaqusExporter("clcircularplatecl_h8_export_$(nperradius)");
    	HEADING(AE, "Clamped square plate with concentrated force");
    	PART(AE, "part1");
    	END_PART(AE);
    	ASSEMBLY(AE, "ASSEM1");
    	INSTANCE(AE, "INSTNC1", "PART1");
    	NODE(AE, fens.xyz);
    	ELEMENT(AE, "c3d8rh", "AllElements", 1, connasarray(femm.integdomain.fes))
    	NSET_NSET(AE, "cnl", cnl)
    	NSET_NSET(AE, "x0nl", x0nl)
    	NSET_NSET(AE, "y0nl", y0nl)
    	NSET_NSET(AE, "cylnl", cylnl)
    	ORIENTATION(AE, "GlobalOrientation", vec([1. 0 0]), vec([0 1. 0]));
    	SOLID_SECTION(AE, "elasticity", "GlobalOrientation", "AllElements", "Hourglassctl");
    	END_INSTANCE(AE);
    	END_ASSEMBLY(AE);
    	MATERIAL(AE, "elasticity")
    	ELASTIC(AE, E, nu)
    	SECTION_CONTROLS(AE, "Hourglassctl", "HOURGLASS=ENHANCED")
    	STEP_PERTURBATION_STATIC(AE)
    	BOUNDARY(AE, "ASSEM1.INSTNC1.x0nl", 1)
    	BOUNDARY(AE, "ASSEM1.INSTNC1.y0nl", 2)
    	BOUNDARY(AE, "ASSEM1.INSTNC1.cylnl", 1)
    	BOUNDARY(AE, "ASSEM1.INSTNC1.cylnl", 2)
    	BOUNDARY(AE, "ASSEM1.INSTNC1.cylnl", 3)
    	CLOAD(AE, "ASSEM1.INSTNC1.cnl", 3, Magnitude/4/length(cnl))
    	END_STEP(AE)
    	close(AE)

    	# File =  "clcircularplatecl_h8_export_$(nperradius).vtk"
    	# vtkexportmesh(File, fens, fes; vectors = [("u", u.values)])
    	# @async run(`"paraview.exe" $File`)

    end

    true
end # bushing_h8_export

function bushing_t10_algo_ms()

	fens, fes  = T10block(ang, L, Re-Ri, nc, nt, nR)
	internal_fenids = selectnode(fens, box = [0 ang 0 L 0 0], inflate = tolerance);

	MR = DeforModelRed3D
	material = MatDeforElastIso(MR, E, nu)
	femm = FEMMDeforLinearMST10(MR, IntegDomain(fes, TetRule(4)), material)
	regions = FDataDict[FDataDict("femm"=>femm)]

	boundaryfes = meshboundary(fes);
	ibcl = selectelem(fens, boundaryfes, box = [0 ang 0 L 0 0], inflate = tolerance);
	ebcl = selectelem(fens, boundaryfes, box = [0 ang 0 L Re-Ri Re-Ri], inflate = tolerance);

	for i in 1:count(fens)
		a = fens.xyz[i,1]; y = fens.xyz[i,2]; r = fens.xyz[i,3]+Ri;
		fens.xyz[i,:] .= (r*sin(a), y, r*cos(a));
	end

	essential_bcs = FDataDict[]
	x0nl = selectnode(fens, box = [0 0 -Inf Inf -Inf Inf], inflate = tolerance);
	push!(essential_bcs, FDataDict("displacement"=> 0.0, "component"=> 1, "node_list"=>x0nl))

	y0nl = selectnode(fens, box = [-Inf Inf 0 0 -Inf Inf], inflate = tolerance);
	push!(essential_bcs, FDataDict("displacement"=> 0.0, "component"=> 2, "node_list"=>y0nl))
	yLnl = selectnode(fens, box = [-Inf Inf L L -Inf Inf], inflate = tolerance);
	push!(essential_bcs, FDataDict("displacement"=> 0.0, "component"=> 2, "node_list"=>yLnl))

	enl = connectednodes(subset(boundaryfes, ebcl))
	push!(essential_bcs, FDataDict("displacement"=> 0.0, "component"=> 1, "node_list"=>enl))
	push!(essential_bcs, FDataDict("displacement"=> 0.0, "component"=> 2, "node_list"=>enl))
	push!(essential_bcs, FDataDict("displacement"=> 0.0, "component"=> 3, "node_list"=>enl))
	inl = connectednodes(subset(boundaryfes, ibcl))
	push!(essential_bcs, FDataDict("displacement"=> 0.0, "component"=> 1, "node_list"=>inl))
	push!(essential_bcs, FDataDict("displacement"=> 0.0, "component"=> 2, "node_list"=>inl))
	push!(essential_bcs, FDataDict("displacement"=> iuz, "component"=> 3, "node_list"=>inl))

	modeldata = FDataDict("fens"=>fens,
	"regions"=>regions,
	"essential_bcs"=>essential_bcs
	)
	modeldata = AlgoDeforLinearModule.linearstatics(modeldata)

	filebase = "bushing_T10_algo_ms"
	modeldata["postprocessing"] = FDataDict("file"=>filebase * "-nodal-", "quantity"=>:pressure, "component"=>collect(1:1), "nodevalmethod"=>:averaging, "reportat"=>:extraptrend)
	modeldata = AlgoDeforLinearModule.exportstress(modeldata)
	modeldata["postprocessing"] = FDataDict("file"=>filebase * "-elwise-", "quantity"=>:pressure, "component"=>collect(1:1), "nodevalmethod"=>:averaging, "reportat"=>:extraptrend)
	modeldata = AlgoDeforLinearModule.exportstresselementwise(modeldata)

    true
end # bushing_t10_algo_ms

function bushing_t10_algo()

	fens, fes  = T10block(ang, L, Re-Ri, nc, nt, nR)
	internal_fenids = selectnode(fens, box = [0 ang 0 L 0 0], inflate = tolerance);

	MR = DeforModelRed3D
	material = MatDeforElastIso(MR, E, nu)
	femm = FEMMDeforLinear(MR, IntegDomain(fes, TetRule(4)), material)
	regions = FDataDict[FDataDict("femm"=>femm)]

	boundaryfes = meshboundary(fes);
	ibcl = selectelem(fens, boundaryfes, box = [0 ang 0 L 0 0], inflate = tolerance);
	ebcl = selectelem(fens, boundaryfes, box = [0 ang 0 L Re-Ri Re-Ri], inflate = tolerance);

	for i in 1:count(fens)
		a = fens.xyz[i,1]; y = fens.xyz[i,2]; r = fens.xyz[i,3]+Ri;
		fens.xyz[i,:] .= (r*sin(a), y, r*cos(a));
	end

	essential_bcs = FDataDict[]
	x0nl = selectnode(fens, box = [0 0 -Inf Inf -Inf Inf], inflate = tolerance);
	push!(essential_bcs, FDataDict("displacement"=> 0.0, "component"=> 1, "node_list"=>x0nl))

	y0nl = selectnode(fens, box = [-Inf Inf 0 0 -Inf Inf], inflate = tolerance);
	push!(essential_bcs, FDataDict("displacement"=> 0.0, "component"=> 2, "node_list"=>y0nl))
	yLnl = selectnode(fens, box = [-Inf Inf L L -Inf Inf], inflate = tolerance);
	push!(essential_bcs, FDataDict("displacement"=> 0.0, "component"=> 2, "node_list"=>yLnl))

	enl = connectednodes(subset(boundaryfes, ebcl))
	push!(essential_bcs, FDataDict("displacement"=> 0.0, "component"=> 1, "node_list"=>enl))
	push!(essential_bcs, FDataDict("displacement"=> 0.0, "component"=> 2, "node_list"=>enl))
	push!(essential_bcs, FDataDict("displacement"=> 0.0, "component"=> 3, "node_list"=>enl))
	inl = connectednodes(subset(boundaryfes, ibcl))
	push!(essential_bcs, FDataDict("displacement"=> 0.0, "component"=> 1, "node_list"=>inl))
	push!(essential_bcs, FDataDict("displacement"=> 0.0, "component"=> 2, "node_list"=>inl))
	push!(essential_bcs, FDataDict("displacement"=> iuz, "component"=> 3, "node_list"=>inl))

	modeldata = FDataDict("fens"=>fens,
	"regions"=>regions,
	"essential_bcs"=>essential_bcs
	)
	modeldata = AlgoDeforLinearModule.linearstatics(modeldata)

	filebase = "bushing_T10_algo"
	modeldata["postprocessing"] = FDataDict("file"=>filebase * "-nodal-", "quantity"=>:pressure, "component"=>collect(1:1), "nodevalmethod"=>:averaging, "reportat"=>:extraptrend)
	modeldata = AlgoDeforLinearModule.exportstress(modeldata)
	modeldata["postprocessing"] = FDataDict("file"=>filebase * "-elwise-", "quantity"=>:pressure, "component"=>collect(1:1), "nodevalmethod"=>:averaging, "reportat"=>:extraptrend)
	modeldata = AlgoDeforLinearModule.exportstresselementwise(modeldata)

    true
end # bushing_t10_algo_ms

#     true
# end # twisted_beam_export

function allrun()
    println("#####################################################")
    println("# bushing_h8_full ")
    bushing_h8_full()
    println("#####################################################")
    println("# bushing_h8_algo_full ")
    bushing_h8_algo_full()
    println("#####################################################")
    println("# bushing_h8_algo_ms ")
    bushing_h8_algo_ms()
    println("#####################################################")
    println("# bushing_h8u_algo_ms ")
    bushing_h8u_algo_ms()
    println("#####################################################")
    println("# bushing_t10_algo_ms")
    bushing_t10_algo_ms()
    return true
end # function allrun

end # module bushing_examples


text:
function M_fast(model::ContactModel, q)
    return model.base.M(q)
end

function B_fast(model::ContactModel, q)
    return model.base.B(q)
end

function A_fast(model::ContactModel, q)
    return model.base.A(q)
end

# function J_fast(model::ContactModel, q)
#     return model.base.J(q)
# end

function C_fast(model::ContactModel, q, q̇)
    return model.base.C(q, q̇)
end

# Dynamics functions
function d_fast(model::ContactModel, h, q0, q1, u1, w1, Λ1, q2)
    return model.dyn.d(h, q0, q1, u1, w1, Λ1, q2)
end


text:
include(joinpath(@__DIR__, "../src/BrainFuck.jl"))
using Test

include("BrainFuck.jl")

max_stars_repo_name: stellartux/Esolangs

max_issues_repo_name: stellartux/Esolangs

max_forks_repo_name: stellartux/Esolangs


text:
# imbibition with wettability change
# solve transport for one species in the aqueous phase
# change the SF accordingly
# plot the production curve
import Plots
include("forced_imbibition_IMPES.jl")
sw0=0.1
sw=0.1
Plots.plot(size=(500,400), xtickfont = Plots.font(10, "Courier"), ytickfont=Plots.font(10, "Courier"),
  xlabel="t [s]", ylabel="Recovery factor", legendfont=Plots.font(10, "Courier"),
  guidefont=Plots.font(12, "Courier"))
t_plot=0.0
R_plot=0.0
for SF0 in 0.1:0.2:0.9
  if sw!=sw0
    sw0=sw.value[2:end-1]
  end
  t,R,sw=forced_imb_impes(SF0, sw0)
  t_plot=t_plot[end]+t
  R_plot=R_plot[end]+R*(1-R_plot[end])
  Plots.plot!(t_plot, R_plot, ylims=(0.0,1.0), linewidth=2,
  label="$SF0 Sal. Factor")
end
Plots.plot!()
Plots.savefig("forced_imb_multiple_sal.png")

max_stars_repo_name: simulkade/peteng

max_issues_repo_name: simulkade/peteng

max_forks_repo_name: simulkade/peteng


text:
nt = 1000
seed = 99
m = ACQR_simula(nt,seed)
y = m["y"]
nx = 4

#######################
print("Testing ACQR_em_run: ")
mest = ACQR_em_run(y,nx,max_iter=10,txo=false)
if sum(mest["A"].^2) > 1e-6
    println("OK")
end


text:
style_graph = create_graph()
v = GraphVertex(1, 1:2)
w = GraphVertex(2, 1:2)

function check_vertex_shape()
    t1 = typeof(@test_throws ArgumentError vertex_shape(:rectangle)) == Test.Pass
    t2 = typeof(@test_throws ArgumentError vertex_shape(:circle)) == Test.Pass

    style1 = vertex_shape(:rectangle; dimensions=(1, 2))
    style2 = vertex_shape(:circle; dimensions=Tuple(1))

    l_v = length(v.opts[:styles])
    l_w = length(w.opts[:styles])

    @add_styles v [style1]
    @add_styles w [style2]

    render(Javis.CURRENT_VIDEO[1];tempdirectory = "images", pathname = "")

    t3 = length(v.opts[:styles]) - l_v == 1 && v.object.opts[:shape] == :rectangle
    t4 = length(w.opts[:styles]) - l_w == 1 && w.object.opts[:shape] == :circle
    return t1 && t2 && t3 && t4
end

function check_vertex_text()
    style1 = vertex_text("text")
    style2 = vertex_text(L"text")

    l_v = length(v.opts[:styles])
    l_w = length(w.opts[:styles])

    @add_styles v [style1]
    @add_styles w [style2]

    render(Javis.CURRENT_VIDEO[1];tempdirectory = "images", pathname = "")

    t1 = length(v.opts[:styles]) - l_v == 1 && haskey(v.object.opts, :text_box)
    t2 = length(w.opts[:styles]) - l_w == 1 && haskey(w.object.opts, :text_box)

    return t1 && t2
end

function check_vertex_fill()
    style1 = vertex_fill(:image, "julia-logo.png")
    style2 = vertex_fill(:color, "red")

    l_v = length(v.opts[:styles])
    l_w = length(w.opts[:styles])

    @add_styles v [style1]
    @add_styles w [style2]

    render(Javis.CURRENT_VIDEO[1];tempdirectory = "images", pathname = "")

    t1 = length(v.opts[:styles]) - l_v == 1 && haskey(v.object.opts, :fill_type)
    t2 = length(w.opts[:styles]) - l_w == 1 && haskey(w.object.opts, :fill_type)

    return t1 && t2
end

function check_vertex_border()
    style = vertex_border()

    l_v = length(v.opts[:styles])

    @add_styles v [style]

    render(Javis.CURRENT_VIDEO[1];tempdirectory = "images", pathname = "")

    t1 = length(v.opts[:styles]) - l_v == 1 && haskey(v.object.opts, :border_width)

    return t1
end

e = GraphEdge(1, 2, 1:2)
f = GraphEdge(2, 1, 1:2)
g = GraphEdge(2, 2, 1:2)

function check_edge_shape()
    style1 = edge_shape(end_offsets=(1, 1))
    style2 = edge_shape(:curved)
    style3 = edge_shape(:curved, end_offsets=(1, 1))

    l_e = length(e.opts[:styles])
    l_f = length(f.opts[:styles])
    l_g = length(g.opts[:styles])

    @add_styles e [style1]
    @add_styles f [style2]
    @add_styles g [style3]

    render(Javis.CURRENT_VIDEO[1];tempdirectory = "images", pathname = "")

    t1 = length(e.opts[:styles]) - l_e == 1 && haskey(e.object.opts, :shape)
    t2 = length(f.opts[:styles]) - l_f == 1 && f.object.opts[:shape] == :curved
    t3 = length(g.opts[:styles]) - l_g == 1 && g.object.opts[:end_offsets] == (1, 1)
    return t1 && t2 && t3
end

function check_edge_style()
    style = edge_style()

    l_e = length(e.opts[:styles])

    @add_styles e [style]

    render(Javis.CURRENT_VIDEO[1];tempdirectory = "images", pathname = "")

    t1 = length(e.opts[:styles]) - l_e == 1 && haskey(e.object.opts, :color)
    return t1
end

function check_edge_arrow()
    style1 = edge_arrow(start=true, finish=true)
    style2 = edge_arrow(finish=true)
    style3 = edge_arrow(start=true)

    l_e = length(e.opts[:styles])
    l_f = length(f.opts[:styles])
    l_g = length(g.opts[:styles])

    @add_styles e [style1]
    @add_styles f [style2]
    @add_styles g [style3]

    render(Javis.CURRENT_VIDEO[1];tempdirectory = "images", pathname = "")

    t1 = length(e.opts[:styles]) - l_e == 1 && haskey(e.object.opts, :arrow_color)
    t2 = length(f.opts[:styles]) - l_f == 1 && f.object.opts[:arrow_finish] == true
    t3 = length(g.opts[:styles]) - l_g == 1 && g.object.opts[:arrow_start] == true
    return t1 && t2
end

function check_edge_text()
    style1 = edge_text("test")

    l_e = length(e.opts[:styles])

    @add_styles e [style1]

    render(Javis.CURRENT_VIDEO[1];tempdirectory = "images", pathname = "")

    t1 = length(e.opts[:styles]) - l_e == 1 && haskey(e.object.opts, :text_content)
    return t1
end


text:
# Archive tests

#
state_history = [UncertainDiscreteState(1.0*idx, 1.0) for idx in 1:10]
times, states, sigmas =
    StateEstimation.vectorize_state_history(state_history)
@test length(times) == 10


text:
module daemon

end

max_stars_repo_name: tmaone/xsh

max_issues_repo_name: turing-complete/xsh

max_forks_repo_name: turing-complete/xsh


text:
"""
Identical to `identity()` used for 1-d distances that do not need to be transformed
before being sorted
"""
function identiy_distance(distances::Array{G, 1}) where G <: Real
    return identity(distances)
end

function identiy_distance(distances::Array{G, 2}) where G <: Real
    return identity(distances)
end

"""
Function to return the permutation of the computed distances that will order from
smallest to largest.  In this case of multiple distance functions a cross dimention
comparison is made by summing the ranks for each dimention and then ordering the
sum from smallest to largest

# Arguments
* `distances::Array{Float64, 1}, SharedArray{Float64, 1}, Array{Array{Float64, 1}, 1}, SharedArray{Array{Float64, 1}, 1}`:
    vector of distances or vector of distance vectors

# Value
permutation of indicies so such that distances are ordered from smallest to largest
"""
function sumrank_distance(distances::Array{G, 2}) where G <: Real
    n, m = size(distances)
    rankArray = Array{Float64}(n, m)
    for jj in 1:m
        rankArray[:, jj] = sortperm(distances[:, jj])
    end
    totals = vec(sum(rankArray, 2))
    return totals
end

"""
Maximum rank seen in each row when ranks are compute for columns
`maxrank_distance(distances)`
A cross dimention comparison is made by ranking each observation across each
dimension and then assigning each observation the maximum across all dimentions.

# Arguments
* `distances::Array{Float64, 2}`: vector of distances or vector of distance vectors

# Value
Vector containing the maximum rank across each row

# Examples
`distances = rand(10,3)`
`maxrank_distance(distances)`
"""
function maxrank_distance(distances::Array{G, 2}) where G <: Real
    n, m = size(distances)
    rankArray = Array{Float64}(n, m)
    for jj in 1:m
        rankArray[:, jj] = sortperm(distances[:, jj])
    end
    totals = vec(maximum(rankArray, 2))
    return totals
end

"""
Normalized compound distance
`normrank_distance(distances)`

# Arguments
* `distances::Array{Float64, 2}`: vector of distances or vector of distance vectors

# Value
Vector containing sum across rows of normalized distances

# Examples
`distances = rand(10,3)`
`maxrank_distance(distances)`
"""
function normrank_distance(distances::Array{G, 2}) where G <: Real
    distmin = minimum(distances, 1)
    diststd = std(distances, 1)
    scaledDistances = (distances .- distmin) ./ diststd
    totals = vec(sum(scaledDistances))
    return totals
end

"""
Efficiently find the points across two distance metrics with exactly N entries <= both
"""
function minarea_distance(x::Array{G, 1}, y::Array{G, 1}, N::Integer) where G <: Real
    permx = sortperm(x)
    rankx = invperm(permx)
    permy = sortperm(y)
    ranky = invperm(permy)
    yranks = falses(y)
    for ii in 1:(N - 1)
        yranks[ranky[permx[ii]]] = true
    end
    ymaxrank = length(yranks)

    xdist = x[permx[N:end]]
    ydist = zeros(xdist)
    #ydist[1] = y[permy[ymaxrank]]
    for (ii, idx) in enumerate(permx[N:end])
        if ranky[idx] < ymaxrank
            yranks[ranky[idx]] = true
            ymaxrank = findprev(yranks, ymaxrank - 1)
        end
        ydist[ii] = y[permy[ymaxrank]]
    end
    return [xdist ydist]
end
#=
x = rand(1000)
y = rand(1000)

permx = sortperm(x)
permy = sortperm(y)

nless = Array{Int64}(1000, 1000)
x10 = Array{Float64}(0)
y10 = Array{Float64}(0)

@time for ii in 1:1000
    for jj in 1:1000
        nless[ii, jj] = sum(((x .<= x[ii]) .* (y .<= y[jj])))
        if nless[ii, jj] == 10
            push!(x10, x[ii])
            push!(y10, y[jj])
        end
    end
end

xyarea = x10 .* y10
extrema(xyarea)

permx[10]
permy[end]
nless[281, 412]

(x[permx[10]], y[permy[end]])

using RCall

R"plot($x, $y)"
R"points($x10, $y10, col = 'blue')"
R"hist($xyarea)"
R"boxplot($xyarea)"
=#


text:
# Analyzing a wooden puzzle I have at home and working generating additional
# puzzles from it

# Each stick has up to 6 groves, 3 on each sides. I number these groves as
#
#  6 5 4
# ‾V‾‾‾V‾
# _Λ_Λ___
#  3 2 1
#
# (U+203E : Overline, U+039B / U+03BB: Lambda)
#
# And then can compactly encode each stick as binary number. This last stick is
# 0b101110 = 46. Since there are some rotations we can do, this same stick can
# be represented by 29, 43, 46 or 63. I take the lowest number (29) to be the
# canonical representation of the stick.

using JuMP, GLPK

function notch_list(repr::Int)::Vector{Int}
    d = digits(repr, base=2, pad=6)
    @assert length(d) == 6 "The integer representation $repr is not well formed."
    d
end

function orient_notch_list(ls, orientation)
    if orientation == 1
        ls
    elseif orientation == 2
        [ls[3], ls[2], ls[1], ls[6], ls[5], ls[4]]
    elseif orientation == 3
        [ls[4], ls[5], ls[6], ls[1], ls[2], ls[3]]
    else
        [ls[6], ls[5], ls[4], ls[3], ls[2], ls[1]]
    end
end

draw_top_notch(isNotch::Int) = isNotch == 1 ? "V" : "‾"
draw_bot_notch(isNotch::Int) = isNotch == 1 ? "Λ" : "_"

"""Use println(draw_stick(29)) to draw a single stick to the console."""
function draw_stick(repr::Int; orientation=1)::String
    d = orient_notch_list(notch_list(repr), orientation)
    """‾$(draw_top_notch(d[6]))‾$(draw_top_notch(d[5]))‾$(draw_top_notch(d[4]))‾
    _$(draw_bot_notch(d[3]))_$(draw_bot_notch(d[2]))_$(draw_bot_notch(d[1]))_"""
end

notch_count(repr::Int)::Int = count_ones(repr)

# List of all the sticks that are possible.
# Note that sticks 0, 1, 2, 3, 5, 7 are terminal.
canonical_sticks = [0, 1, 2, 3, 5, 7, 9, 10, 11, 12, 13, 14, 15, 18, 19, 21, 23, 27, 29, 30, 31, 45, 47, 63]
@assert length(canonical_sticks) == 24

"""Check for horizontal symmetry (i.e. up stays up)"""
function is_h_symmetric(repr::Int)::Bool
    d = notch_list(repr)
    d[1] == d[3] && d[4] == d[6]
end

"""Check for vertical symmetry (i.e. up becomes down)"""
function is_v_symmetric(repr::Int)::Bool
    d = notch_list(repr)
    d[1] == d[4] && d[2] == d[5] && d[3] == d[6]
end

"""Check for rotational symmetry. Note that having any two symmetries implies the third one."""
function is_hv_symmetric(repr::Int)::Bool
    d = notch_list(repr)
    d[1] == d[6] && d[2] == d[5] && d[3] == d[4]
end

@assert is_h_symmetric(45)
@assert is_h_symmetric(21)
@assert !is_h_symmetric(9)

@assert is_v_symmetric(45)
@assert !is_v_symmetric(21)
@assert is_v_symmetric(9)

function reorient_index(orientation::Int, index::Int)::Int
    if orientation == 1
        index
    elseif orientation == 2 # h flip
        Dict(1 => 3, 2 => 2, 3 => 1, 4 => 6, 5 => 5, 6 => 4)[index]
    elseif orientation == 3 # v flip
        Dict(1 => 4, 2 => 5, 3 => 6, 4 => 1, 5 => 2, 6 => 3)[index]
    else # hv flip
        Dict(1 => 6, 2 => 5, 3 => 4, 4 => 3, 5 => 2, 6 => 1)[index]
    end
    end

function notch_at(repr::Int, orientation::Int, index::Int)::Int
    notch_list(repr)[reorient_index(orientation, index)]
end

@assert notch_at(2, 1, 1) == 0
@assert notch_at(23, 3, 6) == 1
@assert notch_at(12, 2, 1) == 1
@assert notch_at(11, 4, 6) == 1

# Here is a simple example I worked out. The marbles in the bottom form the shape
#
# V V ‾ | 1 2 3 | Here each marble is a ‾, as it works as a not-grove.
# V ‾ V | 4 5 6 | Meaning this setup has the representation 0b010101011 = 171
# ‾ V ‾ | 7 8 9 | where representations range from 0 to 511 (inclusive).
#
# On this we will stack two layers, made out of [ 1, 2, 7, 10, 13, 13 ] where the
# 13 is included twice. But we'll place it in two different orientations.

marbles_repr = 171
marbles = digits(marbles_repr, base=2, pad=9)
@assert marbles == [1, 1, 0, 1, 0, 1, 0, 1, 0]

sticks = [1, 2, 7, 10, 13]
stick_count = [1, 1, 1, 1, 2]
layers = 2
@assert length(sticks) == length(stick_count)
@assert sum(stick_count) >= 3 * layers

# Setting up a model to optimize
m = Model(GLPK.Optimizer)

# Variable name is (stick)(position)(orientation=[0, h, v, hv])
@variable(m, spo[1:length(sticks), 1:6, 1:4], Bin)

# for all symmetric sticks, force the symmetry partners to be zero.
for i in 1:length(sticks)
    if is_h_symmetric(sticks[i])
        if is_v_symmetric(sticks[i])
            # Fully symmetric, set all but one orientation to 0
            @constraint(m, spo[i, 1:end, 2:end] .== 0)
        else
            # Only h symmetric, but not v symmetric. This means we eliminate
            # orientations 2 and 4.
            @constraint(m, spo[i, 1:end, 2] .== 0)
            @constraint(m, spo[i, 1:end, 4] .== 0)
        end
    elseif is_v_symmetric(sticks[i])
        # Only v symmetric, but not h symmetric. We eliminate orientations
        # 3 and 4.
        @constraint(m, spo[i, 1:end, 3] .== 0)
        @constraint(m, spo[i, 1:end, 4] .== 0)
    elseif is_hv_symmetric(sticks[i])
        # Only rotation symmetric, not mirror symmetric. We still eliminate
        # 3 and 4 but now they are paired differently than before.
        @constraint(m, spo[i, 1:end, 3] .== 0)
        @constraint(m, spo[i, 1:end, 4] .== 0)
    end
    end

# Constrain that every stick is only used at most `count` times.
@constraint(m, sum(spo, dims=(2, 3)) .<= stick_count )

# Constraint that every position contains at most one stick.
@constraint(m, sum(spo, dims=(1, 3)) .<= 1)

# The objective is to place as many sticks as possible.
@objective(m, Max, sum(spo))

# Now we need add the condition for sticks properly sitting on each other.
# This condition depends on the actual sticks that are in the mix and is quite
# heterogeneous.
# We'll also need an initial and terminal condition.

# sum over all sticks and all orientations, filter by "Has a notch".
function notch_sum(sticks, spo, notch_index, position)
    sum(notch_at.(sticks, transpose(1:4), notch_index) .* spo[1:end, position, 1:4])
end

# Initial layer: Marbles + first layer.

# Positions that are "at the bottom of" the touchpoints.
bottom_positions = nothing
# Positions that are "on top of" the touchpoints.
top_positions = collect(1:3)

@constraint(m, marbles[1] + notch_sum(sticks, spo, 1, top_positions[1]) == 1)
@constraint(m, marbles[2] + notch_sum(sticks, spo, 2, top_positions[1]) == 1)
@constraint(m, marbles[3] + notch_sum(sticks, spo, 3, top_positions[1]) == 1)
@constraint(m, marbles[4] + notch_sum(sticks, spo, 1, top_positions[2]) == 1)
@constraint(m, marbles[5] + notch_sum(sticks, spo, 2, top_positions[2]) == 1)
@constraint(m, marbles[6] + notch_sum(sticks, spo, 3, top_positions[2]) == 1)
@constraint(m, marbles[7] + notch_sum(sticks, spo, 1, top_positions[3]) == 1)
@constraint(m, marbles[8] + notch_sum(sticks, spo, 2, top_positions[3]) == 1)
@constraint(m, marbles[9] + notch_sum(sticks, spo, 3, top_positions[3]) == 1)

# Layer 2: first layer + second layer

bottom_positions = collect(1:3)
top_positions = collect(4:6)

@constraint(m, notch_sum(sticks, spo, 4, bottom_positions[1]) + notch_sum(sticks, spo, 1, top_positions[1]) == 1)
@constraint(m, notch_sum(sticks, spo, 5, bottom_positions[1]) + notch_sum(sticks, spo, 1, top_positions[2]) == 1)
@constraint(m, notch_sum(sticks, spo, 6, bottom_positions[1]) + notch_sum(sticks, spo, 1, top_positions[3]) == 1)
@constraint(m, notch_sum(sticks, spo, 4, bottom_positions[2]) + notch_sum(sticks, spo, 2, top_positions[1]) == 1)
@constraint(m, notch_sum(sticks, spo, 5, bottom_positions[2]) + notch_sum(sticks, spo, 2, top_positions[2]) == 1)
@constraint(m, notch_sum(sticks, spo, 6, bottom_positions[2]) + notch_sum(sticks, spo, 2, top_positions[3]) == 1)
@constraint(m, notch_sum(sticks, spo, 4, bottom_positions[3]) + notch_sum(sticks, spo, 3, top_positions[1]) == 1)
@constraint(m, notch_sum(sticks, spo, 5, bottom_positions[3]) + notch_sum(sticks, spo, 3, top_positions[2]) == 1)
@constraint(m, notch_sum(sticks, spo, 6, bottom_positions[3]) + notch_sum(sticks, spo, 3, top_positions[3]) == 1)

# Layer 3: last layer must be smooth at the top.

bottom_positions = collect(4:6)
top_positions = nothing

@constraint(m, notch_sum(sticks, spo, 4, bottom_positions[1]) == 0)
@constraint(m, notch_sum(sticks, spo, 5, bottom_positions[1]) == 0)
@constraint(m, notch_sum(sticks, spo, 6, bottom_positions[1]) == 0)
@constraint(m, notch_sum(sticks, spo, 4, bottom_positions[2]) == 0)
@constraint(m, notch_sum(sticks, spo, 5, bottom_positions[2]) == 0)
@constraint(m, notch_sum(sticks, spo, 6, bottom_positions[2]) == 0)
@constraint(m, notch_sum(sticks, spo, 4, bottom_positions[3]) == 0)
@constraint(m, notch_sum(sticks, spo, 5, bottom_positions[3]) == 0)
@constraint(m, notch_sum(sticks, spo, 6, bottom_positions[3]) == 0)

# Optimizing the model.
optimize!(m)
objective_value(m)

values = value.(spo)
for p in 1:6
    s, o = Tuple(argmax(values[1:end, p, 1:4]))
    println(draw_stick(sticks[s]; orientation=o))
end

max_stars_repo_name: kreibaum/stick-puzzle

max_issues_repo_name: kreibaum/stick-puzzle

max_forks_repo_name: kreibaum/stick-puzzle


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.6.2
#     language: julia
#     name: julia-1.6
# ---

# %%
using Base.Threads
using ProgressMeter
using Plots
nthreads()
xoxp1(x) = x/(x+1) # {x > 0} -> {0 < y < 1}

# %%
using DynamicalSystems
using StaticArrays

# %%
function f(x, p)
   r, k = p
   a, mu, d = 5.0, 0.5, 0.2
   x1new = x[1]*exp(r*(1-x[1]/k)-x[2]/(a+x[1]^2))
   x2new = x[2]*exp(mu*x[1]/(a+x[1]^2)-d)
   x1new, x2new
end

f(x, p, t) = SVector(f(x, p))

# %%
function traj2d!(tr1, tr2, f, p, x0, niters, nwarmup)
    x = x0
    for _ in 1:nwarmup
        x = f(x, p)
    end
    @inbounds tr1[1], tr2[1] = x
    for i in 2:niters+1
        x = f(x, p)
        @inbounds tr1[i], tr2[i] = x
    end
end

# %%
function traj2d_ds!(tr1, tr2, dds, niters, nwarmup)
    integ = integrator(dds, dds.u0)
    nwarmup != 0 && step!(integ, nwarmup)
    @inbounds tr1[1], tr2[1] = DynamicalSystems.DynamicalSystemsBase.obtain_access(integ.u, nothing)
    for i in 2:niters+1
        step!(integ)
        @inbounds tr1[i], tr2[i] = DynamicalSystems.DynamicalSystemsBase.obtain_access(integ.u, nothing)
    end
end

# %%
niters = 10
nwarmup = 10
tr1 = zeros(Int, niters+1)
tr2 = similar(tr1)
traj2d!(tr1, tr2, ((a, b), p) -> (b, a+b), nothing, (0, 1), niters, nwarmup)
@show tr1 tr2;

# %%
function findperiod(x; tol=1e-3)
    n = length(x)
    @inbounds for k in 2:(n ÷ 2 + 1)
        if abs(x[k] - x[1]) ≤ tol
            all(j -> abs(x[j] - x[j-k+1]) ≤ tol, k:n) && return k - 1
        end
    end
    return n
end

# %%
x = [1, 2, 3, 4, 1, 2, 3]
findperiod(x)

# %%
x = [1, 2, 3, 4, 1, 2, 3, 4]
findperiod(x)

# %%
x = [1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3]
findperiod(x)

# %%
x = @. sinpi(0:0.1:10) + 1e-4randn()
findperiod(x)

# %%
niters = 2000
nwarmup = 50000
tr1_ds = zeros(niters+1)
tr2_ds = similar(tr1)
x0 = (0.4, 0.5)
p = (2.5, 4.0)
dds = DiscreteDynamicalSystem(f, SVector(x0), collect(p))
traj2d_ds!(tr1_ds, tr2_ds, dds, niters, nwarmup)
@show tr1_ds[1:10];

# %%
niters = 2000
nwarmup = 50000
tr1 = zeros(niters+1)
tr2 = similar(tr1)
x0 = (0.4, 0.5)
p = (2.5, 4.0)
traj2d!(tr1, tr2, f, p, x0, niters, nwarmup)
@show tr1[1:10];

# %%
tr1_ds == tr1

# %%
@btime traj2d!($tr1, $tr2, f, p, x0, niters, nwarmup)
@btime traj2d_ds!($tr1_ds, $tr2_ds, dds, niters, nwarmup);

# %%
function calcperiod2d!(tr1, tr2, f, p, x0, niters, nwarmup; tol=1e-3)
    traj2d!(tr1, tr2, f, p, x0, niters, nwarmup)
    findperiod(tr1; tol)
end

# %%
using BenchmarkTools
@btime calcperiod2d!($(zeros(2001)), $(zeros(2001)), f, (2.6, 4.0), (0.4, 0.5), 2000, 50000; tol=1e-3)

# %%
function calcperiod2d(f = f, niters = 2000, nwarmup = 50000, x0 = (0.4, 0.5),
        r = (1:0.05:5)[1:end-1], k = (2:0.05:5)[1:end-1], tol=1e-3
    )
    tr1 = Vector{Float64}(undef, niters+1)
    tr2 = similar(tr1)
    period = Matrix{Int}(undef, length(r), length(k))
    prog = Progress(length(r)*length(k))
    for j in eachindex(k)
        for i in eachindex(r)
            @inbounds period[i, j] = calcperiod2d!(tr1, tr2, f, (r[i], k[j]), x0, niters, nwarmup; tol)
            next!(prog)
        end
    end
    period, r, k
end

period, r, k = @time calcperiod2d()
heatmap(r, k, xoxp1.(period'); xlabel="r", ylabel="k", title="period/(period + 1)")

# %%
function calcperiod2d_threads(f = f, niters = 2000, nwarmup = 50000, x0 = (0.4, 0.5),
        r = (1:0.05:5)[1:end-1], k = (2:0.05:5)[1:end-1], tol=1e-3
    )
    tr1 = [Vector{Float64}(undef, niters+1) for _ in 1:nthreads()]
    tr2 = [Vector{Float64}(undef, niters+1) for _ in 1:nthreads()]
    period = Matrix{Int}(undef, length(r), length(k))
    #prog = Progress(length(r)*length(k))
    @threads for j in eachindex(k)
        tid = threadid()
        for i in eachindex(r)
            @inbounds period[i, j] = calcperiod2d!(tr1[tid], tr2[tid], f, (r[i], k[j]), x0, niters, nwarmup; tol)
            #next!(prog)
        end
    end
    period, r, k
end

# %%
period, r, k_th, r, k = @time calcperiod2d_threads(); flush(stdout)
period, r, k_th, r, k = @time calcperiod2d_threads(); flush(stdout)
period, r, k_th, r, k = @time calcperiod2d_threads(); flush(stdout)
heatmap(r, k, xoxp1.(period_th'); xlabel="r", ylabel="k", title="period/(period + 1)")

# %%
period, r, k_th = @time calcperiod2d_threads(); flush(stdout)
period, r, k_th = @time calcperiod2d_threads(); flush(stdout)
period, r, k_th = @time calcperiod2d_threads(); flush(stdout)
heatmap(r, k, xoxp1.(period_th'); xlabel="r", ylabel="k", title="period/(period + 1)")

# %%
period == period_th

# %%
period'

# %%
@btime calcperiod2d_threads(f; niters = 2000, nwarmup = 50000, x0 = (0.4, 0.5),
    r = (1:0.05:5)[1:end-1], k = (2:0.05:5)[1:end-1], tol=1e-3)

# %%
@code_warntype calcperiod2d_threads(f; niters = 2000, nwarmup = 50000, x0 = (0.4, 0.5),
    r = (1:0.05:5)[1:end-1], k = (2:0.05:5)[1:end-1], tol=1e-3)

# %%

max_stars_repo_name: genkuroki/public

max_issues_repo_name: genkuroki/public

max_forks_repo_name: genkuroki/public


text:
@testset "nl_write" begin
    # Turn on debug mode so files persist
    old_debug = ReSHOP.CONFIG[:debug]
    ReSHOP.setdebug(true)

#    filename = "test"
#    filepath = joinpath(ReSHOP.solverdata_dir, filename)
#    ReSHOP.clean_solverdata()

#    context("all temp files deleted successfully") do
#        @fact length(readdir(ReSHOP.solverdata_dir)) --> 1
#    end

    m = Model(solver=ReSHOP.ReSHOPSolver())
    @variable(m, x >= 0)
    @objective(m, Min, x)
    solve(m)

    # Reset debug mode and clean up
    ReSHOP.setdebug(old_debug)
    ReSHOP.clean_solverdata()
end


text:
#=
Created on Friday 4 December 2020
Last update: Monday 27 September 2021

@author: Bram De Jaegher
bram.de.jaegher@gmail.com

Templates heavily based on the MIT course "Computational Thinking"

https://computationalthinking.mit.edu/
=#

correct(text=rand(yays)) = MD(Admonition("correct", "Got it!", [text]))
keep_working(text=md"The answer is not quite right.") = MD(Admonition("danger", "Keep working on it!", [text]))
not_defined(variable_name) = MD(Admonition("danger", "Oopsie!", [md"Make sure that you define a variable called **$(Markdown.Code(string(variable_name)))**"]))

yays = [md"Great!", md"Yay ❤", md"Great! 🎉", md"Well done!",
        md"Keep it up!", md"Good job!", md"Awesome!", md"You got the right answer!",
				md"Let's move on to the next section."]

hint(text) = Markdown.MD(Markdown.Admonition("hint", "Hint", [text]))
still_missing(text=md"Replace `missing` with your answer.") = MD(Admonition("warning", "Here we go!", [text]))

fyi(text) = Markdown.MD(
	Markdown.Admonition("info",
		"Additional info",
		[fyi_css,
			text
		]
	)
)

max_stars_repo_name: Beramos/PlutoTeaching

max_issues_repo_name: Beramos/PlutoTeaching

max_forks_repo_name: Beramos/PlutoTeaching


text:
# A transformation for testing the local affine mapping
immutable SquareMe <: Transformation; end
(::SquareMe)(x) = x.^2
CoordinateTransformations.transform_deriv(::SquareMe, x0) = diagm(2*x0)

@testset "Common Transformations" begin
    @testset "AffineMap" begin
        @testset "Simple" begin
            M = [1 2; 3 4]
            v = [-1, 1]
            x = [1,0]
            y = [0.0,1.0]
            A = AffineMap(M,v)
            @test A(x) == M*x + v
        end

        @testset "composition " begin
            M1 = [1 2; 3 4]
            v1 = [-1, 1]
            A1 = AffineMap(M1,v1)
            M2 = [0 1; 1 0]
            v2 = [-2, 0]
            A2 = AffineMap(M2,v2)
            x = [1,0]
            y = [0,1]
            @test A1(x) == M1*x + v1
            @test A1(y) == M1*y + v1
            @test (A2∘A1)(x) == M2*(M1*x + v1) + v2
            @test (A2∘A1)(y) == M2*(M1*y + v1) + v2
        end

        @testset "inverse" begin
            M = [1.0 2.0; 3.0 4.0]
            v = [-1.0, 1.0]
            x = [1,0]
            y = [0.0,1.0]
            A = AffineMap(M,v)
            @test inv(A)(A(x)) ≈ x
            @test inv(A)(A(y)) ≈ y
        end

        @testset "Affine approximation" begin
            S = SquareMe()
            x0 = [1,2,3]
            dx = 0.1*[1,-1,1]
            A = AffineMap(S, x0)
            @test isapprox(S(x0 + dx), A(x0 + dx), atol=maximum(2*dx.^2))
        end
    end

    @testset "LinearMap" begin
        M = [1 2; 3 4]
        x = [1,0]
        y = [0,1]
        L = LinearMap(M)
        @test L(x) == M*x
        @test inv(L)(x) == inv(M)*x
        @test inv(L)(y) == inv(M)*y
        @test inv(L)(L(x)) ≈ x
        @test inv(L)(L(y)) ≈ y
        @test (L∘L)(x) == (M*M)*x
    end

    @testset "Translation" begin
        x = SVector(1.0, 2.0)
        trans = Translation(2.0, -1.0)

        # Inverse
        @test inv(trans) == Translation(-2.0, 1.0)

        # Composition
        @test trans ∘ trans == Translation(4.0, -2.0)

        # Transform
        @test trans(x) === SVector(3.0, 1.0)
        @test trans(collect(x)) === SVector(3.0, 1.0)

        # Transform derivative
        m1 = transform_deriv(trans, x)
        m2 = [m1[i,j] for i=1:2, j=1:2] # m1 might be a UniformScaling
        @test m2 == eye(2)

        # Transform parameter derivative
        m1 = transform_deriv_params(trans, x)
        m2 = [m1[i,j] for i=1:2, j=1:2] # m1 might be a UniformScaling
        @test m2 == eye(2)
    end

#=
    @testset "Rotation2D on Polar" begin
        p = Polar(2.0, 1.0)
        trans = Rotation2D(1.0)

        # Inverse
        @test inv(trans) == Rotation2D(-1.0)

        # Composition
        @test trans ∘ trans == Rotation2D(2.0)

        # Transform
        @test trans(p) == Polar(2.0, 2.0)

        # Transform derivative
        @test transform_deriv(trans, p) == [0 0; 0 1]

        # Transform parameter derivative
        @test transform_deriv_params(trans, p) == [0; 1]
    end

    @testset "Rotation2D" begin
        x = SVector(2.0, 0.0)
        x2 = CartesianFromPolar()(Polar(2.0, 1.0))
        trans = Rotation2D(1.0)

        # Constructor
        @test trans.cos == cos(1.0)
        @test trans.sin == sin(1.0)

        # Inverse
        @test inv(trans) == Rotation2D(-1.0)

        # Composition
        @test trans ∘ trans == Rotation2D(2.0)

        # Transform
        @test trans(x) ≈ x2
        @test SVector(trans(Tuple(x))) ≈ x2
        @test trans(collect(x)) ≈ collect(x2)

        # Transform derivative
        x = SVector(2.0,1.0)
        x_gn = SVector(Dual(2.0, (1.0,0.0)), Dual(1.0, (0.0,1.0)))
        x2_gn = trans(x_gn)
        m_gn = @SMatrix [partials(x2_gn[1], 1) partials(x2_gn[1], 2);
                     partials(x2_gn[2], 1) partials(x2_gn[2], 2) ]
        m = transform_deriv(trans, x)
        @test m ≈ m_gn

        # Transform parameter derivative
        trans_gn = Rotation2D(Dual(1.0, (1.0)))
        x = SVector(2.0,1.0)
        x2_gn = trans_gn(x)
        m_gn = Mat(partials(x2_gn[1], 1), partials(x2_gn[2], 1))
        m = transform_deriv_params(trans, x)
        @test m ≈ m_gn
    end

    @testset "Rotation (3D)" begin

        @testset "Rotation{Void} (rotation matrix parameterization)" begin
            θx = 0.1
            θy = 0.2
            θz = 0.3

            sx = sin(θx)
            cx = cos(θx)
            sy = sin(θy)
            cy = cos(θy)
            sz = sin(θz)
            cz = cos(θz)

            Rx = [ cx  sx  0;
                  -sx  cx  0;
                    0   0  1]

            Ry = [1   0  0 ;
                  0  cy -sy;
                  0  sy  cy]

            Rz = [cx  0 -sx;
                  0   1  0 ;
                  sx  0  cx]

            R = Mat{3,3,Float64}(Rx*Ry*Rz)

            x = SVector(1.0, 2.0, 3.0)
            trans = Rotation(R)

            @test inv(trans).matrix ≈ R'
            @test (trans ∘ trans).matrix ≈ R*R

            y = trans(x)
            @test y == R * SVector(1.0, 2.0, 3.0)
            @test trans(Tuple(x)) == Tuple(R * SVector(1.0, 2.0, 3.0))
            @test trans(collect(x)) == R * SVector(1.0, 2.0, 3.0)

            x_gn = SVector(Dual(1.0,(1.,0.,0.)), Dual(2.0,(0.,1.,0.)), Dual(3.0,(0.,0.,1.)))
            y_gn = trans(x_gn)
            M_gn = Mat{3,3,Float64}(vcat(ntuple(i->[partials(y_gn[i], j) for j = 1:3].', 3)...))
            M = transform_deriv(trans, x)
            @test M ≈ M_gn

            g11 = Dual(R[1,1],(1.,0.,0.,0.,0.,0.,0.,0.,0.))
            g12 = Dual(R[1,2],(0.,1.,0.,0.,0.,0.,0.,0.,0.))
            g13 = Dual(R[1,3],(0.,0.,1.,0.,0.,0.,0.,0.,0.))
            g21 = Dual(R[2,1],(0.,0.,0.,1.,0.,0.,0.,0.,0.))
            g22 = Dual(R[2,2],(0.,0.,0.,0.,1.,0.,0.,0.,0.))
            g23 = Dual(R[2,3],(0.,0.,0.,0.,0.,1.,0.,0.,0.))
            g31 = Dual(R[3,1],(0.,0.,0.,0.,0.,0.,1.,0.,0.))
            g32 = Dual(R[3,2],(0.,0.,0.,0.,0.,0.,0.,1.,0.))
            g33 = Dual(R[3,3],(0.,0.,0.,0.,0.,0.,0.,0.,1.))

            G = @SMatrix [g11 g12 g13;
                      g21 g22 g23;
                      g31 g32 g33 ]

            trans_gn = Rotation(G)
            y_gn = trans_gn(x)

            M_gn = Mat{3,9,Float64}(vcat(ntuple(i->[partials(y_gn[i], j) for j = 1:9].', 3)...))
            M = transform_deriv_params(trans, x)
            @test M ≈ M_gn
        end

        @testset "Rotation{Quaternion} (quaternion parameterization)" begin
            v = [0.35, 0.45, 0.25, 0.15]
            v = v / vecnorm(v)
            q = Quaternion(v[1],v[2],v[3],v[4],true)

            trans = Rotation(q)
            x = SVector(1.0, 2.0, 3.0)

            @test inv(trans) ≈ Rotation(inv(Quaternion(v[1],v[2],v[3],v[4])))
            @test trans ∘ trans ≈ Rotation(trans.matrix * trans.matrix)

            y = trans(x)
            @test y ≈ SVector(3.439024390243902,-1.1463414634146332,0.9268292682926829)

            x_gn = SVector(Dual(1.0,(1.,0.,0.)), Dual(2.0,(0.,1.,0.)), Dual(3.0,(0.,0.,1.)))
            y_gn = trans(x_gn)
            M_gn = Mat{3,3,Float64}(vcat(ntuple(i->[partials(y_gn[i], j) for j = 1:3].', 3)...))
            M = transform_deriv(trans, x)
            @test M ≈ M_gn

            v_gn = [Dual(v[1],(1.,0.,0.,0.)), Dual(v[2],(0.,1.,0.,0.)), Dual(v[3],(0.,0.,1.,0.)), Dual(v[4],(0.,0.,0.,1.))]
            q_gn = Quaternion(v_gn[1],v_gn[2],v_gn[3],v_gn[4],true)
            trans_gn = Rotation(q_gn)
            y_gn = trans_gn(x)
            M_gn = Mat{3,4,Float64}(vcat(ntuple(i->[partials(y_gn[i], j) for j = 1:4].', 3)...))
            M = transform_deriv_params(trans, x)
            # Project both to tangent plane of normalized quaternions (there seems to be a change in definition...)
            proj = Mat{4,4,Float64}(eye(4) - v*v')
            @test M*proj ≈ M_gn*proj
        end

        @testset "Rotation{EulerAngles} (Euler angle parameterization)" begin
            θx = 0.1
            θy = 0.2
            θz = 0.3

            trans = Rotation(EulerAngles(θx, θy, θz))
            x = SVector(1.0, 2.0, 3.0)

            @test inv(trans) == Rotation(trans.matrix')
            @test trans ∘ trans == Rotation(trans.matrix * trans.matrix)

            y = trans(x)
            @test y == SVector(0.9984766744283545,2.1054173473736495,2.92750100324502)

            x_gn = SVector(Dual(1.0,(1.,0.,0.)), Dual(2.0,(0.,1.,0.)), Dual(3.0,(0.,0.,1.)))
            y_gn = trans(x_gn)
            M_gn = Mat{3,3,Float64}(vcat(ntuple(i->[partials(y_gn[i], j) for j = 1:3].', 3)...))
            M = transform_deriv(trans, x)
            @test M ≈ M_gn

            # Parameter derivative not defined
            # TODO?
        end

        @testset "RotationXY, RotationYZ and RotationZX" begin
            # RotationXY
            x = SVector(2.0, 0.0, 0.0)
            x2 = CartesianFromSpherical()(Spherical(2.0, 1.0, 0.0))
            trans = RotationXY(1.0)

            # Constructor
            @test trans.cos == cos(1.0)
            @test trans.sin == sin(1.0)

            # Inverse
            @test inv(trans) == RotationXY(-1.0)
            @test RotationYX(1.0) == RotationXY(-1.0)

            # Composition
            @test trans ∘ trans == RotationXY(2.0)

            # Transform
            @test trans(x) ≈ x2
            @test SVector(trans(Tuple(x))) ≈ x2
            @test SVector(trans(collect(x))) ≈ x2

            # Transform derivative
            x = SVector(2.0,1.0,3.0)
            x_gn = SVector(Dual(2.0, (1.0,0.0,0.0)), Dual(1.0, (0.0,1.0,0.0)), Dual(3.0, (0.0,0.0,1.0)))
            x2_gn = trans(x_gn)
            m_gn = @SMatrix [partials(x2_gn[1], 1) partials(x2_gn[1], 2) partials(x2_gn[1], 3);
                         partials(x2_gn[2], 1) partials(x2_gn[2], 2) partials(x2_gn[2], 3);
                         partials(x2_gn[3], 1) partials(x2_gn[3], 2) partials(x2_gn[3], 3) ]
            m = transform_deriv(trans, x)
            @test m ≈ m_gn

            # Transform parameter derivative
            trans_gn = RotationXY(Dual(1.0, (1.0)))
            x = SVector(2.0,1.0,3.0)
            x2_gn = trans_gn(x)
            m_gn = Mat(partials(x2_gn[1], 1), partials(x2_gn[2], 1), partials(x2_gn[3], 1))
            m = transform_deriv_params(trans, x)
            @test m ≈ m_gn

            # RotationYZ
            x = SVector(0.0, 2.0, 0.0)
            x2 = CartesianFromSpherical()(Spherical(2.0, pi/2, 1.0))
            trans = RotationYZ(1.0)

            # Constructor
            @test trans.cos == cos(1.0)
            @test trans.sin == sin(1.0)

            # Inverse
            @test inv(trans) == RotationYZ(-1.0)
            @test RotationZY(1.0) == RotationYZ(-1.0)

            # Composition
            @test trans ∘ trans == RotationYZ(2.0)

            # Transform
            @test trans(x) ≈ x2
            @test SVector(trans(Tuple(x))) ≈ x2
            @test SVector(trans(collect(x))) ≈ x2

            # Transform derivative
            x = SVector(2.0,1.0,3.0)
            x_gn = SVector(Dual(2.0, (1.0,0.0,0.0)), Dual(1.0, (0.0,1.0,0.0)), Dual(3.0, (0.0,0.0,1.0)))
            x2_gn = trans(x_gn)
            m_gn = @SMatrix [partials(x2_gn[1], 1) partials(x2_gn[1], 2) partials(x2_gn[1], 3);
                         partials(x2_gn[2], 1) partials(x2_gn[2], 2) partials(x2_gn[2], 3);
                         partials(x2_gn[3], 1) partials(x2_gn[3], 2) partials(x2_gn[3], 3) ]
            m = transform_deriv(trans, x)
            @test m ≈ m_gn

            # Transform parameter derivative
            trans_gn = RotationYZ(Dual(1.0, (1.0)))
            x = SVector(2.0,1.0,3.0)
            x2_gn = trans_gn(x)
            m_gn = Mat(partials(x2_gn[1], 1), partials(x2_gn[2], 1), partials(x2_gn[3], 1))
            m = transform_deriv_params(trans, x)
            @test m ≈ m_gn

            # RotationZX
            x = SVector(2.0, 0.0, 0.0)
            x2 = CartesianFromSpherical()(Spherical(2.0, 0.0, -1.0))
            trans = RotationZX(1.0)

            # Constructor
            @test trans.cos == cos(1.0)
            @test trans.sin == sin(1.0)

            # Inverse
            @test inv(trans) == RotationZX(-1.0)
            @test RotationXZ(1.0) == RotationZX(-1.0)

            # Composition
            @test trans ∘ trans == RotationZX(2.0)

            # Transform
            @test trans(x) ≈ x2
            @test SVector(trans(Tuple(x))) ≈ x2
            @test SVector(trans(collect(x))) ≈ x2

            # Transform derivative
            x = SVector(2.0,1.0,3.0)
            x_gn = SVector(Dual(2.0, (1.0,0.0,0.0)), Dual(1.0, (0.0,1.0,0.0)), Dual(3.0, (0.0,0.0,1.0)))
            x2_gn = trans(x_gn)
            m_gn = @SMatrix [partials(x2_gn[1], 1) partials(x2_gn[1], 2) partials(x2_gn[1], 3);
                         partials(x2_gn[2], 1) partials(x2_gn[2], 2) partials(x2_gn[2], 3);
                         partials(x2_gn[3], 1) partials(x2_gn[3], 2) partials(x2_gn[3], 3) ]
            m = transform_deriv(trans, x)
            @test m ≈ m_gn

            # Transform parameter derivative
            trans_gn = RotationZX(Dual(1.0, (1.0)))
            x = SVector(2.0,1.0,3.0)
            x2_gn = trans_gn(x)
            m_gn = Mat(partials(x2_gn[1], 1), partials(x2_gn[2], 1), partials(x2_gn[3], 1))
            m = transform_deriv_params(trans, x)
            @test m ≈ m_gn
        end

        @testset "euler_rotation() and composed derivatives" begin
            x = SVector(2.0,1.0,3.0)
            trans = euler_rotation(0.1,0.2,0.3)
            x2 = SVector(2.730537054338937,0.8047190852558106,2.428290466296628)

            #@test trans.t1.t1 == RotationXY(0.1)
            #@test trans.t1.t2 == RotationYZ(0.2)
            #@test trans.t2 == RotationZX(0.3)

            @test inv(trans) ≈ RotationZX(-0.3) ∘ (RotationYZ(-0.2) ∘ RotationXY(-0.1))

            @test trans(x) ≈ x2

            # Transform derivative
            x = SVector(2.0,1.0,3.0)
            x_gn = SVector(Dual(2.0, (1.0,0.0,0.0)), Dual(1.0, (0.0,1.0,0.0)), Dual(3.0, (0.0,0.0,1.0)))
            x2_gn = trans(x_gn)
            m_gn = @SMatrix [partials(x2_gn[1], 1) partials(x2_gn[1], 2) partials(x2_gn[1], 3);
                         partials(x2_gn[2], 1) partials(x2_gn[2], 2) partials(x2_gn[2], 3);
                         partials(x2_gn[3], 1) partials(x2_gn[3], 2) partials(x2_gn[3], 3) ]
            m = transform_deriv(trans, x)
            @test m ≈ m_gn

            # Transform parameter derivative

            #trans_gn = euler_rotation(Dual(0.1, (1.0, 0.0, 0.0)), Dual(0.2, (0.0, 1.0, 0.0)), Dual(0.3, (0.0, 0.0, 1.0)))
            #x = SVector(2.0,1.0,3.0)
            #x2_gn = trans_gn(x)
            #m_gn = @SMatrix [partials(x2_gn[1], 1) partials(x2_gn[1], 2) partials(x2_gn[1], 3);
            #             partials(x2_gn[2], 1) partials(x2_gn[2], 2) partials(x2_gn[2], 3);
            #             partials(x2_gn[3], 1) partials(x2_gn[3], 2) partials(x2_gn[3], 3)]
            #m = transform_deriv_params(trans, x)
            #@test m ≈ m_gn

        end
    end
=#
end


text:
function print_base(energy_totals::Dict)

                println("month", " ", "consumer_energy", " ",
                                "PV_energy", " ",
                                "injection_grid", " ",
                                "grid_energy", " ",
                                "global_allowance", " ",
                                "allowance", " ",
                                "withdrawn_energy", " ",
                                "utility_supplied_energy", " ",
                                "carryover"," ",
                                "peak_power", " ",
                                "peak_demand")

        for i in 1:12

                println(i, " ",round(energy_totals[i]["consumer_energy"], digits=2), " ",
                                round(energy_totals[i]["PV_energy"], digits =2), " ",
                                round(energy_totals[i]["injection_grid"], digits = 2), " ",
                                round(energy_totals[i]["grid_energy"], digits=2), " ",
                                round(energy_totals[i]["global_allowance"], digits=2), " ",
                                round(energy_totals[i]["allowance"],digits=2), " ",
                                round(energy_totals[i]["withdrawn_energy"],digits=2), " ",
                                round(energy_totals[i]["utility_supplied_energy"], digits=2), " ",
                                round(energy_totals[i]["carry_over"], digits=2), " ",
                                round(energy_totals[i]["peak_power"], digits=2), " ",
                                round(energy_totals[i]["peak_demand"], digits=2)
                                )

     end

end

function print_TMT(energy_totals::Dict)

                println("month", " | ","consumer_energy_peak", " ",
                                "consumer_energy_valley", " ",
                                "consumer_energy_night", " ",
                                " | ",
                                "PV_energy", " ",
                                "total_injection", " ",
								"global_allowance", " ",
                                "withdrawn_energy", " ",
                                " | ",
                                "grid_energy_peak", " ",
                                "grid_energy_valley", " ",
                                "grid_energy_night", " ",
                                " | ",
                                "global_generation"," ",
                                "global_withdrawal", " ",
                                "carry_over", " ",
								" | ",
								"peak_power_peak", " ",
                                "peak_power_valley", " ",
                                "peak_power_night", " ",
                                " | ",
                                "peak_demand_peak", " ",
                                "peak_demand_valley", " ",
                                "peak_demand_night", " ",
                                " | ",)

     for i in 1:12
                println(i," | ",round(energy_totals[i]["consumer_energy_peak"], digits=2), " ",
                                round(energy_totals[i]["consumer_energy_valley"], digits=2), " ",
                                round(energy_totals[i]["consumer_energy_night"], digits=2), " ",
                                " | ",
                                round(energy_totals[i]["PV_energy"], digits =2), " ",
                                round(energy_totals[i]["total_injection"], digits = 2), " ",
								round(energy_totals[i]["global_allowance"], digits=2), " ",
                                round(energy_totals[i]["withdrawn_energy"],digits=2), " ",
                                " | ",
                                round(energy_totals[i]["grid_energy_peak"], digits=2), " ",
                                round(energy_totals[i]["grid_energy_valley"], digits=2), " ",
                                round(energy_totals[i]["grid_energy_night"], digits=2), " ",
                                " | ",
                                round(energy_totals[i]["global_generation"], digits=2)," ",
                                round(energy_totals[i]["global_withdrawal"], digits=2), " ",
                                round(energy_totals[i]["carry_over"], digits=2), " ",
								" | ",
								round(energy_totals[i]["peak_power_peak"], digits=2), " ",
                                round(energy_totals[i]["peak_power_valley"], digits=2), " ",
                                round(energy_totals[i]["peak_power_night"], digits=2), " ",
                                " | ",
                                round(energy_totals[i]["peak_demand_peak"], digits=2), " ",
                                round(energy_totals[i]["peak_demand_valley"], digits=2), " ",
                                round(energy_totals[i]["peak_demand_night"], digits=2), " ",)
        end
end

# Billings printout

function print_TMT_bills()

            println(round(bill["grid_energy_cost_peak"], digits=2), " ",
                    round(bill["grid_energy_cost_valley"], digits=2), " ",
                    round(bill["grid_energy_cost_night"], digits=2), " ",
                    " | ",
                    round(bill["grid_demand_cost_peak"], digits=2), " ",
                    round(bill["grid_demand_cost_valley"], digits=2), " ",
                    round(bill["grid_demand_cost_night"], digits=2), " ",
                    " | ",
                    round(bill["withdrawn_energy_cost"], digits =2), " ",
                    " | ",
                    round(bill["savings_energy_cost_peak"], digits=2)," ",
                    round(bill["savings_energy_cost_valley"], digits=2)," ",
                    round(bill["savings_energy_cost_night"], digits=2)," ",
                    " | ",
                    round(bill["savings_demand_cost_peak"], digits=2)," ",
                    round(bill["savings_demand_cost_valley"], digits=2)," ",
                    round(bill["savings_demand_cost_night"], digits=2))

end

stripChar = (s, r) -> replace(s, Regex("[$r]") => " ")

function plot_bill(bill::Dict, fields::Array{String}, digits::Int64 = 1)
	cum_sum = Array{Float64,1}(undef, 12)
    for i in 1:length(fields)
        var = [bill[m]["$(fields[i])"] for m in 1:12]
        bar(collect(1:12), bottom = cum_sum, var, label=stripChar("$(fields[i])", "_"), edgecolor = "black", linewidth = 0.3)
        cum_sum += var
    end
    legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.);
    xticks(collect(1:12));
    y_ul = maximum([round(bill[m]["counterfactual_cost"], RoundUp, sigdigits = digits) for m in 1:12])
	xlabel("Month")
    ylabel("Total Bill [Colones]")
	ylim(0,y_ul)
	title("Consumer Monthly Electric Utility Bill Breakdown")
end

function plot_bill(bill1::Dict, bill2::Dict, fields::Array{String}, digits::Int64 = 1)
    cum_sum1 = Array{Float64,1}(undef, 12)
    cum_sum2 = Array{Float64,1}(undef, 12)

	for i in 1:length(fields)
        var1 = [bill1[m]["$(fields[i])"] for m in 1:12]
        bar(collect(1:12) .- 0.03, bottom = cum_sum1, var1, color = cpalette10[i], edgecolor = "black", linewidth = 0.3, label=stripChar("$(fields[i])", "_"), align="edge", width= -0.4)
        cum_sum1 += var1
        var2 = [bill2[m]["$(fields[i])"] for m in 1:12]
        bar(collect(1:12) .+ 0.03, bottom = cum_sum2, var2, color = cpalette10[i], edgecolor = "black", linewidth = 0.3, align="edge", width= 0.4)
        cum_sum2 += var2
    end
    legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.);
    y_ul1 = maximum([round(bill1[m]["counterfactual_cost"], RoundUp, sigdigits = digits) for m in 1:12])
	y_ul2 = maximum([round(bill2[m]["counterfactual_cost"], RoundUp, sigdigits = digits) for m in 1:12])
	xlabel("Month")
    ylabel("Total Bill [Colones]")
	ylim(0,max(y_ul1,y_ul2))
	title("Consumer's Change in Montly Electric Bill")
end

function plot_savings(bill::Dict, digits::Int64 = 1)
	cum_sum = Array{Float64,1}(undef, 12)
	for i in keys(bill[1]["savings"])
		var = [bill[m]["savings"][i] for m in 1:12]
		bar(collect(1:12), bottom = cum_sum, var, label=stripChar(i, "_"), edgecolor = "black", linewidth = 0.3)
		cum_sum += var
	end
	var1 = [bill[m]["withdrawn_energy_cost"] for m in 1:12]
	bar(collect(1:12), -1*var1, label="withdrawn energy cost", edgecolor = "black", linewidth = 0.3)
	legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.);
	xticks(collect(1:12));
	y_ul = maximum([round(bill[m]["counterfactual_cost"], RoundUp, sigdigits = digits) for m in 1:12])
	y_lb = minimum(-1*var1)
	xlabel("Month")
    ylabel("Total Savings [Colones]")
	ylim(y_lb,y_ul)
	axhline(0, color="black", lw=0.6)
	title("Consumer Monthly Electric Utility Bill Savings")
end

function plot_utility_change(bill::Dict, digits::Int64 = 1)
	cum_sum = Array{Float64,1}(undef, 12)
	for i in keys(bill[1]["utility_change"])
		if !(i in ["consumer_savings", "losses"])
            var = [bill[m]["utility_change"][i] for m in 1:12]
            sum(var .< 0) > 0  ? @info("Losses for energy displacement") : true
			bar(collect(1:12), bottom = cum_sum, var, label=stripChar(i, "_"), edgecolor = "black", linewidth = 0.3)
            cum_sum += var
        end
	end
	var1 = [bill[m]["utility_change"]["consumer_savings"] for m in 1:12]
	bar(collect(1:12), var1, label="consumer savings", edgecolor = "black", linewidth = 0.3)
    var2 = [bill[m]["utility_change"]["losses"] for m in 1:12]
    plot(collect(1:12), var2, label = "net losses", color="red")
	legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.);
	xticks(collect(1:12));
	y_ul = maximum([round(bill[m]["counterfactual_cost"], RoundUp, sigdigits = digits) for m in 1:12])
    y_lb = minimum(var1)
	xlabel("Month")
    ylabel("Total Change for Utility [Colones]")
	ylim(y_lb,y_ul)
	axhline(0, color="black", lw=0.6)
	title("Utility Monthly Revenue Change")
end

max_stars_repo_name: arnavgautam/PVCR

max_issues_repo_name: arnavgautam/PVCR

max_forks_repo_name: arnavgautam/PVCR


text:
## for symbolic functions (dsolve.jl)
(u::SymFunction)(x::Base.Dict) = throw(ArgumentError("IVPsolutions can only be called with symbolic objects"))
(u::SymFunction)(x::Base.Pair) = throw(ArgumentError("IVPsolutions can only be called with symbolic objects"))
function (u::SymFunction)(x)
    if u.n == 0
        PyObject(u)(x) #u.x(PyObject(x))
    else
        __x = Sym("__x")
        diff(u.x(PyObject(__x)), __x, u.n)(__x => x)
    end
end

(u::SymFunction)(x, y...) = u.n== 0 ? u.x(map(PyObject, vcat(x, y...))...) : error("Need to implement derivatives of symbolic functions of two or more variables")

## Call symbolic object with natural syntax
## ex(x=>val)
## how to do from any symbolic object?
(ex::Sym)() = ex
function (ex::Sym)(args...)
    xs = free_symbols(ex)
    if length(xs) >= 1
        subs(ex, collect(zip(xs, args))...)
    else
        if classname(ex) == "Symbol"
            ex
        else
            pycall(PyObject(ex), PyAny, args...)
        end
    end
end
(ex::Sym)(x::Dict) = subs(ex, x)
(ex::Sym)(x::Pair...) = subs(ex, x...)


text:
pystr_fromUTF8(x::Ptr, n::Integer) = pynew(errcheck(C.PyUnicode_DecodeUTF8(x, n, C_NULL)))
pystr_fromUTF8(x) = pystr_fromUTF8(pointer(x), sizeof(x))

"""
    pystr(x)

Convert `x` to a Python `str`.
"""
pystr(x) = pynew(errcheck(@autopy x C.PyObject_Str(getptr(x_))))
pystr(x::String) = pystr_fromUTF8(x)
pystr(x::SubString{String}) = pystr_fromUTF8(x)
pystr(x::Char) = pystr(string(x))
pystr(::Type{String}, x) = (s=pystr(x); ans=pystr_asstring(s); pydel!(s); ans)
export pystr

pystr_asUTF8bytes(x::Py) = GC.@preserve x pynew(errcheck(C.PyUnicode_AsUTF8String(getptr(x))))
pystr_asUTF8vector(x::Py) = (b=pystr_asUTF8bytes(x); ans=pybytes_asvector(b); pydel!(b); ans)
pystr_asstring(x::Py) = (b=pystr_asUTF8bytes(x); ans=pybytes_asUTF8string(b); pydel!(b); ans)

function pystr_intern!(x::Py)
    ptr = Ref(getptr(x))
    C.PyUnicode_InternInPlace(ptr)
    setptr!(x, ptr[])
end

pyconvert_rule_str(::Type{String}, x::Py) = pyconvert_return(pystr_asstring(x))
pyconvert_rule_str(::Type{Symbol}, x::Py) = pyconvert_return(Symbol(pystr_asstring(x)))
pyconvert_rule_str(::Type{Char}, x::Py) = begin
    s = pystr_asstring(x)
    if length(s) == 1
        pyconvert_return(first(s))
    else
        pyconvert_unconverted()
    end
end

pyisstr(x) = pytypecheckfast(x, C.Py_TPFLAGS_UNICODE_SUBCLASS)

pyconvert_rule_fast(::Type{String}, x::Py) = pyisstr(x) ? pyconvert_rule_str(String, x) : pyconvert_unconverted()
pyconvert_rule_fast(::Type{Symbol}, x::Py) = pyisstr(x) ? pyconvert_rule_str(Symbol, x) : pyconvert_unconverted()
pyconvert_rule_fast(::Type{Char}, x::Py) = pyisstr(x) ? pyconvert_rule_str(Char, x) : pyconvert_unconverted()


text:
using MiniFB

# This function displays a time-varying mix of colors on screen. Demonstrates how to create a buffer and render it to a window.
function plasma()
    pallete = zeros(UInt32, 512)
    WIDTH = 320
    HEIGHT = 240
    inc = 90 / 64;
    for c in 1:64
        col = round(Int, ((255 * sin( (c-1) * inc * π / 180)) + 0.5));
        pallete[64*0 + c] = mfb_rgb(col,     0,       0);
        pallete[64*1 + c] = mfb_rgb(255,     col,     0);
        pallete[64*2 + c] = mfb_rgb(255-col, 255,     0);
        pallete[64*3 + c] = mfb_rgb(0,       255,     col);
        pallete[64*4 + c] = mfb_rgb(0,       255-col, 255);
        pallete[64*5 + c] = mfb_rgb(col,     0,       255);
        pallete[64*6 + c] = mfb_rgb(255,     0,       255-col);
        pallete[64*7 + c] = mfb_rgb(255-col, 0,       0);
    end
    window = mfb_open_ex("Plasma Test", WIDTH, HEIGHT, MiniFB.WF_RESIZABLE);
    g_buffer = zeros(UInt32, WIDTH * HEIGHT)
    mfb_set_target_fps(10);
    time=0
    while mfb_wait_sync(window)
        time_x = sin(time * π / 180);
        time_y = cos(time * π / 180);
        i = 1;
        for y in 1:HEIGHT
            dy = cos((y * time_y) * π / 180);
            for x in 1:WIDTH
                dx = sin((x * time_x) * π / 180);
                idx = round(Int, ((2 + dx + dy) * 0.25 * 511) + 1)
                g_buffer[i] = pallete[idx];
                i += 1
            end
        end
        time += 1
        state = mfb_update(window, g_buffer);
        if state != MiniFB.STATE_OK
            break;
        end
    end
    mfb_close(window)
end

# Call the function
plasma()

# ![](../assets/plasma.png)


text:
export VectorField, TensorField

struct VectorField{PD,T,N}
	x::Array{T,N}
	y::Array{T,N}
	z::Array{T,N}
end
function VectorField{T}(yeegrid::YeeGrid{PD,N};init::Function = zeros) where {T,PD,N}
	VectorField{PD,T,N}(init.(T,size(yeegrid))...)
end
struct TensorField{PD,T,N}
	xx::Array{T,N}
	xy::Array{T,N}
	xz::Array{T,N}
	yx::Array{T,N}
	yy::Array{T,N}
	yz::Array{T,N}
	zx::Array{T,N}
	zy::Array{T,N}
	zz::Array{T,N}
end

function TensorField{T}(yeegrid::YeeGrid{PD,N};init::Function = zeros) where {T,PD,N}
	TensorField{PD,T,N}(init.(T,size(yeegrid))...,init.(T,size(yeegrid))...,init.(T,size(yeegrid))...)
end


text:
using Distributions, Plots
using BioLizard.OrdinalRegr: fitordinalsvm, OrdinalSVM, predict, score

Σ = [3.0 1; 1 2]

dist1 = MultivariateNormal([0, 0], Σ)
dist2 = MultivariateNormal([1, 3], Σ)
dist3 = MultivariateNormal([3, 9], Σ)

X = [rand(dist1, 20) rand(dist2, 20) rand(dist3, 20)]'
y = vcat([3 for i in 1:20], [5 for i in 1:20], [8 for i in 1:20])

pcl = scatter(X[1:20,1], X[1:20,2], label="class 1")
scatter!(X[21:40,1], X[21:40,2], label="class 2")
scatter!(X[41:60,1], X[41:60,2], label="class 3")

model = fitordinalsvm(y, X)

phist = histogram(score(model, X[1:20,:]), label="class 1", alpha=0.8)
histogram!(score(model, X[21:40,:]), label="class 2", alpha=0.8)
histogram!(score(model, X[41:60,:]), label="class 3", alpha=0.8)
vline!(model.b, label="b")

plot(pcl, phist, layout=(2,1))


text:
@testset "test constructor" begin
    c = 1.0
    e = @SVector [1.0, 2.0, 3.0]
    r = 0.1

    x1 = RevisedAffineForm(c, e, r)
    @test x1.c == c
    @test x1.ϵ == e
    @test x1.Δ == r

    @test RevisedAffineForm(1.0, e, big(0.1)) isa RevisedAffineForm{3,BigFloat}

    x2 = RevisedAffineForm(1, [1, 2, 3], 0.1)
    @test x1 == x2

    x3 = RevisedAffineForm(1, 3)
    @test x3.c === 1.0
    @test iszero(x3.ϵ)
    @test iszero(x3.Δ)

    c1 = Interval(0.1, 0.3)
    x4 = RevisedAffineForm(c1, 3)
    @test x4.c == mid(c1)
    @test iszero(x4.ϵ)
    @test x4.Δ == radius(c1)

    c2 = 0.1 .. 0.1
    e2 = [0.2 .. 0.2, 0.3 .. 0.3, 0.4 .. 0.4]
    x5 = RevisedAffineForm(c2, e2, 0)
    @test x5.Δ == sup(sum(Interval.(radius.(e2))) + radius(c2))

    x6 = x1 ± 1
    @test x6.c == x1.c
    @test x6.ϵ == x1.ϵ
    @test x6.Δ == 1.1

    @test string(x1) == "1.0 + 1.0ϵ₁ + 2.0ϵ₂ + 3.0ϵ₃ ± 0.1"
end


text:
using Pipe
using Test
_macroexpand(x) = macroexpand(Main, x)

#No change to nonpipes functionality
@test _macroexpand( :(@pipe a) ) == :a #doesn't change single inputs
@test _macroexpand( :(@pipe b(a)) ) == :(b(a)) #doesn't change inputs that a function applications

#Compatable with Julia 0.3 piping functionality
@test _macroexpand( :(@pipe a|>b) ) == :(b(a)) #basic
@test _macroexpand( :(@pipe a|>b|>c) ) == :(c(b(a)))  #Keeps chaining 3
@test _macroexpand( :(@pipe a|>b|>c|>d) ) == :(d(c(b(a)))) #Keeps chaining 4

@test _macroexpand( :(@pipe a|>b(x)) ) == :((b(x))(a))  #applying to function calls returning functions
@test _macroexpand( :(@pipe a(x)|>b ) ) == :(b(a(x)))   #feeding functioncall results on wards

@test _macroexpand(:(@pipe 1|>a)) ==:(a(1)) #Works with literals (int)
@test _macroexpand(:(@pipe "foo"|>a)) == :(a("foo")) #Works with literal (string)
@test _macroexpand( :(@pipe a|>bb[2])) == :((bb[2])(a)) #Should work with RHS that is a array reference

#Marked locations
@test _macroexpand( :(@pipe a|>b(_) ) ) == :(b(a)) #Marked location only
@test _macroexpand( :(@pipe a|>b(x,_) ) ) == :(b(x,a)) # marked 2nd (and last)
@test _macroexpand( :(@pipe a|>b(_,x) ) ) == :(b(a,x)) # marked first
@test _macroexpand( :(@pipe a|>b(_,_) ) ) == :(b(a,a)) # marked double (Not certain if this is a good idea)
@test _macroexpand( :(@pipe a|>bb[2](x,_))) == :((bb[2])(x,a)) #Should work with RHS that is a array reference

#marked Unpacking
@test _macroexpand( :(@pipe a|>b(_...) ) ) == :(b(a...)) # Unpacking
@test _macroexpand( :(@pipe a|>bb[2](_...))) == :((bb[2])(a...)) #Should work with RHS of arry ref and do unpacking

#Mixing modes
@test _macroexpand( :(@pipe a|>b|>c(_) ) ) == :(c(b(a)))
@test _macroexpand( :(@pipe a|>b(x,_)|>c|>d(_,y) ) ) == :(d(c(b(x,a)),y))
@test _macroexpand( :(@pipe a|>b(xb,_)|>c|>d(_,xd)|>e(xe) |>f(xf,_,yf) ) ) == :(f(xf,(e(xe))(d(c(b(xb,a)),xd)),yf)) #Very Complex


text:
# findall

function Base.findall(pred::Function, d::AbstractDictionary)
    out = empty(keys(d))
    @inbounds for (i, v) in pairs(d)
        if pred(v)
            insert!(out, i)
        end
    end
    return out
end

# Note this implements Boolean indexing. Compare:
#     array[array_of_bools]
#     getindices(dict, findall(dict_of_bools))
@inline Base.findall(d::AbstractDictionary) = findall(identity, d)

# findfirst

# TODO what do do about findlast, findnext, findprev?
#  * findnext would require us to be able to restart iteration from an index or token
#    (this might be useful even for a hashmap, e.g. if `findall` allocates too much)
#  * findlast and findprev would require thinking about Iterators.reverse


text:
export Toeplitz, embed
"""
[`Toeplitz` matrix](http://en.wikipedia.org/wiki/Toeplitz_matrix)

Input is a vector of odd length.

```julia
julia> Toeplitz(collect(-4:4))
5x5 Toeplitz{Int64}:
 0  -1  -2  -3  -4
 1   0  -1  -2  -3
 2   1   0  -1  -2
 3   2   1   0  -1
 4   3   2   1   0
```
"""
struct Toeplitz{T} <: AbstractArray{T, 2}
	c :: Vector{T}
end

getindex(T::Toeplitz, i::Int, j::Int) = T.c[i-j+div(length(T.c)+1,2)]
isassigned(T::Toeplitz, i::Int, j::Int) = isassigned(T.c,i-j+div(length(T.c)+1,2))

size(T::Toeplitz, r::Int) = (r==1 || r==2) ? 1 + div(length(T.c),2) :
    throw(ArgumentError("Invalid dimension $r"))
size(T::Toeplitz) = size(T,1), size(T,2)

# Fast matrix x vector multiplication via embedding Toeplitz() into Circulant()
function *(A::Toeplitz{T},x::Vector{T}) where T
    n=length(A.c)
    k=div(n+1,2)
    C=Circulant([A.c[k:n];A.c[1:k-1]])
    (C*[x;zeros(T,k-1)])[1:k]
end

function mul!(y::StridedVector{T},A::Toeplitz{T},x::StridedVector{T}) where T
    n=length(A.c)
    k=div(n+1,2)
    C=Circulant([A.c[k:n];A.c[1:k-1]])
    xx=[x;zeros(T,k-1)]
    yy=mul!(similar(xx),C,xx)
    copyto!(y, 1, yy, 1, length(y))
    return y
end

function Matrix(To::Toeplitz{T}) where T
	n=size(To, 1)
	M=Array{T}(undef,n,n)
	for i=1:n
		M[i:n,i] = To.c[n:2n-i]
		M[1:i-1,i] = To.c[n-i+1:n-1]
	end
	M
end

function embed(To::Toeplitz{T}) where T
    return Circulant([To.c[div(end+1,2):end];To.c[1:div(end-1,2)]])
end


text:
import ..internal
mutable struct NodeBase
    rcl_node::Ref{internal.rcl_node_t}

    NodeBase(node_name, namespace) = begin
        rcl_node = Ref(internal.rcl_get_zero_initialized_node())
        default_options = Ref(internal.rcl_node_get_default_options())
        internal.checkcall(internal.rcl_node_init(rcl_node, node_name, namespace, default_options))

        new(rcl_node)
    end
end

NodeBase(node_name) = NodeBase(node_name, "/")
Base.finalize(rcl_node::Ref{internal.rcl_node_t}) = internal.checkcall(internal.rcl_node_fini(rcl_node); suppress_throw=true)
Base.finalize(nb::NodeBase) = finalize(nb.rcl_node)

mutable struct CallbackGroup
end

function name(nb::NodeBase)
    n = internal.rcl_node_get_name(nb.rcl_node)
    n == C_NULL && error("node name returned NULL. Is the node NULL or un-initialized?")
    unsafe_string(n)
end

function namespace(nb::NodeBase)
    n = internal.rcl_node_get_namespace(nb.rcl_node)
    n == C_NULL && error("node name returned NULL. Is the node NULL or un-initialized?")
    unsafe_string(n)
end

function context(nb::NodeBase) end
# function rcl_node_handle(nb::NodeBase) end
function CallbackGroup(nb::NodeBase) end
function default_callback_group(nb::NodeBase) end
import Base: in, isvalid
function in(cg::CallbackGroup, nb::NodeBase) end
isvalid(nb::NodeBase) = internal.rcl_node_is_valid(nb.rcl_node, C_NULL)
function callback_groups(nb::NodeBase) end
function callback_groups!(cbs, nb::NodeBase) end
function associated_with_executor_atomic(nb::NodeBase) :: Threads.Atomic{Bool} end
function get_notify_guard_condition(nb::NodeBase) end
function acquire_notify_guard_condition_lock(nb::NodeBase) end


text:
using Drawing

with(Axes(centred=true), Font(family="Sans", slant="normal", weight="normal", size=1)) do
    move(-1, -1)
    text("poop")
end


text:
using ProjectPath
using Test

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


text:
module CIFAR100

import ..Datasets.unpack

function getdata(dir::String)
    mkpath(dir)
    url = "https://www.cs.toronto.edu/~kriz/cifar-100-binary.tar.gz"
    println("Downloading $url...")
    path = download(url)
    run(unpack(path,dir,".gz",".tar"))
end

function readdata(data::Vector{UInt8})
    n = Int(length(data)/3074)
    x = Matrix{Float64}(3072, n)
    y = Matrix{Int}(2, n)
    for i = 1:n
        k = (i-1) * 3074 + 1
        y[:,i] = data[k:k+1]
        x[:,i] = data[k+2:k+3073] / 255
    end
    x = reshape(x, 32, 32, 3, n)
    x, y
end

function traindata(dir::String)
    file = joinpath(dir, "cifar-100-binary","train.bin")
    isfile(file) || getdata(dir)
    readdata(open(read,file))
end

function testdata(dir::String)
    file = joinpath(dir, "cifar-100-binary","test.bin")
    isfile(file) || getdata(dir)
    readdata(open(read,file))
end

end