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

using Markdown
using InteractiveUtils

# ╔═╡ 892fb737-6fc0-4438-b586-039a440b01f9
using LinearAlgebra, SparseArrays

# ╔═╡ 069cf190-77cc-4757-9db3-1945ddd07438
using DifferentialEquations, Plots

# ╔═╡ 7bc3e874-b7be-11eb-0792-cf4e21e8c8b2
md"# Heat equation with Proportional control

Consider a one-dimensional rod $\Omega = (0,L)$ and time $t \in (0,T).$
"

# ╔═╡ 1158de2e-b7a6-43f3-a4e6-2ef9598a4ca5
L = 0.1 # Length of rod

# ╔═╡ 30cebec3-9c70-4c6a-94a8-fdcbe6cda783
Tf = 2000.0 # Final simulation time

# ╔═╡ 93fb0f90-c454-410f-aa5e-cd8e6c055dc0
md"The rod has the physical properties
- thermal conductivity $\lambda$,
- specific heat capacity $c$,
- mass density $\rho$.

The diffusivity constant is given by $\alpha = \frac{\lambda}{c \rho}$. The heat equation is noted as

$\dot{\vartheta}(t,x) = \alpha \frac{\partial^2}{\partial x^2} \vartheta(t,x)$

for $(t,x) \in (0,T) \times \Omega$ with initial condition

$\vartheta(0,x) = \vartheta_{0}(x)$

for $x \in \overline{\Omega}$."

# ╔═╡ b294bbfc-45e7-4a15-a839-5ceb56746be7
begin
	λ = 45.0;    # Thermal conductivity
	ρ = 7800.0;  # Mass density
	cap = 480.0; # Specific heat capacitivity
end

# ╔═╡ 543b3f13-4cef-46f1-9141-172b81e7d48b
α = λ/(cap * ρ)   # Diffusivity

# ╔═╡ 1deb117d-21f7-4329-a7db-2be26e43129a
md"## Boundary conditions

Only on the right side of the rod, heat transfer

$\phi_{t}(t,x=L) = - h (\vartheta(t,x=L) - \vartheta_{amb})$

and heat radiation

$\phi_{r}(t,x=L) = - \epsilon \varrho (\vartheta(t,x=L)^4 - \vartheta_{amb}^4)$

is assumed to operate. Parameter $\vartheta_{amb}$ describes the ambient temperature, $h>0$ is called the heat transfer coefficient, $\epsilon \in (0,1)$ is called emissivity and $\varrho\approx 5.67 \cdot 10^{-8}$ is known as Stefan-Boltzmann constant. For simplicity, $k=\epsilon \cdot \varrho$.

The sum of $\phi_{t}$ and $\phi_{r}$

$\phi_{out}(t,L) = \phi_{t}(t,L) + \phi_{r}(t,L) = - h (\vartheta(t,L) - \vartheta_{amb}) - k (\vartheta(t,L)^4 - \vartheta_{amb}^4)$

is the flux from the rod to the environment."

# ╔═╡ 78292b42-a203-496c-a143-def3ff128f69
begin
	h = 10.0; # Heat transfer coefficient
	ϵ = 0.6;  # Emissivity
	sb = 5.67*10^(-8) # Stefan-Boltzmann constant
	k = ϵ * sb; 	  # Radiation coefficient
end

# ╔═╡ 0ecea33b-bf24-48eb-b4b4-bbf288e62e6b
θamb = 298.0 # Ambient temperature in Kelvin

# ╔═╡ 3a059f15-d850-468a-bb1d-f42ef5b24609
ϕout(θ) = -h * (θ - θamb)  - k*(θ^4 - θamb^4)

# ╔═╡ 7202aa44-b393-4e0b-9c34-2a9b9376194c
md"### Induced heat

On the left side of the rod, a heat source is assumed as

$\phi_{in}(t,x=0) = b ~ u(t)$

with constant $b>0$, input signal

$u(t) = K_{p} ~ e(t)$

and error

$e(t) = y_{ref}(t) - y(t).$

The output is considered as the temperature on the right side, e.g.

$y(t) = \vartheta(t, x=L)$

and the reference is as an arbitrary fixed value

$y_{ref}(t) = y_{ref} > 0.$

Thus, the natural (Robin) boundary conditions are noted as

$\lambda ~ \left.\frac{\partial}{\partial x} \vartheta(\cdot, x) \right\rvert_{x = 0} \cdot \vec{n} ~=~ \phi_{in}(t,0) = b ~ u(t)$

and

$\lambda ~  \left. \frac{\partial}{\partial x} \vartheta(\cdot, x) \right\rvert_{x = L} \cdot \vec{n} ~=~ \phi_{out}(t,x) = - h (\vartheta(t,L) - \vartheta_{amb}) - k (\vartheta(t,L)^4 - \vartheta_{amb}^4)$

with outer [normal vector](https://en.wikipedia.org/wiki/Normal_(geometry)) $\vec{n}$ on the left or right boundary. Here, the normal vector $\vec{n} = -1$ on the left boundary at $x=0$, and  $\vec{n} = 1$ on the right boundary at $x=L$."

# ╔═╡ 2775433b-d8b1-48be-9b60-08641db15803
b = 1;

# ╔═╡ 4c59909b-ffb6-4dfc-94ab-142d41deeeab
Kp = 10^3  # Proportional gain

# ╔═╡ 4f0c3671-d3bf-4e36-a1f3-8e8d0ac0e454
yref = 400.0 # Reference temperature

# ╔═╡ 23343fa0-6033-498d-b089-d43e4382e332
u_in(err) = Kp * err # input signals

# ╔═╡ daa17037-fd28-460c-97f5-826742b9a419
md"### Spatial approximation

The one-dimensional rod is discretized as a one-dimensional grid with N points and the finite discretization is noted by

$\Delta x = \frac{L}{N-1}.$

The resulting grid points have the position $x^{0} = 0$, $x^{n} = n~\Delta x$ and $x^{N-1} = L$. The second order derivative is approximated using the [Taylor series](https://en.wikipedia.org/wiki/Taylor_series) as

$\left. \frac{\partial^2 f(x)}{\partial x^2} \right\rvert_{x = \tilde{x}} \approx \frac{1}{\Delta x^2} \left[ f(\tilde{x} - \Delta x) - 2 f(\tilde{x}) + f(\tilde{x} + \Delta x)  \right].$

This [Finite Difference](https://en.wikipedia.org/wiki/Finite_difference_method) scheme is written in matrix-vector form as

$\frac{\partial^2 \vartheta(t,x)}{\partial x^2} \approx \frac{1}{\Delta x^2}
\begin{pmatrix}
-2 & 1 & 0 & \cdots &  & 0 \\
1 & -2 & 1 & 0  & \cdots & \vdots \\
0 & 1 & -2 & 1 & \ddots &  \\
 & & \ddots & \ddots & \ddots &  & \\
& & & 1 & -2 & 1 \\
& & &  & 1 & -2 \\
\end{pmatrix}
~
\begin{pmatrix}
\vartheta(t, x^{0}) \\
\vartheta(t, x^{1}) \\
\vdots \\
\vartheta(t, x^{N-1}) \\
\end{pmatrix}.$
"

# ╔═╡ 0309bf5e-4033-4832-86ed-146aea0863c0
N = 101 # Number of grid elements

# ╔═╡ 5c9ec071-6140-4af6-a418-df79c10daf5e
Δx = L/(N-1) # Finite discretization

# ╔═╡ 905337c5-9eee-4b67-b215-bfd7ca905e4c
# Diffusion matrix
M = spdiagm(-1 => ones(N-1), 0 => -2*ones(N), 1 => ones(N-1));

# ╔═╡ f0e9f197-303b-48a6-a540-9b11513050d4
Matrix(M)[1:5,1:5]

# ╔═╡ 0e155080-2cfe-45c0-8166-ba91f95d6f8b
md"### Approximated boundary conditions

The spatial approximation of the heat equation

$\frac{\partial^2 \vartheta(t,x)}{\partial x^2} \approx \frac{1}{\Delta x^2} \left( \vartheta(\cdot, x^{n-1}) - 2 ~ \vartheta(\cdot, x^{n}) + \vartheta(\cdot, x^{n+1}) \right)$

can not be evaluated directly at the grid points $x^{0}$ and $x^{N-1}$ because they depend on values from the not-existing grid points $x^{-1}$ and $x^{N}$. Therefore, the Neumann boundary condition is used to find $\vartheta(\cdot, x^{-1})$ and $\vartheta(\cdot, x^{N})$.

The Robin boundary condition states that the heat flux at both sides is described by

$\lambda \left. \frac{\partial \vartheta(\cdot,x)}{\partial x}  \right\rvert_{x = 0} \cdot \vec{n} = -1 \cdot \lambda \frac{\vartheta(\cdot, x^{1}) - \vartheta(\cdot, x^{-1})}{2 \Delta x}   = \phi_{in}(t,x^{0})$

on the left boundary and

$\lambda \left. \frac{\partial \vartheta(\cdot,x)}{\partial x}  \right\rvert_{x = L} \cdot \vec{n} = \lambda \frac{\vartheta(\cdot, x^{N}) - \vartheta(\cdot, x^{N-2})}{2 \Delta x}  = \phi_{out}(t,x^{N-1})$

on the right boundary. Reformulating both equations lead to

$\vartheta(\cdot, x^{-1}) = \vartheta(\cdot, x^{1}) + \frac{2 \Delta x}{\lambda} ~ \phi_{in}(t,x^{0})$

and

$\vartheta(\cdot, x^{N}) = \vartheta(\cdot, x^{N-2})  + \frac{2 \Delta x}{\lambda} ~ \phi_{out}(t,x^{N-1}).$

On the left boundary one yields

$\left. \frac{\partial^2 \vartheta(t,x)}{\partial x^2} \right\rvert_{x = 0} \approx \frac{1}{\Delta x^2} \left( - 2 ~ \vartheta(\cdot, x^{0}) + 2 ~ \vartheta(\cdot, x^{1}) \right) + \frac{2}{\lambda ~ \Delta x} \phi_{in}(t,x^{0})$

and analog on the right boundary

$\left. \frac{\partial^2 \vartheta(t,x)}{\partial x^2} \right\rvert_{x = 0} \approx \frac{1}{\Delta x^2} \left(2 ~ \vartheta(\cdot, x^{N-2}) - 2 ~ \vartheta(\cdot, x^{N-1}) \right) + \frac{2}{\lambda ~ \Delta x} \phi_{out}(t,x^{N-1}).$

The induced and emitted flux are noted as

$\Phi_{in}(t) =
\begin{pmatrix}
b  \\
0 \\
\vdots \\
0 \\
\end{pmatrix}
~ u(t)

\quad \text{and} \quad

\Phi_{out}(t) =
\begin{pmatrix}
0 \\
\vdots \\
0 \\
1 \\
\end{pmatrix}
~\phi_{out}(t, x^{N-1})$

and the diffusion matrix is described by

$M =
\begin{pmatrix}
-2 & 2 & 0 & \cdots &  & 0 \\
1 & -2 & 1 & 0  & \cdots & \vdots \\
0 & 1 & -2 & 1 & \ddots &  \\
 & & \ddots & \ddots & \ddots &  & \\
& & & 1 & -2 & 1 \\
& & &  & 2 & -2 \\
\end{pmatrix}.$
"

# ╔═╡ c61962e8-86bd-45ed-a62d-ffc812d8e35e
# First row
M[1,2] = 2;

# ╔═╡ 33baab03-fbc1-4944-bc16-e4022b051c69
# Last row
M[end,end-1] = 2;

# ╔═╡ facdbabf-6e2b-418a-a07e-36188244e585
md"### Heat equation in State-space representation

After the spatial approximation the heat equation has a form of

$\dot{\theta}(t) =  \frac{\alpha}{\Delta x^2} ~ M ~ \theta(t) + \frac{2 \alpha}{\Delta x} ~ \Phi_{in}(t) + \frac{2 \alpha}{\Delta x} ~ \Phi_{out}(t) \tag{1}$

with M as diffusion matrix, temperature vector $\quad \theta(t) = \left( \vartheta(t, x^{0}), \cdots, \vartheta(t, x^{N-1}) \right)^{\top}$.

The original ODE $(1)$ is recasted in [state-space representation](https://en.wikipedia.org/wiki/State-space_representation) as

$\dot{\theta}(t) = A ~ \theta(t) + B ~ u(t) + E ~ w(t)$
$y(t) = C \theta(t)$

with
- system matrix $A =  \frac{\alpha}{\Delta x^2} ~ M$,
- input vector
$B = \frac{2 ~ \alpha}{\lambda ~ \Delta x}
\begin{pmatrix}
b \\
0 \\
\vdots \\
0
\end{pmatrix},$

- disturbance vector
$E = \frac{2 ~ \alpha}{\lambda ~ \Delta x}
\begin{pmatrix}
0 \\
\vdots \\
0 \\
1
\end{pmatrix},$

- emitted flux $w(t) =  \phi_{out}(t,x^{N-1})$ and
- output vector $C = (0, \cdots, 0, 1)$

to describe the input-state-output behaviour including the integral controller. The emitted heat is considered as a disturbance. The PI controller is introduced as

$u(t) = K_{p} ~ e(t)$

with error
$e(t) = y_{ref}(t) - y(t) = y_{ref}(t) - C ~ \theta(t)$

The thermal dynamics with PI control is described by

$\dot{\theta}(t) = A ~ \theta(t) + B ~ K_{p} ~ e(t) + E ~ w(t).$

This is an ordinary differential equation (ODE) that can be solved with common solvers like forward Euler method or Runge-Kutta scheme.

"

# ╔═╡ 5ef85c8f-cb69-42df-a098-ab98903a5f32
A = α/(Δx^2) * M

# ╔═╡ 25bb45ab-55af-4b8b-a11b-05be68631b78
begin
	B = spzeros(N); # Input vector
	B[1] = 2 * (α/(λ*Δx)) * b
	E = spzeros(N); # Disturbance vector
	E[end] = 2 * (α/(λ*Δx)) * b
	C = spzeros(1,N); # Output vector
	C[end] = 1;
end

# ╔═╡ 81a2d204-7386-4ed7-af89-7c076940817b
w(θ) = ϕout(θ)

# ╔═╡ 113d8c87-f4e6-4f1b-b188-879e7ad1e146
function heat_eq!(dx,x,p,t)
    global u_hist

    err = yref - x[end]
    u = u_in(err)

    dx .= A * x + B * u + E * ϕout( x[end] )    # Integration of temperatures

  end

# ╔═╡ d9138271-79fc-4b7c-8501-e8d8f3304cd5
md"### Initial Conditions

For simplicity, the initial data of the heat equation is assumed as $\vartheta_{0}(x)= 1000$ Kelvin for all $x \in \left[0,L\right]$. This function is approximated to gain the initial conditions of the ODE

$\theta(0) = \left( 300, \cdots , 300 \right)^{\top}.$
"

# ╔═╡ db6248c7-24b8-4f44-82d3-4892203297b5
θ₀ = 300.0 * ones(N)

# ╔═╡ 837632eb-b4db-47ec-8e4a-32a0fe81060b
md"## Simulation"

# ╔═╡ 3f382fe5-2aa6-4d82-b939-977985ae8bbb
Δt = 10^(-2) # Sampling time

# ╔═╡ 02ba374a-15b4-4481-baf6-c1a3ce45ed5e
tspan = (0.0, Tf)

# ╔═╡ 0db82493-4dda-4992-aa5b-63734b556a90
prob = ODEProblem( heat_eq!, θ₀, tspan ) # ODE Problem

# ╔═╡ f7dcb8a6-0eb4-4762-942e-53362d4fe1ba
sol = solve(prob,Euler(),dt=Δt,progress=true, saveat=1.0, save_start=true); # Solving the ODE

# ╔═╡ 9c38b2b3-1622-41f3-a2c2-ff79b552e454
# 1-dimensional grid
xspan = 0 : Δx : L

# ╔═╡ ba66c778-6c1c-4e32-934d-e8aca0eb1c66
heatmap(sol.t, xspan, sol[1:end,1:end], xaxis="Time [s]", yaxis="Position [m]", title="Evolution of temperature")

# ╔═╡ 3dcfed19-8235-4309-b1c7-8d4e30ba6701
plot(sol.t,[sol[1,:], sol[end-1,:]], label=["Left" "Right"], title="Temperature at the left/right end", legend=:bottomright)

# ╔═╡ Cell order:
# ╟─7bc3e874-b7be-11eb-0792-cf4e21e8c8b2
# ╠═1158de2e-b7a6-43f3-a4e6-2ef9598a4ca5
# ╠═30cebec3-9c70-4c6a-94a8-fdcbe6cda783
# ╟─93fb0f90-c454-410f-aa5e-cd8e6c055dc0
# ╠═b294bbfc-45e7-4a15-a839-5ceb56746be7
# ╠═543b3f13-4cef-46f1-9141-172b81e7d48b
# ╟─1deb117d-21f7-4329-a7db-2be26e43129a
# ╠═78292b42-a203-496c-a143-def3ff128f69
# ╠═0ecea33b-bf24-48eb-b4b4-bbf288e62e6b
# ╠═3a059f15-d850-468a-bb1d-f42ef5b24609
# ╟─7202aa44-b393-4e0b-9c34-2a9b9376194c
# ╠═2775433b-d8b1-48be-9b60-08641db15803
# ╠═4c59909b-ffb6-4dfc-94ab-142d41deeeab
# ╠═4f0c3671-d3bf-4e36-a1f3-8e8d0ac0e454
# ╠═23343fa0-6033-498d-b089-d43e4382e332
# ╟─daa17037-fd28-460c-97f5-826742b9a419
# ╠═892fb737-6fc0-4438-b586-039a440b01f9
# ╠═0309bf5e-4033-4832-86ed-146aea0863c0
# ╠═5c9ec071-6140-4af6-a418-df79c10daf5e
# ╠═905337c5-9eee-4b67-b215-bfd7ca905e4c
# ╠═f0e9f197-303b-48a6-a540-9b11513050d4
# ╟─0e155080-2cfe-45c0-8166-ba91f95d6f8b
# ╠═c61962e8-86bd-45ed-a62d-ffc812d8e35e
# ╠═33baab03-fbc1-4944-bc16-e4022b051c69
# ╟─facdbabf-6e2b-418a-a07e-36188244e585
# ╠═5ef85c8f-cb69-42df-a098-ab98903a5f32
# ╠═25bb45ab-55af-4b8b-a11b-05be68631b78
# ╠═81a2d204-7386-4ed7-af89-7c076940817b
# ╠═113d8c87-f4e6-4f1b-b188-879e7ad1e146
# ╟─d9138271-79fc-4b7c-8501-e8d8f3304cd5
# ╠═db6248c7-24b8-4f44-82d3-4892203297b5
# ╠═837632eb-b4db-47ec-8e4a-32a0fe81060b
# ╠═069cf190-77cc-4757-9db3-1945ddd07438
# ╠═3f382fe5-2aa6-4d82-b939-977985ae8bbb
# ╠═02ba374a-15b4-4481-baf6-c1a3ce45ed5e
# ╠═0db82493-4dda-4992-aa5b-63734b556a90
# ╠═f7dcb8a6-0eb4-4762-942e-53362d4fe1ba
# ╠═9c38b2b3-1622-41f3-a2c2-ff79b552e454
# ╠═ba66c778-6c1c-4e32-934d-e8aca0eb1c66
# ╠═3dcfed19-8235-4309-b1c7-8d4e30ba6701


text:
# byproducts.jl

"""
    StartWorkers(nwrkrs::Int)

Start workers if needed.
"""
function StartWorkers(nwrkrs::Int)
   set_workers = nwrkrs
   nworkers() < set_workers ? addprocs(set_workers) : nothing
   nworkers()
end

"""
    TaskDriver(indx,fn)

Broacast / distribute task (fn; e.g. loop_task1) over indices (indx; e.g. file indices)

Examples:

```
using CbiomesProcessing, Distributed, SparseArrays
TaskDriver(1,CbiomesProcessing.loop_task1)

StartWorkers(4)
@everywhere using CbiomesProcessing, SparseArrays
TaskDriver(1:4,CbiomesProcessing.loop_task1)
```

Visualize results:

```
using FortranFiles, Plots
k=1
recl=720*360*4
fil="diags_interp/ETAN/ETAN.0000000732.data"
f =  FortranFile(fil,"r",access="direct",recl=recl,convert="big-endian")
tmp=read(f,rec=k,(Float32,(720,360))); close(f)
heatmap(tmp)
```
"""
function TaskDriver(indx::Union{UnitRange{Int},Array{Int,1},Int},fn::Function)
    i=collect(indx)
    length(i)>1 ? i=distribute(i) : nothing
    isa(i,DArray) ? println(i.indices) : nothing
    fn.(i)
end

"""
    MetaFileRead(filIn::String)

Reads a meta file generated by MITgcm
"""
function MetaFileRead(FileName::String)

    MetaFile=FileName[1:end-5]*".meta"
    f = open(MetaFile)
    lines = readlines(f)
    close(f)

    MetaFile=Dict("MetaFile" => MetaFile)
    while !isempty(lines)
        line=popfirst!(lines)
        i0=findfirst(isequal('='), line)
        i1=findfirst(isequal(';'), line)
        !isnothing(i0) ? nam=strip(line[1:i0-1]) : nam=""
        val=nothing
        #show(line)
        if nam=="dimList"
            #special case: dimList
            val=fill(0.,(MetaFile["nDims"],3))
            for ii=1:MetaFile["nDims"]
                line=popfirst!(lines)
                tmp1=split(line,",")
                #tmp1=map(x->(v = tryparse(Int,x); ismissing(v) ? 0.0 : v),tmp1)
                val[ii,:]=parse.(Int,tmp1[1:3])
            end
            line=popfirst!(lines)
        elseif nam=="fldList"
            #special case: fldList
            line=popfirst!(lines)
            tmp1=split(line,"'")
            val=String.(tmp1[2:2:end])
            line=popfirst!(lines)
        elseif nam=="dataprec"
            #sepcial case: dataprec
            tmp1=split(line)
            tmp1[4]=="'float32'" ? val=Float32 : val=Float64
        elseif nam=="nDims"
            #sepcial case: nDims
            tmp1=split(line[i0+1:i1-1])
            val=parse(Int64,tmp1[2])
        end
        #
        if ~isnothing(val)
            tmp2=Dict(nam => val)
            MetaFile=merge(MetaFile,tmp2)
        end
    end

    return MetaFile
end

"""
    MatrixInterp(in::Array{T,N},MTRX,siz) where {T,N}

Interpolate `in` using `MTRX` to grid of size `siz`.
"""
function MatrixInterp(in::Array{T,N},MTRX::SparseMatrixCSC,siz) where {T,N}
    #input
    l=size(in,1)*size(in,2);
    m=size(in,3);
    tmp1=reshape(in,l,m)
    tmp0=Float64.(.!(isnan.(tmp1)))
    tmp1[isnan.(tmp1)].=0.
    siz=siz[1],siz[2],m
    #matrix product
    tmp0=MTRX*tmp0
    tmp1=MTRX*tmp1
    tmp1=tmp1./tmp0
    #this may be redundant:
    tmp1[tmp0 .== 0.] .= NaN
    #output
    out=reshape(tmp1,siz)
    m==1 ? out=dropdims(out,dims=3) : nothing
    return out
end

ext: jl

max_issues_repo_name: darwinproject/CBIOMES-Processing-legacy

max_forks_repo_name: darwinproject/CBIOMES-Processing-legacy


text:
using Test
using QuantumOptics, CollectiveSpins
using Statistics

@testset "symmetric-equations" begin

cs = CollectiveSpins

# System parameters
a = 0.18
γ = 1
e_dipole = [0, 0, 1]
T = [0:0.05:5;]

# =========== Test square configuration =========================

system_square = SpinCollection(cs.geometry.square(a), e_dipole; gammas=γ)
N_square = length(system_square.spins)

# Initial state (Bloch state)
phi0_square = [0. for i=1:N_square]
phi1_square = [0., 0.5*pi, 1.5*pi, 1.0*pi]
phi2_square = [0., 1.0*pi, 1.0*pi, 2.0*pi]
theta_square = ones(Float64, N_square)*pi/2.

# Meanfield
state0_mf = cs.meanfield.blochstate(phi0_square, theta_square)
state1_mf = cs.meanfield.blochstate(phi1_square, theta_square)
state2_mf = cs.meanfield.blochstate(phi2_square, theta_square)

tout, state0_mf_t = cs.meanfield.timeevolution(T, system_square, state0_mf)
tout, state1_mf_t = cs.meanfield.timeevolution(T, system_square, state1_mf)
tout, state2_mf_t = cs.meanfield.timeevolution(T, system_square, state2_mf)

# Symmetric meanfield
state0_sym = cs.meanfield.blochstate(0., pi/2., 1)

Ωeff0_, Γeff0_ = cs.effective_interaction.square_orthogonal(a)
Ωeff0, Γeff0 = cs.effective_interaction_rotated.square_orthogonal(a, 0)
Ωeff1, Γeff1 = cs.effective_interaction_rotated.square_orthogonal(a, 1)
Ωeff2, Γeff2 = cs.effective_interaction_rotated.square_orthogonal(a, 2)
@test (Ωeff0_-Ωeff0) < 1e-12
@test (Γeff0_-Γeff0) < 1e-12

tout, state0_mfsym_t = cs.meanfield.timeevolution_symmetric(T, state0_sym, Ωeff0, Γeff0)
tout, state1_mfsym_t = cs.meanfield.timeevolution_symmetric(T, state0_sym, Ωeff1, Γeff1)
tout, state2_mfsym_t = cs.meanfield.timeevolution_symmetric(T, state0_sym, Ωeff2, Γeff2)

for i=1:length(T)
    state0_mf_rotated = state0_mf_t[i]
    state1_mf_rotated = cs.meanfield.rotate([0.,0.,1.], -phi1_square, state1_mf_t[i])
    state2_mf_rotated = cs.meanfield.rotate([0.,0.,1.], -phi2_square, state2_mf_t[i])
    @test var(cs.meanfield.sx(state0_mf_rotated)) < 1e-12
    @test var(cs.meanfield.sy(state0_mf_rotated)) < 1e-12
    @test var(cs.meanfield.sz(state0_mf_rotated)) < 1e-12
    @test var(cs.meanfield.sx(state1_mf_rotated)) < 1e-12
    @test var(cs.meanfield.sy(state1_mf_rotated)) < 1e-12
    @test var(cs.meanfield.sz(state1_mf_rotated)) < 1e-12
    @test var(cs.meanfield.sx(state2_mf_rotated)) < 1e-12
    @test var(cs.meanfield.sy(state2_mf_rotated)) < 1e-12
    @test var(cs.meanfield.sz(state2_mf_rotated)) < 1e-12

    @test cs.meanfield.sx(state0_mf_rotated)[1]-cs.meanfield.sx(state0_mfsym_t[i])[1] < 1e-12
    @test cs.meanfield.sy(state0_mf_rotated)[1]-cs.meanfield.sy(state0_mfsym_t[i])[1] < 1e-12
    @test cs.meanfield.sz(state0_mf_rotated)[1]-cs.meanfield.sz(state0_mfsym_t[i])[1] < 1e-12
    @test cs.meanfield.sx(state1_mf_rotated)[1]-cs.meanfield.sx(state1_mfsym_t[i])[1] < 1e-12
    @test cs.meanfield.sy(state1_mf_rotated)[1]-cs.meanfield.sy(state1_mfsym_t[i])[1] < 1e-12
    @test cs.meanfield.sz(state1_mf_rotated)[1]-cs.meanfield.sz(state1_mfsym_t[i])[1] < 1e-12
    @test cs.meanfield.sx(state2_mf_rotated)[1]-cs.meanfield.sx(state2_mfsym_t[i])[1] < 1e-12
    @test cs.meanfield.sy(state2_mf_rotated)[1]-cs.meanfield.sy(state2_mfsym_t[i])[1] < 1e-12
    @test cs.meanfield.sz(state2_mf_rotated)[1]-cs.meanfield.sz(state2_mfsym_t[i])[1] < 1e-12
end

# =========== Test cube configuration ===========================

system_cube = SpinCollection(cs.geometry.cube(a), e_dipole; gammas=γ)
N_cube = length(system_cube.spins)

# Initial state (Bloch state)

dphi1_cube = pi
phi0_cube = [0. for i=1:N_cube]
phi1_cube = [0., 0., 0., 0., dphi1_cube, dphi1_cube, dphi1_cube, dphi1_cube]
theta_cube = ones(Float64, N_cube)*pi/2.

# Meanfield
state0_mf = cs.meanfield.blochstate(phi0_cube, theta_cube)
state1_mf = cs.meanfield.blochstate(phi1_cube, theta_cube)

tout, state0_mf_t = cs.meanfield.timeevolution(T, system_cube, state0_mf)
tout, state1_mf_t = cs.meanfield.timeevolution(T, system_cube, state1_mf)

# Symmetric meanfield
state0_sym = cs.meanfield.blochstate(0., pi/2., 1)

Ωeff0_, Γeff0_ = cs.effective_interaction.cube_orthogonal(a)
Ωeff0, Γeff0 = cs.effective_interaction_rotated.cube_orthogonal(a, 0.)
Ωeff1, Γeff1 = cs.effective_interaction_rotated.cube_orthogonal(a, dphi1_cube)
@test (Ωeff0_-Ωeff0) < 1e-12
@test (Γeff0_-Γeff0) < 1e-12

tout, state0_mfsym_t = cs.meanfield.timeevolution_symmetric(T, state0_sym, Ωeff0, Γeff0)
tout, state1_mfsym_t = cs.meanfield.timeevolution_symmetric(T, state0_sym, Ωeff1, Γeff1)

for i=1:length(T)
    state0_mf_rotated = state0_mf_t[i]
    state1_mf_rotated = cs.meanfield.rotate([0.,0.,1.], -phi1_cube, state1_mf_t[i])

    @test var(cs.meanfield.sx(state0_mf_rotated)) < 1e-8
    @test var(cs.meanfield.sy(state0_mf_rotated)) < 1e-8
    @test var(cs.meanfield.sz(state0_mf_rotated)) < 1e-8
    @test var(cs.meanfield.sx(state1_mf_rotated)) < 1e-8
    @test var(cs.meanfield.sy(state1_mf_rotated)) < 1e-8
    @test var(cs.meanfield.sz(state1_mf_rotated)) < 1e-8

    @test cs.meanfield.sx(state0_mf_rotated)[1]-cs.meanfield.sx(state0_mfsym_t[i])[1] < 1e-8
    @test cs.meanfield.sy(state0_mf_rotated)[1]-cs.meanfield.sy(state0_mfsym_t[i])[1] < 1e-8
    @test cs.meanfield.sz(state0_mf_rotated)[1]-cs.meanfield.sz(state0_mfsym_t[i])[1] < 1e-8
    @test cs.meanfield.sx(state1_mf_rotated)[1]-cs.meanfield.sx(state1_mfsym_t[i])[1] < 1e-8
    @test cs.meanfield.sy(state1_mf_rotated)[1]-cs.meanfield.sy(state1_mfsym_t[i])[1] < 1e-8
    @test cs.meanfield.sz(state1_mf_rotated)[1]-cs.meanfield.sz(state1_mfsym_t[i])[1] < 1e-8
end

end # testset


text:
include("euler/euler.jl")
using .Grid: max_path_sum
using BenchmarkTools
BenchmarkTools.DEFAULT_PARAMETERS.samples = 100

compute(triangle::String)::Int = max_path_sum([parse.(Int, split(line, " ")) for line ∈ split(triangle, "\n")])

small_triangle = """3
7 4
2 4 6
8 5 9 3"""

big_triangle = """75
95 64
17 47 82
18 35 87 10
20 04 82 47 65
19 01 23 75 03 34
88 02 77 73 07 63 67
99 65 04 28 06 16 70 92
41 41 26 56 83 40 80 70 33
41 48 72 33 47 32 37 16 94 29
53 71 44 65 25 43 91 52 97 51 14
70 11 33 28 77 73 17 78 39 68 17 57
91 71 52 38 17 14 91 43 58 50 27 29 48
63 66 04 68 89 53 67 30 73 16 69 87 40 31
04 62 98 27 23 09 70 98 73 93 38 53 60 04 23"""

compute(small_triangle)

compute(big_triangle)

@benchmark compute(big_triangle)

max_stars_repo_name: Dynortice/Project-Euler

max_issues_repo_name: Dynortice/Project-Euler

max_forks_repo_name: Dynortice/Project-Euler


text:
using LinearAlgebra
import LinearAlgebra:dot, norm

#=
Note that all of the inner product operations defined here are based on a
domain filled with dual cells, in which the dual nodes are at their centers, and surrounded
by a layer of padding (i.e. ghost) cells.
=#

"""
    dot(p1::Nodes{Dual},p2::Nodes{Dual}) -> Real

Computes the inner product between two sets of dual node data on the same grid.
"""
function dot(p1::Nodes{Dual,NX,NY},p2::Nodes{Dual,NX,NY}) where {NX,NY}
  # remember that sizes NX and NY include the ghost cells
  dims = node_inds(Dual,(NX,NY))
  return dot(p1[2:dims[1]-1,2:dims[2]-1],p2[2:dims[1]-1,2:dims[2]-1])/((NX-2)*(NY-2))
end

"""
    dot(p1::Nodes{Primal},p2::Nodes{Primal}) -> Real

Computes the inner product between two sets of primal node data on the same grid.
"""
function dot(p1::Nodes{Primal,NX,NY},p2::Nodes{Primal,NX,NY}) where {NX,NY}

  dims = node_inds(Primal,(NX,NY))

  # interior
  tmp = dot(p1[2:dims[1]-1,2:dims[2]-1],p2[2:dims[1]-1,2:dims[2]-1])

  # boundaries
  tmp += 0.5*dot(p1[1,2:dims[2]-1],      p2[1,2:dims[2]-1])
  tmp += 0.5*dot(p1[dims[1],2:dims[2]-1],p2[dims[1],2:dims[2]-1])
  tmp += 0.5*dot(p1[2:dims[1]-1,1],      p2[2:dims[1]-1,1])
  tmp += 0.5*dot(p1[2:dims[1]-1,dims[2]],p2[2:dims[1]-1,dims[2]])

  # corners -- use dot to ensure it works for complex types
  tmp += 0.25*dot([p1[1,1],p1[dims[1],1],p1[1,dims[2]],p1[dims[1],dims[2]]],
      [p2[1,1],p2[dims[1],1],p2[1,dims[2]],p2[dims[1],dims[2]]])
  #tmp += 0.25*(p1[1,1]*p2[1,1]             + p1[dims[1],1]*p2[dims[1],1] +
  #             p1[1,dims[2]]*p2[1,dims[2]] + p1[dims[1],dims[2]]*p2[dims[1],dims[2]])

  return tmp/((NX-2)*(NY-2))
end

"""
    dot(p1::XEdges{Dual},p2::XEdges{Dual}) -> Real

Computes the inner product between two sets of dual x-edge component data on the same grid.
"""
function dot(p1::XEdges{Dual,NX,NY},p2::XEdges{Dual,NX,NY}) where {NX,NY}

  udims = xedge_inds(Dual,(NX,NY))

  # interior
  tmp = dot(p1[2:udims[1]-1,2:udims[2]-1],p2[2:udims[1]-1,2:udims[2]-1])

  # boundaries
  tmp += 0.5*dot(p1[1,       2:udims[2]-1],p2[1,       2:udims[2]-1])
  tmp += 0.5*dot(p1[udims[1],2:udims[2]-1],p2[udims[1],2:udims[2]-1])

  return tmp/((NX-2)*(NY-2))
end

"""
    dot(p1::YEdges{Dual},p2::YEdges{Dual}) -> Real

Computes the inner product between two sets of dual y-edge component data on the same grid.
"""
function dot(p1::YEdges{Dual,NX,NY},p2::YEdges{Dual,NX,NY}) where {NX,NY}

  vdims = yedge_inds(Dual,(NX,NY))

  # interior
  tmp = dot(p1[2:vdims[1]-1,2:vdims[2]-1],p2[2:vdims[1]-1,2:vdims[2]-1])

  # boundaries
  tmp += 0.5*dot(p1[2:vdims[1]-1,1],       p2[2:vdims[1]-1,1])
  tmp += 0.5*dot(p1[2:vdims[1]-1,vdims[2]],p2[2:vdims[1]-1,vdims[2]])

  return tmp/((NX-2)*(NY-2))
end

"""
    dot(p1::XEdges{Primal},p2::XEdges{Primal}) -> Real

Computes the inner product between two sets of primal x-edge component data on the same grid.
"""
function dot(p1::XEdges{Primal,NX,NY},p2::XEdges{Primal,NX,NY}) where {NX,NY}

  udims = xedge_inds(Primal,(NX,NY))

  # interior
  tmp = dot(p1[2:udims[1]-1,2:udims[2]-1],p2[2:udims[1]-1,2:udims[2]-1])

  # boundaries
  tmp += 0.5*dot(p1[2:udims[1]-1,1],       p2[2:udims[1]-1,1])
  tmp += 0.5*dot(p1[2:udims[1]-1,udims[2]],p2[2:udims[1]-1,udims[2]])

  return tmp/((NX-2)*(NY-2))
end

"""
    dot(p1::YEdges{Primal},p2::YEdges{Primal}) -> Real

Computes the inner product between two sets of primal y-edge component data on the same grid.
"""
function dot(p1::YEdges{Primal,NX,NY},p2::YEdges{Primal,NX,NY}) where {NX,NY}

  vdims = yedge_inds(Primal,(NX,NY))

  # interior
  tmp = dot(p1[2:vdims[1]-1,2:vdims[2]-1],p2[2:vdims[1]-1,2:vdims[2]-1])

  # boundaries
  tmp += 0.5*dot(p1[1,       2:vdims[2]-1],p2[1,       2:vdims[2]-1])
  tmp += 0.5*dot(p1[vdims[1],2:vdims[2]-1],p2[vdims[1],2:vdims[2]-1])

  return tmp/((NX-2)*(NY-2))
end

"""
    dot(p1::Edges{Dual/Primal},p2::Edges{Dual/Primal}) -> Real

Computes the inner product between two sets of dual edge data on the same grid.
"""
dot(p1::Edges{C,NX,NY},p2::Edges{C,NX,NY}) where {C<:CellType,NX,NY} =
      dot(p1.u,p2.u) + dot(p1.v,p2.v)

######## NORMS ###########

"""
    norm(p::GridData) -> Real

Computes the L2 norm of data on a grid.
"""
norm(p::GridData) = sqrt(dot(p,p))

######## INTEGRALS ###########

# This function computes an integral by just taking the inner product with
# another set of cell data uniformly equal to 1
"""
    integrate(p::GridData) -> Real

Computes a numerical quadrature of node data.
"""
function integrate(p::T) where {T<:GridData}
  p2 = T()
  fill!(p2,1) # fill it with ones
  return dot(p2,p) # order matters for complex types
end


text:
"""
$(TYPEDEF)

Direction observation information of a `Point2` variable.
"""
mutable struct PriorPoint2{T} <: IncrementalInference.FunctorSingleton where {T <: IIF.SamplableBelief}
  Z::T
  # W::Array{Float64,1} # TODO, deprecate the weight parameter
  PriorPoint2{T}() where T = new()
  PriorPoint2{T}(dist::T) where {T <: IIF.SamplableBelief} = new{T}(dist)
end
PriorPoint2(z::T) where {T <: IIF.SamplableBelief} = PriorPoint2{T}(z)
function PriorPoint2D(mu, cov, W)
  @warn "PriorPoint2D(mu, cov, W) is deprecated, use PriorPoint{T}(T(...)) instead -- e.g. PriorPoint2{MvNormal}(MvNormal(...) or any other Distributions.Distribution type instead."
  PriorPoint2{MvNormal{Float64}}(MvNormal(mu, cov))
end
function getSample(p2::PriorPoint2, N::Int=1)
  return (rand(p2.Z, N),)
end

"""
$(TYPEDEF)
"""
mutable struct Point2Point2{D <: IIF.SamplableBelief} <: FunctorPairwise
    Zij::D
    Point2Point2{T}() where T = new{T}()
    Point2Point2{T}(x::T) where {T <: IIF.SamplableBelief} = new{T}(x)
end
Point2Point2(x::T) where {T <: IIF.SamplableBelief} = Point2Point2{T}(x)
function getSample(pp2::Point2Point2, N::Int=1)
  return (rand(pp2.Zij,N),  )
end
function (pp2r::Point2Point2{T})(
            res::Array{Float64},
            userdata::FactorMetadata,
            idx::Int,
            meas::Tuple,
            xi::Array{Float64,2},
            xj::Array{Float64,2} ) where T
  #
  res[1]  = meas[1][1,idx] - (xj[1,idx] - xi[1,idx])
  res[2]  = meas[1][2,idx] - (xj[2,idx] - xi[2,idx])
  nothing
end

"""
$(TYPEDEF)
"""
mutable struct Point2Point2WorldBearing{T} <: IncrementalInference.FunctorPairwise where {T <: IIF.SamplableBelief}
    Z::T
    rangemodel::Rayleigh
    # zDim::Tuple{Int, Int}
    Point2Point2WorldBearing{T}() where T = new{T}()
    Point2Point2WorldBearing{T}(x::T) where {T <: IIF.SamplableBelief} = new{T}(x, Rayleigh(100))
end
Point2Point2WorldBearing(x::T) where {T <: IIF.SamplableBelief} = Point2Point2WorldBearing{T}(x)
function getSample(pp2::Point2Point2WorldBearing, N::Int=1)
  sp = Array{Float64,2}(undef, 2,N)
  sp[1,:] = rand(pp2.Z,N)
  sp[2,:] = rand(pp2.rangemodel,N)
  return (sp, )
end
function (pp2r::Point2Point2WorldBearing)(
          res::Array{Float64},
          userdata::FactorMetadata,
          idx::Int,
          meas::Tuple,
          pi::Array{Float64,2},
          pj::Array{Float64,2} )
  #
  # noisy bearing measurement
  z, r = meas[1][1,idx], meas[1][2,idx]
  dx, dy = pj[1,idx]-pi[1,idx], pj[2,idx]-pi[2,idx]
  res[1] = z - atan(dy,dx)
  res[2] = r - norm([dx; dy])
  nothing
end

"""
$(TYPEDEF)

Will be deprecated, use `addFactor!(.., nullhypo=)` instead (work in progress)
"""
mutable struct PriorPoint2DensityNH <: IncrementalInference.FunctorSingletonNH
  belief::BallTreeDensity
  nullhypothesis::Distributions.Categorical
  PriorPoint2DensityNH() = new()
  PriorPoint2DensityNH(belief, p::Distributions.Categorical) = new(belief, p)
  PriorPoint2DensityNH(belief, p::Vector{Float64}) = new(belief, Distributions.Categorical(p))
end
function getSample(p2::PriorPoint2DensityNH, N::Int=1)
  return (rand(p2.belief, N), )
end

"""
$(TYPEDEF)

Will be deprecated, use `addFactor!(.., nullhypo=)` instead (work in progress)
"""
mutable struct PackedPriorPoint2DensityNH <: IncrementalInference.PackedInferenceType
    rpts::Vector{Float64} # 0rotations, 1translation in each column
    rbw::Vector{Float64}
    dims::Int
    nh::Vector{Float64}
    PackedPriorPoint2DensityNH() = new()
    PackedPriorPoint2DensityNH(x1,x2,x3, x4) = new(x1, x2, x3, x4)
end
function convert(::Type{PriorPoint2DensityNH}, d::PackedPriorPoint2DensityNH)
  return PriorPoint2DensityNH(
            kde!(EasyMessage( reshapeVec2Mat(d.rpts, d.dims), d.rbw, (:Euclid, :Euclid))),
            Distributions.Categorical(d.nh)  )
end
function convert(::Type{PackedPriorPoint2DensityNH}, d::PriorPoint2DensityNH)
  return PackedPriorPoint2DensityNH( getPoints(d.belief)[:], getBW(d.belief)[:,1], Ndim(d.belief), d.nullhypothesis.p )
end

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

"""
$(TYPEDEF)

Serialization type for `PriorPoint2`.
"""
mutable struct PackedPriorPoint2  <: IncrementalInference.PackedInferenceType
    str::String
    PackedPriorPoint2() = new()
    PackedPriorPoint2(x::String) = new(x)
end

function convert(::Type{PriorPoint2}, d::PackedPriorPoint2)
  # Cov = reshapeVec2Mat(d.vecCov, d.dimc)
  distr = extractdistribution(d.str)
  return PriorPoint2{typeof(distr)}(distr)
end
function convert(::Type{PackedPriorPoint2}, d::PriorPoint2)
  # v2 = d.mv.Σ.mat[:];
  return PackedPriorPoint2(string(d.Z))
end

"""
$(TYPEDEF)

Serialization type for `Point2Point2WorldBearing`.
"""
mutable struct PackedPoint2Point2WorldBearing  <: IncrementalInference.PackedInferenceType
    str::String
    # NOTE Not storing rangemodel which may cause inconsistencies if the implementation parameters change
    PackedPoint2Point2WorldBearing() = new()
    PackedPoint2Point2WorldBearing(x::String) = new(x)
end
function convert(::Type{PackedPoint2Point2WorldBearing}, d::Point2Point2WorldBearing)
  return PackedPoint2Point2WorldBearing( string(d.Z) )
end
function convert(::Type{Point2Point2WorldBearing}, d::PackedPoint2Point2WorldBearing)
  return Point2Point2WorldBearing( extractdistribution(d.str) )
end

"""
$(TYPEDEF)

Serialization type for `Point2Point2`.
"""
mutable struct PackedPoint2Point2 <: IncrementalInference.PackedInferenceType
    # mu::Vector{Float64}
    # sigma::Vector{Float64}
    # sdim::Int
    str::String
    PackedPoint2Point2() = new()
    PackedPoint2Point2(s::AS) where {AS <: AbstractString} = new(s)
end
function convert(::Type{Point2Point2}, d::PackedPoint2Point2)
  return Point2Point2( extractdistribution(d.str) )
end
function convert(::Type{PackedPoint2Point2}, d::Point2Point2)
  return PackedPoint2Point2( string(d.Zij) )
end


text:
using JuMP, SatisfiabilityInterface

@testset "JuMP" begin
    model = JuMP.Model(SatisfiabilityInterface.Optimizer{Bool})
    @variable(model, x, Bin)
    @variable(model, y, Bin)
    @variable(model, z, Bin)
    @constraint(model, [x, y, z] in SatisfiabilityInterface.CNF([true, true, true]))
    @constraint(model, [x, y, z] in SatisfiabilityInterface.CNF([false, true, true]))
    @constraint(model, [x, y]    in SatisfiabilityInterface.CNF([false, false]))
    @constraint(model, [y, z]    in SatisfiabilityInterface.CNF([false, false]))
    @constraint(model, [x, z]    in SatisfiabilityInterface.CNF([true, false]))
    optimize!(model)
    @test raw_status(model) == "sat"
    @test termination_status(model) == MOI.OPTIMAL
    @test primal_status(model) == MOI.FEASIBLE_POINT
    @test value(x) == 0.0
    @test value(y) == 1.0
    @test value(z) == 0.0
end


text:
include("../NativeInterface/native.jl")

function call_basic_triangulation(flat_vertices::Vector{Cdouble}, flat_vertices_map::Vector{Cint})
    return NativeInterface.basic_triangulation(flat_vertices, flat_vertices_map)
end

function call_constrained_triangulation(flat_vertices::Vector{Cdouble}, flat_vertices_map::Vector{Cint}, flat_vertices_edge::Vector{Cint})
    return NativeInterface.constrained_triangulation(flat_vertices, flat_vertices_map, flat_vertices_edge)
end

function call_constrained_triangulation_bounded(flat_vertices::Vector{Cdouble}, flat_vertices_map::Vector{Cint}, flat_vertices_edge::Vector{Cint}, flat_boundary_edges::Vector{Cint})
    options = NativeInterface.TriangulateOptions()
    options.pslg = true

    return NativeInterface.constrained_triangulation_bounded(flat_vertices, flat_vertices_map, flat_vertices_edge, flat_boundary_edges, options)
end

function call_constrained_triangulation_bounded(flat_vertices::Vector{Cdouble}, flat_vertices_map::Vector{Cint}, flat_vertices_edge::Vector{Cint}, flat_boundary_edges::Vector{Cint}, flat_holes::Vector{Cdouble})
    options = NativeInterface.TriangulateOptions()
    options.pslg = true

    return NativeInterface.constrained_triangulation_bounded(flat_vertices, flat_vertices_map, flat_vertices_edge, flat_boundary_edges, flat_holes, options)
end


text:
# This file was a part of Julia.
# License is MIT: http://julialang.org/license
if VERSION >= v"0.5-"
    using BaseTestDeprecated
    const Test = BaseTestDeprecated
else
    using Base.Test
end

# test file to test testing

# Test @test
@test true
@test 1 == 1
@test 1 != 2
@test strip("\t  hi   \n") == "hi"
@test strip("\t  this should fail   \n") != "hi"

scary = Test.Error("hi",DimensionMismatch,[])
@test sprint(showerror,scary) == "test error in expression: hi\nDimensionMismatch"

a = Array(Float64, 2, 2, 2, 2, 2)
a[1,1,1,1,1] = 10
@test a[1,1,1,1,1] == 10
@test a[1,1,1,1,1] != 2

@test rand() != rand()

# Test with_handler
successflag = false
failureflag = false
errorflag = false
test_handler(r::Test.Success) = !successflag
test_handler(r::Test.Failure) = !failureflag
test_handler(r::Test.Error) = !errorflag

Test.with_handler(test_handler) do
    @test true
    @test successflag
    @test !failureflag
    @test !errorflag
    successflag = false
    @test false
    @test !successflag
    @test failureflag
    @test !errorflag
    failureflag = false
    @test error("throw error")
    @test !successflag
    @test !failureflag
    @test errorflag
end

# Test evaluation of comparison tests
i7586_1() = 1
i7586_2() = 7
i7586_3() = 9

comparison_flags_s = [false,false,false]
comparison_flags_f = [false,false,false]
function test_handler2(r::Test.Success)
    comparison_flags_s[1] = (r.resultexpr.args[1] == 1)
    comparison_flags_s[2] = (r.resultexpr.args[3] == 7)
    comparison_flags_s[3] = (r.resultexpr.args[5] == 9)
end

function test_handler2(r::Test.Failure)
    comparison_flags_f[1] = (r.resultexpr.args[1] == 1)
    comparison_flags_f[2] = (r.resultexpr.args[3] == 7)
    comparison_flags_f[3] = (r.resultexpr.args[5] == 10)
end

Test.with_handler(test_handler2) do
    @test i7586_1() <= i7586_2() <= i7586_3()
    @test i7586_1() >= i7586_2() >= 10
end
@test all(comparison_flags_s)
@test all(comparison_flags_f)

# Test @test_throws
domainerror_thrower() = throw(DomainError())
boundserror_thrower() = throw(BoundsError())
error_thrower() = error("An error happened")
@test_throws DomainError domainerror_thrower()
@test_throws BoundsError boundserror_thrower()

failureflag = false
successflag = false
Test.with_handler(test_handler) do
    @test_throws DomainError boundserror_thrower()
    @test failureflag
    @test_throws DomainError domainerror_thrower()
    @test successflag
end

# Test @test_approx_eq
# TODO
@test isapprox(.1+.1+.1, .3)
@test !isapprox(.1+.1+.1, .4)

@test_throws ErrorException Test.test_approx_eq(ones(10),ones(11),1e-8,"a","b")
@test_throws ErrorException Test.test_approx_eq(ones(10),zeros(10),1e-8,"a","b")

# Test @test_approx_eq_eps
# TODO


text:
@testset "normal_sphere_x_axis" begin
    s = RayTracer.sphere()
    n = RayTracer.normal_at(s, RayTracer.point(1,0,0))
    @test n ≈ RayTracer.vector(1,0,0)
end

@testset "normal_sphere_y_axis" begin
    s = RayTracer.sphere()
    n = RayTracer.normal_at(s, RayTracer.point(0,1,0))
    @test n ≈ RayTracer.vector(0,1,0)
end

@testset "normal_sphere_z_axis" begin
    s = RayTracer.sphere()
    n = RayTracer.normal_at(s, RayTracer.point(0,0,1))
    @test n ≈ RayTracer.vector(0,0,1)
end

@testset "normal_sphere_random_point" begin
    s = RayTracer.sphere()
    n = RayTracer.normal_at(s, RayTracer.point(√3/3,√3/3,√3/3))
    @test n ≈ RayTracer.vector(√3/3,√3/3,√3/3)
end

@testset "normal_is_normalized" begin
    s = RayTracer.sphere()
    n = RayTracer.normal_at(s, RayTracer.point(√3/3,√3/3,√3/3))
    @test n ≈ RayTracer.normalize(n)
end

@testset "normal_translated_sphere" begin
    s = RayTracer.sphere()
    RayTracer.set_transform!(s, RayTracer.translation(0,1,0))
    n = RayTracer.normal_at(s,RayTracer.point(0,1+√2/2,-√2/2))
    @test n ≈ RayTracer.vector(0,√2/2,-√2/2)
end

@testset "normal_transformed_sphere" begin
    s = RayTracer.sphere()
    m = RayTracer.scaling(1,0.5,1)*RayTracer.rotation_z(π/5)
    RayTracer.set_transform!(s, m)
    n = RayTracer.normal_at(s,RayTracer.point(0,√2/2,-√2/2))
    @test ≈(n,RayTracer.vector(0,0.97014,-0.24254),1e-4)
end

@testset "default_material_sphere" begin
    s = RayTracer.sphere()
    m = s.material
    @test m ≈ RayTracer.material()
end

@testset "assigning_material_to_sphere" begin
    s = RayTracer.sphere()
    m = RayTracer.material()
    m.ambient = 1
    s.material = m
    @test s.material ≈ m
end


text:
using PolyhedralRelaxations
import Memento
using JuMP

const PR = PolyhedralRelaxations

PR.logger_config!("debug")

using Test
using Cbc
using Ipopt

cbc_optimizer = JuMP.optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
ipopt_optimizer =
    JuMP.optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0, "sb" => "yes")

@testset "PolyhedralRelaxations" begin

    include("api_tests.jl")
    include("error_tests.jl")
    include("objective_tests.jl")
    include("tolerance_tests.jl")

end


text:
"""
    localImprove{Tv,Ti}(G::SparseMatrixCSC{Tv,Ti}, A::Array{Int64,1}; epsSigma=-1.0, err=1e-10, maxSize = max(G.n, G.m)

The LocalImprove function, from the Orrechia-Zhu paper. Given a graph and an initial set, finds a set of smaller conductance
based on the starting set using a localized version of max-flow.

Small discussion: When adding in the neighbors of the initial component, if the resulting  conductance is worse than the initial one,
the algorithm will add more and more vertices until hitting a better conductance. However, if we fix a certain  maximum size for our component,
it might be the case that this new conductance will always be worse than what we had initially. Thus, if we run the algorithm with a small maxSize,
our initial conductance might be the best solution we can raech.

* G is the given graph, A is the initial set
* epsSigma is a measure of the quality of the returning set (the smaller the better). It's defaulted to volume(A) / volume(V - A)
* err is the numerical error considered throughout the algorithm. It's defaulted to 1e-10
* maxSize is the maximum allowed size for the flow graph at any iteration of the algorithm. It's defaulted to |V|
"""
function localImprove(G::SparseMatrixCSC{Tv,Ti}, A::Array{Int64,1}; epsSigma=-1.0, err=1e-10, maxSize = max(G.n, G.m)) where {Tv,Ti}
  #=
    Notes: err smaller than 1e-13 might give weird behavior
  =#

  n = max(G.n, G.m)

  # set epsSigma if it was not set
  minEpsSigma = getVolume(G, A) / getVolume(G, setdiff(collect(1:n), A))
  if epsSigma == -1
    epsSigma = minEpsSigma
  end

  # show warnings for epsSigma
  if epsSigma < minEpsSigma
    print_with_color(:red, "eps-sigma should be greater than ")
    println(minEpsSigma)
  end

  # show warnings for maxSize
  nbrs = Set{Int64}()
  for i in 1:length(A)
    push!(nbrs, A[i])
    for j in 1:deg(G,A[i])
      push!(nbrs, nbri(G,A[i],j))
    end
  end
  if length(nbrs) > maxSize
    print_with_color(:red, string("maxSize should be at least ", length(nbrs), " .....finishing execution"))
    return [],0
  end

  # the actual algorithm
  alphaMin = 0.0
  alphaMax = 1.0
  while alphaMax - err > alphaMin
    alpha = (alphaMin + alphaMax) / 2

    if abs(localFlow(G, A, alpha, epsSigma, maxSize)[2] - getVolume(G, A)) < err
      alphaMin = alpha
    else
      alphaMax = alpha
    end
  end

  # check if the initial conductance is better than what we got after running flow
  if compConductance(G, A) < compConductance(G, localFlow(G, A, alphaMax, epsSigma, maxSize)[1])
    return A, compConductance(G, A)
  else
    return localFlow(G, A, alphaMax, epsSigma, maxSize)
  end

end # localImprove

" The LocalFlow function, from the Orecchia-Zhu paper "
function localFlow(G::SparseMatrixCSC{Tv,Ti}, A::Array{Int64,1}, alpha::Float64, epsSigma::Float64, maxSize = max(G.n, G.m)) where {Tv,Ti}

  # compute the number of vertices
  n = max(G.n, G.m)

  # compute the number of Iterations for which we'll do flow
  volA = getVolume(G, A)
  I = 5 / alpha * log(3 * volA / alpha)

  #=
    Initially: {A} + {Neighbors(A)}
    newID will represent the index in GPrime of a certain vertex in G
    oldID will represent the index in G of a vertex in GPrime
  =#

  newID = Dict{Int64,Int64}()
  oldID = Dict{Int64,Int64}()
  s = maxSize + 1
  t = maxSize + 2

  # initialize GPrime with the input set, A. The current index for new vertices will be "current"
  GPrime = initGPrime(G, A, newID, oldID, alpha, maxSize)
  current = length(A)

  # "considered" will keep track of the vertices interesting for the flow graph.
  considered = Set{Int64}()
  push!(considered, s)
  push!(considered, t)
  for i in 1:length(A)
    push!(considered, A[i])
  end
  for i in 1:length(A)
    u = A[i]
    for j in 1:deg(G,u)
      v = nbri(G,u,j)

      if (v in considered) == false
        push!(considered, v)

        # add this new vertex to GPrime
        current = current + 1

        addToGPrime(G, GPrime, newID, oldID, v, current, alpha, epsSigma, maxSize)
      end
    end
  end

  # do the iterative part of the algorithm: lines 4-9
  totalflow = 0
  for iter in 1:floor(Int64, I)
    saturated, flowInc = localBlockFlow(GPrime, s, t)
    saturated = [oldID[u] for u in saturated]
    if flowInc == 0
      break
    end

    # update the total flow
    totalflow = totalflow + flowInc

    # update the vertices considered for the flow
    newbatch = Set{Int64}()
    for u in saturated
      for i in 1:deg(G, u)
        v = nbri(G, u, i)

        if (v in considered) == false
          push!(newbatch, v)
          push!(considered, v)
        end
      end
    end

    # see if the size of the graph exceeds the maximum size imposed by the user
    if (length(considered) > maxSize)
      continue
    end

    # update GPrime
    for u in newbatch
      current = current + 1
      addToGPrime(G, GPrime, newID, oldID, u, current, alpha, epsSigma, maxSize)
    end

  end

  return [oldID[u] for u in getCutSet(GPrime, s, t)], totalflow

end # localFlow

" Initialize GPrime with the set A and edges of type s->u"
function initGPrime(G::SparseMatrixCSC{Tv,Ti}, A::Array{Int64,1}, newID::Dict{Int64,Int64}, oldID::Dict{Int64,Int64}, alpha::Float64, maxSize::Int64) where {Tv,Ti}
  GPrime = [Tuple{Int64,Float64}[] for i in 1:(maxSize + 2)]
  s = maxSize + 1

  for i in 1:length(A)
    newID[A[i]] = i
    oldID[i] = A[i]

    # add a directed edges from s to i
    push!(GPrime[s], (newID[A[i]], wdeg(G, A[i])))
    push!(GPrime[newID[A[i]]], (s, 0.0))
  end

  # add edges between the vertices in A
  for u in A
    for i in 1:deg(G,u)
      v = nbri(G,u,i)
      if haskey(newID, v)
        push!(GPrime[newID[u]], (newID[v], weighti(G,u,i) / alpha))
      end
    end
  end

  return GPrime
end # initGPrime

" Add a new vertex to GPrime "
function addToGPrime(G::SparseMatrixCSC{Tv,Ti}, GPrime::Array{Array{Tuple{Int64,Float64},1},1},
                     newID::Dict{Int64,Int64}, oldID::Dict{Int64,Int64},
                     u::Int64, newU::Int64, alpha::Float64, epsSigma::Float64, maxSize::Int64) where {Tv,Ti}
  t = maxSize + 2

  # 1. store indices for the new vertex
  newID[u] = newU
  oldID[newU] = u

  # 2. add an edge from newU to t (directed edge)
  push!(GPrime[newU], (t, epsSigma * wdeg(G, u)))
  push!(GPrime[t], (newU, 0))

  # 3. add weights between newU and his neighbors already evaluated
  for i in 1:deg(G, u)
    v = nbri(G, u, i)
    if haskey(newID, v)
      newV = newID[v]
      push!(GPrime[newU], (newV, weighti(G, u, i) / alpha))
      push!(GPrime[newV], (newU, weighti(G, u, i) / alpha))
    end
  end
end

" Compute block flow between s and t"
function localBlockFlow(G::Array{Array{Tuple{Int64,Float64},1},1}, s::Int64, t::Int64)

  backInd = backIndices(G)

  n = length(G)
  inQ = zeros(Bool, n)
  prev = [(0,0) for i in 1:n]
  Q = zeros(Int64, n)

  # do one iteration of block-flow
  left = 0
  right = 1
  Q[right] = s
  inQ[s] = 1

  while left < right
    left = left + 1
    u = Q[left]
    if (u == t)
      continue
    end

    for i in 1:length(G[u])
      v,w = G[u][i]

      if (inQ[v] == false) && (w > 0)
        right = right + 1
        Q[right] = v
        inQ[v] = true
        prev[v] = (u, i)
      end
    end
  end

  if inQ[t] == false
    return [], 0
  else
    totalflow = 0
    saturated = []

    for i in 1:length(G[t])
      # ignore the current vertex if flow hasn't passed through it
      u,w = G[t][i]
      if inQ[u] == false
        continue
      end

      # see the index of t in u's list
      indexT = 1
      for j in 1:length(G[u])
        if G[u][j][1] == t
          indexT = j
          break
        end
      end

      # indexT = backInd[t][i]
      prev[t] = (u, indexT)
      u = t

      # compute the augmentation flow
      currentflow = typemax(Float64)
      while u != s
        prevU = prev[u][1]
        indexPrevU = prev[u][2]
        currentflow = min(currentflow, G[prevU][indexPrevU][2])
        u = prevU
      end

      if currentflow != 0
        # update edges with augmentation flow
        u = t
        while u != s
          prevU = prev[u][1]
          # the index of u in prevU's neighbor list
          indexPrevU = prev[u][2]

          # the index of prevU in u's neighbor list
          indexU = backInd[prevU][indexPrevU]

          # between prevU and u
          G[prevU][indexPrevU] = (G[prevU][indexPrevU][1], G[prevU][indexPrevU][2] - currentflow)

          # between u and prevU
          G[u][indexU] = (G[u][indexU][1], G[u][indexU][2] + currentflow)

          u = prevU
        end
      end

      # check if the vertex to the sink is saturated
      u = G[t][i][1]
      if G[u][indexT][2] == 0
        push!(saturated, u)
      end

      totalflow = totalflow + currentflow
    end

    return saturated, totalflow
  end

end # localBlockFlow

" Get the min cut from the source - return all vertices in the cut besides the source "
function getCutSet(G::Array{Array{Tuple{Int64,Float64},1},1}, s::Int64, t::Int64)

  n = length(G)
  Q = zeros(Int64, n)
  inQ = zeros(Bool, n)

  left = 0
  right = 1
  Q[right] = s
  inQ[s] = true

  while left < right
    left = left + 1
    u = Q[left]

    for i in 1:length(G[u])
      v,w = G[u][i]

      if inQ[v] == 0 && w > 0
        right = right + 1
        Q[right] = v
        inQ[v] = 1
      end
    end
  end

  # consider only vertices in the initial graph into the getCutSet
  cut = (Int64)[]
  for i in 1:right
    if Q[i] != s && Q[i] != t
      push!(cut, Q[i])
    end
  end

  return cut

end # getCutSet

"""
  prn{Tv, Ti}(G::SparseMatrixCSC{Tv,Ti}, s::Array{Int64,1}, phi::Float64, b::Int64)

The PageRank-Nibble cutting algorithm from the Anderson/Chung/Lang paper\n
s is a set of starting vertices, phi is a constant in (0, 1], and b is an integer in [1, [log m]]

phi is a bound on the quality of the conductance of the cut - the smaller the phi, the higher the quality.
b is used to handle precision throughout the algorithm - the higher the b, the greater the precision.
"""
function prn(G::SparseMatrixCSC{Tv,Ti}, s::Array{Int64,1}, phi::Float64, b::Int64) where {Tv,Ti}

  m = div(nnz(G), 2)

  # show warnings
  if (!(0 < phi && phi <= 1))
    print_with_color(:red, "phi should be in (0, 1]\n")
  end
  if (!(1 <= b && b <= round(Int64, log(m))))
    print_with_color(:red, string("b should be in [1, ", string(round(Int64, log(m))), "]\n"))
  end

  alpha = phi * phi / (225 * log(100 * sqrt(m)))
  eps = 1 / (2^b * 48 * log(m))

  p = apr(G, s, alpha, eps)

  # iterate through support set and find a cut-set S that fits our constraints
  S = Set(Int64[])
  volS = 0
  volSc = 2 * m
  obound = 0

  while !isempty(p)
    u = DataStructures.dequeue!(p)
    push!(S, u)

    # update volumes and number of connecting edges
    volS = volS + deg(G, u)
    volSc = volSc - deg(G, u)
    for j in 1:deg(G, u)
      v = nbri(G, u, j)

      if v in S
        obound = obound - 1
      else
        obound = obound + 1
      end
    end

    # next are ignoring conditions

    # ignore from conductance
    conductance = obound / (min(volS, volSc))
    if (conductance >= phi)
      continue
    end

    # ignore from set size
    if !(2^b < volS && volS < 2.0/3.0 * (2 * m))
      continue
    end

    return collect(S)
  end

  return Int64[]

end # prn_local

"""
Computes an approximate page rank vector from a starting set s, an alpha and an epsilon
The algorithm follows the Anderson,Chung,Lang paper and Dan Spielman's lecture notes
"""
function apr(G::SparseMatrixCSC{Tv,Ti}, s::Array{Int64,1}, alpha::Float64, eps::Float64) where {Tv,Ti}

  # p is ordered by p[i] / deg(i) (ordering required for prn)
  # r is oredered by r[i] - eps * deg(G,u). r is initially mimicking the unit vector
  p = DataStructures.PriorityQueue{Int64,Float64,Base.Order.ReverseOrdering{Base.Order.ForwardOrdering}}(Base.Order.Reverse)
  pelems = Set(s)

  r = DataStructures.PriorityQueue{Int64,Float64,Base.Order.ReverseOrdering{Base.Order.ForwardOrdering}}(Base.Order.Reverse)
  relems = Set(s)

  for u in s
      p[u] = 0
      r[u] = 1 / length(s)^2 - eps * deg(G,u)
  end

  # we are moving mass from a node u only if it has more mass than eps * deg(G,u)
  # making the inequality strictly greater than 0 deals with u being an isolated node
  while DataStructures.peek(r)[2] > 0

    u,ru = DataStructures.peek(r)
    ru = ru + eps * deg(G,u)

    # check if u is in the priority queue for p
    if (u in pelems) == false
      p[u] = 0
      push!(pelems, u)
    end

    # update p[u] & r[u]
    p[u] = p[u] + alpha * ru / deg(G,u)
    r[u] = -eps * deg(G,u) # this means it's set to 0

    for i in 1:deg(G,u)
      v = nbri(G,u,i)

        # check if v is in the priority queue for r
      if (v in relems) == false
        r[v] = -eps * deg(G,v) # this means it's set to 0
        push!(relems, v)
      end

      # update u's neighbors in r
        r[v] = r[v] + (1 - alpha) * ru / deg(G, u)

    end

  end

  return p

end # apr_local


text:
using Test
using LogicCircuits

@testset "Smoothing tests" begin

    c1 = load_logical_circuit(zoo_psdd_file("plants.psdd"))[end]
    c2 = load_logical_circuit(zoo_sdd_file("random.sdd"))[end]
    c3 = smooth(c1)
    c4 = smooth(c2)

    @test !is_smooth(c1)
    @test !is_smooth(c2)
    @test is_smooth(c3)
    @test is_smooth(c4)

    @test c1 !== c3
    @test c2 !== c4

    @test smooth(c3) === c3
    @test smooth(c4) === c4

end


text:
"""
# COMBinatorial Optimization in Julia

The functions listed below are exported, see help on each of these using the ? operator.\n

### blackbox combinatorial optization

`cakewalk` - main entry point for the package. \n

### clique finding in undirected graphs

`clique_cw` - uses cakewalk to find cliques. \n
`is_local_max_clique` - validate that a clique is a locally maximal. \n

### k-medoids

Two standard algorithms for the k-medoids problem are supplied. A third algorithm
that utilizes a greedy search, and whose starting point is optimized by cakewalk
will be added sometime soon. \n

`kmedoids_voronoi` \n
`kmedoids_pam` \n

"""
module Combo

using StatsBase
using Iterators

const Float = Float64

import Base: get
include("utils.jl")
include("Dist.jl")
include("cakewalk.jl")
export cakewalk
include("clique.jl")
export clique_cw, is_local_max_clique
include("kmedoids.jl")
export kmedoids_voronoi, kmedoids_pam

end


text:
# Julia wrapper for header: /Users/jgoldfar/Documents/misc/julia/Hypre/deps/usr/include/HYPRE_DistributedMatrixPilutSolver_protos.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

function HYPRE_NewDistributedMatrixPilutSolver()
    ccall((:HYPRE_NewDistributedMatrixPilutSolver,libHypre),Cint,())
end

function HYPRE_FreeDistributedMatrixPilutSolver()
    ccall((:HYPRE_FreeDistributedMatrixPilutSolver,libHypre),Cint,())
end

function HYPRE_DistributedMatrixPilutSolverInitialize()
    ccall((:HYPRE_DistributedMatrixPilutSolverInitialize,libHypre),Cint,())
end

function HYPRE_DistributedMatrixPilutSolverSetMatrix()
    ccall((:HYPRE_DistributedMatrixPilutSolverSetMatrix,libHypre),Cint,())
end

function HYPRE_DistributedMatrixPilutSolverGetMatrix()
    ccall((:HYPRE_DistributedMatrixPilutSolverGetMatrix,libHypre),Cint,())
end

function HYPRE_DistributedMatrixPilutSolverSetNumLocalRow()
    ccall((:HYPRE_DistributedMatrixPilutSolverSetNumLocalRow,libHypre),Cint,())
end

function HYPRE_DistributedMatrixPilutSolverSetFactorRowSize()
    ccall((:HYPRE_DistributedMatrixPilutSolverSetFactorRowSize,libHypre),Cint,())
end

function HYPRE_DistributedMatrixPilutSolverSetDropTolerance()
    ccall((:HYPRE_DistributedMatrixPilutSolverSetDropTolerance,libHypre),Cint,())
end

function HYPRE_DistributedMatrixPilutSolverSetMaxIts()
    ccall((:HYPRE_DistributedMatrixPilutSolverSetMaxIts,libHypre),Cint,())
end

function HYPRE_DistributedMatrixPilutSolverSetup()
    ccall((:HYPRE_DistributedMatrixPilutSolverSetup,libHypre),Cint,())
end

function HYPRE_DistributedMatrixPilutSolverSolve()
    ccall((:HYPRE_DistributedMatrixPilutSolverSolve,libHypre),Cint,())
end


text:
module Display

using Base.Random: UUID
using Pkg3.Types

export print_project_diff, print_manifest_diff

const colors = Dict(
    ' ' => :white,
    '+' => :light_green,
    '-' => :light_red,
    '↑' => :light_yellow,
    '~' => :light_yellow,
    '↓' => :light_magenta,
    '?' => :red,
)
const color_dark = :light_black

function status(env::EnvCache, mode::Symbol)
    project₀ = project₁ = env.project
    manifest₀ = manifest₁ = env.manifest
    if env.git != nothing
        git_path = LibGit2.path(env.git)
        project_path = relpath(env.project_file, git_path)
        manifest_path = relpath(env.manifest_file, git_path)
        project₀ = read_project(git_file_stream(env.git, "HEAD:$project_path", fakeit=true))
        manifest₀ = read_manifest(git_file_stream(env.git, "HEAD:$manifest_path", fakeit=true))
    end
    if mode == :project || mode == :combined
        # TODO: handle project deps missing from manifest
        m₀ = filter_manifest(in_project(project₀["deps"]), manifest₀)
        m₁ = filter_manifest(in_project(project₁["deps"]), manifest₁)
        info("Status ", pathrepr(env, env.project_file))
        print_diff(manifest_diff(m₀, m₁))
    end
    if mode == :manifest
        info("Status ", pathrepr(env, env.manifest_file))
        print_diff(manifest_diff(manifest₀, manifest₁))
    elseif mode == :combined
        p = !in_project(merge(project₀["deps"], project₁["deps"]))
        m₀ = filter_manifest(p, manifest₀)
        m₁ = filter_manifest(p, manifest₁)
        diff = filter!(x->x.old != x.new, manifest_diff(m₀, m₁))
        if !isempty(diff)
            info("Status ", pathrepr(env, env.manifest_file))
            print_diff(diff)
        end
    end
end

function print_project_diff(env₀::EnvCache, env₁::EnvCache)
    pm₀ = filter_manifest(in_project(env₀.project["deps"]), env₀.manifest)
    pm₁ = filter_manifest(in_project(env₁.project["deps"]), env₁.manifest)
    diff = filter!(x->x.old != x.new, manifest_diff(pm₀, pm₁))
    if isempty(diff)
        print_with_color(color_dark, " [no changes]\n")
    else
        print_diff(diff)
    end
end

function print_manifest_diff(env₀::EnvCache, env₁::EnvCache)
    diff = manifest_diff(env₀.manifest, env₁.manifest)
    diff = filter!(x->x.old != x.new, diff)
    if isempty(diff)
        print_with_color(color_dark, " [no changes]\n")
    else
        print_diff(diff)
    end
end

struct VerInfo
    hash::SHA1
    ver::Union{VersionNumber,Void}
end

vstring(a::VerInfo) =
    a.ver == nothing ? "[$(string(a.hash)[1:16])]" : "v$(a.ver)"

Base.:(==)(a::VerInfo, b::VerInfo) =
    a.hash == b.hash && a.ver == b.ver

≈(a::VerInfo, b::VerInfo) = a.hash == b.hash &&
    (a.ver == nothing || b.ver == nothing || a.ver == b.ver)

struct DiffEntry
    uuid::UUID
    name::String
    old::Union{VerInfo,Void}
    new::Union{VerInfo,Void}
end

function print_diff(io::IO, diff::Vector{DiffEntry})
    same = all(x.old == x.new for x in diff)
    for x in diff
        warnings = String[]
        if x.old != nothing && x.new != nothing
            if x.old ≈ x.new
                verb = ' '
                vstr = vstring(x.new)
            else
                if x.old.hash != x.new.hash && x.old.ver != x.new.ver
                    verb = x.old.ver == nothing || x.new.ver == nothing ||
                           x.old.ver == x.new.ver ? '~' :
                           x.old.ver < x.new.ver  ? '↑' : '↓'
                else
                    verb = '?'
                    msg = x.old.hash == x.new.hash ?
                        "hashes match but versions don't: $(x.old.ver) ≠ $(x.new.ver)" :
                        "versions match but hashes don't: $(x.old.hash) ≠ $(x.new.hash)"
                    push!(warnings, msg)
                end
                vstr = x.old.ver == x.new.ver ? vstring(x.new) :
                    vstring(x.old) * " ⇒ " * vstring(x.new)
            end
        elseif x.new != nothing
            verb = '+'
            vstr = vstring(x.new)
        elseif x.old != nothing
            verb = '-'
            vstr = vstring(x.old)
        else
            verb = '?'
            vstr = "[unknown]"
        end
        v = same ? "" : " $verb"
        print_with_color(color_dark, " [$(string(x.uuid)[1:8])]")
        print_with_color(colors[verb], "$v $(x.name) $vstr\n")
    end
end
print_diff(diff::Vector{DiffEntry}) = print_diff(STDOUT, diff)

function manifest_by_uuid(manifest::Dict)
    entries = Dict{UUID,Dict}()
    for (name, infos) in manifest, info in infos
        uuid = UUID(info["uuid"])
        haskey(entries, uuid) && warn("Duplicate UUID in manifest: $uuid")
        entries[uuid] = merge(info, Dict("name" => name))
    end
    return entries
end

function name_ver_info(info::Dict)
    name = info["name"]
    hash = haskey(info, "hash-sha1") ? SHA1(info["hash-sha1"]) : nothing
    ver = haskey(info, "version") ? VersionNumber(info["version"]) : nothing
    name, VerInfo(hash, ver)
end

function manifest_diff(manifest₀::Dict, manifest₁::Dict)
    diff = DiffEntry[]
    entries₀ = manifest_by_uuid(manifest₀)
    entries₁ = manifest_by_uuid(manifest₁)
    for uuid in union(keys(entries₀), keys(entries₁))
        name₀ = name₁ = v₀ = v₁ = nothing
        haskey(entries₀, uuid) && ((name₀, v₀) = name_ver_info(entries₀[uuid]))
        haskey(entries₁, uuid) && ((name₁, v₁) = name_ver_info(entries₁[uuid]))
        name₀ == nothing && (name₀ = name₁)
        name₁ == nothing && (name₁ = name₀)
        if name₀ == name₁
            push!(diff, DiffEntry(uuid, name₀, v₀, v₁))
        else
            push!(diff, DiffEntry(uuid, name₀, v₀, nothing))
            push!(diff, DiffEntry(uuid, name₁, nothing, v₁))
        end
    end
    sort!(diff, by=x->(x.name, x.uuid))
end

function filter_manifest!(predicate::Function, manifest::Dict)
    empty = String[]
    for (name, infos) in manifest
        filter!(infos) do info
            predicate(name, info)
        end
        isempty(infos) && push!(empty, name)
    end
    for name in empty
        pop!(manifest, name)
    end
    return manifest
end
filter_manifest(predicate::Function, manifest::Dict) =
    filter_manifest!(predicate, deepcopy(manifest))

in_project(deps::Dict) = (name::String, info::Dict) ->
    haskey(deps, name) && haskey(info, "uuid") && deps[name] == info["uuid"]

end # module


text:
module MinimaxAdaptiveControl

using LinearAlgebra
using JuMP

include("mat_comp.jl")

"""
    Candidate{T<:Number}

# Fields:
`A::AbstractMatrix{T}` System matrix

`B::AbstractMatrix{T}` Input Gain matrix

`K::AbstractMatrix{T}` ``H_\\infty``` feedback gain

`P::AbstractMatrix{T}` Stationary solution to the Riccati equation

`hist::Base.RefValue{<:Real}` History, ``\\sum _{k=0}^N \\lambda^{N-k}|x_{k+1} - Ax_k -Bu_k|^2``

`lam::T` Forgetting factor ``\\lambda ``
"""
struct Candidate{T<:Number}
    A::AbstractMatrix{T}
    B::AbstractMatrix{T}
    K::AbstractMatrix{T}
    P::AbstractMatrix{T}
    hist::Base.RefValue{<:Real}
    lam::T
end

"""
    MAController{T<:Number}

# Fields
`candidates::AbstractArray{Candidate{T},1}` Set of system candidates

`xprev::AbstractArray{T,1}` Previous state

`argmin::Base.RefValue{Int64}` Index of minimizing controller
"""

struct MAController{T<:Number}
    candidates::AbstractArray{Candidate{T},1}
    xprev::AbstractArray{T,1}
    γ::Real
    Q::AbstractMatrix{T}
    R::AbstractMatrix{T}
    argmin::Base.RefValue{Int64}
end

# Constructors
"""
    mac = MAController(As, Bs, Q, R, γ, x0, lam)

Create a MAController (minimax adaptive controller) `mac::MAController{T}`
with matrices containing elements of type T, solving
\\min\\_u\\max\\_w,i\\sum (|x_t|^2\\_Q + |u_t|^2\\_R) - \\gamma^2 \\sum |w_t|^2
...
# Arguments
- `As::AbstractVector{<:AbstractMatrix{T}}` A vector of state matrices
- `Bs::AbstractVector{<:AbstractMatrix{T}}` A vector of input gains
- `Q::AbstractMatrix,`                      State penalty matrix |x|_Q^2
- `R::AbstractMatrix,`                      Control penalty matrix |u|_R^2
- `γ::Real,`                                Disturbance penalty, \\H_\\infty gain
- `x0::AbstractVector{<:Number}`            Initial state
- `lam::Real=1`                             Forgetting factor
...
"""
function MAController(
    As::AbstractVector{<:AbstractMatrix{T}}, Bs::AbstractVector{<:AbstractMatrix{T}},
    Q::AbstractMatrix, R::AbstractMatrix, γ::Real,
    x0::AbstractVector{<:Number}, lam=1) where T<:Number
    N = length(As)

    @assert length(As) == length(Bs) "As and Bs" must be of the same length
    @assert length(As) > 0 "As and Bs cannot be nonempty"
    @assert γ > 0 "γ must be strictly positive"
    @assert issemiposdef(Q) "Q must be positive semidefinite"
    @assert isposdef(R) "R nust be positive definite"

    n = size(As[1])[1]  # Number of states
    m = size(Bs[1])[2]  # Number of inputs

    # Extend R to put on standard DARE form
    Rext = [R zeros(m,n)
            zeros(n,m) -γ^2*I(n)]

    candidates = Vector{Candidate{T}}(undef, N)
    for i = 1:N
        A = copy(As[i])
        B = copy(Bs[i])
        # Extend B to put on standard DARE form
        Bext = [B I(n)]
        P = dare(A, Bext, Matrix{Float64}(Q), Rext)
        Kext = (Rext + Bext'*P*Bext)\(Bext'*P*A)
        K = Kext[1:m, 1:n]
        candidates[i] = Candidate(A, B, K,P, Base.RefValue{Float64}(0),T(lam))
    end
    return MAController{T}(candidates,copy(x0),γ,Matrix{T}(Q),Matrix{T}(R), Base.RefValue{Int64}(1))
end

# Controls
"""
    update!(mac,x, u)

Update the internal states of the controller based
on current state `x` and control signal `u`.
...
# Arguments:
- `mac::MAController`       MinimaxAdaptiveControl controller object
- `x::AbstractArray{T,1}`   The current state
- `u::Abstractarray{T,1}`   The previous control signal
...

"""
function update!(mac::MAController,
    x::AbstractArray{T,1}, u::AbstractArray{T,1}) where T<:Number
    lowestcost = Inf
    argmin = 1
    k = 1
    for c in mac.candidates
        c.hist[] += c.lam*(c.A*mac.xprev + c.B*u - x)'*(c.A*mac.xprev + c.B*u - x)
        if c.hist[] < lowestcost
            lowestcost = c.hist[]
            argmin = k
        end
        k +=1
    end
    mac.argmin[] = argmin
    mac.xprev[:] = x
end

"""
    K(mac::MAController)

Select the feedback gain such that u = -Kx

"""
function K(mac::MAController)
    return mac.candidates[mac.argmin[]].K
end

"""
    Vbar(mac, T, x)

Compute the current upper bound of the value function.

...
# Arguments
- `mac::MAController`: Controller object
- `T::AbstractMatrix`: T - matrix, can be synthesized using Tval(...)
- `x::AbstractArray`: Next state
...
"""
function Vbar(mac::MAController, T::AbstractMatrix, x::AbstractArray)
    N = length(mac.candidates)
    maxval = -Inf
    minhist = Inf
    secminhist = Inf
    for i = 1:N
        c = x'*mac.candidates[i].P*x -mac.γ^2*mac.candidates[i].hist[]
        maxval < c ? maxval = c : nothing
        hist = mac.candidates[i].hist[]
        hist < secminhist ? secminhist = hist : nothing
        hist <= minhist ? (minhist, secminhist) = (hist, minhist) : nothing
    end
    d =
    return max(maxval, x'*T*x - (mac.γ^2)/2*(minhist + secminhist))
end

"""
    T(mac, model, tol = 0.01)
synthesize a common T using convex programming such that
inequalities (19) and (20) are fulfilled.
...
# Arguments:
- `mac::MAController` MinimaxAdaptiveControl controller object
- `model::JuMP.Model` A user supplied JuMP model. Currently the solvers Mosek and Hypatia works.
- `tol::Real = 0.01`  Forcing (19) and (20) to hold with margin, i.e. T - tol*I >= ...
...
"""
function Tsyn(mac::MAController{M}, model::JuMP.Model, tol::Real = 0.01) where M<:Number
    N_systems = length(mac.candidates)
    n = size(mac.Q)[1]
    @variable(model, T[1:n, 1:n] in PSDCone())
    @variable(model,t)
    @constraint(model, Symmetric(mac.γ^2*I(n)-T) in PSDCone())
    for i = 1:N_systems
        @constraint(model, Symmetric(T - mac.candidates[i].P) in PSDCone())
        for j = 1:N_systems
            @constraint(model, X(mac, T, i,j, tol) in PSDCone())
            if i!=j
                for k = 1:N_systems
                    @constraint(model, Z(mac, T, i,j,k, tol) in PSDCone())
                end
            end
        end
    end

    @constraint(model, Symmetric(t*I(n) - T) in PSDCone()) # Regularize the operator norm of T
    @objective(model, Min, t)
    optimize!(model)
    T = value.(T)
    stat = termination_status(model)
    return T, stat
end

"""
    X(mac,T,i, k, tol) where P<:Number

Synthesize a symmetric matrix X such that inequality (19) is satisfied iff
X is positive semidefinite.
...
# Arguments:
- `mac::MAController`   MinimaxAdaptiveControl controller object
- `T::AbstractMatrix`   Either a matrix or a JuMP variable
- `i::Integer`          Index variable,
- `k::Integer`          Index Variable
- `tol::Real = 0.01`    Forcing (19) to hold with margin, i.e. T - tol*I >= ...
...
"""
function X(mac::MAController{P},
    T::AbstractMatrix,
    i::Integer, k::Integer, tol::Real = 0.01) where P<:Number
    Q = mac.Q
    R = mac.R
    Kk = mac.candidates[k].K
    Ai = mac.candidates[i].A
    Bi = mac.candidates[i].B
    Pi = mac.candidates[i].P

    Aik = Ai - Bi*Kk

    Xik = [(T - Q - Kk'*R*Kk)-tol*I(size(Pi)[1]) Aik'
            Aik (inv(Pi) - mac.γ^(-2)*I(size(Pi)[1]))]
    return Xik
end

"""
    Z(mac, T, i, j, k)

Synthesize a symmetric matrix Z such that inequality (20) is satisfied iff
Z is positive semidefinite.

...
# Arguments:
- `mac::MAController`   MinimaxAdaptiveControl controller object
- `T::AbstractMatrix`   Either a matrix or a JuMP variable
- `i::Integer`          Index variable
- `j::Integer`          Index variable
- `k::Integer`          Index Variable
- `tol::Real = 0.01`    Forcing (20) to hold with margin, i.e. T - tol*I >= ...
...
"""
function Z(mac::MAController{P},
    T::AbstractMatrix,
    i::Integer, j::Integer, k::Integer, tol::Real = 0.01) where P<:Number

    Q = mac.Q
    R = mac.R
    γ = mac.γ
    Kk = mac.candidates[k].K
    Ai = mac.candidates[i].A
    Bi = mac.candidates[i].B
    Pi = mac.candidates[i].P

    Aj = mac.candidates[j].A
    Bj = mac.candidates[j].B
    Pj = mac.candidates[j].P

    Aik = Ai - Bi*Kk
    Ajk = Aj - Bj*Kk

    Zijk11 = T - Q - Kk'*R*Kk + γ^2/2*(Aik'*Aik + Ajk'*Ajk) -tol*I(size(Pi)[1])
    Zijk21 = (Aik + Ajk)
    Zijk22 = 4/γ^4*(γ^2*I(size(Pi)[1]) - T)
    Zijk = [Zijk11 Zijk21'
            Zijk21 Zijk22]
    return Zijk
end

export Candidate, MAController, update!, K, Tsyn, Vbar
end


text:
module adapters_postgresql_multiple_databases_test

using Test # @test
using Octo.Adapters.SQLite # Repo Schema from SELECT FROM WHERE

import SQLite # pathof(SQLite)
dbfile = joinpath(dirname(pathof(SQLite)), "..", "test", "Chinook_Sqlite.sqlite")
dbfile2 = joinpath(@__DIR__, "test.sqlite")
cp(dbfile, dbfile2; force=true)
chmod(dbfile2, 0o666)

Repo.debug_sql()

sqlite_connector = Repo.connect(
    adapter = Octo.Adapters.SQLite,
    dbfile = dbfile2,
    multiple = true,
)

using Octo.Adapters.PostgreSQL # Repo Schema Raw

include("options.jl")

pg_connector = Repo.connect(;
    adapter = Octo.Adapters.PostgreSQL,
    multiple = true,
    Options.for_postgresql...
)

struct Employee
end
Schema.model(Employee, table_name="Employee", primary_key="ID")

for c in (sqlite_connector, pg_connector)
    Repo.execute([DROP TABLE IF EXISTS Employee], db=c)
end
Repo.execute(Raw("""
    CREATE TABLE Employee (
        ID INTEGER PRIMARY KEY,
        Name TEXT NOT NULL,
        Salary NUMERIC NOT NULL
    )"""), db=sqlite_connector)
Repo.execute(Raw("""
    CREATE TABLE Employee (
        ID SERIAL,
        Name VARCHAR(255),
        Salary FLOAT(8),
        PRIMARY KEY (ID)
    )"""), db=pg_connector)

result = Repo.insert!(Employee, (Name = "Jessica", Salary= 70000.50); db=sqlite_connector)
@test keys(result) == (:id, :num_affected_rows)
@test result.num_affected_rows == 1
result = Repo.insert!(Employee, (Name = "Jessica", Salary= 70000.50); db=pg_connector)
@test keys(result) == (:id, :num_affected_rows)
@test result.num_affected_rows == 1

Repo.query(Employee, db=sqlite_connector)
Repo.query(Employee, db=pg_connector)

Repo.disconnect(db=sqlite_connector)
Repo.disconnect(db=pg_connector)

end # module adapters_postgresql_multiple_databases_test


text:
function neural_ode(model,x,tspan, solver,
                    args...;kwargs...)
  p = destructure(model)
  dudt_(u::TrackedArray,p,t) = restructure(model,p)(u)
  dudt_(u::AbstractArray,p,t) = Flux.data(restructure(model,p)(u))
  prob = ODEProblem(dudt_,x,tspan,p)
  return diffeq_adjoint(p,prob,solver,args...;u0=x,kwargs...)
end

function neural_ode_rd(model,x,tspan,
                       args...;kwargs...)
  dudt_(u,p,t) = model(u)
  _x = x isa TrackedArray ? x : param(x)
  prob = ODEProblem(dudt_,_x,tspan)
  # TODO could probably use vcat rather than collect here
  solve(prob, args...; kwargs...) |> Tracker.collect
end

function neural_dmsde(model,x,mp,tspan,
                      args...;kwargs...)
  dudt_(u,p,t) = model(u)
  g(u,p,t) = mp.*u
  prob = SDEProblem(dudt_,g,param(x),tspan,nothing)
  # TODO could probably use vcat rather than collect here
  solve(prob, args...; kwargs...) |> Tracker.collect
end

# Flux Layer Interface
struct NeuralODE
    model
    tspan
    solver
    args
    kwargs
end

# Optional args and kwargs
NeuralODE(model,tspan,solver,kwargs)= NeuralODE(model,tspan,solver,(),kwargs)
NeuralODE(model,tspan,solver)= NeuralODE(model,tspan,solver,(),())

# Play nice with Flux
Flux.@treelike NeuralODE
Flux.params(n::NeuralODE) = params(n.model)

function (n::NeuralODE)(x)
    return neural_ode(n.model,x,n.tspan,n.solver,n.args...;n.kwargs...)
end


text:
using Statistics
#using Plots
include("ColorCode.jl")

function simulate(str,error_rate)
  assignment = Dict([(k,1) for k in keys(keyboardStrings)])

  lmModel = LM.getModel()
  lmState = LM.getStartState(lmModel)
  prior = getPrior(lmModel,lmState)
  belief = Belief(prior,9,1)
  history = BeliefHistory()
  certaintyThreshold = 0.95
  changeAssignment(belief,assignment)

  letters = Stack{Symbol}()
  foreach(c->push!(letters,c==' ' ? :SPACE : Symbol(c)), reverse(str))
  commString = ""
  totalClicks = 0.0
  print("starting simulate\n")
  while !isempty(letters)
    letter = pop!(letters)
    clickCount = 0.0
    while true
      if rand() < error_rate
        color = assignment[letter]==1 ? 2 : 1
      else
        color = assignment[letter]
      end
      clickCount += 1.0
      updateBelief(belief,color,assignment)
      newCommString = chooseLetter(belief,commString,certaintyThreshold,history,lmModel,lmState)
      changeAssignment(belief,assignment)
      if length(newCommString) != length(commString)
        new_letter = nothing
        if length(newCommString) > length(commString)
          new_letter = newCommString[end] == ' ' ? :SPACE : Symbol(newCommString[end])
          #print("selected: $(new_letter), click count: $(clickCount)\n")
        end
        #length(newCommString) < length(commString) ? print("selected: $(letter), click count: $(clickCount)\n") : nothing
        commString = newCommString
        if !isnothing(new_letter) && new_letter != letter
          push!(letters,letter)
          push!(letters,:UNDO)
        end
        break
      end
    end
    totalClicks += clickCount
  end
  LM.releaseModel(lmModel)
  LM.releaseState(lmState)
  clicksPerChar = totalClicks/length(str)
  print("error_rate=$(error_rate), average clicks per letter: $(clicksPerChar)\n")
  return clicksPerChar
end

function simulate(str)
  h2(x) = x*log(2,1/x) + (1-x)*log(2,1/(1-x))
  error_rates = 0.01:0.01:0.45
  base = simulate(str,0)
  #p1 = plot(x->1/((1-h2(x))/base),xlim=[0 0.5],legend=false,xticks=error_rates,yticks=0:0.5:35,minorgrid=true)
  #gui()
  results = []
  for r in error_rates
    res = simulate(str,r)
    push!(results,res)
    #plot!(p1,[r],[res],seriestype=:scatter)
    #gui()
  end
  print("results: \n")
  print(results)

  #p2 = plot(x->1-h2(x),xlim=[0 0.5],xticks=0:0.05:0.5,yticks=0:0.1:1.0,minorgrid=true)
  #i = base ./ results
  #plot!(p2,error_rates,i,seriestype=:scatter)
  #gui(p2)
  #wait()
end

if length(ARGS) != 2 && length(ARGS) != 1
    error("usage: julia simulator.jl <infile>.txt <error_rate>")
end

inputfilename = ARGS[1]
error_rate = nothing
if length(ARGS) == 2
  error_rate = parse(Float64,ARGS[2])
end

open(inputfilename) do f
  str = chomp(uppercase(read(f,String)))
  if !isnothing(error_rate)
    simulate(str,error_rate)
  else
    simulate(str)
  end
end

max_stars_repo_name: mrdaly/ColorCode

max_issues_repo_name: mrdaly/ColorCode

max_forks_repo_name: mrdaly/ColorCode


text:
# Simulate MUSIC with a linear array

# before running:
#  include("../src/DirectionFinding.jl")
#  using <favorite plotting program>

using DirectionFinding

begin
    fc = 1e9
    λ = 299.79e6/fc
    N₀ = 0.01
    a1 = 5; a2 = 1.5;
    la = lineararray(11, 2λ, σ²=0.001)
    s1 = SignalSource(a1, halfsine_train(1e6))
    s2 = SignalSource(a2, halfsine_train(1e6))
    sim = Simulation(la, fc, N₀, s1, s2)
    mps = MUSIC(sim, 2.00001e6, 1000)
    ϕ = -π/2:0.01:π/2
    peaks = [p<0 ? p+2π : p for p in findpeaks(mps, ϕ, n=2)]
    println("Angles: $a1 -- $a2")
    println("Found: $(peaks[1]) -- $(peaks[2])")
    plot(ϕ, mps.(ϕ))
end


text:
#===
・サブ配列の転置は，相手側を転置して，積を再転置する方が早い
・similar()は，計算に用いらない場合に使う（値の置換など）
===#

#===
ファイルの読み込み
===#
# include("seqential.jl")

export Affine,LSTM,Sigmoid_with_loss,Dropout

#===
共通型，関数
・全てのレイヤはLayersの子型になる
・forward!(this,x,learn_flg)
・backward!(this,din)
・gpu_flg は add_layer!()で処理してる
===#

#===
計算
・計算式をつなげる場合は@.を付ける
====#

#===
cacheは各structを作成し，保持
===#

#=====Affine==============#
mutable struct Affine <: Layers
    params #W,b
    grads #dW,db
    in
    gpu_flg
    function Affine(input,output)
        W = randn((input,output))
        b = zeros((1,output))
        dW = zeros(size(W))
        db = zeros(size(b))
        new([W,b],[dW,db],nothing,false)
    end
end
function forward!(this::Affine,in,learn_flg)
    W,b = this.params
    this.in = in
    in * W .+ b
end
function backward!(this::Affine,din)
    # mul!(this.grads[1],(this.in)',din)
    this.grads[1] = (din' * this.in)'
    this.grads[2] = sum(din,dims=1)
    din * this.params[1]'
end
function reset!(this::Affine)
    nothing
end
function convert_to_cu!(this::Affine)
    this.params = cu.(this.params)
    this.grads = cu.(this.grads)
    this.in = nothing
    this.gpu_flg = true
    nothing
end
function convert_to_array!(this::Affine)
    this.params = Array.(this.params)
    this.grads = Array.(this.grads)
    this.in = nothing
    this.gpu_flg = false
    nothing
end

#====uni_LSTM=============#
mutable struct uni_LSTM <: Layers
    params #W,b
    grads #dW,db
    x
    h_prev
    c_prev
    i
    f
    g
    o
    c_next
    gpu_flg
    cache
    function uni_LSTM(params,grads,cache,gpu_flg)
        new(params,grads,nothing,nothing,nothing,nothing,nothing,
        nothing,nothing,nothing,gpu_flg,cache)
    end
end
function forward!(this::uni_LSTM,x,h_prev,c_prev)
    #x{N,D}, h_prev{N,H}, c_prev{N,H}
    Wx, Wh, b = this.params #{D*4H,H*4H,1*4H}
    gpu_flg = this.gpu_flg
    N = size(x,1)
    H = size(h_prev,2)

    @unpack A,x_Wx,h_Wh = this.cache

    #N,4H
    if A === nothing || size(A) != (N,4H)
        A = gpu_flg ? CUDA.zeros(N,4H) : zeros(N,4H)
        x_Wx = gpu_flg ? CUDA.zeros(N,4H) : zeros(N,4H)
        h_Wh = gpu_flg ? CUDA.zeros(N,4H) : zeros(N,4H)
        this.cache.A = A
        this.cache.x_Wx = x_Wx
        this.cache.h_Wh = h_Wh
    end

    mul!(x_Wx,x,Wx)
    mul!(h_Wh,h_prev,Wh)
    A .= x_Wx .+ h_Wh .+ b

    #それぞれ{N,H}
    this.f = 1 ./ (1 .+ exp.(-(view(A,:,1:H)))) #Sigmoid
    this.g = tanh.(view(A,:,H+1:2*H))
    this.i = 1 ./ (1 .+ exp.(-(view(A,:,2*H+1:3*H))))
    this.o = 1 ./ (1 .+ exp.(-(view(A,:,3*H+1:size(A,2)))))

    c_next = this.f .* c_prev .+ this.g .* this.i #N,H

    h_next = this.o .* tanh.(c_next) #N,H

    this.x = x
    this.h_prev = h_prev #N,H
    this.c_prev = c_prev #N,H
    this.c_next = c_next

    return h_next, c_next
end
function backward!(this::uni_LSTM,dh_next,dc_next)
    @unpack A = this.cache #dA

    gpu_flg = this.gpu_flg
    Wx, Wh, b = this.params
    i = this.i
    f = this.f
    g = this.g
    o = this.o

    tanh_c_next = tanh.(this.c_next)
    ds = dc_next .+ (dh_next .* o) .* (1 .- tanh_c_next.^2)

    dc_prev = ds .* f

    N,H = size(i)
    #N,4H
    if A === nothing || size(A) != (N,4H)
        A = gpu_flg ? CUDA.zeros(N,4H) : zeros(N,4H)
        this.cache.A = A
    end

    A[:,2H+1:3H] = ds .* g .* i .* (1 .- i)
    A[:,1:H] = ds .* this.c_prev .* f .* (1 .- f)
    A[:,3H+1:4H] = dh_next .* tanh_c_next .* o .* (1 .- o)
    A[:,H+1:2H] = ds .* i .* (1 .- g.^2)

    # this.A .= hcat(df,dg,di,d_o)

    this.grads[1] = (A' * this.x)'
    this.grads[2] = this.h_prev' * A

    this.grads[3] .= sum(A,dims=1)

    dx = A * Wx'
    dh_prev = A * Wh'

    return dx, dh_prev, dc_prev
end
function reset!(this::uni_LSTM)
    nothing
end

#====LSTM===================#
#===
Dropoutを引数に渡すと，時系列方向にも適応できる(p.265)

===#
mutable struct LSTM_cache
    A
    x_Wx
    h_Wh
    function LSTM_cache()
        new(nothing,nothing,nothing)
    end
end
mutable struct LSTM
    params #Wx,Wh,b
    grads #dWx,dWh,db
    layers #複数のLSTMレイヤを管理
    h #次のエポックに引き継ぐ
    c
    stateful
    out_sequence
    dxs
    dropout
    gpu_flg
    cache
    function LSTM(input,output;stateful=true,out_sequence=true,dropout=nothing)
        Wx = randn(input,4*output)/sqrt(input)
        Wh = randn(output,4*output)/sqrt(output)
        b =  zeros(1,4*output)

        dWx = zeros(size(Wx))
        dWh = zeros(size(Wh))
        db =  zeros(size(b))
        new(
            [Wx,Wh,b],
            [dWx,dWh,db],
            nothing,
            nothing,
            nothing,
            stateful,
            out_sequence,
            nothing,
            dropout,
            false,
            LSTM_cache()
        )
    end
end
function forward!(this::LSTM,xs,learn_flg)
    stateful = this.stateful
    gpu_flg = this.gpu_flg
    #バッチサイズ，ブロック数，入力データ数
    N, T, D = size(xs)
    H = size(this.params[2],1) #Wh{H,4H}

    hs = similar(xs,N,T,H)

    if !stateful || this.h === nothing || typeof(this.h)
        this.h = gpu_flg ? CUDA.zeros(N,H) : zeros(N,H)
    end
    if !stateful || this.c === nothing
        this.c = gpu_flg ? CUDA.zeros(N,H) : zeros(N,H)
    end

    #uni_LSTMの初期化
    this.layers = Array{uni_LSTM}(undef,T)

    #Dropoutのマスクのシャッフル(生成)
    if this.dropout !== nothing
        make_mask!(this.dropout,(N,H))
    end

    @inbounds for t in 1:T
        this.layers[t] = uni_LSTM(this.params,this.grads.*0,this.cache,gpu_flg)
        this.h, this.c = forward!(this.layers[t],view(xs,:,t,:),this.h,this.c)
        hs[:,t,:] = this.h #次レイヤへの伝播用
        if this.dropout !== nothing
            #dropoutをかけるのはuni_LSTMに渡す値のみ(次レイヤへは別のdropoutを使う)
            this.h = forward!(this.dropout,this.h,learn_flg,false)
            this.c = forward!(this.dropout,this.c,learn_flg,false)
        end
    end

    return this.out_sequence ? hs : hs[:,end,:] #hs{N,T,H}
end
function backward!(this::LSTM,dhs)
    Wx, Wh, b = this.params
    D = size(Wx,1) #{D,4H}
    N = size(dhs,1) #dhs{N,H} || {N,T,H}
    T = length(this.layers)
    this.dxs = similar(dhs,(N,T,D))
    this.grads *= 0 #勾配初期化

    if this.out_sequence
        backward_dims_3!(this,dhs,T) #dhs{N,T,H}
    else
        backward_dims_2!(this,dhs,T) #dhs{N,H}
    end

    return this.dxs
end
function backward_dims_2!(this::LSTM,dh,T)
    dc = this.gpu_flg ? CUDA.zeros(size(dh)) : zeros(size(dh))
    @inbounds for t in T:-1:1
        uni_layer = this.layers[t]
        this.dxs[:,t,:], dh, dc = backward!(uni_layer,dh,dc)
        #勾配合算
        this.grads += uni_layer.grads
        #dropoutの誤差逆伝播
        if this.dropout !== nothing
            dh = backward!(this.dropout,dh)
            dc = backward!(this.dropout,dc)
        end
    end
    nothing
end
function backward_dims_3!(this::LSTM,dhs,T)
    dh = this.gpu_flg ? CUDA.zeros(size(dhs[:,1,:])) : zeros(size(dhs[:,1,:]))
    dc = this.gpu_flg ? CUDA.zeros(size(dhs[:,1,:])) : zeros(size(dhs[:,1,:]))
    @inbounds @simd for t in T:-1:1
        uni_layer = this.layers[t]
        this.dxs[:,t,:], dh, dc = backward!(uni_layer,view(dhs,:,t,:)+dh, dc)
        #勾配合算
        this.grads += uni_layer.grads
        #dropoutの誤差逆伝播
        if this.dropout !== nothing
            dh = backward!(this.dropout,dh)
            dc = backward!(this.dropout,dc)
        end
    end
    nothing
end
function reset!(this::LSTM)
    this.c = nothing
    this.h = nothing
end
function convert_to_cu!(this::LSTM)
    this.params = cu.(this.params)
    this.grads = cu.(this.grads)
    this.gpu_flg = true
    this.layers = nothing
    this.dxs = nothing
    #==cache==#
    if this.cache.A !== nothing
        this.cache.A = cu(this.cache.A)
        this.cache.x_Wx = cu(this.cache.x_Wx)
        this.cache.h_Wh = cu(this.cache.h_Wh)
    end

    if this.c !== nothing
        this.c = cu(this.c)
        this.h = cu(this.h)
    end
    if this.dropout !== nothing
        convert_to_cu!(this.dropout)
    end
    nothing
end
function convert_to_array!(this::LSTM)
    this.params = Array.(this.params)
    this.grads = Array.(this.grads)
    this.gpu_flg = false
    this.layers = nothing
    this.dxs = nothing
    #==cache==#
    if this.cache.A !== nothing
        this.cache.A = Array(this.cache.A)
        this.cache.x_Wx = Array(this.cache.x_Wx)
        this.cache.h_Wh = Array(this.cache.h_Wh)
    end
    if this.c !== nothing
        this.c = Array(this.c)
        this.h = Array(this.h)
    end
    if this.dropout !== nothing
        convert_to_array!(this.dropout)
    end
    nothing
end

#====Sigmoid_with_loss=================#
mutable struct Sigmoid_with_loss
    params
    grads
    s #スコア
    t
    gpu_flg
    Sigmoid_with_loss() = new(nothing,nothing,nothing,nothing,false)
end
function forward!(this::Sigmoid_with_loss,in,learn_flg)
    s = 1 ./ (1 .+ exp.(-in))
    delta = 1e-7
    this.s = s
    if !learn_flg
        #ただの推論
        return s
    end
    l = -this.t .* log.(delta.+s) .- (1 .-this.t) .* log.(1+delta .-s)
    return l
end
function backward!(this::Sigmoid_with_loss,din)
    return this.s .- this.t
end
function reset!(this::Sigmoid_with_loss)
    this.t = nothing
end
function convert_to_cu!(this::Sigmoid_with_loss)
    this.gpu_flg = true
    this.s = nothing
    this.t = nothing
    nothing
end
function convert_to_array!(this::Sigmoid_with_loss)
    this.gpu_flg = false
    this.s = nothing
    this.t = nothing
    nothing
end

#====dropout===========================#
#===
・LSTMに適応する場合
    通常：LSTMの出力に，異なるmaskをかける
        LSTM
        ↑
        Dropout(0.5) dims=3
        ↑
        LSTM

    変分：LSTM内部に共通のmask，外部も共通mask(生成したmaskの[:,1,:]を使う)
    (内部でdropoutを使用したLSTMの入出力で使われる)
        LSTM(Dropout(0.5)) dims=2
        ↑
        Dropout(0.5,true) dims=2,各時刻に同じmask
        ↑
        LSTM(Dropout(0.5)) dims=2

===#
mutable struct Dropout
    params
    grads
    mask
    ratio
    gpu_flg
    variational_flg #変分=>T方向へ
    function Dropout(ratio,variational_flg=false)
        new(nothing,nothing,nothing,ratio,false,variational_flg)
    end
end
function forward!(this::Dropout,x,learn_flg,shuffle_flg=true)
    if !learn_flg
        return  x .*= (1.0 - this.ratio) #推論
    end

    #LSTM内部のT方向の処理はshuffleしない（代わりに計算前に行う）
    if shuffle_flg || size(this.mask) != size(x)
        make_mask!(this,size(x))
    end

    return x .*= this.mask #学習
end
function backward!(this::Dropout,dx)
    return dx .*= this.mask #通常
end
function reset!(this::Dropout)
    nothing
end
function convert_to_cu!(this::Dropout)
    this.gpu_flg = true
    if this.mask !== nothing
        this.mask = cu(this.mask)
    end
    nothing
end
function convert_to_array!(this::Dropout)
    this.gpu_flg = false
    if this.mask !== nothing
        this.mask = Array(this.mask)
    end
    nothing
end
function make_mask!(this::Dropout,s)
    #新しく作る
    if this.mask === nothing || s != size(this.mask)
        if this.variational_flg
            this.mask = rand(Float64,s[1],1,s[3]) #変分
        else
            this.mask = rand(Float64,s) #通常
        end
        this.mask = map(x->(x.>this.ratio ? 1.0 : 0.0),this.mask)
        if this.gpu_flg
            this.mask = cu(this.mask)
        end
    else
    #shuffle
        if this.gpu_flg
            this.mask = Array(this.mask)
        end
        this.mask = shuffle(this.mask)
        if this.gpu_flg
            this.mask = cu(this.mask)
        end
    end
end

max_stars_repo_name: catmanium/Spike

max_issues_repo_name: catmanium/Spike

max_forks_repo_name: catmanium/Spike


text:
# Solution to Q4e (implement a CrossEntropyLog)
import Base

type CrossEntropyLog{T}
    # Buffer of the targets, needed during fprop and bprop.
    t::Matrix{T}

    CrossEntropyLog() = new(zeros(T, 0, 0))
end

target!{T}(l::CrossEntropyLog{T}, t::Matrix{T}) = l.t = t

function fprop!{T}(l::CrossEntropyLog{T}, x::Matrix{T})
    # No need to remember anything here.
    return -sum(l.t .* x, 2)
end

function bprop!{T}(l::CrossEntropyLog{T}, ∂y::Matrix{T})
    # This expression comes from Q4d. Notice the absence of division,
    # which means there's no potential for instability!
    return ∂y .* -l.t
end

params(::CrossEntropyLog) = ()
grads(::CrossEntropyLog) = ()

Base.show(io::IO, l::CrossEntropyLog) = print(io, "CrossEntropy(Log)")


text:
module SamplePaths

using SampleSimpleFail

end # module


text:
# ------------------------------------------------------------------
# Copyright (c) 2017, Júlio Hoffimann Mendes <juliohm@stanford.edu>
# Licensed under the ISC License. See LICENCE in the project root.
# ------------------------------------------------------------------

function digest(solution::SimulationSolution{<:RegularGrid})
  # solution variables
  variables = collect(keys(solution.realizations))

  # get the size of the grid
  sz = size(solution.domain)

  # build dictionary pairs
  pairs = []
  for var in variables
    reals = map(r -> reshape(r, sz), solution.realizations[var])
    push!(pairs, var => reals)
  end

  # output dictionary
  Dict(pairs)
end