text:
#=
lm-aa
Core methods for LinearMapAA objects.
2019-08-07 Jeff Fessler, University of Michigan
=#

using LinearMaps: LinearMap
using LinearAlgebra: UniformScaling, I
import LinearAlgebra: issymmetric, ishermitian, isposdef
#import LinearAlgebra: mul! #, lmul!, rmul!
import SparseArrays: sparse

export redim, undim

# copy
Base.copy(A::LinearMapAX{T,Do,Di}) where {T,Do,Di} =
	LinearMapAA(A._lmap ; prop=A._prop, T=T, idim=A._idim, odim=A._odim)

# Matrix
Base.Matrix(A::LinearMapAX) = Matrix(A._lmap)

# ndims
# Base.ndims(A::LinearMapAX) = ndims(A._lmap) # 2 for AbstractMatrix
Base.ndims(A::LinearMapAO) = ndims(A._lmap) # 2

"""
    show(io::IO, A::LinearMapAX)
    show(io::IO, ::MIME"text/plain", A::LinearMapAX)
Pretty printing for `display`
"""
function Base.show(io::IO, A::LinearMapAX) # short version
	print(io, isa(A, LinearMapAM) ? "LinearMapAM" : "LinearMapAO")
	print(io, ": $(size(A,1)) × $(size(A,2))")
end

# multi-line version:
Base.show(io::IO, ::MIME"text/plain", A::LinearMapAX{T,Do,Di}) where {T,Do,Di} =
	begin
		show(io, A)
		(A._prop != NamedTuple()) && print(io, "\n$(A._prop)")
		print(io, "\nodim=$(A._odim) idim=$(A._idim) T=$T Do=$Do Di=$Di")
	#	print(io, "\n$(A._lmap)\n") # todo: hide until "show" fixed for LM
		print(io, "\n$(typeof(A._lmap)) ...\n")
	#	tmp = "$(A._lmap)"[1:77]
	#	print(io, " $tmp ..")
	end

# size
Base.size(A::LinearMapAX) = size(A._lmap)
Base.size(A::LinearMapAX, d::Int) = size(A._lmap, d)

"""
    redim(A::LinearMapAX ; idim::Dims=A._idim, odim::Dims=A._odim)

"Reinterpret" the `idim` and `odim` of `A`
"""
function redim(A::LinearMapAX{T} ;
	idim::Dims=A._idim, odim::Dims=A._odim) where {T}

	prod(idim) == prod(A._idim) || throw("incompatible idim")
	prod(odim) == prod(A._odim) || throw("incompatible odim")
	return LinearMapAA(A._lmap ; prop=A._prop, T=T, idim=idim, odim=odim)
end

"""
    undim(A::LinearMapAX)

"Reinterpret" the `idim` and `odim` of `A` to be of AM type
"""
undim(A::LinearMapAX{T}) where {T} =
	LinearMapAA(A._lmap ; prop=A._prop, T=T)

# adjoint
Base.adjoint(A::LinearMapAX) = LinearMapAA(adjoint(A._lmap) ;
    prop=A._prop, idim=A._odim, odim=A._idim, operator=isa(A,LinearMapAO))

# transpose
Base.transpose(A::LinearMapAX) = LinearMapAA(transpose(A._lmap) ;
    prop=A._prop, idim=A._odim, odim=A._idim, operator=isa(A,LinearMapAO))

# eltype
Base.eltype(A::LinearMapAX) = eltype(A._lmap)

# LinearMap algebraic properties
issymmetric(A::LinearMapAX) = issymmetric(A._lmap)
#ishermitian(A::LinearMapAX{<:Real}) = issymmetric(A._lmap)
ishermitian(A::LinearMapAX) = ishermitian(A._lmap)
isposdef(A::LinearMapAX) = isposdef(A._lmap)

# comparison of LinearMapAX objects, sufficient but not necessary
Base.:(==)(A::LinearMapAX, B::LinearMapAX) =
	eltype(A) == eltype(B) &&
		A._lmap == B._lmap && A._prop == B._prop &&
		A._idim == B._idim && A._odim == B._odim

# convert to sparse
sparse(A::LinearMapAX) = sparse(A._lmap)

# add or subtract objects (with compatible idim,odim)
function Base.:(+)(A::LinearMapAX, B::LinearMapAX)
	(A._idim != B._idim) && throw("idim mismatch in +")
	(A._odim != B._odim) && throw("odim mismatch in +")
	LinearMapAA(A._lmap + B._lmap ;
		idim=A._idim, odim=A._odim,
		prop = (sum=nothing,props=(A._prop,B._prop)),
		operator = isa(A, LinearMapAO) & isa(B, LinearMapAO),
	)
end

# Allow LMAA + AM only if Do=Di=1
function Base.:(+)(A::LinearMapAX, B::AbstractMatrix)
    (length(A._idim) != 1 || length(A._odim) != 1) && throw("use redim")
    LinearMapAA(A._lmap + LinearMap(B) ; prop=A._prop,
        operator = isa(A, LinearMapAO),
    )
end

# But allow LMAA + I for any Do,Di
Base.:(+)(A::LinearMapAX, B::UniformScaling) = # A + I -> A + I(N)
    LinearMapAA(A._lmap + B(size(A,2)) ;
        prop = (Aprop=A._prop, Iscale=B(size(A,2))[1]),
        idim = A._idim,
        odim = A._odim,
        operator = isa(A, LinearMapAO),
    )
Base.:(+)(A::AbstractMatrix, B::LinearMapAX) = B + A

Base.:(-)(A::LinearMapAX, B::LinearMapAX) = A + (-1)*B
Base.:(-)(A::LinearMapAX, B::AbstractMatrix) = A + (-1)*B
Base.:(-)(A::AbstractMatrix, B::LinearMapAX) = A + (-1)*B

# A.?
Base.getproperty(A::LinearMapAX, s::Symbol) =
	(s in LMAAkeys) ? getfield(A, s) :
#	s == :m ? size(A._lmap, 1) :
	haskey(A._prop, s) ? getfield(A._prop, s) :
		throw("unknown key $s")

Base.propertynames(A::LinearMapAX) = (propertynames(A._prop)..., LMAAkeys...)


text:
"""
Support for nanomsg library
"""

function nn_strerror(errno::Int32)
    err = ccall((:nn_strerror, "libnanomsg"), Cstring, (Int32,), errno)
    if err == C_NULL
        error("Error nn_strerror errno: ", errno)
    end
    return unsafe_string(err)
end
nn_strerror(errno::Int64) = nn_strerror(Int32(errno))

nn_errno() = ccall((:nn_errno, "libnanomsg"), Int32, ())

function check_err(func_name, rv)
    if rv < 0
        err = nn_strerror(nn_errno())
        error(func_name, err)
    end
end

"""
    NNSocket(protocol; raw=false)

Create an SP socket of type `protocol`.

* `raw`: Set to true to omit end-to-end functionality during creation.
  Used to implement intermediary devices in SP topologies. (default: `false`)
"""
function nn_socket(protocol; raw=false)
    domain = !raw ? 1 : 2 # AF_SP=1 AF_SP_RAW = 2
    fd = ccall((:nn_socket, "libnanomsg"), Int32, (Int32, Int32), domain, protocol)
    return fd
end

"""
    nn_bind(s, addr)

Bind local endpoint to socket s.

Returns positive endpoint ID, else -1 for failure.
"""
nn_bind(s, addr) = ccall((:nn_bind, "libnanomsg"), Int32, (Int32, Ptr{UInt8}), s, addr)

"""
Unlike with traditional BSD sockets, this function operates asynchronously,
and returns to the caller before the operation is complete.

Returns positive endpoint ID, else -1 for failure.
"""
nn_connect(sock, addr) = ccall((:nn_connect, "libnanomsg"), Int32, (Int32, Ptr{UInt8}), sock, addr)

function nn_recv(s ; blocking=true)
    # define NN_MSG ((size_t) -1)
    flags = blocking ? 0 : 1 # NN_DONTWAIT=1
    buf = Vector{Ptr{Cchar}}(undef,1)
    bytes = ccall(
        (:nn_recv, "libnanomsg"),
        Int32,
        (Cint, Ptr{Cchar}, Csize_t, Cint),
        s, pointer(buf), typemax(Csize_t), flags)
    msg = unsafe_string(buf[1])
    ccall((:nn_freemsg, "libnanomsg"), Cint, (Ptr{Cchar},), buf[1])
    return msg
end

function nn_send(sock, msg; blocking=true)
    sz_msg = length(msg)
    flags = blocking ? 0 : 1 # NN_DONTWAIT=1
    bytes = ccall((:nn_send, "libnanomsg"), Int32, (Int32, Ptr{UInt8}, Csize_t, Int32), sock, msg, sz_msg, flags)
    return bytes
end

nn_shutdown(s, endpoint_id) = ccall((:nn_shutdown, "libnanomsg"), Int32, (Int32, Int32), s, endpoint_id)

# Pipeline example

const NN_PROTO_PIPELINE = 5
const NN_PUSH = NN_PROTO_PIPELINE * 16 + 0
const NN_PULL = NN_PROTO_PIPELINE * 16 + 1

function node0(url)
    sock = nn_socket(NN_PULL)
    check_err("nn_socket", sock)
    endpoint = nn_bind(s, url)
    check_err("nn_bind", endpoint)
    while true
        msg = nn_recv(s)
        check_err("nn_recv", msg)
        println("NODE0: RECEIVED ", msg)
    end
end

function node1(url, msg)
    sock = nn_socket(NN_PUSH)
    check_err("nn_socket", sock)
    rv = nn_connect(sock, url)
    check_err("nn_connect", rv)
    bytes = nn_send(sock, msg)
    check_err("nn_send", bytes)
    sleep(1) #  wait for messages to flush before shutting down
    rv = nn_shutdown(s, 0)
    check_err("nn_shutdown", rv)
end

node0("ipc:///tmp/pipeline.ipc")
node1("ipc:///tmp/pipeline.ipc", "hello, world its me")

ext: jl


text:
# This example mirrors rest.jl, but via the CRUD API.

module CRUD

using Discord

# Set this environment variable or replace with your own guild ID.
const GUILD = parse(Discord.Snowflake, get(ENV, "DISCORD_GUILD_ID", "1234567890"))

function main()
    c = Client(ENV["DISCORD_TOKEN"])

    # Send a request.
    future = retrieve(c, Guild, GUILD)
    # Then await it.
    resp = fetch(future)
    # Check for success.
    if !resp.ok
        # You can access the HTTP response if it went through.
        if resp.http_response !== nothing
            println("HTTP status: ", resp.http_response.status)
        end
        # You can also access caught exceptions.
        if resp.exception !== nothing
            println("Request error: ", sprint(showerror, resp.exception))
        end
        error("Retrieving guild failed")
    end

    # Get the actual Guild.
    guild = resp.val
    # If we're really confident, we can skip the error checks.
    guild = fetchval(future)

    println("Guild name: $(guild.name)")

    # Create a channel.
    channel = fetchval(create(c, DiscordChannel, guild; name="foo"))

    # Send a message to the channel.
    message = fetchval(create(c, Message, channel; content="Hello, world!"))
    # React to it.
    create(c, Reaction, message, '👍')
    # Edit it
    update(c, message; content="Goodbye, world!")
    # Delete it.
    delete(c, message)

    # Delete the channel.
    delete(c, channel)

    # Specific to CRUD: We can omit fetch/fetchval with macros.
    @fetch retrieve begin  # Here we're only wrapping calls to retrieve in fetch.
        resp = retrieve(c, Guild, GUILD)  # We don't have to call fetch.
        println("Response success: $(resp.ok)")
        future = create(c, DiscordChannel, resp.val; name="foo")  # Still returns a Future.
    end
    @fetchval retrieve create begin  # Wrapping calls to both retrieve and create.
        guild = retrieve(c, Guild, GUILD)  # We get the guild directly.
        println("Guild name: $(guild.name)")
        channel = create(c, DiscordChannel, guild; name="foo")  # We get the channel.
        future = update(c, channel; name="bar")  # Still returns a Future.
    end
    # If we don't specify any functions to wrap, they all get wrapped.
    @fetchval begin
        guild = retrieve(c, Guild, GUILD)
        channel = create(c, DiscordChannel, guild; name="foo")
        channel = update(c, channel; name="bar")
        delete(c, channel)
    end
end

end

if abspath(PROGRAM_FILE) == @__FILE__
    CRUD.main()
end


text:
struct BoundingBox{N,T}
    intervals::NTuple{N,ClosedInterval{T}}
end

BoundingBox(left::NTuple{N,T}, right::NTuple{N,T}) where {N,T<:Real} =
    BoundingBox(map((l, r) -> l..r, left, right))

BoundingBox(left::CartesianIndex, right::CartesianIndex) =
    BoundingBox(left.I, right.I)

Base.:(|)(bb::BoundingBox, I::CartesianIndex) =
    BoundingBox(map((iv, i) -> min(minimum(iv), i)..max(maximum(iv), i),
                    bb.intervals, I.I))

Statistics.middle(bb::BoundingBox) = map(iv->(minimum(iv)+maximum(iv))/2, bb.intervals)

buffer(bb::BoundingBox, l) = BoundingBox(map(iv -> minimum(iv)-l .. maximum(iv)+l, bb.intervals))


text:
export PlotCamvidLossUnconstrained, PlotCamvidLossConstrained, PlotDataLabelPredictionCamvid, PlotHyperspectralLossConstrained, PlotDataLabelPredictionHyperspectral

function PlotHyperspectralLossConstrained(logs,eval_every)

  iter_ax = range(0,step=eval_every,length=length(logs.dc2_train))

  figure(figsize=(5,4));
  semilogy(iter_ax,logs.dc2_train);title("Squared distance to set");xlabel("Iteration")
  tight_layout()
  savefig("dc2_hyperspectral.png")
end

function PlotDataLabelPredictionHyperspectral(plt_ind::Int,data,label,HN,active_channels,active_z_slice,pos_inds_select,neg_inds_select,tag::String)

  #predict
  p1, prediction, ~ = HN.forward(data[plt_ind], data[plt_ind])
  prediction = prediction |> cpu
  prediction[:,:,active_z_slice,active_channels,1].=softmax(prediction[:,:,active_z_slice,active_channels,1],dims=3);

  prediction = prediction |> cpu

  close("all")
  vmi = 0.0
  vma = 1.0
  figure(figsize=(6,4))
  subplot(2,2,1);imshow(Array(prediction)[3:end-2,3:end-2,active_z_slice,1,1],vmin=vmi,vmax=vma);PyPlot.title("Prediction - class 1");colorbar()
  subplot(2,2,2);imshow(Array(prediction)[3:end-2,3:end-2,active_z_slice,2,1],vmin=vmi,vmax=vma);PyPlot.title("Prediction - class 2");colorbar()
  subplot(2,2,3);imshow(label[plt_ind][3:end-2,3:end-2,active_z_slice,1,1],vmin=vmi,vmax=vma);PyPlot.title("Label - class 1");colorbar()
  subplot(2,2,4);imshow(label[plt_ind][3:end-2,3:end-2,active_z_slice,2,1],vmin=vmi,vmax=vma);PyPlot.title("Label - class 2");colorbar()
  savefig(string(tag,"_prediction_channels.png"))

  #thresholded prediction
  pred_thres = zeros(Int,size(prediction)[1:2])
  pos_inds   = findall(prediction[:,:,active_z_slice,1,1] .> prediction[:,:,active_z_slice,2,1])
  pred_thres[pos_inds] .= 1
  pred_thres2 = zeros(Int,size(prediction)[1:2])
  pos_inds    = findall(prediction[:,:,active_z_slice,1,1] .< prediction[:,:,active_z_slice,2,1])
  pred_thres2[pos_inds] .= 1
  figure(figsize=(5,4));
  imshow(Array(pred_thres)[3:end-2,3:end-2],vmin=vmi,vmax=vma);PyPlot.title(string("Prediction"));#xlabel("x");ylabel("y")
  tight_layout()
  savefig(string(tag,"_prediction.png"))#,bbox="tight")

  println(length(findall(pred_thres.>0))/prod(size(pred_thres)[1:2]))

  figure(figsize=(5,4));
  imshow(label[plt_ind][:,:,active_z_slice,1,1],vmin=vmi,vmax=vma);PyPlot.title(string("Labels"));#xlabel("x");ylabel("y")
  for i=1:length(pos_inds_select)
      PyPlot.scatter(pos_inds_select[i][2],pos_inds_select[i][1],c="r")
  end
  if isempty(neg_inds_select)==false
    for i=1:length(neg_inds_select)
        PyPlot.scatter(neg_inds_select[i][2],neg_inds_select[i][1],c="w")
    end
  end
  tight_layout()
  savefig(string(tag,"_labels.png"))#,bbox="tight")

  figure(figsize=(5,4));
  imshow(Array(pred_thres)[3:end-2,3:end-2]-label[plt_ind][:,:,active_z_slice,1,1][3:end-2,3:end-2],vmin=-1,vmax=1);PyPlot.title(string("Difference"));#xlabel("x");ylabel("y")
  tight_layout()
  savefig(string(tag,"_error.png"))#,bbox="tight")

  figure(figsize=(5,4));
  imshow(Array(pred_thres2)[3:end-2,3:end-2]-label[plt_ind][:,:,active_z_slice,1,1][3:end-2,3:end-2],vmin=-1,vmax=1);PyPlot.title(string("Difference"));#xlabel("x");ylabel("y")
  tight_layout()
  savefig(string(tag,"_error2.png"))#,bbox="tight")

  #Plot data+prediction
  figure(figsize=(9,5))
  subplot(1,2,1);
  imshow(data[plt_ind][3:end-2,3:end-2,40,1,1],cmap="Greys")
  imshow(Array(pred_thres)[3:end-2,3:end-2],alpha=0.30)
  title(string("Data T1 + Prediction - ",tag))
  subplot(1,2,2);
  imshow(data[plt_ind][3:end-2,3:end-2,40,2,1],cmap="Greys")
  imshow(Array(pred_thres)[3:end-2,3:end-2],alpha=0.30)
  title(string("Data T2 + Prediction - ",tag))
  tight_layout()
  savefig(string(tag,"data_plus_pred.png"))

  figure();
  imshow(data[plt_ind][:,:,40,1,1]);title("Slice of T1 data");savefig("data_surface_T1.png")
  imshow(data[plt_ind][:,:,40,2,1]);title("Slice of T2 data");savefig("data_surface_T2.png")

  #plot data cube
  #-----------------------
  close("all")
  index_cube = CartesianIndices(size(data[plt_ind])[1:3]);
  slice1 = vec(index_cube[:,:,end])
  slice1x = zeros(length(slice1)); slice1y = zeros(length(slice1));slice1z = zeros(length(slice1))
  slice1_values1 = zeros(length(slice1));
  slice1_values2 = zeros(length(slice1));
  for i=1:length(slice1)
    slice1x[i] = slice1[i][1]
    slice1y[i] = slice1[i][2]
    slice1z[i] = slice1[i][3]
    slice1_values1[i] = data[plt_ind][slice1[i][1],slice1[i][2],slice1[i][3],1,1]
    slice1_values2[i] = data[plt_ind][slice1[i][1],slice1[i][2],slice1[i][3],2,1]
  end
  figure()
  scatter3D(slice1x,slice1y,slice1z,c=slice1_values1,vmin=-0.2,vmax=0.5);title("Time 1")
  #scatter3D(slice1x,slice1y,slice1z,c=slice1_values2,vmin=-0.2,vmax=0.5);title("Time 2")

  slice1 = vec(index_cube[:,1,:])
  slice1x = zeros(length(slice1)); slice1y = zeros(length(slice1));slice1z = zeros(length(slice1))
  slice1_values1 = zeros(length(slice1));
  slice1_values2 = zeros(length(slice1));
  for i=1:length(slice1)
    slice1x[i] = slice1[i][1]
    slice1y[i] = slice1[i][2]
    slice1z[i] = slice1[i][3]
    slice1_values1[i] = data[plt_ind][slice1[i][1],slice1[i][2],slice1[i][3],1,1]
    slice1_values2[i] = data[plt_ind][slice1[i][1],slice1[i][2],slice1[i][3],2,1]
  end
  scatter3D(slice1x,slice1y,slice1z,c=slice1_values1,vmin=-0.2,vmax=0.5)
  #subplot(1,2,2);scatter3D(slice1x,slice1y,slice1z,c=slice1_values2,vmin=-0.2,vmax=0.5);

  slice1 = vec(index_cube[end,:,:])
  slice1x = zeros(length(slice1)); slice1y = zeros(length(slice1));slice1z = zeros(length(slice1))
  slice1_values1 = zeros(length(slice1));
  slice1_values2 = zeros(length(slice1));
  for i=1:length(slice1)
    slice1x[i] = slice1[i][1]
    slice1y[i] = slice1[i][2]
    slice1z[i] = slice1[i][3]
    slice1_values1[i] = data[plt_ind][slice1[i][1],slice1[i][2],slice1[i][3],1,1]
    slice1_values2[i] = data[plt_ind][slice1[i][1],slice1[i][2],slice1[i][3],2,1]
  end
  scatter3D(slice1x,slice1y,slice1z,c=slice1_values1,vmin=-0.2,vmax=0.5)
  #subplot(1,2,2);scatter3D(slice1x,slice1y,slice1z,c=slice1_values2,vmin=-0.2,vmax=0.5)
  tight_layout()
  savefig("data3D_T1.png",bbox="tight")
  #-----------------------

  #-----------------------
  close("all")
  index_cube = CartesianIndices(size(data[plt_ind])[1:3]);
  slice1 = vec(index_cube[:,:,end])
  slice1x = zeros(length(slice1)); slice1y = zeros(length(slice1));slice1z = zeros(length(slice1))
  slice1_values1 = zeros(length(slice1));
  slice1_values2 = zeros(length(slice1));
  for i=1:length(slice1)
    slice1x[i] = slice1[i][1]
    slice1y[i] = slice1[i][2]
    slice1z[i] = slice1[i][3]
    slice1_values1[i] = data[plt_ind][slice1[i][1],slice1[i][2],slice1[i][3],1,1]
    slice1_values2[i] = data[plt_ind][slice1[i][1],slice1[i][2],slice1[i][3],2,1]
  end
  figure()
  scatter3D(slice1x,slice1y,slice1z,c=slice1_values2,vmin=-0.2,vmax=0.5);title("Time 2")
  #scatter3D(slice1x,slice1y,slice1z,c=slice1_values2,vmin=-0.2,vmax=0.5);title("Time 2")

  slice1 = vec(index_cube[:,1,:])
  slice1x = zeros(length(slice1)); slice1y = zeros(length(slice1));slice1z = zeros(length(slice1))
  slice1_values1 = zeros(length(slice1));
  slice1_values2 = zeros(length(slice1));
  for i=1:length(slice1)
    slice1x[i] = slice1[i][1]
    slice1y[i] = slice1[i][2]
    slice1z[i] = slice1[i][3]
    slice1_values1[i] = data[plt_ind][slice1[i][1],slice1[i][2],slice1[i][3],1,1]
    slice1_values2[i] = data[plt_ind][slice1[i][1],slice1[i][2],slice1[i][3],2,1]
  end
  scatter3D(slice1x,slice1y,slice1z,c=slice1_values2,vmin=-0.2,vmax=0.5)
  #subplot(1,2,2);scatter3D(slice1x,slice1y,slice1z,c=slice1_values2,vmin=-0.2,vmax=0.5);

  slice1 = vec(index_cube[end,:,:])
  slice1x = zeros(length(slice1)); slice1y = zeros(length(slice1));slice1z = zeros(length(slice1))
  slice1_values1 = zeros(length(slice1));
  slice1_values2 = zeros(length(slice1));
  for i=1:length(slice1)
    slice1x[i] = slice1[i][1]
    slice1y[i] = slice1[i][2]
    slice1z[i] = slice1[i][3]
    slice1_values1[i] = data[plt_ind][slice1[i][1],slice1[i][2],slice1[i][3],1,1]
    slice1_values2[i] = data[plt_ind][slice1[i][1],slice1[i][2],slice1[i][3],2,1]
  end
  scatter3D(slice1x,slice1y,slice1z,c=slice1_values2,vmin=-0.2,vmax=0.5)
  #subplot(1,2,2);scatter3D(slice1x,slice1y,slice1z,c=slice1_values2,vmin=-0.2,vmax=0.5)
  tight_layout()
  savefig("data3D_T2.png")
  #-----------------------

end

function PlotCamvidLossUnconstrained(logs,eval_every)

  iter_ax = range(0,step=eval_every,length=length(logs.train))

  figure(figsize=(13,4));
  subplot(1,2,1);plot(iter_ax,logs.train,label="train");title("Training loss labels");
  subplot(1,2,1);plot(iter_ax,logs.val,"--",label="validation");legend();title("Training loss labels");
  xlabel("Iteration number");
  subplot(1,2,2);plot(iter_ax,logs.IoU_train[:,1],label="class1 - train");title("IoU");
  subplot(1,2,2);plot(iter_ax,logs.IoU_train[:,2],label="class2 - train");legend();title("IoU");
  subplot(1,2,2);plot(iter_ax,logs.IoU_val[:,1],"--",label="class1 - validation");title("IoU");
  subplot(1,2,2);plot(iter_ax,logs.IoU_val[:,2],"--",label="class2 - validation");legend();title("IoU");
  xlabel("Iteration number");
  tight_layout()
  savefig("loss_unconstrained.png")

end

function PlotCamvidLossConstrained(logs,eval_every)

  iter_ax = range(0,step=eval_every,length=length(logs.dc2_train))
  figure(figsize=(13,4));
  subplot(1,3,1);semilogy(iter_ax,logs.dc2_train,label="train");legend();title("Squared distance-to-set");
  xlabel("Iteration number");
  subplot(1,3,2);plot(iter_ax,logs.train,label="train");title("Training loss labels");
  subplot(1,3,2);plot(iter_ax,logs.val,"--",label="validation");legend();title("Training loss labels");
  xlabel("Iteration number");
  subplot(1,3,3);plot(iter_ax,logs.IoU_train[:,1],label="class1 - train");title("IoU");
  subplot(1,3,3);plot(iter_ax,logs.IoU_train[:,2],label="class2 - train");legend();title("IoU");
  subplot(1,3,3);plot(iter_ax,logs.IoU_val[:,1],"--",label="class1 - validation");title("IoU");
  subplot(1,3,3);plot(iter_ax,logs.IoU_val[:,2],"--",label="class2 - validation");legend();title("IoU");
  xlabel("Iteration number");
  tight_layout()
  savefig("loss_constrained.png")

end

function PlotDataLabelPredictionCamvid(plt_ind::Int,data,label,HN,active_channels,tag::String)

  #Plot data
  figure(figsize=(5,4))
  imshow(data[plt_ind][:,:,1:3,1])
  title(string("Data - ",tag))
  tight_layout()
  savefig(string(tag,"data.png"),bbox="tight")

  #predict
  p1, prediction, ~ = HN.forward(data[plt_ind], data[plt_ind])
  prediction[:,:,active_channels,1].=softmax(prediction[:,:,active_channels,1],dims=3);

  close("all")
  vmi = 0.0
  vma = 1.0
  figure(figsize=(6,4))
  subplot(2,2,1);imshow(Array(prediction)[3:end-2,3:end-2,1,1],vmin=vmi,vmax=vma);PyPlot.title("Prediction - class 1");colorbar()
  subplot(2,2,2);imshow(Array(prediction)[3:end-2,3:end-2,2,1],vmin=vmi,vmax=vma);PyPlot.title("Prediction - class 2");colorbar()
  subplot(2,2,3);imshow(label[plt_ind][3:end-2,3:end-2,1,1],vmin=vmi,vmax=vma);PyPlot.title("Label - class 1");colorbar()
  subplot(2,2,4);imshow(label[plt_ind][3:end-2,3:end-2,2,1],vmin=vmi,vmax=vma);PyPlot.title("Label - class 2");colorbar()
  tight_layout()
  savefig(string(tag,"_prediction_channels.png"))

  #thresholded prediction
  pred_thres = zeros(Int,size(prediction)[1:2])
  pos_inds = findall(prediction[:,:,1,1] .> 0.65)#prediction[:,:,2,1])
  pred_thres[pos_inds] .= 1
  figure(figsize=(5,4));
  imshow(Array(pred_thres)[3:end-2,3:end-2],vmin=vmi,vmax=vma);PyPlot.title(string("Prediction - ",tag));#xlabel("x");ylabel("y")
  tight_layout()
  savefig(string(tag,"_prediction.png"))#,bbox="tight")

  #Plot data+prediction
  figure(figsize=(5,4))
  imshow(data[plt_ind][3:end-2,3:end-2,1:3,1])
  imshow(Array(pred_thres)[3:end-2,3:end-2],alpha=0.45)
  title(string("Image+Prediction - ",tag))
  tight_layout()
  savefig(string(tag,"data_plus_pred.png"))
  return
end


text:
# ======================================================
# AbstractValGraph structure
# ======================================================

"""
    AbstractValGraph{V, V_VALS, E_VALS, G_VALS} <: AbstractGraph{V}

Abstract value graph with vertex type `V`, vertex values of types `V_VALS`,
edge values of types `E_VALS` and graph values of types `G_VALS`.

### See also
[`AbstractGraph`](@ref), [`AbstractEdgeValGraph`](@ref)
"""
abstract type AbstractValGraph{V<:Integer, V_VALS, E_VALS, G_VALS} <: AbstractGraph{V} end

#  ------------------------------------------------------
#  specialised types
#  ------------------------------------------------------

OneEdgeValGraph{V, V_VALS, E_VAL} = AbstractValGraph{V, V_VALS, <: AbstractNTuple{1, E_VAL}}

ZeroEdgeValGraph{V, V_VALS} = AbstractValGraph{V, V_VALS, <: AbstractNTuple{0}}

OneVertexValGraph{V, V_VAL, E_VALS} = AbstractValGraph{V, <: AbstractNTuple{1, V_VAL}, E_VALS}

ZeroVertexValGraph{V, E_VALS} = AbstractValGraph{V, <: AbstractNTuple{0}, E_VALS}

OneGraphValGraph{V, V_VALS, E_VALS, G_VAL} = AbstractValGraph{V, V_VALS, E_VALS, <: AbstractNTuple{1, G_VAL}}

ZeroGraphValGraph{V, V_VALS, E_VALS} = AbstractValGraph{V, V_VALS, E_VALS, <: AbstractNTuple{0}}

# ======================================================
# Type parameter information
# ======================================================

#  ------------------------------------------------------
#  eltype
#  ------------------------------------------------------

Graphs.eltype(::Type{<:AbstractValGraph{V}}) where {V} = V

# This should not be necessary, as Base implements `eltype(x) = eltype(typeof(x))`
# but unfortunately Graphs.jl redefines `eltype(::AbstractGraph)` as not defined
Graphs.eltype(g::AbstractValGraph) = eltype(typeof(g))

#  ------------------------------------------------------
#  vertexvals_type
#  ------------------------------------------------------

# TODO might also implement this for SimpleGraph
"""
    vertexvals_type(g::AbstractValGraph)
    vertexvals_type(::Type{<:AbstractValGraph})

Return the types of the vertex values of a graph `g`.
"""
vertexvals_type(g::AbstractValGraph) = vertexvals_type(typeof(g))

vertexvals_type(::Type{<:AbstractValGraph{V, V_VALS}}) where {V, V_VALS} = V_VALS

#  ------------------------------------------------------
#  edgevals_type
#  ------------------------------------------------------

# TODO might also implement this for SimpleGraph
"""
    edgevals_type(g::AbstractValGraph)
    edgevals_type(::Type{<:AbstractValGraph})

Return the types of the edge values of a graph `g`.
"""
edgevals_type(g::AbstractValGraph) = edgevals_type(typeof(g))

edgevals_type(::Type{<:AbstractValGraph{V, V_VALS, E_VALS}}) where {V, V_VALS, E_VALS} = E_VALS

# TODO test
"""
    edgevals_type(g::AbstractValGraph, key)
    edgevals_type(::Type{<:AbstractValGraph}, key)

Return the type of the edge values for a specific `key` of a graph `g`.
"""
edgevals_type(g::AbstractValGraph, key) = edgevals_type(typeof(g), key)

function edgevals_type(G::Type{<:AbstractValGraph}, key::Integer)

    return fieldtype(edgevals_type(G), Int(key))
end

function edgevals_type(G::Type{<:AbstractValGraph}, key::Symbol)

    return fieldtype(edgevals_type(G), key)
end

#  ------------------------------------------------------
#  graphvals_type
#  ------------------------------------------------------

graphvals_type(g::AbstractValGraph) = graphvals_type(typeof(g))

graphvals_type(::Type{<:AbstractValGraph{V, V_VALS, E_VALS, G_VALS}}) where {V, V_VALS, E_VALS, G_VALS} = G_VALS

#  ------------------------------------------------------
#  hasedgekey & hasedgekey_or_throw
#  ------------------------------------------------------

"""
    hasedgekey(g::AbstractValGraph, key)
    hasedgekey(::Type{<:AbstractValGraph}, key)

Return true if `key` is an edge value key for this graph.
"""
hasedgekey(g::AbstractValGraph, key) = hasedgekey(typeof(g), key)

function hasedgekey(G::Type{<:AbstractValGraph}, key::Symbol)

    return hasfield(edgevals_type(G), key)
end

function hasedgekey(G::Type{<:AbstractValGraph}, key::Integer)

    E_VALS = edgevals_type(G)
    return key in OneTo(length(E_VALS.types))
end

function hasedgekey_or_throw(G::Type{<:AbstractValGraph}, key)
    hasedgekey(G, key) && return nothing

    error("$key is not a valid edge key for this graph.")
end

hasedgekey_or_throw(g::AbstractValGraph, key) = hasedgekey_or_throw(typeof(g), key)

#  ------------------------------------------------------
#  hasvertexkey & hasvertexkey_or_throw
#  ------------------------------------------------------

"""
    hasvertexkey(g::AbstractValGraph, key)
    hasvertexkey(::Type{<:AbstractValGraph}, key)

Return true if `key` is a vertex value key for this graph.
"""
function hasvertexkey(G::Type{<:AbstractValGraph}, key::Symbol)

    return hasfield(vertexvals_type(G), key)
end

function hasvertexkey(G::Type{<:AbstractValGraph}, key::Integer)

    V_VALS = vertexvals_type(G)
    return key in OneTo(length(V_VALS.types))
end

hasvertexkey(g::AbstractValGraph, key) = hasvertexkey(typeof(g), key)

function hasvertexkey_or_throw(G::Type{<:AbstractValGraph}, key)
    hasvertexkey(G, key) && return nothing

    error("$key is not a valid vertex key for this graph.")
end

hasvertexkey_or_throw(g::AbstractValGraph, key) = hasvertexkey_or_throw(typeof(g), key)

#  ------------------------------------------------------
#  hasgraphkey & hasgraphkey_or_throw
#  ------------------------------------------------------

function hasgraphkey(G::Type{<:AbstractValGraph}, key::Symbol)

    return hasfield(graphvals_type(G), key)
end

function hasgraphkey(G::Type{<:AbstractValGraph}, key::Integer)

    G_VALS = graphvals_type(G)
    return key ∈ OneTo(length(G_VALS.types))
end

hasgraphkey(g::AbstractValGraph, key) = hasgraphkey(typeof(g), key)

function hasgraphkey_or_throw(G::Type{<:AbstractValGraph}, key)
    hasgraphkey(G, key) && return nothing

    error("$key is not a valid graph key for this graph.")
end

hasgraphkey_or_throw(g::AbstractValGraph, key) = hasgraphkey_or_throw(typeof(g), key)

# ======================================================
#  show
# ======================================================

function Base.show(io::IO, ::MIME"text/plain", g::AbstractValGraph)
    print(io, "{$(nv(g)), $(ne(g))}")
    print(io, " ", is_directed(g) ? "directed" : "undirected")
    print(io, " $(typeof(g).name.name) with")
    println(io)
    print(io, "              eltype: $(eltype(g))")
    println(io)
    print(io, "  vertex value types: $(typetuple(vertexvals_type(g)))")
    println(io)
    print(io, "    edge value types: $(typetuple(edgevals_type(g)))")
    println(io)
    print(io, "   graph value types: $(typetuple(graphvals_type(g)))")
end

# ======================================================
# Partial default implementation of Graphs.jl interface
# ======================================================

Graphs.vertices(g::AbstractValGraph) = OneTo{eltype(g)}(nv(g))

Graphs.has_vertex(g::AbstractValGraph, v) = v ∈ vertices(g)

"""
    edges(g::AbstractValGraph[, key])

Return the edges of `g`. By default add no edge values
but when `key=:` then add edge values.
"""
Graphs.edges(g::AbstractValGraph, key=nothing) = ValEdgeIter(g, key)

Graphs.edgetype(g::AbstractValGraph) = eltype(edges(g))

Graphs.ne(g::AbstractValGraph) = length(edges(g))

# TODO a Base.Generator would be better here, but it causes problems with sort. Maybe
# add a custom iterator
Graphs.outneighbors(g::AbstractValGraph, u) = [v for v ∈ vertices(g) if has_edge(g, u, v)]

Graphs.inneighbors(g::AbstractValGraph, v) = is_directed(g) ? [u for u ∈ vertices(g) if has_edge(g, u, v)] : outneighbors(g, v)

# ======================================================
# Modifiers
# ======================================================

#  -----------------------------------------------------
#  add_edge!
#  -----------------------------------------------------

"""
    add_edge!(g::AbstractValGraph, s, d, [values])
    add_edge!(g::AbstractValGraph, s, d; val=value)

Add an edge `e = (s, d, values)` to a graph `g` and set the edge values.

If `g` does not have and edge values, `values` can be omitted.

For graphs with a single edge value, that value can also be specified with
the `val` keyword argument.

Return `true` if the edge was added successfully, otherwise return `false`.
If the edge already exists, return `false` but still change the edge values.
"""
function add_edge! end

Graphs.add_edge!(g::ZeroEdgeValGraph, s, d) = add_edge!(g, s, d, ())

Graphs.add_edge!(g::OneEdgeValGraph, s, d; val) = add_edge!(g, s, d, tuple(val))

#  -----------------------------------------------------
#  add_vertex!
#  -----------------------------------------------------

"""
    add_vertex!(g::AbstractValGraph[, values])
    add_vertex!(g::AbstractValGraph; val=value)

Add an vertex to a graph `g` and set the vertex values.

If `g` does not have vertex values, `values` can be omitted.

For graphs with a single vertex value, that value can also be specified with
the `val` keyword argument.

Return `true` if the vertex was added successfully, otherwise return `false`.
"""
function add_vertex! end

Graphs.add_vertex!(g::ZeroVertexValGraph) = add_vertex!(g, ())

Graphs.add_vertex!(g::OneVertexValGraph; val) = add_vertex!(g, tuple(val))

# ======================================================
# Accessors
# ======================================================

#  -----------------------------------------------------
#  get_vertexval
#  -----------------------------------------------------

"""
    get_vertexval(g::AbstractValGraph, v, key)

Return the vertex value for vertex `v` and key `key`. If `g` has a single
vertex value, `key` can be omitted.
"""
function get_vertexval end

get_vertexval(g::AbstractValGraph, v, key::Symbol) =
    get_vertexval(g, v, Base.fieldindex(vertexvals_type(g), key))

get_vertexval(g::OneVertexValGraph, v) = get_vertexval(g, v, 1)

"""
    get_vertexval(g::AbstractValGraph, v, :)

Return all vertex values for vertex `v`.
"""
function get_vertexval(g::AbstractValGraph, v, ::Colon)

    V_VALS = vertexvals_type(g)
    return V_VALS(get_vertexval(g, v, i) for i in OneTo(length(V_VALS.types)))
end

#  -----------------------------------------------------
#  get_edgeval
#  -----------------------------------------------------

# TODO unfortunately I have not found a way yet to make the generic
# get_edgeval(g, s, d, ::Symbol) function typesafe yet so that concrete
# graph types might need to implement this function for extra performance

"""
    get_edgeval(g::AbstractValGraph, s, d, key)

Return the value associated with the edge `s -> d` for the key `key` in `g`.

Throw an exception if the graph does not contain such an edge or if the key is not a valid edge key.

For graphs that only have one value per edge, `key` can be omitted.

### See also
[`get_edgeval_or`](@ref), [`set_edgeval!`](@ref)

### Examples
```jldoctest
julia> gv = ValDiGraph(path_digraph(3), edgeval_types=(a=Float64,), edgeval_init=(s, d) -> (rand(MersenneTwister(0)),))
{3, 2} directed ValDiGraph with
              eltype: Int64
  vertex value types: ()
    edge value types: (a = Float64,)
   graph value types: ()

julia> get_edgeval(gv, 1, 2, :a)
0.8236475079774124

julia> get_edgeval(gv, 1, 2, 1)
0.8236475079774124

julia> get_edgeval(gv, 1, 2)
0.8236475079774124

julia> get_edgeval(gv, 1, 3)
ERROR: No such edge

julia> get_edgeval(gv, 1, 2, :b)
ERROR: b is not a valid edge key for this graph.
```
"""
get_edgeval(g::AbstractValGraph, s, d, key::Symbol) =
    get_edgeval(g, s, d, Base.fieldindex(edgevals_type(g), key))

get_edgeval(g::OneEdgeValGraph, s, d) = get_edgeval(g, s, d, 1)

"""
    get_edgeval(g::AbstractValGraph, s, d, :)

Return all values associated with the edge `s -> d` in `g`.

Throw an exception if the graph does not contain such an edge.

### See also
[`get_edgeval_or`](@ref), [`set_edgeval!`](@ref)

### Examples
```jldoctest
julia> gv = ValDiGraph(path_digraph(3), edgeval_types=(a=Float64, b=Int), edgeval_init=(s, d) -> (rand(MersenneTwister(0)), 10))
{3, 2} directed ValDiGraph with
              eltype: Int64
  vertex value types: ()
    edge value types: (a = Float64, b = Int64)
   graph value types: ()

julia> get_edgeval(gv, 1, 2, :)
(a = 0.8236475079774124, b = 10)

julia> get_edgeval(gv, 1, 3, :)
ERROR: Values not found
```
"""
function get_edgeval(g::AbstractValGraph, s, d, ::Colon)

    E_VALS = edgevals_type(g)
    return E_VALS(get_edgeval(g, s, d, i) for i in OneTo(length(E_VALS.types)))
end

#  -----------------------------------------------------
#  get_edgeval_or
#  -----------------------------------------------------

"""
    get_edgeval_or(g::ValGraph, s, d, alternative)
    get_edgeval_or(g::ValDiGraph, s, d, alternative)
    get_edgeval_or(g::ValOutDiGraph, s, d, alternative)

Return all values associated with the edge `s -> d` in `g`.

If there is no such edge return `alternative`.

### See also
[`get_edgeval`](@ref), [`set_edgeval!`](@ref)

### Examples
```jldoctest
julia> gv = ValDiGraph(path_digraph(3), edgeval_types=(a=Float64, b=Int), edgeval_init=(s, d) -> (rand(MersenneTwister(0)), 10))
{3, 2} directed ValDiGraph with
              eltype: Int64
  vertex value types: ()
    edge value types: (a = Float64, b = Int64)
   graph value types: ()

julia> get_edgeval_or(gv, 1, 2, :, missing)
(a = 0.8236475079774124, b = 10)

julia> get_edgeval_or(gv, 1, 3, :, missing)
missing
```
"""
get_edgeval_or(g::AbstractValGraph, s, d, key, alternative) =
    has_edge(g, s, d) ? get_edgeval(g, s, d, key) : alternative

get_edgeval_or(g::OneEdgeValGraph, s, d, alternative) = get_edgeval_or(g, s, d, 1, alternative)

#  -----------------------------------------------------
#  get_graphval
#  -----------------------------------------------------

"""
    get_graphval(g::AbstractValGraph, key)

Return the graph value for the key `key`. If `g` has a single
graph value, `key` can be omitted.
"""
function get_graphval end

get_graphval(g::AbstractValGraph, key::Symbol) =
    get_graphval(g, Base.fieldindex(graphvals_type(g), key))

get_graphval(g::OneGraphValGraph) = get_graphval(g, 1)

"""
    get_graphval(g::AbstractValGraph, :)

Return all graph values.
"""
function get_graphval(g::AbstractValGraph, ::Colon)

    G_VALS = graphvals_type(g)
    return G_VALS(get_graphval(g, i) for i in OneTo(length(G_VALS.types)))
end

#  -----------------------------------------------------
#  set_edgeval!
#  -----------------------------------------------------

set_edgeval!(g::OneEdgeValGraph, s, d, value) = set_edgeval!(g, s, d, 1, value)

set_edgeval!(g::AbstractValGraph, s, d, key::Symbol, value) =
    set_edgeval!(g, s, d, Base.fieldindex(edgevals_type(g), key), value)

"""
    set_edgeval!(g::AbstractValGraph, s, d, :, values)

Set the values of the edge `e: s -> d` to `values`. Return `true` if such an edge exists and
`false` otherwise.
"""
function set_edgeval!(g::AbstractValGraph, s, d, ::Colon, values)

    has_edge(g, s, d) || return false

    E_VALS = edgevals_type(g)

    # TODO currently we cannot convert from tuples to named tuples or vice versa
    values = convert_to_tuple(E_VALS, values)
    for (i, value) in enumerate(values)
        set_edgeval!(g, s, d, i, value)
    end

    return true
end

#  ------------------------------------------------------
#  set_vertexval!
#  ------------------------------------------------------

set_vertexval!(g::OneVertexValGraph, v, value) = set_vertexval!(g, v, 1, value)

set_vertexval!(g::AbstractValGraph, v, key::Symbol, value) =
    set_vertexval!(g, v, Base.fieldindex(vertexvals_type(g), key), value)

"""
    set_vertexval!(g::AbstractValGraph, v, :, values)

Set the values of the vertex `v` to `values`. Return `true` if such an vertex exists
and `false` otherwise.
"""
function set_vertexval!(g::AbstractValGraph, v, ::Colon, values)

    has_vertex(g, v) || return false

    V_VALS = vertexvals_type(g)

    # TODO currently we cannot convert from tuples to named tuples or vice versa
    values = convert_to_tuple(V_VALS, values)
    for (i, value) in enumerate(values)
        set_vertexval!(g, v, i, value)
    end

    return true
end

#  ------------------------------------------------------
#  set_graphval!
#  ------------------------------------------------------

"""
    set_graphval!(g::AbstractValGraph, key, value)

Set the graph value specified by `g` to `value`. If `g` has only a single graph
value, `key` can be omitted.

### See also
[`get_graphval`](@ref), [`set_vertexval!`](@ref), [`set_edgeval!`](@ref)
"""
function set_graphval! end

set_graphval!(g::OneGraphValGraph, value) = set_graphval!(g, 1, value)

set_graphval!(g::AbstractValGraph, key::Symbol, value) =
    set_graphval!(g, Base.fieldindex(graphvals_type(g), key), value)

"""
    set_graphval!(g::AbstractValGraph, :, values)

Set all graph values of `g` to `values`.
### See also
[`get_graphval`](@ref), [`set_vertexval!`](@ref), [`set_edgeval!`](@ref)
"""
function set_graphval!(g::AbstractValGraph, ::Colon, values)

    G_VALS = graphvals_type(g)

    # TODO currently we cannot convert from tuples to named tuples or vice versa
    values = convert_to_tuple(G_VALS, values)
    for (key, value) in enumerate(values)
        set_graphval!(g, key, value)
    end

    return true
end

#  ------------------------------------------------------
#  outedgevals
#  ------------------------------------------------------

"""
    outedgevals(g::AbstractValGraph, v [, key])

Return an iterator of edge values of outgoing edges from `v` to its neighbors.

If `g` has multiple edge values, the key cannot be omitted.
The order of the neighbors is the same as for `outneighbors(g, v)`.
"""
function outedgevals end

outedgevals(g::OneEdgeValGraph, v) = outedgevals(g, v, 1)

outedgevals(g::AbstractValGraph, v, key::Symbol) =
    outedgevals(g, v, Base.fieldindex(edgevals_type(g), key))

outedgevals(g::AbstractValGraph, u, key::Integer) =
    [get_edgeval(g, u, v, key) for v in outneighbors(g, u)]

"""
    outedgevals(g::AbstractValGraph, v, :)

Return an iterator of all edge values of outgoing edges from `v` to its neighbors.
"""
outedgevals(g::AbstractValGraph, u, ::Colon) =
    [get_edgeval(g, u, v, :) for v in outneighbors(g, u)]

#  ------------------------------------------------------
#  inedgevals
#  ------------------------------------------------------

"""
    inedgevals(g::AbstractValGraph, v [, key])

Return an iterator of edge values of ingoing edges from neighbors of `v`.

If `g` has multiple edge values, the key cannot be omitted.
The order of the neighbors is the same as for `inneighbors(g, v)`.
"""
function inedgevals end

inedgevals(g::OneEdgeValGraph, v) = inedgevals(g, v, 1)

inedgevals(g::AbstractValGraph, v, key::Symbol) =
    inedgevals(g, v, Base.fieldindex(edgevals_type(g), key))

inedgevals(g::AbstractValGraph, v, key::Integer) =
    [get_edgeval(g, u, v, key) for u in inneighbors(g, v)]

"""
    inedgevals(g::AbstractValGraph, v, :)

Return an iterator of all edge values of ingoing edges from neighbors of `v`.
"""
inedgevals(g::AbstractValGraph, v, ::Colon) =
    [get_edgeval(g, u, v, :) for u in inneighbors(g, v)]

# ======================================================
# Edge Iterator
# ======================================================

# TODO support also Int as key
struct ValEdgeIter{G<:AbstractValGraph, key} <: AbstractEdgeIter
    graph::G

    function ValEdgeIter{G}(g::G, key::Union{Colon, Nothing}) where {G}

        return new{G, key}(g)
    end
end

ValEdgeIter(g, key) = ValEdgeIter{typeof(g)}(g, key)

Base.length(iter::ValEdgeIter) = count(_ -> true, iter)

function Base.eltype(::Type{<:ValEdgeIter{G, key}}) where {G, key}

    E = is_directed(G) ? ValDiEdge : ValEdge
    V = eltype(G)
    E_VALS = edgevals_type(G)

    if key == nothing
        # TODO it might better to return an empty named tuple type in case
        # E_VALS is a named tuple, but then we need to adjust the iterators
        return E{V, Tuple{}}
    end

    return E{V, E_VALS}
end
Base.eltype(iter::ValEdgeIter) = eltype(typeof(iter))

function Base.iterate(iter::ValEdgeIter)

    g = iter.graph

    vertices_iter = vertices(g)
    vertices_next = iterate(vertices_iter)
    vertices_next == nothing && return nothing

    u, _ = vertices_next
    outneighbors_iter = outneighbors(g, u)
    outneighbors_next = iterate(outneighbors_iter)

    iterate(iter, (vertices_iter, vertices_next, outneighbors_iter, outneighbors_next))
end

# TODO this function might be faster if we could ensure
# that none of the fields of state can be nothing
# TODO might replace get_edgeval with outedgevals
function Base.iterate(iter::ValEdgeIter{G, key}, state) where {G, key}

    g = iter.graph

    vertices_iter, vertices_next, outneighbors_iter, outneighbors_next = state

    # at this point vertices_next != nothing
    u, vertices_state = vertices_next

    while outneighbors_next == nothing
        vertices_next = iterate(vertices_iter, vertices_state)
        vertices_next == nothing && return nothing
        u, vertices_state = vertices_next
        outneighbors_iter = outneighbors(g, u)
        outneighbors_next = iterate(outneighbors_iter)
        if !is_directed(g)
            while outneighbors_next != nothing && outneighbors_next[1] < u
                outneighbors_next = iterate(outneighbors_iter, outneighbors_next[2])
            end
        end
    end

    v, outneighbors_state = outneighbors_next
    edge = if key == nothing
                eltype(iter)(u, v, ())
            else
                eltype(iter)(u, v, get_edgeval(g, u, v, :))
            end

    outneighbors_next = iterate(outneighbors_iter, outneighbors_state)

    return edge, (vertices_iter, vertices_next, outneighbors_iter, outneighbors_next)

end


text:
function main()

    add::Int64 = 0
    st = time_ns();

    for index in 0:100000000
        add = add + index
    end

    println("==========================")
    println("julia 테스트")
    println(add)
    println(time_ns() - st)
end;

main()

max_stars_repo_name: bluetsys/language-benchmarking

max_issues_repo_name: bluetsys/language-benchmarking

max_forks_repo_name: bluetsys/language-benchmarking


text:
using Makie, GLMakie, GeometryBasics, Rotations

"""
    conformal_cubed_sphere_mapping(x, y)

Conformal mapping of a cube onto a sphere. Maps `(x, y)` on the north-pole face of a cube
to (X, Y, Z) coordinates in physical space. The face is oriented normal to Z-axis with
X and Y increasing with x and y.

The input coordinates must lie within the range -1 <= x <= 1,  -1 <= y <= 1.

This numerical conformal mapping is described by Rančić et al. (1996).

This is a Julia translation of MATLAB code from MITgcm [1] that is based on
Fortran 77 code from Jim Purser & Misha Rančić.

[1] http://wwwcvs.mitgcm.org/viewvc/MITgcm/MITgcm_contrib/high_res_cube/matlab-grid-generator/map_xy2xyz.m?view=markup
Author of conformal_cubed_sphere_mapping: ["ali.hh.ramadhan@gmail.com <ali.hh.ramadhan@gmail.com>"]
"""
A_Rancic = [
    +0.00000000000000,
    +1.47713062600964,
    -0.38183510510174,
    -0.05573058001191,
    -0.00895883606818,
    -0.00791315785221,
    -0.00486625437708,
    -0.00329251751279,
    -0.00235481488325,
    -0.00175870527475,
    -0.00135681133278,
    -0.00107459847699,
    -0.00086944475948,
    -0.00071607115121,
    -0.00059867100093,
    -0.00050699063239,
    -0.00043415191279,
    -0.00037541003286,
    -0.00032741060100,
    -0.00028773091482,
    -0.00025458777519,
    -0.00022664642371,
    -0.00020289261022,
    -0.00018254510830,
    -0.00016499474461,
    -0.00014976117168,
    -0.00013646173946,
    -0.00012478875823,
    -0.00011449267279,
    -0.00010536946150,
    -0.00009725109376
]

function conformal_cubed_sphere_mapping(x, y)
    X = xᶜ = abs(x)
    Y = yᶜ = abs(y)

    kxy = yᶜ > xᶜ

    xᶜ = 1 - xᶜ
    yᶜ = 1 - yᶜ

    kxy && (xᶜ = 1 - Y)
    kxy && (yᶜ = 1 - X)

    Z = ((xᶜ + im * yᶜ) / 2)^4
    W = W_Rancic(Z)

    im³ = im^(1/3)
    ra = √3 - 1
    cb = -1 + im
    cc = ra * cb / 2

    W = im³ * (W * im)^(1/3)
    W = (W - ra) / (cb + cc * W)
    X, Y = reim(W)

    H = 2 / (1 + X^2 + Y^2)
    X = X * H
    Y = Y * H
    Z = H - 1

    if kxy
        X, Y = Y, X
    end

    y < 0 && (Y = -Y)
    x < 0 && (X = -X)

    # Fix truncation for x = 0 or y = 0.
    x == 0 && (X = 0)
    y == 0 && (Y = 0)

    return X, Y, Z
end

W_Rancic(Z) = sum(A_Rancic[k] * Z^(k-1) for k in 1:length(A_Rancic))

function cubed_sphere_warp(
    a,
    b,
    c,
    R = max(abs(a), abs(b), abs(c)),
)

    # @show (a,b,c)
    fdim = argmax(abs.((a, b, c)))
    # @show fdim
    M = max(abs.((a, b, c))...)
    if fdim == 1 && a < 0
        # left face
        x1, x2, x3 = conformal_cubed_sphere_mapping(-b / M, c / M)
        x1, x2, x3 = RotX(π/2) * RotY(-π/2) * [x1, x2, x3]
    elseif fdim == 2 && b < 0
        # front face
        x1, x2, x3 = conformal_cubed_sphere_mapping(a / M, c / M)
        x1, x2, x3 = RotX(π/2) * [x1, x2, x3]
    elseif fdim == 1 && a > 0
        # right face
        x1, x2, x3 = conformal_cubed_sphere_mapping(b / M, c / M)
        x1, x2, x3 = RotX(π/2) * RotY(π/2) * [x1, x2, x3]
    elseif fdim == 2 && b > 0
        # back face
        x1, x2, x3 = conformal_cubed_sphere_mapping(a / M, -c / M)
        x1, x2, x3 = RotX(-π/2) * [x1, x2, x3]
    elseif fdim == 3 && c > 0
        # top face
        x1, x2, x3 = conformal_cubed_sphere_mapping(a / M, b / M)
    elseif fdim == 3 && c < 0
        # bottom face
        x1, x2, x3 = conformal_cubed_sphere_mapping(a / M, -b / M)
        x1, x2, x3 = RotX(π) * [x1, x2, x3]
    else
        error("invalid case for cubed_sphere_warp(::ConformalCubedSphere): $a, $b, $c")
    end

    return x1 * R, x2 * R, x3 * R

end

n = 6

r = range(-1,1,length=n+1)

xa = zeros(n+1,n+1)
xb = zeros(n+1,n+1)
xc = zeros(n+1,n+1)

# top face coordinates
xa .= r
xb .= r'
xc .= 1

a = [cubed_sphere_warp(a,b,c)[1] for (a,b,c) in zip(xa,xb,xc)]
b = [cubed_sphere_warp(a,b,c)[2] for (a,b,c) in zip(xa,xb,xc)]
c = [cubed_sphere_warp(a,b,c)[3] for (a,b,c) in zip(xa,xb,xc)]

sc = Scene(limits = HyperRectangle(Vec3f0(-2), Vec3f0(2)))

wireframe!(sc,
           a,b,c,
           show_axis=false,
           linewidth=1.2)
wireframe!(sc,
           a,b,-c,
           show_axis=false,
           linewidth=1.2)
wireframe!(sc,
           c,a,b,
           show_axis=false,
           linewidth=1.2)
wireframe!(sc,
           -c,a,b,
           show_axis=false,
           linewidth=1.2)
wireframe!(sc,
           b,c,a,
           show_axis=false,
           linewidth=1.2)
wireframe!(sc,
           b,-c,a,
           show_axis=false,
           linewidth=1.2)

sc


text:
using JuMP, Gurobi, StatsFuns

# problem parameters
NumAssets = [1000] # number of variables minus one
epsilon = logspace(-4.5,-0.3011,50)
epsilon = reverse(epsilon) # constraint violation tolerance

# iterate over the number of assets
for ass = 1:size(NumAssets,1)

	# random variable distribution parameters
	mean = Array{Float64}(NumAssets[ass])
	stdev = Array{Float64}(NumAssets[ass])
	for i = 1:NumAssets[ass]
		mean[i] = 1.05 + 0.3*(NumAssets[ass] - i)/(NumAssets[ass]-1)
		stdev[i] = (0.05 + 0.6*(NumAssets[ass] - i)/(NumAssets[ass]-1))/3
	end

	# matrix with stdevs on the diagonal
	sigma = zeros(Float64, NumAssets[ass], NumAssets[ass])
	for i = 1:NumAssets[ass]
		sigma[i,i] = stdev[i]
	end

	# iterate over the epsilons
	for iter_e = 1:size(epsilon,1)

		println("Iter: (",ass,",",iter_e,")  ","NumAssets: ",NumAssets[ass],"  epsilon: ",epsilon[iter_e])

		# create folder for storing results
		dirName = "C:/Users/rkannan/Desktop/SA for CCP/portfolio_analytical/" *
					string(NumAssets[ass]) * "/" * string(iter_e) * "/"
		mkpath(dirName)

		# write details of case study to text file
		details_file = dirName * "portfolio1.txt"
		open(details_file, "w") do f
			write(f,"Model: portfolio1 \n")
			write(f,"Method: EXACT \n")
			write(f,"Number of assets: $(NumAssets[ass]) \n")
			write(f,"epsilon: $(epsilon[iter_e]) \n")
		end
		# write epsilon to text file
		epsilon_file = dirName * "epsilon.txt"
		open(epsilon_file, "w") do f
			write(f,"$(epsilon[iter_e]) \n")
		end

		# pass params as keyword arguments to GurobiSolver
		m = Model(solver=GurobiSolver(Presolve=0,OutputFlag=0))

		@variable(m, objvar)
		@variable(m, 0 <= y[1:NumAssets[ass]] <= 1)

		@objective(m, Max, objvar)
		@constraint(m, sum(y[i] for i = 1:NumAssets[ass]) == 1)
		@constraint(m, soc, sum(mean[i]*y[i] for i = 1:NumAssets[ass]) >= objvar + norminvcdf(1-epsilon[iter_e])*norm(sigma*y))

		status = solve(m)

		objvar_opt = getvalue(objvar)
		time_to_opt = getsolvetime(m)

		y_opt = getvalue(y)
		soln_file = dirName * "soln.txt"
		open(soln_file, "w") do f
			for i = 1:NumAssets[ass]
				write(f,"$(y_opt[i]) ")
			end
		end

		# write objective
		obj_file = dirName * "objvals.txt"
		writeMode = "w"
		open(obj_file, writeMode) do f
			write(f,"$objvar_opt \n")
		end
		# write solution time
		time_file = dirName * "solntime.txt"
		open(time_file, writeMode) do f
			write(f,"$time_to_opt \n")
		end

		println("Objective: ",objvar_opt,"  Time: ",time_to_opt)
		println("")

		sleep(2)

	end

	sleep(5)
end

max_stars_repo_name: rohitkannan/SA-for-CCP

max_issues_repo_name: rohitkannan/SA-for-CCP

max_forks_repo_name: rohitkannan/SA-for-CCP


text:
module SAMTools

## Modules Used
using CFTime
using Crayons, Crayons.Box
using Dates
using DelimitedFiles
using Glob
using JLD2
using NCDatasets
using Printf
using Statistics

## Exporting the following functions:
export
        templaterun, templateexp, templatehpc, overwritetemplate,
        sammakefile, sammodules, samscratch

function __init__()
    samfol = joinpath(DEPOT_PATH[1],"files","SAMTools")
    if !isdir(samfol); mkpath(samfol) end
    overwritetemplate()
end

## Including other files in the module
include("setup.jl")

# include("analysis/topofatmosbalance.jl")

end


text:
"""
    BTEMobjective( a, X )
Returns the scalar BTEM objective function obtained from the linear combination vector `a` and loadings `X`.
*Note: This is not the function used in the original paper. This will be updated... it was written from memory.*
"""
function BTEMobjective( a, X )
    LinComb = a' * X
    Deriv = Scale1Norm( ScaleMinMax( FirstDerivative( LinComb ) ) )
    H = entropy( Deriv )
    #These aren't the penalties used in the original method. I don't have the paper handy...
    Negatives = LinComb .< 1e-6
    Penalty = sum( abs.( LinComb[ Negatives ] ) )
    return H + Penalty
end

"""
    BTEM(X, bands = nothing; Factors = 3, particles = 50, maxiters = 1000)
Returns a single recovered spectra from a 2-Array `X`, the selected `bands`, number of `Factors`, using a Particle Swarm Optimizer.
*Note: This is not the function used in the original paper. This will be updated... it was written from memory. Also the original method uses Simulated Annealing not PSO.*
Band-Target Entropy Minimization (BTEM):  An Advanced Method for Recovering Unknown Pure Component Spectra. Application to the FTIR Spectra of Unstable Organometallic Mixtures. Wee Chew,Effendi Widjaja, and, and Marc Garland. Organometallics 2002 21 (9), 1982-1990. DOI: 10.1021/om0108752
"""
function BTEM(X, bands = nothing; Factors = 3, particles = 50, maxiters = 1000)
    ( Obs, Vars ) = size( X )
    if isa( bands, Nothing ); bands = 1 : Vars; end
    pca = PCA( X )
    (Obs, Vars) = size(Mixture)
    if isa( bands, Nothing ); bands = 1 : Vars; end
    pca = PCA(Mixture; Factors = Factors);

    (a, score, otherparticles) = PSO(objective, Bounds(-10, 10, Factors), Bounds(-0.1, 0.1, Factors), particles;
                    tolerance = 1e-6, maxiters = maxiters,
                    InertialDecay = 0.25, PersonalWeight = 0.5, GlobalWeight = 0.5, InternalParams = pca.Loadings[:,bands]  )
    #PSO is pretty nice... It returns N other particles,
    #But whatever... Let's just return the main one... It's also a lot faster then MO-SA and easier to write...
    return a
end

"""
    NMF(X; Factors = 1, tolerance = 1e-7, maxiters = 200)
Performs a variation of non-negative matrix factorization on Array `X` and returns the a 2-Tuple of (Concentration Profile, Spectra)
*Note: This is not a coordinate descent based NMF. This is a simple fast version which works well enough for chemical signals*

Algorithms for non-negative matrix factorization. Daniel D. Lee. H. Sebastian Seung. NIPS'00 Proceedings of the 13th International Conference on Neural Information Processing Systems. 535-54
"""
function NMF(X; Factors = 1, tolerance = 1e-7, maxiters = 200)
    (Obs, Vars) = size(X)
    W = abs.( randn( Obs , Factors ) )
    H = abs.( randn( Factors, Vars ) )

    Last = zeros(Obs,Vars)
    tolerancesq = tolerance ^ 2
    iters = 0
    #Monitor change in F norm...
    while (sum( ( Last .- W * H ) .^ 2)  > tolerancesq) && (iters < maxiters)
        Last = W * H
        H .*= ( W' * X ) ./ ( W' * W * H )
        W .*= ( X * H' ) ./ ( W * H * H')
        iters += 1
    end
    return (W, H)
end

"""
    SIMPLISMA(X; Factors = 1, alpha = 0.05, includedvars = 1:size(X)[2], SecondDeriv = true)
Performs SIMPLISMA on Array `X` using either the raw spectra or the Second Derivative spectra.
alpha can be set to reduce contributions of baseline, and a list of included variables in the determination
of pure variables may also be provided.
Returns a tuple of the following form: (Concentraion Profile, Pure Spectral Estimates, Pure Variables)

W. Windig, Spectral Data Files for Self-Modeling Curve Resolution with Examples Using the SIMPLISMA Approach, Chemometrics and Intelligent Laboratory Systems, 36, 1997, 3-16.
"""
function SIMPLISMA(X; Factors = 1, alpha = 0.03, includedvars = 1:size(X)[2],
                    SecondDeriv = true)
    @warn("SIMPLISMA has not been completely tested for correctness.")
    Xcpy = deepcopy(X)
    X = X[:,includedvars]
    if SecondDeriv
        X = map( x -> max( x, 0.0 ), -SecondDerivative( X ) )
    end
    purvarindex = []
    (obs, vars) = size(X)
    pureX = zeros( Factors, vars )
    weights = zeros( vars )

    Col_Std = Statistics.std(X, dims = 1) .* sqrt( (obs - 1) / obs);
    Col_Mu = Statistics.mean(X, dims = 1);
    Robust_Col_Mu = Col_Mu .+ (alpha .* reduce(max, Col_Mu) );
    NormFactor = sqrt.( ((Col_Std .+ Robust_Col_Mu).^ 2) .+ (Col_Mu .^ 2) )
    Xp = X ./ NormFactor

    purity = Col_Std ./ Robust_Col_Mu

    for i in 1 : (Factors)
        for j in 1 : vars
           purvarmatrix = Xp[ : , vcat( purvarindex, j) ] ;
           O = (purvarmatrix' * purvarmatrix) ./ obs
           weights[j] = det( O );
       end
       purity_Spec = weights .* purity'
       push!(purvarindex, argmax(purity_Spec)[1])
       pureX[i,:] = purity_Spec;
    end
    pureinX = Xcpy[ : , includedvars[purvarindex] ]
    purespectra = pureinX \ Xcpy
    pureabundance = Xcpy / purespectra

    scale = LinearAlgebra.Diagonal(1.0 ./ sum(pureabundance, dims = 2))
    pureabundance = pureabundance * scale
    purespectra = Base.inv( scale ) * purespectra
    return (pureabundance, purespectra, includedvars[purvarindex])
end

"""
    FNNLS( A, b; LHS = false, maxiters = 500 )
Uses an implementation of Bro et. al's Fast Non-Negative Least Squares on the matrix `A` and vector `b`.
Returns regression coefficients in the form of a vector.

Bro, R., de Jong, S. (1997) A fast non-negativity-constrained least squares algorithm. Journal of Chemometrics, 11, 393-401.
"""
function FNNLS(A, b; maxiters = 500)
    ATA = A' * A
    ATb = A' * b
    X = zeros( size( ATA )[ 2 ] )
    tolerance = eps(Float16) * reduce(max,sum(abs.(ATA), dims = 1)) * prod(size(ATb))
    P = zeros( size( ATA )[ 2 ] ) .|> Int
    R = collect( 1 : length( ATb ) ) .|> Int
    r = zeros( size( ATA )[ 2 ] )
    Inds = collect( 1 : length( ATb ) ) .|> Int
    Rinds = collect( 1 : length( ATb ) ) .|> Int

    W = ATb
    inneriters = 0
    zeroint = 0 |> Int
    while any( R .> zeroint ) && any(W[Rinds] .> tolerance)
        Rinds = Inds[ R .> zeroint ]
        j = Rinds[argmax( W[Rinds] )]
        P[j] = j ; R[j] = zeroint
        Pinds = Inds[ P .> zeroint ] ; Rinds = Inds[ R .> zeroint ]
        #Update R & P
        r[Pinds] = Base.inv( ATA[ Pinds, Pinds ]) * ATb[ Pinds ]
        r[Rinds] .= 0.0
        while any(r[Pinds] .<= tolerance) && (inneriters < maxiters)
            Select = (r .<= tolerance) .& (P .> zeroint)
            alpha = reduce( min, X[Select] ./ ( X[Select] .- r[Select] ) )
            X .+= alpha .* ( r .- X )
            Select = Inds[(abs.(X) .< tolerance) .& ( P .!= zeroint )]
            R[Select] .= Select
            P[Select] .= zeroint
            Pinds = Inds[ P .> zeroint ] ; Rinds = Inds[ R .> zeroint ]
            r[Pinds] = Base.inv( ATA[ Pinds, Pinds ]) * ATb[ Pinds ]
            r[Rinds] .= 0.0
            inneriters += 1
        end
        X = r
        W = ATb .- (ATA * X)
    end

    return X
end

"""
    UnimodalFixedUpdate(x)
This function performs a unimodal least squares update at a fixed maximum for a vector x.
This is faster then UnimodalUpdate() but, is less accurate.

Bro R., Sidiropoulos N. D.. Least Squares Algorithms Under Unimodality and Non-Negativity Constraints
"""
function UnimodalFixedUpdate(x::Array{Float64,1})
    bins = length(x)
    if bins == 1
        return x
    end
    maxindx = argmax(x)[1]
    #Handle edge cases
    if (maxindx == 1)
        return reverse(MonotoneRegression(reverse(x)))
    elseif (maxindx == bins)
        return MonotoneRegression(x)
    end
    bLeft = MonotoneRegression(x[1:(maxindx-1)])
    bRight = reverse(MonotoneRegression(reverse(x[(maxindx+1):end])))
    if x[maxindx] > max( bLeft[maxindx-1], bRight[1] )
       return [bLeft; x[maxindx]; bRight];
    end
    return [bLeft; x[maxindx]; bRight]
end

"""
    UnimodalUpdate(x)

This function performs a unimodal least squares update for a vector x.
This is slower then UnimodalUpdate() but, is more accurate.

Bro R., Sidiropoulos N. D.. Least Squares Algorithms Under Unimodality and Non-Negativity Constraints. Journal of Chemometrics, June 3, 1997
"""
function UnimodalUpdate(x::Array{Float64,1})
    bins = length(x)
    if bins == 1
        return x
    end
    maxindx = argmax(x)[1]
    bLeft = MonotoneRegression(x)
    bRight = reverse(MonotoneRegression(reverse(x)))
    (SSE, BestSSE) = (Inf, Inf)
    (b, bestB) = (zeros( bins ),zeros( bins )) #Dummy initialization
    for (indx, value) in enumerate( bRight .+ bLeft )
        if indx == 1
            bRightFine = reverse(MonotoneRegression(reverse(x[(indx+1):end])))
            b = [ x[indx]; bRightFine ]
        elseif indx == bins
            bLeftFine = MonotoneRegression(x[1:(indx-1)])
            b = [ bLeftFine; x[indx] ]
        elseif (x[indx] >= (value / 2.0))
            bLeftFine = MonotoneRegression(x[1:(indx-1)])
            bRightFine = reverse(MonotoneRegression(reverse(x[(indx+1):end])))
            b = [ bLeftFine; x[indx]; bRightFine ]
        end
        SSE = sum(abs2.(b .- x))
        if SSE < BestSSE
            bestB = copy(b)
            BestSSE = SSE
        end
    end
    return bestB
end

"""
    UnimodalLeastSquares(x)

This function performs a unimodal least squares regression for a matrix A and b (X and Y).

Bro R., Sidiropoulos N. D.. Least Squares Algorithms Under Unimodality and Non-Negativity Constraints.Journal of Chemometrics, June 3, 1997
"""
function UnimodalLeastSquares(  A::Union{Array{Float64,1}, Array{Float64,2}},
                                b::Union{Array{Float64,1}, Array{Float64,2}};
                                maxiters::Int = 500, fixed::Bool = false)
    (obs, vars) = size( A )
    obspreds = size( b )
    (obs, preds) = (length(obspreds) > 1) ? obspreds : (obspreds, 1)
    B = randn( vars, preds )
    #Obs x Vars * Vars x Preds = Obs x Preds
    iter = 0
    LastB = ones( vars, preds ) .* Inf
    while (sum(abs2.(LastB .- B) ) / sum(abs2.( B )) > 1e-11) && (iter < maxiters)
        LastB[:,:] = B
        for var in 1:vars
            Cols = vcat(collect.([1:(var-1), (var+1):vars])...)
            beta = LinearAlgebra.pinv(A[ :, var ]) * (b - ( A[:,Cols] * B[Cols,:] ))
            if length(beta) == 1
                B[var, :] = beta
            else
                B[var, :] = (fixed) ? UnimodalFixedUpdate( beta' ) : UnimodalUpdate( beta' )
            end
        end
        iter += 1
    end
    return B
end

"""
    MCRALS(X, C, S = nothing; norm = (false, false), Factors = 1, maxiters = 20, constraintiters = 500, nonnegative = (false, false), unimodalS = false, fixedunimodal = false  )
Performs Multivariate Curve Resolution using Alternating Least Squares on `X` taking initial estimates for either `S` or `C`.
The number of maximum iterations for the ALS updates can be set `maxiters`.
S or C can be constrained by their `norm`(true/false,true/false), or by nonnegativity by using `nonnegative` arguments (true/false,true/false).
S can also be constrained by unimodality(`unimodalS`). Two unimodal algorithms exist the `fixedunimodal`(true), and the quadratic (false).
The number of resolved `Factors` can also be set.
The number of maximum iterations for constraints can be set by `constraintiters`.

Tauler, R. Izquierdo-Ridorsa, A. Casassas, E. Simultaneous analysis of several spectroscopic titrations with self-modelling curve resolution.Chemometrics and Intelligent Laboratory Systems. 18, 3, (1993), 293-300.
"""
function MCRALS(X::Union{Array{Float64,1}, Array{Float64,2}},
                C::Union{Nothing, Array{Float64,1}, Array{Float64,2}},
                S::Union{Nothing, Array{Float64,1}, Array{Float64,2}} = nothing;
                Factors::Int = 1, maxiters::Int = 20, constraintiters::Int = 50,
                norm::Tuple{Bool,Bool} = (false, false),
                nonnegative::Tuple{Bool,Bool} = (false, false),
                unimodalS::Bool = false, fixedunimodal::Bool = false )
    @assert all( isa.( [ C , S ], Nothing ) ) == false
    lowestErr = Inf
    elements = prod( size(X) )
    err = zeros(maxiters)
    D = X
    isC = isa(C, Nothing)
    isS = isa(S, Nothing)
    C = isC ? zeros(size(X)[1], Factors) : C[:, 1:Factors]
    S = isS ? zeros(Factors, size(X)[2]) : S[1:Factors, :]
    C ./= norm[1] ? sum(C, dims = 2) : 1.0
    S ./= norm[2] ? sum(S, dims = 1) : 1.0
    output = (C, S, err)
    for iter in 1 : maxiters
        if !isS
            if nonnegative[2]
                for obs in 1 : size( X )[ 1 ]
                    C[obs,:] = FNNLS(S', X[obs,:]; maxiters = constraintiters)
                end
            else
                C = X * LinearAlgebra.pinv(S)
            end
            if norm[1]
                C ./= sum(C, dims = 2)
            end
            isC = false
            D = C * S
        end
        if !isC
            if unimodalS
                S = UnimodalLeastSquares(C, X; maxiters = constraintiters, fixed = fixedunimodal)
                if nonnegative[1]
                    S[S .< 0.0] .= 0.0
                end
            elseif nonnegative[1]
                for var in 1:size(X)[2]
                    S[:,var] = FNNLS(C, vec(X[:,var]); maxiters = constraintiters)
                end
            else
                S = LinearAlgebra.pinv(C) * X
            end
            if norm[2]
                S ./= sum(S, dims = 1)
            end
            isS = false
            D = C * S
        end
        err[iter] = sum(abs2.(X .- D)) / elements
        if err[iter] < lowestErr
            lowestErr = err[iter]
            output = (C, S, err)
        end
    end
    return output
end

"""
    ITTFA(X; Factors = 1, maxiters = 500, threshold = 1e-8,
                    nonnegativity = true)

This function performs a Iterative Target Transform Factor Analysis (ITTFA).
It's untested. Use at your own risk.

Chemometric characterization of batch reactions. ISA Transactions 38(3):211-216. July 1999.
DOI: 10.1016/S0019-0578(99)00022-1
"""
function ITTFA(X; Factors = 1, Components = Factors, maxiters = 500,
                threshold = 1e-8, nonnegativity = true)
    @warn("ITTFA has not been tested for correctness.")
    rows, vars = size( X );
    Result = zeros( Components, rows )
    selectedneedles = zeros( Components )
    PCA = LinearAlgebra.svd( X );
    T = PCA.U
    needlematrix = LinearAlgebra.Diagonal( ones( rows ) );
    estimates = zeros( rows, Factors )
    #Find the best locations for needles based on their projected scores
    possiblecomponents = zeros( Factors, Components )
    possiblelocations = zeros( Factors, Components )
    for f in 1 : Factors
        projscore = zeros( rows )
        for i in 1 : rows #Maybe they limit these to the window span?
            needlevector = needlematrix[ :, i ]
            needlevector_old = copy( needlevector );
            needlevector = T[:,f] * T[:,f]' * needlevector
            projscore[ i ] = sqrt( sum( (needlevector .- needlevector_old) .^ 2 ) )
        end
        #Get the smallest projections of the needles!
        possiblelocations[f, : ] = sortperm( projscore )[ 1 : Components ]
        possiblecomponents[f, : ] = projscore[ Int.( possiblelocations[f, : ] ) ]
    end
    #Let's find the best needles! by consuming the smallest projected norms
    bestlocations = zeros(Components)
    for c in 1 : Components
        bestcomponents, idx = findmin( possiblecomponents )
        possiblecomponents[idx] = Inf
        bestlocations[c] = possiblelocations[ idx ]
        #Note this allows for duplicates... can fix...
        #toberemoved = findall( possiblelocations .== possiblelocations[ idx ] )
        #possiblecomponents[toberemoved] .= Inf
    end
    #Now let's assume we have good needles and go forward with curve resolution via ITTFA
    for c in 1 : Components
        #Initialize optimized needle for each component
        needlevector_old = zeros( rows )
        needlevector = needlematrix[ :, Int( bestlocations[c] ) ]
        iter = 0
        while ( iter < maxiters ) && ( sum( ( needlevector - needlevector_old ) .^ 2 ) > 1e-6 )
            needlevector_old = copy( needlevector );
            needlevector = T[:,1:c] * T[:,1:c]' * needlevector
            if nonnegativity
                needlevector[ needlevector .< 0.0 ] .= 0.0
            end
            iter += 1
    	end
        Result[c, : ] = copy(needlevector)
    end
    #Project back to X using ALS
    purespectra = LinearAlgebra.pinv( Result )'  * X
    pureabundance = X * LinearAlgebra.pinv( purespectra )
    return purespectra, pureabundance, bestlocations
end

#I believe the SIMPLISMA implementation below has errors. It's a super neat algorithm, but it does
#not display expected behaviour...

# """
#     SIMPLISMA(X; Factors = 1)
# Performs SIMPLISMA on Array `X`.
# Returns a tuple of the following form: (Concentraion Profile, Pure Spectral Estimates, Pure Variables)
# *Note: This is not the traditional SIMPLISMA algorithm presented by Willem Windig.*
# REAL-TIME WAVELET COMPRESSION AND SELF-MODELING CURVE RESOLUTION FOR ION MOBILITY SPECTROMETRY. PhD. Dissertation. 2003. Guoxiang Chen.
# """
# function SIMPLISMA(X; Factors = 1)
#     (obs, vars) = size(X)
#     PurestVar = ones(Factors) .|> Int
#     Ortho = zeros(obs, Factors)
#     SSE = StatsBase.sum(X .^ 2, dims = 1)
#     e = (SSE .- StatsBase.sum(X, dims = 1).^2) / obs #RSE/SSE
#
#     PurestVar[1] = argmax(vec(e))
#     Intensity = X[:, PurestVar[1]]
#     Ortho[:,1] = Intensity ./ sqrt( Intensity' * Intensity )#2-norm
#     for F in 2 : Factors
#         proj = sum( (Ortho[:,1:(F-1)]' * X) .^ 2, dims = 1)
#         p = vec(e .* ( 1.0 .- ( proj ./ SSE ) ) )
#         p[PurestVar[1:F]] .= -Inf
#         PurestVar[F] = argmax( p )
#         Intensity = X[ : , PurestVar[ F ] ]
#         OrthTmp = Intensity .- sum(proj * ( proj' * Intensity' ) )
#         Ortho[:,F] = OrthTmp ./ sqrt.( OrthTmp' * OrthTmp )
#     end
#
#     C = X[:, PurestVar]
#     S = LinearAlgebra.pinv(C) * X
#     magnitude = vec( sum( S .^ 2, dims = 2 ) )
#     for F in 1:Factors
#         S[F,:] = (magnitude[F] <= 1e-8) ? (S[F,:] .* 0.0) : (S[F,:] ./ sqrt(magnitude[F]))
#     end
#     S .+= StatsBase.mean(X, dims = 1)
#     return (C, S, PurestVar)
# end


text:
__precompile__()

"""
Module to prepare convolutive mixtures.
And ICA is used for Blind Source Separation (BSS) of these mixtures.
"""
module ConvMix

import JuMIT.Conv
import JuMIT.CICA
import JuMIT.Grid
import JuMIT.Misfits
using DataFrames
using RecipesBase

mutable struct Param
	ns::Int # number of independent components
	obs::Vector{Conv.Param{Float64,2}} # observed independent components
	nomix::Vector{Conv.Param{Float64,2}} # unmixed data for simple xcorr
	mix::Vector{Conv.Param} # mixture, results of simple xcorr stored here
	"ica model"
	ica::CICA.ICA # results after performing ICA
	icaunmix::Vector{Conv.Param} # data after ica unmixing, this can be used for further processing
	"error in ge compared to gb"
	err::DataFrames.DataFrame
end

@userplot Plot

@recipe function f(p::Plot, rvec=nothing)
	pa=p.args[1]
end

"""
Greens functions are input to this constructor.
Source signal models are
* `ntgf` : samples in Green function
* `nt` : samples in the long time series

* `ssparsepvec` : spartisty

"""
function Param(;ntgf=1, # instantaneos
	       nt=500,
	       nr=2,
	       ns=2,
	       gfobs=[randn(ntgf,nr) for is in 1:ns],
	       wavobs=[randn(nt) for is in 1:ns],
	       ica_func::Function=(x,ns)->CICA.ICA(x,ns),
	)

	obs=[Conv.Param(ntwav=nt, ntd=nt, ntgf=ntgf, dims=(nr,), wavlags=[nt-1, 0]) for is in 1:ns]
	nomix=deepcopy(obs);
	mix=deepcopy(obs);
	icaunmix=deepcopy(obs);

	for is in 1:ns
		# copy observed Greens functions
		copy!(obs[is].gf,gfobs[is])
		# copy oroginal source wavelets
		for ir in 1:nr
			for i in 1:nt
				obs[is].wav[i,ir]=wavobs[is][i]
			end
		end
	end

	# model observed data for each source
	for is in 1:ns
		Conv.mod!(obs[is], :d)
	end

	# copy d and wav to nomix
	for is in 1:ns
		copy!(nomix[is].wav, obs[is].wav)
		copy!(nomix[is].d, obs[is].d)
	end

	# estimated gf without mixing and xcorr
	for is in 1:ns
		Conv.mod!(nomix[is], :gf)
	end

	#
	# copy sum of d and wav to mix
	for is in 1:ns
		for iss in 1:ns
			for i in eachindex(mix[is].d)
				mix[is].d[i] += obs[iss].d[i]
			end
		end
		copy!(mix[is].wav, obs[is].wav)
	end

	# copy wavelets to icaunmix, replace this step with blind deconvolution?
	for is in 1:ns
		copy!(icaunmix[is].wav, obs[is].wav)
	end

	# estimated gf after mixing and xcorr
	for is in 1:ns
		Conv.mod!(mix[is], :gf)
	end

	x=transpose(rfft(mix[1].d, [1]))
	ica=ica_func(x, ns)

	err=DataFrame(rec_src=vec([(ir,is) for ir in 1:nr, is in 1:ns]), gf_no_mix_xcorr=zeros(nr*ns), gf_mix_xcorr=zeros(nr*ns), d_mix_ica=zeros(nr*ns),
	       gf_mix_ica=zeros(nr*ns))

	pa=Param(ns, obs, nomix,  mix, ica, icaunmix, err)

	update_err(pa)

	return pa

end

function update_err(pa)
	nr=size(pa.obs[1].gf,2)
	ns=pa.ns
	for is in 1:pa.ns
		for ir in 1:nr
			d1=view(pa.nomix[is].gf, :, ir)
			d2=view(pa.mix[is].gf, :, ir)
			d0=view(pa.obs[is].gf, :, ir)
			pa.err[:gf_no_mix_xcorr][ir+(is-1)*nr]=Misfits.error_after_scaling(d1, d0)[1]
			pa.err[:gf_mix_xcorr][ir+(is-1)*nr]=Misfits.error_after_scaling(d2, d0)[1]
		end
	end
end
function Base.print(pa::Param)
	println("Convolutive Mixture\t")
	println("===================")
	println(pa.err)
end

function update_ica(cmix)
	x=rfft(cmix.d, [2])
	cmix.ica=ica_func(x, 2)
end

"""
Do deblending for every receiver individually.
And then update `icaunmix.d`
Then, `err_ica` is also updated with scaled least-squares error.
"""
function fastica!(pa::Param)
	nr=size(pa.obs[1].gf,2)
	ns=pa.ns

	for ir in 1:nr
		pa.ica.magic_recv=ir
		CICA.fastica!(pa.ica)

		s=transpose((pa.ica.s))
		s=irfft(s,pa.obs[1].ntd,[1])

		for is in 1:ns
			ss=view(s,:,is)

			err=zeros(ns)
			for iss=1:ns
				d0=view(pa.obs[iss].d,:,ir)
				err[iss]=Misfits.error_after_scaling(ss, d0)[1]
			end
			pa.err[:d_mix_ica][ir+(is-1)*nr]=minimum(err)

			# fill to the closest source, only for synthetic!
			d=view(pa.icaunmix[indmin(err)].d,:,ir)
			copy!(d,ss)
		end
	end

	# use xcorr to estimate green functions
	for is in 1:ns
		Conv.mod!(pa.icaunmix[is], :gf)
	end

	# compute error in estimated gf assuming source signatures are known
	for is in 1:pa.ns
		for ir in 1:nr
			d=view(pa.icaunmix[is].gf, :, ir)
			d0=view(pa.obs[is].gf, :, ir)
			pa.err[:gf_mix_ica][ir+(is-1)*nr]=Misfits.error_after_scaling(d, d0)[1]
		end
	end

end

end


text:
using Revise # reduce recompilation time
using Plots
# using Documenter
using LinearAlgebra
using SparseArrays
using BenchmarkTools
using UnPack
using ToeplitzMatrices

push!(LOAD_PATH, "./src")
using CommonUtils
using Basis1D
using Basis2DQuad
using UniformQuadMesh

using SetupDG

# TODO: refactor the code
function fS(uL, uR)
    return 1/6*(uL^2+uR^2+uL*uR)
end

function burgers_exact_sol_2D(u0,x,y,T,dt)
    Nsteps = ceil(Int,T/dt)
    dt = T/Nsteps
    u = u0(x,y) # computed at input points
    for i = 1:Nsteps
        t = i*dt
        u .= @. u0(x-u*t,y) # evolve solution at quadrature points using characteristics
    end
    return u
end

S_N(x) = @. sin(pi*x/h)/(2*pi/h)/tan(x/2)
"""
Vandermonde matrix of sinc basis functions determined by h,
evaluated at r
"""
function vandermonde_Sinc(h,r)
    N = convert(Int, 2*pi/h)
    V = zeros(length(r),N)
    for n = 1:N
        V[:,n] = S_N(r.-n*h)
    end
    V[1,1] = 1
    return V
end

"Constants"
const sp_tol = 1e-12
# u0(x,y) = @. -sin(pi*x)*sin(y)
u0(x,y) = @. exp(-1*((y-pi)^2+2*pi*x^2))
# u0(x,y) = @. exp(-1*((y-pi)^2))
# u0(x,y) = @. exp(-10*(x^2))

"Program parameters"
compute_L2_err = false

"Approximation Parameters"
N_P = 8;    # The order of approximation in polynomial dimension
Np_P = N_P+1;
Np_F = 16;    # The order of approximation in Fourier dimension
CFL  = 0.2;
T    = 3.0;  # End time

"Time integration Parameters"
rk4a,rk4b,rk4c = rk45_coeffs()
CN = (max(N_P,Np_F)+1)*(max(N_P,Np_F)+2)/2  # estimated trace constant for CFL
dt = CFL / CN
Nsteps = convert(Int,ceil(T/dt))
dt = T/Nsteps

"Initialize Reference Element in Fourier dimension"
h = 2*pi/Np_F
column = [0; .5*(-1).^(1:Np_F-1).*cot.((1:Np_F-1)*h/2)]
column2 = [-pi^2/3/h^2-1/6; -((-1).^(1:Np_F-1)./(2*(sin.((1:Np_F-1)*h/2)).^2))]
Ds = Toeplitz(column,column[[1;Np_F:-1:2]])
D2s = Toeplitz(column2,column2[[1;Np_F:-1:2]])
s = LinRange(h,2*pi,Np_F)

"Initialize Reference Element in polynomial dimension"
rd = init_reference_interval(N_P);
@unpack r,rq,wq,rf,Dr,VDM,Vq,Vf,Vp,M,Pq,LIFT,nrJ = rd
wf = [1;1]
Nq_P = length(rq)
rp = LinRange(-1,1,100)
Vp = vandermonde_1D(N_P,rp)/VDM # Redefine Vp
Lq = LIFT

"Initialize Mesh"
s,r = meshgrid(s,r)
Nfp = 2*Np_F

"scale by Fourier dimension"
M = h*M
wq = h*wq
wf = h*wf

"Hybridized operators"
Ef = Vf*Pq
Br = diagm(wf.*nrJ)
Qr = Pq'*M*Dr*Pq
Qs = Ds
Drh = [Vq*Dr*Pq-1/2*Vq*Lq*diagm(nrJ)*Vf*Pq  1/2*Vq*Lq*diagm(nrJ);
       -1/2*diagm(nrJ)*Vf*Pq                1/2*diagm(nrJ)]
Qrh = 1/2*[Qr-Qr' Ef'*Br;
           -Br*Ef Br]
Qrh_skew = 1/2*(Qrh-Qrh')
Qsh = Qs # Not the SBP operator, weighted when flux differencing

ops = (Vq,Vf,wq,wf,Pq,Lq,Qrh,Qsh,D2s,nrJ,h)
function rhs(u,ops,compute_rhstest)
    Vq,Vf,wq,wf,Pq,Lq,Qrh,Qsh,D2s,nrJ,h = ops
    # Volume term
    uq,uf = (A->A*u).((Vq,Vf))
    uh = [uq;uf]
    w = [wq;wf]
    ∇fh = zeros(size(uh))

    # Flux differencing in x direction
    for k = 1:size(uh,2)
        ∇fh[:,k] += 2*Qrh.*[fS(uL,uR) for uL in uh[:,k], uR in uh[:,k]]*ones(size(uh,1),1)
    end

    # Flux differencing in y direction
    for k = 1:size(uh,1)
        ∇fh[k,:] += 2*pi*[wq;0;0][k]*Qsh.*[fS(uL,uR) for uL in uh[k,:], uR in uh[k,:]]*ones(size(uh,2),1)
    end

    # Spatial term
    ∇f = [Pq Lq]*diagm(1 ./ w)*∇fh

    # Flux term
    uM = Vf*u
    uP = [uM[2,:] uM[1,:]]'
    LF = @. max(abs(uP),abs(uM))*(uP-uM)
    uflux = diagm(nrJ)*(@. fS(uM,uP)-uM*uM/2)-LF

    # Artificial Viscosity
    # ϵ = 0.05
    ϵ = 0.0
    Δf = ϵ*u*D2s'
    rhsu = -(∇f+Lq*uflux-Δf)

    rhstest = 0
    if compute_rhstest
        rhstest += sum(diagm(wq)*uq.*(Vq*rhsu))
    end

    return rhsu,rhstest
end

u = @. u0(r,s)
# Test single element
resu = zeros(size(u))
for i = 1:Nsteps
    rhstest = 0
    for INTRK = 1:5
        compute_rhstest = INTRK==5
        rhsu,rhstest = rhs(u,ops,compute_rhstest)
        @. resu = rk4a[INTRK]*resu + dt*rhsu
        @. u += rk4b[INTRK]*resu
    end

    if i%10==0 || i==Nsteps
        println("Time step: $i out of $Nsteps with rhstest = $rhstest")
    end
end

# Plotting
gr(aspect_ratio=1/pi, legend=false,
   markerstrokewidth=0, markersize=2.5,
   xlims=(-1,1),ylims=(h,2*pi),showaxis=false)
plot()

sp = LinRange(h,2*pi,100)
VDM_F = vandermonde_Sinc(h,sp)
V2 = vandermonde_1D(1,sp)/vandermonde_1D(1,LinRange(h,2*pi,Np_F))

plt = scatter(Vp*r*V2',Vp*s*V2',Vp*u*VDM_F',zcolor=Vp*u*VDM_F',camera=(0,90))
display(plt)
# exact_sol = burgers_exact_sol_2D(u0,Vp*r*V2',Vp*s*V2',T,dt/100)
# scatter(Vp*r*V2',Vp*s*V2',exact_sol,zcolor=exact_sol,camera=(0,90))

# L2 error
if compute_L2_err
    rq2,wq2 = gauss_quad(0,0,N_P+2)
    wq2 = wq2*h
    Vq2 = vandermonde_1D(N_P,rq2)/VDM
    xq2 = Vq2*r
    yq2 = Vq2*s

    exact_sol = burgers_exact_sol_2D(u0,xq2,yq2,T,dt/100)
    L2_err = sum(diagm(wq2)*(Vq2*u-exact_sol).^2)
    @show L2_err
end

max_stars_repo_name: yiminllin/ESDG-Fourier

max_issues_repo_name: yiminllin/ESDG-Fourier

max_forks_repo_name: yiminllin/ESDG-Fourier


text:
#!/bin/bash
# -*- mode: julia -*-
#=
JULIA="${JULIA:-julia --color=yes --startup-file=no}"
export JULIA_PROJECT="$(dirname ${BASH_SOURCE[0]})"

set -ex
${JULIA} -e 'using Pkg; Pkg.instantiate();
             Pkg.develop(PackageSpec(path=pwd()))'
exec ${JULIA} "${BASH_SOURCE[0]}" "$@"
=#

include("build.jl")
include("deploy.jl")


text:
using CholmodSolve2
using Base.Test

m = 400
n = 500
A = randn(m, n)
Q = sparse([I A'; A -I])
y = randn(m+n)
x = similar(y)
F = ldltfact(Q)

precompile(A_ldiv_B!,typeof.((x, F, y)))
val, t, bytes, gctime, memallocs = @timed A_ldiv_B!(x, F, y)
@test bytes < 50000# should be roughly 900*(1+1+4)*8=43200
@test Q*x ≈ y
# Solve again, without allocs
val, t, bytes, gctime, memallocs = @timed A_ldiv_B!(x, F, y)
@test bytes < 400 # should be 192
@test Q*x ≈ y

# Compile for array of dimension 2 (still same size)
y2 = Array{Float64,2}(m+n,1); y2 .= x
x2 = similar(y2)
val, t, bytes, gctime, memallocs = @timed A_ldiv_B!(x2, F, y2)

y2 = [y 2y]
x2 = similar(y2)
# Test allocation with double size
# We will probably not allocate too much here (only X_Handle, and B_Handle)
val, t, bytes, gctime, memallocs = @timed A_ldiv_B!(x2, F, y2)
@test bytes < 30000# should be roughly 900*2*2*8=28800
@test Q*x2 ≈ y2

# Go back to 1 dimension and make sure no allocations
val, t, bytes, gctime, memallocs = @timed A_ldiv_B!(x, F, y)
@test bytes < 400 # should be 192
@test Q*x ≈ y

# Try with many dimensions (since CHOLMOD solves with 4 columns at a time)
y2 = [y 2y 3y 4y]
x2 = similar(y2)
# We will probably not allocate too much here (only X_Handle, and B_Handle)
val, t, bytes, gctime, memallocs = @timed A_ldiv_B!(x2, F, y2)
@test bytes < 60000 # should be roughly 900*2*4*8=57600
@test Q*x2 ≈ y2

y2 = [y2 y2]
x2 = similar(y2)
# We also allocate new Y_Handle (not only X_Handle, and B_Handle)
val, t, bytes, gctime, memallocs = @timed A_ldiv_B!(x2, F, y2)
@test bytes < 180000 # should be roughly 900*3*8*8=172800
@test Q*x2 ≈ y2

y2 = 2y2
x2 = similar(y2)
# Test with views
x22 = view(x2, :, 1:4)
y22 = view(y2, :, 1:4)
precompile(A_ldiv_B!,typeof.((x22, F, y22)))
val, t, bytes, gctime, memallocs = @timed A_ldiv_B!(x22, F, y22)
@test bytes < 400 # should be roughly 900*3*8*8=172800
@test Q*x22 ≈ y22
@test Q*x2[1:size(y2,1), 1:4] ≈ y2[1:size(y2,1), 1:4]

# Test with non-contiguous
x22 = view(x2, 1:2:2size(x2,1))
y22 = view(y2, 1:2:2size(y2,1))
y2 .= randn(size(y2))
precompile(A_ldiv_B!,typeof.((x22, F, y22)))
val, t, bytes, gctime, memallocs = @timed A_ldiv_B!(x22, F, y22)
@test bytes < 400 # should be roughly 900*3*8*8=172800
@test Q*x22 ≈ y22

# Change size of F
m = 200;  n = 300;
A2 = randn(m, n);
Q2 = sparse([I A2'; A2 -I]);
y3 = randn(m+n); x3 = randn(m+n);
F2 = ldltfact(Q2)

val, t, bytes, gctime, memallocs = @timed A_ldiv_B!(x3, F2, y3)
# We don't update size of Y here, so CHOLMOD will reallocate it
#@test bytes < 400
@test Q2*x3 ≈ y3

# Just test that print_common doesn't fail
ok = CholmodSolve2.print_common()
@test ok == 1
# Expected result:
# memory blocks in use:          39
# memory in use (MB):           8.9
# peak memory usage (MB):      16.6

F = 1.0
F2 = 1.0
gc()

#Should have less allocated now
ok = CholmodSolve2.print_common()
@test ok == 1
# memory blocks in use:          10
# memory in use (MB):           0.2
# peak memory usage (MB):      16.6


text:
import RData
using RNGPool
using CoupledConditionalSMC
using BenchmarkTools

# data generated with A = 0.95 (not 0.9)
dataJLS = RData.load("demo/ar1data.RData")
ys = vec(dataJLS["observations"])

include("../demo/lgModel.jl")

# they start at time 0 with observations from time 1,2,...
# so we have an initial distribution of N(0,1+0.9^2)
theta = LGTheta(0.9, 1.0, 1.0, 1.0, 0.0, 1.81)
model = makeLGModel(theta, ys)
lM = LinearGaussian.makelM(theta)

function foo(model, lM, ccsmcio, algorithm)
  for i in 1:ccsmcio.n
    ccsmcio.ref1[i].x = -1.0
    ccsmcio.ref2[i].x = 1.0
  end
  ccXpf!(model, lM, ccsmcio, algorithm, true)
end

n = 500
N = 512
ccsmcio = CCSMCIO{model.particle, model.pScratch}(N, n)

algorithm = :AT
initializeCCSMC(model, lM, ccsmcio, algorithm, true)
@btime foo(model, lM, ccsmcio, algorithm)

algorithm = :AS
initializeCCSMC(model, lM, ccsmcio, algorithm, true)
@btime foo(model, lM, ccsmcio, algorithm)

algorithm = :BS
initializeCCSMC(model, lM, ccsmcio, algorithm, true)
@btime foo(model, lM, ccsmcio, algorithm)

n = 1000
N = 1024
ccsmcio = CCSMCIO{model.particle, model.pScratch}(N, n)

algorithm = :AT
initializeCCSMC(model, lM, ccsmcio, algorithm, true)
@btime foo(model, lM, ccsmcio, algorithm)

algorithm = :AS
initializeCCSMC(model, lM, ccsmcio, algorithm, true)
@btime foo(model, lM, ccsmcio, algorithm)

algorithm = :BS
initializeCCSMC(model, lM, ccsmcio, algorithm, true)
@btime foo(model, lM, ccsmcio, algorithm)


text:
# Gradient regression
function get_gain_gpu(::L, ∑δ::AbstractVector{T}, λ::T) where {L <: GradientRegression,T <: AbstractFloat}
    gain = ∑δ[1]^2 / (∑δ[2] + λ * ∑δ[3]) / 2
    return gain
end

# Gaussian regression
function get_gain_gpu(::L, ∑δ::AbstractVector{T}, λ::T) where {L <: GaussianRegression,T <: AbstractFloat}
    gain = ∑δ[1]^2 / (∑δ[3] + λ * ∑δ[5]) / 2 + ∑δ[2]^2 / (∑δ[4] + λ * ∑δ[5]) / 2
    return gain
end

#####################
# linear
#####################
function kernel_linear_δ!(δ::CuDeviceMatrix{T}, p::CuDeviceMatrix{T}, y::CuDeviceVector{T}) where {T <: AbstractFloat}
    i = threadIdx().x + (blockIdx().x - 1) * blockDim().x
    ϵ = Cfloat(2e-7)
    if i <= length(y)
        @inbounds δ[i,1] = 2 * (p[i] - y[i]) * δ[i,3]
        @inbounds δ[i,2] =  max(ϵ, 2 * δ[i,3])
    end
    return
end

# base approach - block built along the cols first, the rows (limit collisions)
function update_grads_gpu!(loss::Linear, δ::CuMatrix{T}, p::CuMatrix{T}, y::CuVector{T}; MAX_THREADS=1024) where {T <: AbstractFloat}
    thread_i = min(MAX_THREADS, length(y))
    threads = (thread_i)
    blocks = ceil.(Int, (length(y)) ./ threads)
    @cuda blocks = blocks threads = threads kernel_linear_δ!(δ, p, y)
    return
end

#####################
# Logistic
#####################
function kernel_logistic_δ!(δ::CuDeviceMatrix{T}, p::CuDeviceMatrix{T}, y::CuDeviceVector{T}) where {T <: AbstractFloat}
    i = threadIdx().x + (blockIdx().x - 1) * blockDim().x
    ϵ = Cfloat(2e-7)
    if i <= length(y)
        @inbounds δ[i,1] = (sigmoid(p[i]) * (1 - y[i]) - (1 - sigmoid(p[i])) * y[i]) * δ[i,3]
        @inbounds δ[i,2] = max(ϵ, sigmoid(p[i]) * (1 - sigmoid(p[i])) * δ[i,3])
    end
    return
end

# base approach - block built along the cols first, the rows (limit collisions)
function update_grads_gpu!(loss::Logistic, δ::CuMatrix{T}, p::CuMatrix{T}, y::CuVector{T}; MAX_THREADS=1024) where {T <: AbstractFloat}
    thread_i = min(MAX_THREADS, length(y))
    threads = (thread_i)
    blocks = ceil.(Int, (length(y)) ./ threads)
    @cuda blocks = blocks threads = threads kernel_logistic_δ!(δ, p, y)
    return
end

################################################################################
# Gaussian - http://jrmeyer.github.io/machinelearning/2017/08/18/mle.html
# pred[i][1] = μ
# pred[i][2] = log(σ)
################################################################################
function kernel_gauss_δ!(δ::CuDeviceMatrix{T}, p::CuDeviceMatrix{T}, y::CuDeviceVector{T}) where {T <: AbstractFloat}
    i = threadIdx().x + (blockIdx().x - 1) * blockDim().x
    ϵ = Cfloat(2e-7)
    if i <= length(y)

        # first order gradients
        δ[i,1] = (p[i,1] - y[i]) / max(ϵ, exp(2f0 * p[i,2])) *  δ[i,5]
        δ[i,2] = (1f0 - (p[i,1] - y[i])^2 / max(ϵ, exp(2f0 * p[i,2]))) *  δ[i,5]

        # second order gradients
        δ[i,3] =  max(ϵ, δ[i,5] / max(ϵ, exp(2 * p[i,2])))
        δ[i,4] = max(ϵ, 2 * δ[i,5] / max(ϵ, exp(2 * p[i,2])) * (p[i,1] - y[i])^2)
    end
    return
end

# base approach - block built along the cols first, the rows (limit collisions)
function update_grads_gpu!(loss::Gaussian, δ::CuMatrix{T}, p::CuMatrix{T}, y::CuVector{T}; MAX_THREADS=1024) where {T <: AbstractFloat}
    thread_i = min(MAX_THREADS, length(y))
    threads = (thread_i)
    blocks = ceil.(Int, (length(y)) ./ threads)
    @cuda blocks = blocks threads = threads kernel_gauss_δ!(δ, p, y)
    return
end


text:
@info "Test composite pattern..."

pic = Picture()
add!(pic, Line(5.0))
add!(pic, Text_("ABC"))
add!(pic, Circle(5.0))
draw(pic)

pic2 = Picture()
add!(pic2, Text_("DEF"))
add!(pic2, pic)
draw(pic2)


text:
"""
update!(sim::PlanktonSimulaiton)
update the `PlanktonSimulaiton` for `sim.iterations` time steps.
"""
function update!(sim::PlanktonSimulation)
if sim.input.vels ≠ (;)
    for t in 1:sim.iterations
        t_vel = Int(sim.model.t÷sim.input.ΔT_vel)+1 # starting from 1
        vel_copy!(sim.model.timestepper.vel₁, sim.input.vels.u[:,:,:,t_vel],
                  sim.input.vels.v[:,:,:,t_vel], sim.input.vels.w[:,:,:,t_vel], sim.model.grid)

        # tc = sim.model.t % 86400 ÷ sim.ΔT + 1 # time of day in ΔT as Int index for temp and PAR
        t_par = Int(sim.model.t÷sim.input.ΔT_PAR)+1 # starting from 1
        copyto!(sim.model.timestepper.PARF, sim.input.PARF[:,:,t_par])

        t_temp = Int(sim.model.t÷sim.input.ΔT_temp)+1 # starting from 1
        copy_interior!(sim.model.timestepper.temp, sim.input.temp[:,:,:,t_temp], sim.model.grid)

        TimeStep!(sim.model, sim.ΔT, sim.diags)

        write_output!(sim.output_writer, sim.model, sim.diags, sim.ΔT)
    end
else
    for t in 1:sim.iterations
        t_par = Int(sim.model.t÷sim.input.ΔT_PAR)+1 # starting from 1
        copyto!(sim.model.timestepper.PARF, sim.input.PARF[:,:,t_par])

        t_temp = Int(sim.model.t÷sim.input.ΔT_temp)+1 # starting from 1
        copy_interior!(sim.model.timestepper.temp, sim.input.temp[:,:,:,t_temp], sim.model.grid)

        TimeStep!(sim.model, sim.ΔT, sim.diags)

        write_output!(sim.output_writer, sim.model, sim.diags, sim.ΔT)
    end
end
end


text:
module FreeEnergyMultimodeQED

ENV["CPLEX_STUDIO_BINARIES"] = "/opt/ibm/ILOG/CPLEX_Studio1210/cplex/bin/x86-64_linux"

using Random                     # utils
using LinearAlgebra              # model
using CPLEX, Dates, Ipopt, JuMP  # optimiziers
using DelimitedFiles             # simulations, parsing
using DataFrames, JLD, StatsBase # parsing

export coprimeFractionsGenerator
export gradMagnetizationEntropy, gradModelFreeEnergy, magnetizationEntropy, modelDimension, modelEntropy, modelFreeEnergy, modelInteraction
export modelFindMinimum
export simulateGS, simulateMinima
export analyzeGS, analyzeMinima, checkMinimaFound, getTotalRuns, parseRawData, preprocessRawDataframe

include("utils.jl")             # utility functions
include("model.jl")             # model-related functions
include("optimizers.jl")        # optimizers to find minima of the free energy of the model
include("simulations.jl")       # utlities for simulations
include("parsing.jl")           # utilities for data processing

end