text:
module Lux

# Accelerator Support
using CUDA
using CUDA.CUDNN
# Neural Network Backend
using NNlib
import NNlibCUDA: batchnorm, ∇batchnorm, CUDNNFloat
# Julia StdLibs
using Random, Statistics, LinearAlgebra, SparseArrays
# Parameter Manipulation
using Functors, Setfield
import Adapt: adapt, adapt_storage
# Arrays
using FillArrays, ComponentArrays
# Automatic Differentiation
using ChainRulesCore, Zygote
# Optimization
using Optimisers
# Optional Dependency
using Requires

const use_cuda = Ref{Union{Nothing, Bool}}(nothing)

# Data Transfer Utilities
include("adapt.jl")
# Utilities
include("utils.jl")
# Core
include("core.jl")
# Layer Implementations
include("layers/basic.jl")
include("layers/normalize.jl")
include("layers/conv.jl")
include("layers/dropout.jl")
include("layers/recurrent.jl")
# Neural Network Backend
include("nnlib.jl")
# Pretty Printing
include("layers/display.jl")
# AutoDiff
include("autodiff.jl")
# Flux to Lux
function __init__()
    @require Flux="587475ba-b771-5e3f-ad9e-33799f191a9c" include("transform.jl")
end

# Data Transfer
export cpu, gpu
# Layers
export Chain, Parallel, SkipConnection, PairwiseFusion, BranchLayer
export Dense, Scale
export Conv, MaxPool, MeanPool, GlobalMaxPool, GlobalMeanPool, AdaptiveMaxPool,
       AdaptiveMeanPool, Upsample
export Dropout, VariationalHiddenDropout
export BatchNorm, GroupNorm
export WeightNorm
export NoOpLayer, ReshapeLayer, SelectDim, FlattenLayer, WrappedFunction, ActivationFunction
export RNNCell, LSTMCell, GRUCell

end


text:
using AutoHOOT
using Test

const ad = AutoHOOT.autodiff
const gops = AutoHOOT.graphops

@testset "test rewrite expression" begin
    a1 = ad.Variable(name = "a1", shape = [3, 2])
    a2 = ad.Variable(name = "a2", shape = [2, 3])
    x = ad.einsum("ik,kj->ij", a1, a2)
    y = ad.einsum("sm,ml->sl", a1, a2)
    gops.rewrite_einsum_expr(x)
    gops.rewrite_einsum_expr(y)
    @test x.einsum_subscripts == y.einsum_subscripts
end

@testset "test optimizer and simplify" begin
    A = ad.Variable(name = "A", shape = [2, 2])
    out = ad.einsum("ab,bc->ac", A, ad.einsum("ab,bc->ac", ad.identity(2), ad.identity(2)))
    newout_optimize = gops.optimize(out)
    newout_simplify = gops.simplify(out)
    @test newout_optimize == A
    @test newout_simplify == A
end

ext: jl


text:
###################
## Packages ########
###################

using CSV
using DataFrames
using JuMP
using Gurobi
using ProgressMeter
using ElectronDisplay
using Dates

# Loading the project module, containing essential functions and structs
include("model/src/Sesam.jl")
using .Sesam

###################
# SETS, PARAMS ####
###################

# Case
case = "ieee_rts_24_sesam_2019"

# Time frame
t_winter = collect(1:168)
t_spring = collect(3025:3192)
t_summer = collect(5593:5760)
t_autumn = collect(7225:7392)

time = t_summer	# Select time frame

# Read nodes and load from .csv
nodes_df = CSV.read(string("data/", case, "/nodes.csv"))
load_df = CSV.read(string("data/", case, "/load.csv"))

# Create set of nodes
nodes = Nodes(nodes_df, load_df)

# Read power plants from .csv
powerplants_df = CSV.read(string("data/", case, "/powerplants.csv"))

# Read fuel costs from .csv
fuelcost_df = CSV.read(string("data/", case, "/fuelcost.csv"))

# Create set of power plants
pp = PowerPlants(powerplants_df, fuelcost_df)

# Read renewables and availabilities from .csv
renewables_df = CSV.read(string("data/", case, "/renewables.csv"))
avail_solar_df = CSV.read(string("data/", case, "/avail_solar.csv"))
avail_wind_df = CSV.read(string("data/", case, "/avail_wind.csv"))

# Create set of renewables
res = Renewables(renewables_df, avail_solar_df, avail_wind_df)

# Read lines from .csv
lines_df = CSV.read(string("data/", case, "/lines.csv"))
lines = Lines(lines_df)

### With renewables
###################
# Run DCOPF #######
###################

DCOPF_mod, DCOPF_P_opt, DCOPF_P_R_opt, DCOPF_Θ_opt, DCOPF_P_inj_opt,
	DCOPF_P_flow_opt = DCOPF(time, nodes, lines, pp, res)

# Optimisation status
JuMP.termination_status(DCOPF_mod)
JuMP.primal_status(DCOPF_mod)
JuMP.dual_status(DCOPF_mod)

# Results
JuMP.objective_value(DCOPF_mod)
DCOPF_P_opt
DCOPF_P_R_opt
DCOPF_Θ_opt
DCOPF_P_inj_opt
DCOPF_P_flow_opt

###################
# Run ED ##########
###################

ED_mod, ED_P_opt, ED_P_R_opt, ED_P_R_opt_dist, ED_price = EconomicDispatch(time, nodes, pp, res)

# Optimisation status
JuMP.termination_status(ED_mod)
JuMP.primal_status(ED_mod)
JuMP.dual_status(ED_mod)

# Results
JuMP.objective_value(ED_mod)
ED_P_opt
ED_P_R_opt
ED_P_R_opt_dist
ED_price

###################
# Run CM ##########
###################

CM_mod, CM_ΔP_up_opt, CM_ΔP_dn_opt, CM_ΔP_R_up_opt, CM_ΔP_R_dn_opt, CM_Θ_opt,
	CM_P_inj_opt, CM_P_flow_opt, CM_price =
	CongestionManagement(time, nodes, lines, ED_price, ED_P_opt, pp,
	res, ED_P_R_opt_dist)

# Optimisation status
JuMP.termination_status(CM_mod)
JuMP.primal_status(CM_mod)
JuMP.dual_status(CM_mod)

# Results
JuMP.objective_value(CM_mod)
sum(CM_ΔP_up_opt)
sum(CM_ΔP_dn_opt)
sum(CM_ΔP_R_up_opt)
sum(CM_ΔP_R_dn_opt)
CM_Θ_opt
CM_P_inj_opt
CM_P_flow_opt

###################
# Run CM + PtG ####
###################

CM_PtG_mod, CM_PtG_ΔP_up_opt, CM_PtG_ΔP_dn_opt, CM_PtG_ΔP_R_up_opt,
	CM_PtG_ΔP_R_dn_opt, CM_PtG_Θ_opt, CM_PtG_P_inj_opt, CM_PtG_P_flow_opt,
	CM_PtG_P_syn_opt, CM_PtG_D_PtG_opt, CM_PtG_L_syn_opt =
	CongestionManagementPtG(time, nodes, lines, ED_price, ED_P_opt, pp,
	res, ED_P_R_opt_dist)

# Optimisation status
JuMP.termination_status(CM_PtG_mod)
JuMP.primal_status(CM_PtG_mod)
JuMP.dual_status(CM_PtG_mod)

# Results
JuMP.objective_value(CM_PtG_mod)

sum(CM_PtG_ΔP_up_opt)
sum(CM_PtG_ΔP_dn_opt)

sum(CM_PtG_ΔP_R_up_opt)
sum(CM_PtG_ΔP_R_dn_opt)

sum(CM_PtG_P_syn_opt)
sum(CM_PtG_D_PtG_opt)
sum(CM_PtG_L_syn_opt)


text:
#nearest neighbors tour implementation

function nearest_neighbor{T<:Real}(distmat::Matrix{T}, firstcity::Int = 0)
	numCities = size(distmat, 1)
	# put first city on path
	path = Vector{Int}()
	if firstcity == 0
		firstcity = rand(1:numCities)
	end
	push!(path, firstcity)
	# cities to visit
	citiesToVisit = setdiff(1:numCities,path)
	# nearest neighbor loop
	while !isempty(citiesToVisit)
		curCity = path[end]
    dists = distmat[citiesToVisit,curCity]
		_, nextInd = findmin(dists)
		nextCity = citiesToVisit[nextInd]
		push!(path, nextCity)
		deleteat!(citiesToVisit, nextInd)
	end

	return path
end


text:
###
# Authors: Daniel Thuerck
# NEC Laboratories Europe GmbH, Copyright (c) 2021, All rights reserved.
#
# THIS HEADER MAY NOT BE EXTRACTED OR MODIFIED IN ANY WAY.
#
# SOFTWARE LICENSE AGREEMENT
# ACADEMIC OR NON-PROFIT ORGANIZATION NONCOMMERCIAL RESEARCH USE ONLY
# BY USING OR DOWNLOADING THE SOFTWARE, YOU ARE AGREEING TO THE TERMS
# OF THIS LICENSE AGREEMENT. IF YOU DO NOT AGREE WITH THESE TERMS, YOU MAY
# NOT USE OR DOWNLOAD THE SOFTWARE.
#
# This is a license agreement ("Agreement") between your academic institution or non-profit
# organization or self (called "Licensee" or "You" in this Agreement) and NEC Laboratories Europe
# GmbH (called "Licensor" in this Agreement). All rights not specifically granted to you in this
# Agreement are reserved for Licensor.
#
# RESERVATION OF OWNERSHIP AND GRANT OF LICENSE: Licensor retains exclusive
# ownership of any copy of the Software (as defined below) licensed under this Agreement and hereby
# grants to Licensee a personal, non-exclusive, non-transferable license to use the Software for
# noncommercial research purposes, without the right to sublicense, pursuant to the terms and
# conditions of this Agreement. NO EXPRESS OR IMPLIED LICENSES TO ANY OF LICENSOR’S
# PATENT RIGHTS ARE GRANTED BY THIS LICENSE. As used in this Agreement, the term
# "Software" means (i) the actual copy of all or any portion of code for program routines made
# accessible to Licensee by Licensor pursuant to this Agreement, inclusive of backups, updates, and/or
# merged copies permitted hereunder or subsequently supplied by Licensor, including all or any file
# structures, programming instructions, user interfaces and screen formats and sequences as well as any
# and all documentation and instructions related to it, and (ii) all or any derivatives and/or modifications
# created or made by You to any of the items specified in (i).
#
# CONFIDENTIALITY/PUBLICATIONS: Licensee acknowledges that the Software is proprietary to
# Licensor, and as such, Licensee agrees to receive all such materials and to use the Software only in
# accordance with the terms of this Agreement. Licensee agrees to use reasonable effort to protect the
# Software from unauthorized use, reproduction, distribution, or publication. All publication materials
# mentioning features or use of this software must explicitly include an acknowledgement the software
# was developed by NEC Laboratories Europe GmbH.
#
# COPYRIGHT: The Software is owned by Licensor.
#
# PERMITTED USES: The Software may be used for your own noncommercial internal research
# purposes. You understand and agree that Licensor is not obligated to implement any suggestionsand/or
# feedback you might provide regarding the Software, but to the extent Licensor does so, you are
# not entitled to any compensation related thereto.
#
# DERIVATIVES: You may create derivatives of or make modifications to the Software, however, You
# agree that all and any such derivatives and modifications will be owned by Licensor and become a part
# of the Software licensed to You under this Agreement. You may only use such derivatives and
# modifications for your own noncommercial internal research purposes, and you may not otherwise
# use, distribute or copy such derivatives and modifications in violation of this Agreement.
#
# BACKUPS: If Licensee is an organization, it may make that number of copies of the Software
# necessary for internal noncommercial use at a single site within its organization provided that all
# information appearing in or on the original labels, including the copyright and trademark notices are
# copied onto the labels of the copies.
#
# USES NOT PERMITTED: You may not distribute, copy or use the Software except as explicitly
# permitted herein. Licensee has not been granted any trademark license as part of this Agreement.
# Neither the name of NEC Laboratories Europe GmbH nor the names of its contributors may be used to
# endorse or promote products derived from this Software without specific prior written permission.
# You may not sell, rent, lease, sublicense, lend, time-share or transfer, in whole or in part, or provide
# third parties access to prior or present versions (or any parts thereof) of the Software.
#
# ASSIGNMENT: You may not assign this Agreement or your rights hereunder without the prior
# written consent of Licensor. Any attempted assignment without such consent shall be null and void.
#
# TERM: The term of the license granted by this Agreement is from Licensee's acceptance of this
# Agreement by downloading the Software or by using the Software until terminated as provided below.
# The Agreement automatically terminates without notice if you fail to comply with any provision of
# this Agreement. Licensee may terminate this Agreement by ceasing using the Software. Upon any
# termination of this Agreement, Licensee will delete any and all copies of the Software. You agree that
# all provisions which operate to protect the proprietary rights of Licensor shall remain in force should
# breach occur and that the obligation of confidentiality described in this Agreement is binding in
# perpetuity and, as such, survives the term of the Agreement.
#
# FEE: Provided Licensee abides completely by the terms and conditions of this Agreement, there is no
# fee due to Licensor for Licensee's use of the Software in accordance with this Agreement.
#
# DISCLAIMER OF WARRANTIES: THE SOFTWARE IS PROVIDED "AS-IS" WITHOUT
# WARRANTY OF ANY KIND INCLUDING ANY WARRANTIES OF PERFORMANCE OR
# MERCHANTABILITY OR FITNESS FOR A PARTICULAR USE OR PURPOSE OR OF NON-
# INFRINGEMENT. LICENSEE BEARS ALL RISK RELATING TO QUALITY AND
# PERFORMANCE OF THE SOFTWARE AND RELATED MATERIALS.
#
# SUPPORT AND MAINTENANCE: No Software support or training by the Licensor is provided as
# part of this Agreement.
#
# EXCLUSIVE REMEDY AND LIMITATION OF LIABILITY: To the maximum extent permitted
# under applicable law, Licensor shall not be liable for direct, indirect, special, incidental, or
# consequential damages or lost profits related to Licensee's use of and/or inability to use the Software,
# even if Licensor is advised of the possibility of such damage.
#
# EXPORT REGULATION: Licensee agrees to comply with any and all applicable export control laws,
# regulations, and/or other laws related to embargoes and sanction programs administered by law.
#
# SEVERABILITY: If any provision(s) of this Agreement shall be held to be invalid, illegal, or
# unenforceable by a court or other tribunal of competent jurisdiction, the validity, legality and
# enforceability of the remaining provisions shall not in any way be affected or impaired thereby.
#
# NO IMPLIED WAIVERS: No failure or delay by Licensor in enforcing any right or remedy under this
# Agreement shall be construed as a waiver of any future or other exercise of such right or remedy by
# Licensor.
#
# GOVERNING LAW: This Agreement shall be construed and enforced in accordance with the laws of
# Germany without reference to conflict of laws principles. You consent to the personal jurisdiction of
# the courts of this country and waive their rights to venue outside of Germany.
#
# ENTIRE AGREEMENT AND AMENDMENTS: This Agreement constitutes the sole and entire
# agreement between Licensee and Licensor as to the matter set forth herein and supersedes any
# previous agreements, understandings, and arrangements between the parties relating hereto.
#
# THIS HEADER MAY NOT BE EXTRACTED OR MODIFIED IN ANY WAY.
###

include("typestructs.jl")

# register composite types, implicit and explicit
mutable struct type_definition_registry
    defs :: Dict
    hash_name_map
    found_order
    anonymous_counter

    function type_definition_registry()
        new(Dict(), Dict(), [], 1)
    end
end

function import_struct_def(type_reg :: type_definition_registry, id::identifier_t)

    # create type from name
    name_type = eval(id.fqname)

    # create code for struct definition
    member_expr = []
    for (fname, ftype) in zip(fieldnames(name_type), name_type.types)
        ptype_expr = Meta.parse(string(ftype))
        push!(
            member_expr,
            Expr(
                Symbol("::"),
                fname,
                ptype_expr
            )
        )
    end

    struct_expr = Expr(
        Symbol("struct"),
        true,
        id.fqname,
        Expr(
            Symbol("block"),
            member_expr...
        )
    )

    # dump(struct_expr)
    p_struct = param_from_expr!(type_reg, struct_expr, id)
end

function prepend_id!(id_s)
    # only modify short name
    id_s.name = Symbol("jl_" * String(id_s.name))
end

# parsing from expr
function param_from_expr!(type_reg :: type_definition_registry, expr, p_id)

    # parse type
    p_name = nothing
    p_type = nothing
    p_managed = false
    p_composite = false

    if isa(expr, Symbol)
        sym_string = string(expr)

        if expr in supported_continuous_type_symbols

            # primitive type
            p_type = continuous_t(
                parse(Int, sym_string[6:length(sym_string)]),
                false
            )

        elseif expr in supported_int_type_symbols
            # primitive type

            if sym_string == "Bool"
                p_type = integer_t(1, true)
            elseif occursin("UInt", sym_string)
                p_type = integer_t(parse(Int, sym_string[5:length(sym_string)]), true)
            else
                p_type = integer_t(parse(Int, sym_string[4:length(sym_string)]), false)
            end

        elseif expr in [:Any, :Nothing]
            p_type = void_t()

        elseif expr in keys(type_reg.hash_name_map)
            # must be a reference to a custom (but explicit) struct whose
            # definition we must already have encountered
            p_managed = false
            p_type = struct_ref_t(type_reg.hash_name_map[expr])
        end
    else
        # not a symbol, so must be an expression (pointer or array / vector)
        if expr.head == Symbol("curly")

            if expr.args[1] == Symbol("Ptr")
                pointee_param = param_from_expr!(type_reg, expr.args[2], nothing)

                p_type = pointer_t(pointee_param.type)

            elseif expr.args[1] == Symbol("Vector")
                scalar_param = param_from_expr!(type_reg, expr.args[2], nothing)

                p_type = array_t(scalar_param.type, 1)
                p_composite = true
                p_managed = true

            elseif expr.args[1] == Symbol("Matrix")
                scalar_param = param_from_expr!(type_reg, expr.args[2], nothing)

                p_type = array_t(scalar_param.type, 2)
                p_composite = true
                p_managed = true

            elseif expr.args[1] == Symbol("Array")
                scalar_param = param_from_expr!(type_reg, expr.args[2], nothing)
                num_dims = expr.args[3]

                p_type = array_t(scalar_param.type, num_dims)
                p_composite = true
                p_managed = true

            elseif expr.args[1] == Symbol("Tuple")
                members = []

                for (ix, m_expr) in enumerate(expr.args)
                    if ix > 1
                        m_param = param_from_expr!(type_reg, m_expr, nothing)
                        push!(members, m_param)
                    end
                end

                p_type = tuple_t(members)
                p_composite = true
                p_managed = true

            end

        elseif expr.head == Symbol(".")

            # seems to be a name of some sorts -> without any other type,
            # must be a custom struct

            p_managed = false
            struct_id = parse_id(expr)

            if !(struct_id.name in keys(type_reg.hash_name_map))
                # need to retrieve this struct from the runtime
                import_struct_def(type_reg, struct_id)
            end

            p_id = type_reg.defs[type_reg.hash_name_map[struct_id.name]].name
            p_type = struct_ref_t(p_id)

        elseif expr.head == Symbol("::")
            # named member of struct
            p_name = expr.args[1]
            named_param = param_from_expr!(type_reg, expr.args[2], nothing)
            p_type = named_param.type
            p_managed = named_param.managed

        elseif expr.head == Symbol("block")
            # struct definition
            members = []

            for (ix, m_expr) in enumerate(expr.args)
                m_param = param_from_expr!(type_reg, m_expr, nothing)
                push!(members, m_param)
            end

            p_type = struct_def(p_id, members)
            p_composite = true
            p_managed = true

        elseif expr.head == Symbol("struct")
            # reference to struct
            struct_id = parse_id(expr.args[2])
            struct_name = struct_id.name

            if !(struct_name in keys(type_reg.hash_name_map))
                pparam = param_from_expr!(type_reg, expr.args[3], struct_id)
            end

            p_managed = false
            p_id = type_reg.defs[type_reg.hash_name_map[struct_name]].name
            p_type = struct_ref_t(p_id)
        end
    end

    if p_type !== nothing
        if p_composite
            hash_name = Symbol((p_id === nothing) ? hash_type(p_type) : p_id.name)

            if !(hash_name in keys(type_reg.hash_name_map))
                alias_name = hash_name
                if p_id === nothing
                    alias_name = Symbol("jl_implicit" * string(type_reg.anonymous_counter))
                    type_reg.anonymous_counter += 1
                else
                    # custom struct -> prepend "jl_"
                    alias_name = Symbol("jl_" * string(alias_name))

                    prepend_id!(p_id)
                end

                # register this structure
                # println("NEW: ", hash_name, " -> ", alias_name, " -> ", p_id)

                type_reg.defs[alias_name] = p_type
                type_reg.hash_name_map[hash_name] = alias_name
                push!(type_reg.found_order, alias_name)

            end
        end

        return param(p_name, p_type, p_managed, nothing)
    end

    error("Unsupported type:", expr)
end

function create_composite_struct_def!(type_reg :: type_definition_registry, composite_tpl)
    composite_param = param_from_expr!(type_reg, composite_tpl[2],
        (composite_tpl[1] !== nothing) ? composite_tpl[1] : nothing)

    return composite_param
end

function embed_replacements!(type_reg :: type_definition_registry)
    for key in type_reg.found_order
        def = type_reg.defs[key]

        # println(def)

        # process membersc
        if typeof(def) in [struct_def, tuple_t]
            for ix in 1:length(def.members)
                if def.members[ix].managed
                    hash_name = Symbol(hash_type(def.members[ix]))
                    new_param_name = type_reg.hash_name_map[hash_name]

                    def.members[ix].managed = true
                    def.members[ix].replacement = struct_ref_t(
                        identifier_t(new_param_name, [new_param_name], new_param_name))

                elseif isa(def.members[ix].type, struct_ref_t)
                    def.members[ix].replacement = def.members[ix].type
                end
            end
        end

        # generate new structure definitions for implicit types (replace
        # by structure defs)
        if typeof(def) == array_t
            replacement = struct_def(
                identifier_t(key, [key], key),
                [
                    param(
                        Symbol("vals"),
                        ptr_t(def.pointee_t),
                        false,
                        nothing)
                ])

            for ix in 1:def.dims
                push!(
                    replacement.members,
                    param(
                        Symbol("dim" * string(ix)),
                        integer_t(64, true),
                        false,
                        nothing
                    ))
            end

            type_reg.defs[key] = replacement

        elseif typeof(def) == tuple_t
            replacement = struct_def(
                identifier_t(key, [key], key),
                [])

            for (ix, member) in enumerate(def.members)
                push!(
                    replacement.members,
                    param(
                        Symbol("elem" * string(ix)),
                        member.type,
                        member.managed,
                        member.replacement
                    )
                )
            end

            type_reg.defs[key] = replacement
        end
    end
end

function replace_function_param!(fun_param, type_reg :: type_definition_registry)

    # use a pointer for composite types
    if isa(fun_param.type, composite_t)
        composite_hash = Symbol(hash_type(fun_param.type))

        if composite_hash in keys(type_reg.hash_name_map)
            hname = type_reg.hash_name_map[composite_hash]
            fun_param.replacement = struct_ref_t(
                identifier_t(hname, [hname], hname))

        elseif composite_hash in keys(type_reg.defs)
            # covers user-define structs
            fun_param.replacement = fun_param.type
        end
    end
end

max_stars_repo_name: nec-research/jCall

max_issues_repo_name: nec-research/jCall

max_forks_repo_name: nec-research/jCall


text:
module EmpiricalRisks

using ArrayViews

import Base.LinAlg: BlasReal
import Base.LinAlg.BLAS: axpy!, gemv!, gemm!

export

    ## common.jl
    no_op,
    shrink,

    ## loss.jl

    Loss,
    UnivariateLoss,
    MultivariateLoss,

    AbsLoss,
    SqrLoss,
    QuantileLoss,
    EpsilonInsLoss,
    HuberLoss,
    HingeLoss,
    SqrHingeLoss,
    SmoothedHingeLoss,
    LogisticLoss,
    SumLoss,
    SumSqrLoss,
    MultiLogisticLoss,

    value,
    value!,
    deriv,
    grad,
    grad!,
    value_and_deriv,
    value_and_grad,
    value_and_grad!,

    ## prediction.jl

    PredictionModel,
    UnivariatePredictionModel,
    MultivariatePredictionModel,
    LinearPred,
    MvLinearPred,
    AffinePred,
    MvAffinePred,

    inputlen,
    inputsize,
    outputlen,
    outputsize,
    paramlen,
    paramsize,
    ninputs,
    isvalidparam,
    predict,
    predict!,

    # risks

    RiskModel,
    SupervisedRiskModel,
    riskmodel,
    value_and_addgrad!,

    # regularizers

    Regularizer,
    ZeroReg,
    SqrL2Reg,
    L1Reg,
    ElasticReg,
    NonNegReg,
    SimplexReg,
    L1BallReg,

    prox, prox!

# source files
include("common.jl")
include("prediction.jl")
include("loss.jl")
include("risks.jl")
include("regularizer.jl")

end # module


text:
############################################################################################
"""
$(SIGNATURES)
Sort `val` elements at last column according to `buffer` indices, using `ancestors` as buffer.
Resort all previous columns accordingly.

# Examples
```julia
```

"""
function shuffle!(
    val::AbstractMatrix{P}, buffer::AbstractArray{P}, ancestor::AbstractVector{I}
) where {P,I<:Integer}
    ## Assign temporary variables
    Nparticles = size(ancestor, 1)
    Ndata = size(val, 2)
    ArgCheck.@argcheck Nparticles == size(val, 1) == length(buffer)
    ## Reshuffle val - can be done inplace via buffer
    @inbounds for t in Ndata:-1:1
        for idx in Base.OneTo(Nparticles)
            #!NOTE: Ancestor is only a buffer, so we do not need change order of it after shuffling.
            buffer[idx] = val[ancestor[idx], t]
        end
        for idx in Base.OneTo(Nparticles)
            val[idx, t] = buffer[idx]
        end
    end
    return nothing
end

############################################################################################
"""
$(SIGNATURES)
Sort `val` elements from `iterₘₐₓ` up until 'lookback' according to `ancestor` indices at `iterₘₐₓ`.
`ancestor` at `iterₘₐₓ - lookback' will be adjusted accordingly.

# Examples
```julia
```

"""
function shuffle_forward!(
    val::AbstractMatrix{P},
    buffer_val::AbstractArray{P},
    ancestor::AbstractMatrix{I},
    buffer_ancestor::AbstractVector{I},
    lookback::Integer,
    iterₘₐₓ::Integer,
) where {P,I<:Integer}
    ## Assign temporary variables and check for sizes
    Nparticles = size(buffer_val, 1)
    ArgCheck.@argcheck Nparticles == size(val, 1) == size(ancestor, 1)
    ArgCheck.@argcheck 0 < iterₘₐₓ <= size(val, 2) && size(val, 2) == size(ancestor, 2)
    ArgCheck.@argcheck lookback >= 0

    ## 1 Resample current particles
    @inbounds for idx in Base.OneTo(Nparticles)
        #!NOTE: buffer_ancestor assigned for higher order memory case below
        buffer_ancestor[idx] = ancestor[idx, iterₘₐₓ]
        buffer_val[idx] = val[buffer_ancestor[idx], iterₘₐₓ]
    end
    ## Assign resampled particle
    @inbounds for idx in Base.OneTo(Nparticles)
        val[idx, iterₘₐₓ] = buffer_val[idx]
    end
    ## 2 If higher order memory, resample relevant previous particles, and adjust ancestors.
    if lookback > 1
        ## Assign last possible time point for lookback,
        #!NOTE: i.e. lookback == 2 means current particles AND particles at iterₘₐₓ-1 are resampled.
        #!NOTE: Minimum iter₀ is at index 1
        iter₀ = max(1, iterₘₐₓ - lookback + 1)
        @inbounds for t in (iterₘₐₓ - 1):-1:iter₀
            ## Assign buffer for resampled particle from ancestor at iterₘₐₓ
            for idx in Base.OneTo(Nparticles)
                buffer_val[idx] = val[buffer_ancestor[idx], t]
            end
            ## Assign resampled particle and set current ancestor as non-resampled
            for idx in Base.OneTo(Nparticles)
                ## All but the last ancestors are set back to default order 1,2,3,4,..
                ancestor[idx, t] = idx
                val[idx, t] = buffer_val[idx]
            end
        end
        ## last ancestor in lookback inherits the ancestors of the current iteration
        @inbounds for idx in Base.OneTo(Nparticles)
            ancestor[idx, iter₀] = ancestor[idx, iterₘₐₓ]
            #!NOTE: Set back ancestors at iterₘₐₓ to default order as particles before already resampled
            ancestor[idx, iterₘₐₓ] = idx
        end
    end
    return nothing
end

"""
$(SIGNATURES)
Similar to `suffle`, sort `val` elements at last column according to `buffer` indices, using `ancestors` as buffer.
`buffer` indices will be different for each column, and will be adjusted accordingly.

# Examples
```julia
```

"""
function shuffle_backward!(
    val::AbstractMatrix{P},
    buffer_val::AbstractArray{P},
    ancestor::AbstractMatrix{I},
    buffer_ancestor::AbstractVector{I},
) where {P,I<:Integer}
    ## Assign relevant looping variables
    Nparticles = size(buffer_val, 1)
    iter₀ = 1
    iterₘₐₓ = size(ancestor, 2)
    ## 1 Assign Ancestors theoretical iterₘₐₓ+1 = [1,2,3,...]
    @inbounds for iter in Base.OneTo(Nparticles)
        buffer_ancestor[iter] = iter
    end
    ## 2 Resample particles in correct order
    @inbounds for t in (iterₘₐₓ - 1):-1:iter₀
        ## Reassign current ancestor
        for idx in Base.OneTo(Nparticles)
            buffer_ancestor[idx] = ancestor[buffer_ancestor[idx], t + 1]
            buffer_val[idx] = val[buffer_ancestor[idx], t]
        end
        ## Assign resampled particle and set current ancestor as non-resampled
        for idx in Base.OneTo(Nparticles)
            #!NOTE: If this line is used, we can recall shuffle_backward!() as many times as we want and particles are in correct order
            ancestor[idx, t + 1] = idx
            val[idx, t] = buffer_val[idx]
        end
    end
    return nothing
end

############################################################################################
#export
export shuffle!, shuffle_forward!, shuffle_backward!


text:
#!/usr/bin/env julia

using WriteVTK
import Compat.UTF8String
typealias FloatType Float32
const vtk_filename_noext = "collection"

function main()
    # Define grid.
    const Ni, Nj, Nk, Nt = 20, 30, 40, 4

    x = zeros(FloatType, Ni)
    y = zeros(FloatType, Nj)
    z = zeros(FloatType, Nk)

    [x[i] = i*i/Ni/Ni for i = 1:Ni]
    [y[j] = sqrt(j/Nj) for j = 1:Nj]
    [z[k] = k/Nk for k = 1:Nk]

    # Create some scalar and vectorial data that decays exponentially in
    # time
    p = zeros(FloatType, Ni, Nj, Nk, Nt)
    q = zeros(FloatType, Ni, Nj, Nk, Nt)
    vec = zeros(FloatType, 3, Ni, Nj, Nk, Nt)

    for t = 1:Nt, k = 1:Nk, j = 1:Nj, i = 1:Ni
        p[i, j, k, t] = exp(-t)*(i*i + k)
        q[i, j, k, t] = exp(-t)*k*sqrt(j)
        vec[1, i, j, k] = exp(-t)*i
        vec[2, i, j, k] = exp(-t)*j
        vec[3, i, j, k] = exp(-t)*k
    end

    # Create some scalar data at grid cells.
    # Note that in structured grids, the cells are the hexahedra formed between
    # grid points.
    cdata = zeros(FloatType, Ni-1, Nj-1, Nk-1, Nt)
    for t = 1:Nt, k = 1:Nk-1, j = 1:Nj-1, i = 1:Ni-1
        cdata[i, j, k,t] = exp(-t)*(2i + 3k * sin(3*pi * (j-1) / (Nj-2)))
    end

    # Test extents (this is optional!!)
    ext = [0, Ni-1, 0, Nj-1, 0, Nk-1] + 42

    # Initialise pvd container file
    pvd = paraview_collection(vtk_filename_noext)

    # Create files for each time-step and add them to the collection
    vtk = []
    for it = 0:Nt-1
      # vtk = vtk_grid(string(vtk_filename_noext,@sprintf("_%02i",it)),x,y,z;extent=ext)
      vtk = vtk_grid(@sprintf("%s_%02i", vtk_filename_noext, it), x, y, z;
                     extent=ext)
      # Add data for current time-step
      vtk_point_data(vtk, p[:,:,:,it+1], "p_values")
      vtk_point_data(vtk, q[:,:,:,it+1], "q_values")
      vtk_point_data(vtk, vec[:,:,:,:,it+1], "myVector")
      vtk_cell_data(vtk, cdata[:,:,:,it+1], "myCellData")
      vtk_save(vtk)
      collection_add_timestep(pvd, vtk, Float64(it+1))
    end

    # Save and close vtk file.
    outfiles = vtk_save(pvd)
    println("Saved:", [" "^3 * s for s in outfiles]...)

    return outfiles::Vector{UTF8String}
end

main()


text:
# This file was generated, do not modify it. # hide
rank([W V])

max_stars_repo_name: cscherrer/InformativePrior

max_issues_repo_name: cscherrer/InformativePrior

max_forks_repo_name: cscherrer/InformativePrior


text:
function SparseFieldMatrix(arg0::Field)
    return SparseFieldMatrix((Field,), arg0)
end

function SparseFieldMatrix(arg0::Field, arg1::jint, arg2::jint)
    return SparseFieldMatrix((Field, jint, jint), arg0, arg1, arg2)
end

function SparseFieldMatrix(arg0::FieldMatrix)
    return SparseFieldMatrix((FieldMatrix,), arg0)
end

function SparseFieldMatrix(arg0::SparseFieldMatrix)
    return SparseFieldMatrix((SparseFieldMatrix,), arg0)
end

function add_to_entry(obj::SparseFieldMatrix, arg0::jint, arg1::jint, arg2::FieldElement)
    return jcall(obj, "addToEntry", void, (jint, jint, FieldElement), arg0, arg1, arg2)
end

function copy(obj::SparseFieldMatrix)
    return jcall(obj, "copy", FieldMatrix, ())
end

function create_matrix(obj::SparseFieldMatrix, arg0::jint, arg1::jint)
    return jcall(obj, "createMatrix", FieldMatrix, (jint, jint), arg0, arg1)
end

function get_column_dimension(obj::SparseFieldMatrix)
    return jcall(obj, "getColumnDimension", jint, ())
end

function get_entry(obj::SparseFieldMatrix, arg0::jint, arg1::jint)
    return jcall(obj, "getEntry", FieldElement, (jint, jint), arg0, arg1)
end

function get_row_dimension(obj::SparseFieldMatrix)
    return jcall(obj, "getRowDimension", jint, ())
end

function multiply_entry(obj::SparseFieldMatrix, arg0::jint, arg1::jint, arg2::FieldElement)
    return jcall(obj, "multiplyEntry", void, (jint, jint, FieldElement), arg0, arg1, arg2)
end

function multiply_transposed(obj::SparseFieldMatrix, arg0::FieldMatrix)
    return jcall(obj, "multiplyTransposed", FieldMatrix, (FieldMatrix,), arg0)
end

function set_entry(obj::SparseFieldMatrix, arg0::jint, arg1::jint, arg2::FieldElement)
    return jcall(obj, "setEntry", void, (jint, jint, FieldElement), arg0, arg1, arg2)
end

function transpose_multiply(obj::SparseFieldMatrix, arg0::FieldMatrix)
    return jcall(obj, "transposeMultiply", FieldMatrix, (FieldMatrix,), arg0)
end


text:
# Wrapper for fitting Count Data Models

"""
    fit(y, model, MLEControl, printResults = true, initiate = "first")
Function to fit model to count time series.

* `y`: Time series
* `model`: Model specification
* `MLEControl`: Estimation settings
* `printResults`: Print results to console?
* `initiate`: How is conditional mean recursion initiated? "first", "marginal" or "intercept"

If MLEControl is not put in, default settings are used. Caution: This might lead to
poor initial values.
"""

# Function if MLEControl is given
import StatsBase.fit

function fit(y::Array{T, 1} where T<:Integer,
    model::T where T<:INGARCH,
    MLEControl::MLEControl;
    printResults::Bool = true,
    initiate::String = "first")

    T = length(y)

    if typeof(model) == INGARCHModel
        q = length(model.pastMean)
        Q = maximum(model.pastMean)
    else
        q = 0
        Q = 0
    end

    if typeof(model) != IIDModel
        p = length(model.pastObs)
        P = maximum(model.pastObs)
    else
        p = 0
        P = 0
    end

    r = length(model.external)
    M = maximum([P, Q])

    logl = (model.link == "Log")
    lin = !logl

    nb = model.distr == "NegativeBinomial"

    zi = model.zi

    nObs = T - M
    nPar = 1 + p + q + r + nb + zi

    optimFun = BFGS()
    if MLEControl.optimizer == "LBFGS"
        optimFun = LBFGS()
    end

    if MLEControl.optimizer == "NelderMead"
        optimFun = NelderMead()
    end

    if ll(y, model, MLEControl.init, initiate = initiate)[1] == -Inf
        error("Infinite likelihood for initial parameters. Provide other values.")
    end

    initVec = par2θ(MLEControl.init, model)
    optRaw = optimize(vars -> -ll(y, model, vars, initiate = initiate)[1],
                      initVec,
                      optimFun,
                      Optim.Options(iterations =  MLEControl.maxEval))
    cvg = Optim.converged(optRaw)
    estsVec = optRaw.minimizer
    ests = θ2par(estsVec, model)
    λ = LinPred(y, model, ests)

    temp = ll(y, model, ests, initiate = initiate)
    LLmax = temp[1]
    LLs = temp[2]

    resi = y .- λ

    se = zeros(nPar)
    CI = zeros(2, nPar)

    ciProb = false
    if nb
        if ests.ϕ[1] > 10000
            ciProb = true
            println("No confidence intervals computed due to large estimate of ϕ.")
        end
    end

    if MLEControl.ci & !ciProb
        H = Calculus.hessian(vars -> ll(y, model, vars, initiate = initiate)[1], estsVec)
        se = sqrt.(-diag(inv(H)))
        quan = quantile(Normal(), 0.975)
        for i = 1:nPar
            CI[:, i] = estsVec[i] .+ [-1, 1].*se[i].*quan
        end
    end

    res = INGARCHresults(y, estsVec, ests, λ, resi, LLmax, LLs, nPar, nObs, se, CI, model, cvg, MLEControl)

    if printResults
        show(res)
    end

    return res
end

function fit(y::Array{T, 1} where T<:Integer,
    model::T where T<:INGARCH;
    printResults::Bool = true,
    initiate::String = "first")

    MLEControl = MLESettings(y, model)

    fit(y, model, MLEControl, printResults = printResults, initiate = initiate)
end

function fit(y::Array{T, 1} where T<:Integer,
    model::T where T<:INARMA,
    MLEControl::MLEControl;
    printResults::Bool = true)

    T = length(y)

    if typeof(model) in [INARMAModel, INMAModel]
        q = length(model.pastMean)
        Q = maximum(model.pastMean)

        if Q > 2
            error("INARMA(p, q) / INMA(q) with q > 2 not yet supported.")
        end
    else
        q = 0
        Q = 0
    end

    if typeof(model) != INMAModel
        p = length(model.pastObs)
        P = maximum(model.pastObs)
    else
        p = 0
        P = 0
    end

    r = length(model.external)
    M = maximum([P, Q])

    logl1 = (model.link[1] == "Log")
    logl2 = (model.link[2] == "Log")

    nb1 = model.distr[1] == "NegativeBinomial"
    nb2 = model.distr[2] == "NegativeBinomial"
    if r == 0
        nb2 = false
    else
        if sum(model.external) == 0
            nb2 = false
        end
    end

    zi = model.zi

    nObs = T - M
    nPar = 1 + p + q + r + nb1 + nb2 + zi

    optimFun = BFGS()
    if MLEControl.optimizer == "LBFGS"
        optimFun = LBFGS()
    end

    if MLEControl.optimizer == "NelderMead"
        optimFun = NelderMead()
    end

    if ll(y, model, MLEControl.init)[1] == -Inf
        error("Infinite likelihood for initial parameters. Provide other values.")
    end

    initVec = par2θ(MLEControl.init, model)

    optRaw = optimize(vars -> -ll(y, model, vars)[1],
                      initVec,
                      optimFun,
                      Optim.Options(iterations =  MLEControl.maxEval))

    cvg = Optim.converged(optRaw)
    estsVec = optRaw.minimizer
    ests = θ2par(estsVec, model)

    temp = ll(y, model, ests)
    LLmax = temp[1]
    LLs = temp[2]

    se = zeros(nPar)
    CI = zeros(2, nPar)

    ciProb = false
    if nb1 | nb2
        if any(ests.ϕ .> 10000)
            ciProb = true
            println("No confidence intervals computed due to large estimate of ϕ.")
        end
    end

    if MLEControl.ci & !ciProb
        H = Calculus.hessian(vars -> ll(y, model, vars)[1], estsVec)
        se = sqrt.(-diag(inv(H)))
        quan = quantile(Normal(), 0.975)
        for i = 1:nPar
            CI[:, i] = estsVec[i] .+ [-1, 1].*se[i].*quan
        end
    end

    res = INARMAresults(y, estsVec, ests, LLmax, LLs, nPar, nObs, se, CI, model, cvg, MLEControl)

    if printResults
        show(res)
    end

    return res
end

function fit(y::Array{T, 1} where T<:Integer,
    model::T where T<:INARMA;
    printResults::Bool = true)

    MLEControl = MLESettings(y, model)

    fit(y, model, MLEControl, printResults = printResults)
end


text:
export dense,
       dims,
       dim,
       mindim,
       diaglength

# dim and dims are used in the Tensor interface, overload
# base Dims here
dims(ds::Dims) = ds

dims(::Tuple{}) = ()

dense(ds::Dims) = ds

dense(::Type{DimsT}) where {DimsT<:Dims} = DimsT

dim(ds::Dims) = prod(ds)

dim(ds::Dims,i::Int) = dims(ds)[i]

mindim(inds::Tuple) = minimum(dims(inds))

diaglength(inds::Tuple) = mindim(inds)

"""
strides(ds::Dims)

Get the strides of the dimensions.

This is unexported, call with NDTensors.strides.
"""
strides(ds::Dims) = Base.size_to_strides(1, dims(ds)...)

"""
stride(ds::Dims, k::Int)

"""
stride(ds::Dims, k::Int) = strides(ds)[k]

# This is to help with some generic programming in the Tensor
# code (it helps to construct a Tuple(::NTuple{N,Int}) where the
# only known thing for dispatch is a concrete type such
# as Dims{4})
similar_type(::Type{<:Dims},
             ::Type{Val{N}}) where {N} = Dims{N}

# This is to help with ITensor compatibility
dir(::Int) = 0

# This is to help with ITensor compatibility
dag(i::Int) = i


text:
# model parameters
A0 = 0.03
Lambda = 2000f0                 # brightness falls with exp(-0.5*Lambda*r^2)
# empirical lambda_x = 1f-6, lambda_a = 1f-2
Tau = Float32[0,0,1f-2,1f-4]  # how much x,y,z,a change between timesteps

N = 3                           # number of neurons
X = 512                         # image width, assuming square X,X images
Z = 20                          # number of layers
ZX = 1.5                        # one z layer is ZX pixels in distance
Z0 = 17                         # Z0+1:Z0+Z is one block
Z1 = Z0+div(Z,2)                # middle of the block
X1drift = 2.5                   # image drift in pixels for first dim
X2drift = 1.0                   # image drift in pixels for second dim

# we map 1:X to 0:1 in the xy directions.
X1 = Array(Float32,X,X); for i=1:X,j=1:X; X1[i,j]=i/X; end
X2 = Array(Float32,X,X); for i=1:X,j=1:X; X2[i,j]=j/X; end

# we map 1:Z to 0:ZX*Z/X in the z direction. some conversions:
# p: pixels, i: images, z: z levels, x: x units
p2x(p::Int)=Float32(p/X)
x2p{T}(x::T)=round(Int,x*X)
z2x(z::Int)=Float32(z*ZX/X)
x2z{T}(x::T)=round(Int,x*X/ZX)
i2z(i::Int)=mod1(i-Z0,Z)
i2x(i::Int)=z2x(i2z(i))

# alloc some temporary arrays
M = [ Array(Float32,X,X) for i=1:20 ]
W = [ Array(Float32,4,N) for i=1:10 ]

function zforw{T}(par::Matrix{T},out::Matrix{T}=M[1]; img=Z1,
                  m1::Matrix{T}=M[2], m2::Matrix{T}=M[3])
    fill!(out, A0)                              # A0 is the background brightness
    X3 = i2x(img)                             # X1,X2,X3 represent pixel positions
    for k=1:size(par,2)
        (x1,x2,x3,a) = par[:,k]                 # x1,x2,x3 represent center position, a is activation
        broadcast!(-, m1, X1, x1)
        broadcast!(*, m1, m1, m1)               # (X1-x1)^2
        broadcast!(-, m2, X2, x2)
        broadcast!(*, m2, m2, m2)               # (X2-x2)^2
        broadcast!(+, m2, m1, m2)               # +(X1-x1)^2
        broadcast!(+, m2, m2, (x3-X3)^2) 	# +(X3-x3)^2
        broadcast!(*, m2, m2, -Lambda/2)
        broadcast!(exp, m2, m2)                 # exp(-Lambda/2 * d2)
        broadcast!(*, m2, m2, a)                # ai * exp(-Lambda/2 * d2)
        broadcast!(+, out, out, m2)             # a0 + sum ai * exp(-Lambda/2 * d2)
    end
    return out
end

function zloss{T}(curr::Matrix{T},prev::Matrix{T},pred::Matrix{T},gold::Matrix{T}; img=Z1,
                  xdiff::Matrix{T}=W[1], ydiff::Matrix{T}=M[10])
    broadcast!(-, xdiff, curr, prev)
    broadcast!(*, xdiff, xdiff, xdiff)
    broadcast!(*, xdiff, xdiff, Tau)
    broadcast!(-, ydiff, pred, gold)
    broadcast!(*, ydiff, ydiff, ydiff)
    return (mean(ydiff) + sum(xdiff))/2
end

# Backward gradient:
function zback{T}(curr::Matrix{T},prev::Matrix{T},pred::Matrix{T},gold::Matrix{T},diff::Matrix{T}=M[4]; img=Z1,
                  delta::Matrix{T}=M[5], dx1::Matrix{T}=M[6], dx2::Matrix{T}=M[7], m1::Matrix{T}=M[8], m2::Matrix{T}=M[9])
    X3 = i2x(img)
    broadcast!(-, diff, curr, prev)
    broadcast!(*, diff, diff, Tau)
    broadcast!(-, delta, pred, gold)
    for k=1:size(curr,2)
        (x1,x2,x3,a) = curr[:,k]
        broadcast!(-, dx1, X1, x1)
        broadcast!(-, dx2, X2, x2)
        dx3 = X3 - x3
        broadcast!(*, m1, dx1, dx1)
        broadcast!(*, m2, dx2, dx2)
        broadcast!(+, m2, m2, m1)
        broadcast!(+, m2, m2, dx3*dx3)
        broadcast!(*, m2, m2, -Lambda/2)
        broadcast!(exp, m2, m2)
        broadcast!(*, m2, m2, delta)
        diff[4,k] += mean(m2)                 # dJ/da = sum delta*exp(-Lambda*dist2)
        broadcast!(*, m2, m2, a*Lambda)
        diff[1,k] += mean(broadcast!(*, dx1, dx1, m2))
        diff[2,k] += mean(broadcast!(*, dx2, dx2, m2))
        diff[3,k] += mean(broadcast!(*, m2, m2, dx3))
    end
    return diff
end

# Gradient check:
function zcheck{T}(curr::Matrix{T},prev::Matrix{T},gold::Matrix{T}; img=Z1, eps=T(1e-4),
                   pred::Matrix{T}=M[11], diff::Matrix{T}=W[2])
    zforw(curr, pred, img=img)
    zback(curr, prev, pred, gold, diff, img=img)
    for i=1:length(curr)
        xi = curr[i]
        curr[i] = xi - eps
        zforw(curr, pred, img=img)
        f1 = zloss(curr, prev, pred, gold)
        curr[i] = xi + eps
        zforw(curr, pred, img=img)
        f2 = zloss(curr, prev, pred, gold)
        curr[i] = xi
        df = (f2-f1)/(2*eps)
        println((i,df,diff[i]))
    end
end

# Training:
function ztrain{T}(curr::Matrix{T}, prev::Matrix{T}, gold::Matrix{T}; img=Z1, iter=50, gmin=1f-5, gmax=1f-4, lr=0f0,
                   pred::Matrix{T}=M[12], disp::Matrix{T}=M[13], diff::Matrix{T}=W[3], prn=false)
    zdraw(curr, gold, img=img)
    for i=1:iter
        zforw(curr, pred, img=img)
        zback(curr, prev, pred, gold, diff, img=img)
        gnorm = vecnorm(diff)
        gnorm > gmax && scale!(gmax/gnorm, diff)
        lr == 0 && (@show lr = 5f-3/gnorm)
        Base.LinAlg.axpy!(-lr, diff, curr)
        zdraw(curr, gold, img=img)
        vecnorm(diff) < gmin && break
        prn && println((:img,img,:loss,zloss(curr, prev, zforw(curr,pred,img=img), gold, img=img),:gnorm,vecnorm(diff),:x,round(Int64,X*curr[:])))
    end
    loss = zloss(curr, prev, zforw(curr,pred,img=img), gold, img=img)
    println((img,i2z(img),xx2pp(curr)))
    curr
end

function xx2pp{T}(x::Matrix{T})
    a = cell(length(x))
    for i=1:4:length(x)
        a[i] = x2p(x[i])
        a[i+1] = x2p(x[i+1])
        a[i+2] = x2z(x[i+2])
        a[i+3] = round(Int,x[i+3]*0xff)
    end
    return a
end

# This loses the centers going from one block to the next.
# So we need to fix the drift using zfix.
function ztrack{T}(y::Array{T}, i1=Z1, i2=size(y,3); lr=100)
    fixed = zfix(y[:,:,i1], i1)
    prev = zinit(fixed, img=i1)
    curr = copy(prev)
    for i=i1:i2
        fixed = zfix(y[:,:,i], i)
        ztrain(curr, prev, fixed; img=i, lr=lr)
        copy!(prev, curr)
    end
end

function zfix{T}(img::Array{T,2}, i::Int; dx=X1drift, dy=X2drift, out::Array{T,2}=M[15])
    copy!(out,img)
    X,Y = size(img)
    z = i2z(i)
    x = round(Int,dx * z); y = round(Int,dy * z)
    ox1 = max(1,1-x); ox2 = min(X,X-x); oy1 = max(1,1-y); oy2 = min(Y,Y-y)
    ix1 = max(1,1+x); ix2 = min(X,X+x); iy1 = max(1,1+y); iy2 = min(Y,Y+y)
    copy!(sub(out,ox1:ox2,oy1:oy2), sub(img,ix1:ix2,iy1:iy2))
    return out
end

# Tracking
# This version starts from the middle of each block, goes back and forth
function ztrack2{T}(y::Array{T}, img1=1, img2=size(y,3);
                    curr::Matrix{T}=W[6], prev::Matrix{T}=W[7])
    while img1 < img2
        i1 = img1; while mod1(i1-Z0,Z) > 1; i1-=1; end
        i3 = img1; while mod1(i3-Z0,Z) < Z; i3+=1; end
        i2 = i1 + 10
        cent = zinit(y[:,:,i2], img=i2)
        copy!(prev, cent); copy!(curr, cent)
        for i=i2:-1:max(1,i1)
            ztrain(curr, prev, y[:,:,i], img=i)
            copy!(prev, curr)
            i==i2 && copy!(cent, curr)
        end
        copy!(prev, cent); copy!(curr, cent)
        for i=i2+1:i3
            ztrain(curr, prev, y[:,:,i], img=i)
            copy!(prev, curr)
        end
        img1 = i3+1
    end
end

function zinit{T}(y::Matrix{T};img=Z1)
    x = Array(T,4,N)
    c = zcenters(y)
    for n=1:N
        x[1,n] = c[n][1]/size(y,1)
        x[2,n] = c[n][2]/size(y,2)
        x[3,n] = i2x(img)
        x[4,n] = mean(y[c[n][1]-2:c[n][1]+2,c[n][2]-2:c[n][2]+2])-A0
    end
    return x
end

function zcenters{T}(x::Matrix{T}; r=30)
    c = Any[]
    for i=1:size(x,1)
        for j=1:size(x,2)
            zlocalmax(x,i,j,r) && push!(c, (i,j,x[i,j]))
        end
    end
    return sort(c, by=(x->x[3]), rev=true)
end

function zlocalmax{T}(x::Matrix{T}, i, j, r)
    for ii=max(1,i-r):min(size(x,1),i+r)
        for jj=max(1,j-r):min(size(x,2),j+r)
            if x[i,j] < x[ii,jj]
                return false
            end
        end
    end
    return true
end

using ImageView

function zview(y;cmax=maximum(y))
    global _zview
    try
        view(_zview,y)
    catch e
        (_zview,_) = view(y,clim=[0,cmax],name="zview",xy=["y","x"])
    end
end

using Images, Colors
using ImageView: annotate!, AnnotationText

function zdraw{T}(curr::Matrix{T}, gold::Matrix{T}; img=Z1, temp::Matrix{T}=M[14])
    global _zdraw
    copy!(temp, gold)
    X3 = i2x(img)                             # X1,X2,X3 represent pixel positions
    # f = a0 + ai * exp(-lambda/2 * d2)
    dmin = -log(0.3)/(Lambda/2)
    dmax = 1f0 # -log(0.05)/(Lambda/2)
    for k=1:size(curr,2)
        (x1,x2,x3,a) = curr[:,k]                 # x1,x2,x3 represent center position, a is activation
        X1a = max(1/X, x1-20/X); X1b = min(1, x1+20/X)
        X2a = max(1/X, x2-20/X); X2b = min(1, x2+20/X)
        for X1=X1a:1/X:X1b, X2=X2a:1/X:X2b
            d2 = (x1-X1)^2 + (x2-X2)^2 + (x3-X3)^2
            dmin <= d2 <= dmax && (temp[x2p(X1),x2p(X2)]=0.0)
        end
    end
    annot = sprint() do s
        println(s, "$img ($(i2z(img)))")
        writedlm(s, reshape(append!(Any["x","y","z","a"], xx2pp(curr)),(4,4)))
    end
    try
        a = first(values(_zdraw[1].annotations))
        a.data.string = annot
        view(_zdraw[1],temp)
    catch e
        _zdraw = view(temp, clim=[0,maximum(temp)], name="zdraw", xy=["y","x"])
        annotate!(_zdraw[1], _zdraw[2], AnnotationText(10, 10, annot, halign="left", valign="top", color=RGB(1,1,1)))
        view(_zdraw[1],temp)
    end
end

max_stars_repo_name: JuliaTagBot/NeuronTracking

max_issues_repo_name: JuliaTagBot/NeuronTracking

max_forks_repo_name: JuliaTagBot/NeuronTracking


text:
export ARC_Sec

function ARC_Sec(h :: AbstractNLPModel,
					  nlpstop :: AbstractStopping;
				 	  kwargs...)

	OK, stop = ARC_generic(h, nlpstop; direction = :Sec, kwargs...)
	return OK, stop
end


text:
"""
```
function scalar_reduce(args...)
```
Each individual iteration returns n scalars. The output is reduced to n vectors,
where the i-th vector contains all of the i-th scalars from each iteration.

The return type of reduce functions must be the same type as the tuple of
arguments passed in. If args is a tuple of Vector{Float64}, then the return
argument will be a Vector{Float64}.

e.g.
a, b = @parallel (scalar_reduce) for i in 1:10000
           [[1], [2]]
       end
a = [1, 1, 1, ...]
b = [2, 2, 2, ...]

Input/Output type: Vector{Vector{Float64}}
"""
function scalar_reduce(args...)
    return_arg = args[1]
    for (i, arg) in enumerate(args[2:end])
        for (j, el) in enumerate(arg)
            append!(return_arg[j], el)
        end
    end
    return return_arg
end

"""
```
function vector_reduce(args...)
```
Each individual iteration returns n Vector types; we vector-reduce to n matrices, where
the i-th column of that matrix corresponds to the i-th vector from an individual iteration.

Input/Output type: Vector{Matrix{Float64}}
"""
function vector_reduce(args...)
    nargs1 = length(args)    # The number of times the loop is run
    nargs2 = length(args[1]) # The number of variables output by a single run

    return_arg = args[1]
    for i in 1:nargs2
        for j in 2:nargs1
            return_arg[i] = hcat(return_arg[i], args[j][i])
        end
    end
    return return_arg
end

"""
```
function scalar_reshape(args...)
```
Function ensures type conformity of the return arguments.
"""
function scalar_reshape(args...)
    n_args = length(args)
    return_arg = Vector{Vector{Float64}}(undef, n_args)
    for i in 1:n_args
        arg = typeof(args[i]) <: Vector ? args[i] : [args[i]]
        return_arg[i] = arg
    end
    return return_arg
end

"""
```
function vector_reshape(args...)
```
Function ensures type conformity of the return arguments.
"""
function vector_reshape(args...)
    n_args = length(args)
    return_arg = Vector{Matrix{Float64}}(undef, n_args)
    for i in 1:n_args
        arg = typeof(args[i]) <: Vector ? args[i] : [args[i]]
        return_arg[i] = reshape(arg, length(arg), 1)
    end
    return return_arg
end

"""
```
sendto(p::Int; args...)
```
Function to send data from master process to particular worker, p.
Code from ChrisRackauckas, avavailable at:
 https://github.com/ChrisRackauckas/ParallelDataTransfer.jl/blob/master/src/ParallelDataTransfer.jl.
"""
function sendto(p::Int; args...)
    for (nm, val) in args
        @spawnat(p, Core.eval(Main, Expr(:(=), nm, val)))
    end
end

"""
```
sendto(ps::AbstractVector{Int}; args...)
```
Function to send data from master process to list of workers.
Code from ChrisRackauckas, available at:
https://github.com/ChrisRackauckas/ParallelDataTransfer.jl/blob/master/src/ParallelDataTransfer.jl.
"""
function sendto(ps::AbstractVector{Int}; args...)
    for p in ps
        sendto(p; args...)
    end
end

function get_cloud(filepath::String)
    return load(filepath, "cloud")
end

function init_stage_print(cloud::Cloud, para_symbols::Vector{Symbol};
                          verbose::Symbol=:low, use_fixed_schedule::Bool = true)
    if VERBOSITY[verbose] >= VERBOSITY[:low]
        if use_fixed_schedule
            println("--------------------------")
            println("Iteration = $(cloud.stage_index) / $(cloud.n_Φ)")
        else
            println("--------------------------")
            println("Iteration = $(cloud.stage_index)")
        end
	    println("--------------------------")
        println("phi = $(cloud.tempering_schedule[cloud.stage_index])")
	    println("--------------------------")
        println("c = $(cloud.c)")
        println("ESS = $(cloud.ESS[cloud.stage_index])   ($(cloud.resamples) total resamples.)")
	    println("--------------------------")
    end
    if VERBOSITY[verbose] >= VERBOSITY[:high]
        μ = weighted_mean(cloud)
        σ = weighted_std(cloud)
        for n=1:length(para_symbols)
            println("$(para_symbols[n]) = $(round(μ[n], digits = 5)), $(round(σ[n], digits = 5))")
	    end
    end
end

function end_stage_print(cloud::Cloud, para_symbols::Vector{Symbol};
                         verbose::Symbol=:low, use_fixed_schedule::Bool = true)
    if VERBOSITY[verbose] >= VERBOSITY[:low]
        total_sampling_time_minutes = cloud.total_sampling_time/60
        if use_fixed_schedule
            expected_time_remaining_sec = (cloud.total_sampling_time/cloud.stage_index) *
                (cloud.n_Φ - cloud.stage_index)
            expected_time_remaining_minutes = expected_time_remaining_sec / 60
        end

        println("--------------------------")
        if use_fixed_schedule
            println("Iteration = $(cloud.stage_index) / $(cloud.n_Φ)")
            println("time elapsed: $(round(total_sampling_time_minutes, digits = 4)) minutes")
            println("estimated time remaining: " *
                    "$(round(expected_time_remaining_minutes, digits = 4)) minutes")
        else
            println("Iteration = $(cloud.stage_index)")
            println("time elapsed: $(round(total_sampling_time_minutes, digits = 4)) minutes")
        end
        println("--------------------------")
        println("phi = $(cloud.tempering_schedule[cloud.stage_index])")
        println("--------------------------")
        println("c = $(cloud.c)")
        println("accept = $(cloud.accept)")
        println("ESS = $(cloud.ESS[cloud.stage_index])   ($(cloud.resamples) total resamples.)")
        println("--------------------------")
    end
    if VERBOSITY[verbose] >= VERBOSITY[:high]
        μ = weighted_mean(cloud)
        σ = weighted_std(cloud)
        for n=1:length(para_symbols)
            println("$(para_symbols[n]) = $(round(μ[n], digits = 5)), $(round(σ[n], digits = 5))")
        end
    end
end

"""
Sparse identity matrix - since deprecated in 0.7
"""
function speye(n::Integer)
    return SparseMatrixCSC{Float64}(I, n, n)
end

"""
Sparse identity matrix - since deprecated in 0.7
"""
function speye(T::Type, n::Integer)
    return SparseMatrixCSC{T}(I, n, n)
end

"""
    <(a::Complex, b::Complex)

Compare real values of complex numbers.
"""
function <(a::Complex, b::Complex)
    return a.re < b.re
end

"""
    <(a::Real, b::Complex)

Compare real values of complex numbers.
"""
function <(a::Real, b::Complex)
    return a < b.re
end

"""
    <(a::Complex, b::Real)

Compare real values of complex numbers.
"""
function <(a::Complex, b::Real)
    return a.re < b
end

function min(a::Complex, b::Real)
    return min(a.re, b)
end

function min(a::Complex, b::Complex)
    return min(a.re, b.re)
end

function min(a::Real, b::Complex)
    return min(a, b.re)
end

function max(a::Complex, b::Real)
    return max(a.re, b)
end

function max(a::Complex, b::Complex)
    return max(a.re, b.re)
end

function max(a::Real, b::Complex)
    return max(a, b.re)
end

function isempty(c::Cloud)
    isempty(c.particles)
end


text:
abstract type AbstractWakeDeflectionModel end

"""
    GaussYawDeflection(horizontal_spread_rate, vertical_spread_rate, alpha_star, beta_star)

Container for parameters related to the Gaussian deflection model presented by Bastankhah and Porte-Agel 2016

# Arguments
- `horizontal_spread_rate::Float`: parameter controlling the horizontal spread of the deficit model. Default value is 0.022.
- `vertical_spread_rate::Float`: parameter controlling the vertical spread of the deficit model. Default value is 0.022.
- `alpha_star::Float`: parameter controlling the impact of turbulence intensity on the length of the near wake. Default value is 2.32.
- `beta_star::Float`: parameter controlling the impact of the thrust coefficient on the length of the near wake. Default value is 0.154.
"""
struct GaussYawDeflection{TF} <: AbstractWakeDeflectionModel
    horizontal_spread_rate::TF
    vertical_spread_rate::TF
    alpha_star::TF
    beta_star::TF
end
GaussYawDeflection() = GaussYawDeflection(0.022, 0.022, 2.32, 0.154)

"""
    GaussYawDeflectionVariableSpread(alpha_star, beta_star, k1, k2, wec_factor)

Container for parameters related to the Gaussian deflection model with yaw presented by Bastankhah and Porte-Agel 2016

# Arguments
- `alpha_star::Float`: parameter controlling the impact of turbulence intensity on the length of the near wake. Default value is 2.32.
- `beta_star::Float`: parameter controlling the impact of the thrust coefficient on the length of the near wake. Default value is 0.154.
- `k1::Float`: first parameter tuning wake spread as based on turbulence intensity
- `k2::Float`: second parameter tuning wake spread as based on turbulence intensity
"""
struct GaussYawVariableSpreadDeflection{TF} <: AbstractWakeDeflectionModel
    alpha_star::TF
    beta_star::TF
    k1::TF
    k2::TF
end
GaussYawVariableSpreadDeflection() = GaussYawVariableSpreadDeflection(2.32, 0.154, 0.3837, 0.003678)
GaussYawVariableSpreadDeflection(x, y) = GaussYawVariableSpreadDeflection(x, y, 0.3837, 0.003678)

"""
    JiminezYawDeflection(horizontal_spread_rate)

Container for parameters related to the Jiminez deflection model

# Arguments
- `horizontal_spread_rate::Float`: parameter controlling the wake spreading rate and deficit decay. Default value is 0.1
"""
struct JiminezYawDeflection{TF} <: AbstractWakeDeflectionModel
    horizontal_spread_rate::TF
end
JiminezYawDeflection() = JiminezYawDeflection(0.1)

"""
    MultizoneDeflection(horizontal_spread_rate)

Container for parameters related to the Jiminez deflection model

# Arguments
- `horizontal_spread_rate::Float`: parameter controlling the wake spreading rate and deficit decay. Default value is 0.1
- `ad::Float`:Helps define the horizontal deflection of the wake at 0 deg yaw
- `bd::Float`:Helps define the horizontal deflection of the wake due to downwind distance at 0 deg yaw
"""
struct MultizoneDeflection{TF} <: AbstractWakeDeflectionModel
    horizontal_spread_rate::TF
    ad::TF
    bd::TF
end
MultizoneDeflection() = MultizoneDeflection(0.15, -4.5, -0.01)

"""
    wake_deflection_model(locx, locy, locz, turbine_id, turbine_definition::TurbineDefinition, model::JiminezYawDeflection, windfarmstate::SingleWindFarmState)

    Calculates the horizontal deflection of the wind turbine wake

    Based on:
    [1] Jiminez 2010 "Wake defl ection of a wind turbine in yaw"
    [2] Gebraad 2014 "Wind plant optimization by yaw control using a parametric wake model"
    this version ignores the corrections made to the yaw model for rotor rotation as described in [2] and
    [3] Thomas 2017 "Improving the FLORIS wind plant model for compatibility with gradient-based optimization"
"""
function wake_deflection_model(locx, locy, locz, turbine_x, turbine_yaw, turbine_ct, turbine_id, rotor_diameter, turbine_local_ti, model::JiminezYawDeflection)

    dx = locx-turbine_x[turbine_id]
    yaw = -turbine_yaw[turbine_id] # Jiminez used opposite rotation convention, hence (-) sign
    ct = turbine_ct[turbine_id]
    diam = rotor_diameter[turbine_id]

    kd = model.horizontal_spread_rate

    initial_yaw_angle = 0.5*((cos(yaw))^2)*sin(yaw)*ct  # [1] eq. 20, [2] eq. 8

    # [2] eq. 10
    a = 2.0*kd*dx/diam + 1.0
    b = initial_yaw_angle*(15.0*a^4+initial_yaw_angle^2)
    c = (30.0*kd/diam)*a^5
    d = initial_yaw_angle*diam*(15.0 + initial_yaw_angle^2)
    e = 30.0*kd

    y_deflection = b/c - d/e

    return y_deflection

end

"""
    wake_deflection_model(locx, locy, locz, turbine_id, turbine_definition::TurbineDefinition, model::MultizoneDeflection, windfarmstate::SingleWindFarmState)

    Calculates the horizontal deflection of the wind turbine wake accounting for both yaw and rotational deflection

    Based on:
    [1] Jiminez 2010 "Wake defl ection of a wind turbine in yaw"
    [2] Gebraad 2014 "Wind plant optimization by yaw control using a parametric wake model"
    this version ignores the corrections made to the yaw model for rotor rotation as described in [2] and
    [3] Thomas 2017 "Improving the FLORIS wind plant model for compatibility with gradient-based optimization"
"""

function wake_deflection_model(locx, locy, locz, turbine_x, turbine_yaw, turbine_ct, turbine_id, rotor_diameter, turbine_local_ti, model::MultizoneDeflection)

    dx = locx-turbine_x[turbine_id]
    yaw = -turbine_yaw[turbine_id] # Jiminez used opposite rotation convention, hence (-) sign
    ct = turbine_ct[turbine_id]
    diam = rotor_diameter[turbine_id]

    kd = model.horizontal_spread_rate
    ad = model.ad
    bd = model.bd

    initial_yaw_angle = 0.5*((cos(yaw))^2)*sin(yaw)*ct  # [1] eq. 20, [2] eq. 8

    # [2] eq. 10
    a = 2.0*kd*dx/diam + 1.0
    b = initial_yaw_angle*(15.0*a^4+initial_yaw_angle^2)
    c = (30.0*kd/diam)*a^5
    d = initial_yaw_angle*diam*(15.0 + initial_yaw_angle^2)
    e = 30.0*kd

    # [2] eq. 10, 11, and 12 define deflection
    yaw_deflection = b/c - d/e
    rotation_deflection = ad + bd*(dx)
    y_deflection = -1*(yaw_deflection + rotation_deflection)

    return y_deflection

end

function _bpa_theta_0(yaw, ct)

    theta0 = (0.3*yaw/cos(yaw))*(1.0-sqrt(1.0-ct*cos(yaw)))

    return theta0
end

function _bpa_deflection(diam, ct, yaw, ky, kz, sigmay, sigmaz, theta0, x0)
    a = theta0*x0/diam
    b = (theta0/14.7)*sqrt(cos(yaw)/(ky*kz*ct))*(2.9-1.3*sqrt(1.0-ct)-ct)
    c = (1.6+sqrt(ct))*(1.6*sqrt(8.0*sigmay*sigmaz/(cos(yaw)*diam^2))-ct)
    d = (1.6-sqrt(ct))*(1.6*sqrt(8.0*sigmay*sigmaz/(cos(yaw)*diam^2))+ct)
    y_deflection = diam*(a+b*log(c/d))
    return y_deflection
end

"""
    wake_deflection_model(locx, locy, locz, turbine_x, turbine_yaw, turbine_ct, turbine_id, rotor_diameter, turbine_local_ti, model::GaussYawDeflection)

    Calculates the horizontal deflection of the wind turbine wake

    Based on:
    [1] Bastankhah and Porte-Agel 2016 "Experimental and theoretical study of
    wind turbine wakes in yawed conditions"
"""
function wake_deflection_model(locx, locy, locz, turbine_x, turbine_yaw, turbine_ct, turbine_id, rotor_diameter, turbine_local_ti, model::GaussYawDeflection)

    dx = locx-turbine_x[turbine_id]
    yaw = turbine_yaw[turbine_id]
    ct = turbine_ct[turbine_id]
    diam = rotor_diameter[turbine_id]
    ti = turbine_local_ti[turbine_id]

    as = model.alpha_star
    bs = model.beta_star
    ky = model.horizontal_spread_rate
    kz = model.vertical_spread_rate

    # [1] eqn 6.12
    theta0 = _bpa_theta_0(yaw, ct)

    # [1] eqn 7.4
    x0 = _gauss_yaw_potential_core(diam, yaw, ct, as, ti, bs)

    # calculate the discontinuity point of the gauss yaw model
    xd = _gauss_yaw_discontinuity(diam, x0, ky, kz, yaw, ct)

    # calculate horizontal wake spread (paper eq: 7.2)
    sigmay = _gauss_yaw_spread_interpolated(diam, ky, dx, x0, yaw, xd)

    # calculate vertical wake spread (paper eq: 7.2)
    sigmaz = _gauss_yaw_spread_interpolated(diam, kz, dx, x0, 0.0, xd)

    y_deflection = _bpa_deflection(diam, ct, yaw, ky, kz, sigmay, sigmaz, theta0, x0)

    return y_deflection
end

"""
    wake_deflection_model(oc, turbine_x, turbine_yaw, turbine_ct, turbine_id, rotor_diameter, turbine_local_ti, model::GaussYawVariableSpreadDeflection)

    Calculates the horizontal deflection of the wind turbine wake. Varies based on local turbulence intensity.

    Based on:
    [1] Bastankhah and Porte-Agel 2016 "Experimental and theoretical study of
    wind turbine wakes in yawed conditions"
    [2] Niayifar and Porte-Agel 2016 "Analytical Modeling of Wind Farms:
    A New Approach for Power Prediction"
"""
function wake_deflection_model(locx, locy, locz, turbine_x, turbine_yaw, turbine_ct, turbine_id, rotor_diameter, turbine_local_ti, model::GaussYawVariableSpreadDeflection)

    dx = locx-turbine_x[turbine_id]
    yaw = turbine_yaw[turbine_id]
    ct = turbine_ct[turbine_id]
    dt = rotor_diameter[turbine_id]
    ti = turbine_local_ti[turbine_id]

    as = model.alpha_star
    bs = model.beta_star

    # [2] calculate wake spread based on local turbulence intensity
    ky = kz = _k_star_func(ti, model.k1, model.k2)

    # [1] eqn 6.12 initial wake angle
    theta0 = _bpa_theta_0(yaw, ct)

    # [1] eqn 7.4
    x0 = _gauss_yaw_potential_core(dt, yaw, ct, as, ti, bs)

    # calculate the discontinuity point of the gauss yaw model
    xd = _gauss_yaw_discontinuity(dt, x0, ky, kz, yaw, ct)

    # calculate horizontal wake spread (paper eq: 7.2)
    sigma_y = _gauss_yaw_spread_interpolated(dt, ky, dx, x0, yaw, xd)

    # calculate vertical wake spread (paper eq: 7.2)
    sigma_z = _gauss_yaw_spread_interpolated(dt, kz, dx, x0, 0.0, xd)

    # finally, calculate deflection
    y_deflection = _bpa_deflection(dt, ct, yaw, ky, kz, sigma_y, sigma_z, theta0, x0)

    return y_deflection
end


text:
using OpenCV
using ClTools
using Base.Test
using ImageCore
using MicroLogging

limit_logging(MicroLogging.Info);

# Visualization.
img = ones(UInt8, 100, 100);
ClTools.visualization.apply_colormap(img; vmin=0, vmax=18, cmap=22);
ClTools.visualization.apply_colormap(img; vmin=0, vmax=18, cmap=:viridis);
img = ones(Float32, 100, 100);
ClTools.visualization.apply_colormap(img; vmin=0, vmax=18, cmap=22);
ClTools.visualization.apply_colormap(img; vmin=0, vmax=18, cmap=:viridis);
duration = @elapsed ClTools.visualization.apply_colormap(img;
                                                         vmin=0, vmax=18,
                                                         cmap=22);
@info "Visualizing a 100x100 image took $(duration)s."

# Configuration.
ncpus = ClTools.config.available_cpu_count()
@info "This machine has $(ncpus) CPUs."
npus = ClTools.config.available_pu_count()
@info "This machine has $(npus) PUs."


text:
using Tar, CodecZlib, GoogleCodeSearch, Sockets, JSON

# A simple search server for demo.
#
# - Fetches source code for Julia packages from github releases
# - Extracts them
# - Indexes them using GoogleCodeSearch.jl
# - Serves a REST API that provides results using the index
#
# To simplify things and highlight the important parts, we use
# - only certain pre-downloaded sources, from the file system
# - a single index, no incremental re-indexing
# - the search server provided by GoogleCodeSearch.jl

const paths = (
    inputs = joinpath(workdir, "inputs"),           # the pre-fetched source bundles
    downloaded = joinpath(workdir, "downloaded"),   # just copy those here and pretend we fetched from github
    extracted = joinpath(workdir, "extracted"),     # extract the tar bundles here
    index = joinpath(workdir, "index"),             # and create the index here
)

function fetch_sources()
    mkpath(paths.downloaded)                        # create working path if it does not exist
    for name in readdir(paths.downloaded)           # clean up any old downloaded sources
        rm(joinpath(paths.downloaded, name))
    end
    for name in readdir(paths.inputs)               # fetch all inputs into downloaded folder
        src = joinpath(paths.inputs, name)
        dest = joinpath(paths.downloaded, name)
        cp(src, dest)
    end
end

function extract_sources()
    mkpath(paths.extracted)                         # create working path if it does not exist
    for name in readdir(paths.extracted)            # clean up any old extracted sources
        rm(joinpath(paths.extracted, name), recursive=true)
    end
    for name in readdir(paths.downloaded)
        src = joinpath(paths.downloaded, name)
        dest = joinpath(paths.extracted, name)

        tar_gz = open(src)
        tar = GzipDecompressorStream(tar_gz)
        Tar.extract(tar, dest)                      # extract all downloaded bundles
    end
end

function index_sources()
    ctx = Ctx(store=paths.index)                    # initialize code search
    for name in readdir(paths.index)                # clean up the old index
        rm(joinpath(paths.index, name))
    end
    for path in readdir(paths.extracted)            # create a new index
        index(ctx, joinpath(paths.extracted, path))
    end
end

function run_search_server()
    ctx = Ctx(store=paths.index)                    # initialize context pointing to index
    run_http(ctx; host=ip"0.0.0.0", port=5555)      # start the search server
end

max_stars_repo_name: tanmaykm/KubernetesIntegratedJulia

max_issues_repo_name: tanmaykm/KubernetesIntegratedJulia

max_forks_repo_name: tanmaykm/KubernetesIntegratedJulia


text:
# code for simulate polr trait from hapmap3 data

using CSV, DataFrames, PolrModels, SnpArrays

# read in Plink fam file
plkfam = CSV.read("hapmap3.fam";
    header=["famid", "perid", "faid", "moid", "sex", "trait"],
    types=[String, String, String, String, Int, Float64],
    delim=' ',
    datarow=1)

# read in Plink bim file
hapmap = SnpArray("hapmap3")
# set SNP 2, 3, 4, 5 as causal snps with effect size βgeno
Xgeno = convert(Matrix{Float64}, hapmap[:, 2:5])
# set sex effect to be 0.5
X = [plkfam[:sex] Xgeno]

# simular ordered categorical trait
srand(123)
J = 4
θ = collect(1.0:J-1)
β = [0.5; 1.25ones(size(Xgeno, 2))]
link = LogitLink()
Y = rpolr(X, β, θ, link)

# write csv file
covdf = [plkfam[1:5] DataFrame(trait=Y)]
CSV.write("covariate.txt", covdf)

max_stars_repo_name: Hua-Zhou/PolrGWAS

max_issues_repo_name: Hua-Zhou/PolrGWAS


text:
esc_nonconstant(x::Number) = x
esc_nonconstant(x) = esc(x)
complements_error(args, str) = error("In @complements($(join(args,","))): ", str)

# for backward compatibility with 0.4
function comparison_to_call(ex)
    if isexpr(ex,:comparison) && length(ex.args) == 3
        return Expr(:call, ex.args[2], ex.args[1], ex.args[3])
    else
        return ex
    end
end

function smooth(c1, c2)
    c1s = Expr(:call, esc(:(^)), c1, 2)
    c2s = Expr(:call, esc(:(^)), c2, 2)
    cs = Expr(:call, esc(:(+)), c1s, c2s, :($mpec_tol))
    csqrt = Expr(:call, esc(:sqrt), cs)

    csum = Expr(:call, esc(:+), c1, c2)
    cc = Expr(:call, esc(:-), csqrt, csum)

    return Expr(:call, esc(:(==)), cc, 0  )
end

function get_complementarity(c1, c2, method)

    if method == :smooth
        expr = smooth(c1, c2)
    elseif method == :simple
        cc = Expr(:call, esc(:(*)), c1, c2)
        expr = Expr(:call, esc(:(<=)), cc, :($mpec_tol)  )
    else
        expr = smooth(c1, c2)
    end

    return expr

end

macro complements(args...)

    # Currently supports
    #  F ⟂ lb <= x <= ub
    #  or
    #  0 <= F ⟂ x >= 0
    #  0 >= F ⟂ x >= 0
    #  0 <= F ⟂ x <= 0
    #  0 >= F ⟂ x <= 0
    #  Note: the bound on F must be on the left hand side, and
    #        the obund on x must be on the right hand side.

    # Many parts of this macro were copied from @variable of JuMP.jl
    # Thanks!

    m = esc(args[1])
    F = args[2]
    x = args[3]

    method = :smooth
    mpec_tol = :(1e-8)
    if length(args) > 3
        method = args[4]
    end

    ub_F = Inf
    lb_F = -Inf
    ub_x = Inf
    lb_x = -Inf

    xhaslb = false
    xhasub = false
    Fhaslb = false
    Fhasub = false

    # Identify the variable bounds. Five (legal) possibilities are "x >= lb",
    # "x <= ub", "lb <= x <= ub", "x == val", or just plain "x"
    if VERSION < v"0.5.0-dev+3231"
        x = comparison_to_call(x)
        F = comparison_to_call(F)
    end

    ############################### x      #####################################
    if isexpr(x, :comparison) # two-sided

        if x.args[2] == :>= || x.args[2] == :≥
            # ub >= x >= lb
            x.args[4] == :>= || x.args[4] == :≥ || complements_error(args, "Invalid variable bounds")
            var = x.args[3]
            lb_x = esc_nonconstant(x.args[5])
            ub_x = esc_nonconstant(x.args[1])

            xhaslb = true
            xhasub = true
        elseif x.args[2] == :<= || x.args[2] == :≤
            # lb <= x <= ub
            var = x.args[3]
            (x.args[4] != :<= && x.args[4] != :≤) &&
                complements_error(args, "Expected <= operator after variable name.")

            lb_x = esc_nonconstant(x.args[1])
            ub_x = esc_nonconstant(x.args[5])

            xhaslb = true
            xhasub = true
        else
            complements_error(args, "Use the form lb <= ... <= ub.")
        end

    elseif isexpr(x, :call)
        if x.args[1] == :>= || x.args[1] == :≥
            # x >= lb
            var = x.args[2]
            @assert length(x.args) == 3
            lb_x = esc_nonconstant(x.args[3])
            ub_x = Inf
            xhaslb = true
            xhasub = false

            # # May also be ub >= x
            # if isa(eval(var), Number)
            #     var = x.args[3]
            #     @assert length(x.args) == 3
            #     ub_x = esc_nonconstant(x.args[2])
            #     lb_x = -Inf
            #     xhaslb = false
            #     xhasub = true
            # end

        elseif x.args[1] == :<= || x.args[1] == :≤
            # x <= ub
            var = x.args[2]
            @assert length(x.args) == 3
            ub_x = esc_nonconstant(x.args[3])
            lb_x = -Inf
            xhaslb = false
            xhasub = true

            # # May also be lb <= x
            # if isa(eval(var), Number)
            #     var = x.args[3]
            #     @assert length(x.args) == 3
            #     ub_x = Inf
            #     lb_x = esc_nonconstant(x.args[2])
            #     xhaslb = true
            #     xhasub = false
            # end

        else
            # # Its a comparsion, but not using <= ... <=
            # complements_error(args, "Unexpected syntax $(string(x)).")
            # No bounds provided - free variable
            # If it isn't, e.g. something odd like f(x), we'll handle later
            var = x
            lb_x = -Inf
            ub_x = Inf
        end

    else
        # No bounds provided - free variable
        # If it isn't, e.g. something odd like f(x), we'll handle later
        var = x
        lb_x = -Inf
        ub_x = Inf
    end
    ############################### x      #####################################

    ############################### F      #####################################
    if isexpr(F, :comparison) # two-sided

        if F.args[2] == :>= || F.args[2] == :≥
            # ub >= x >= lb
            F.args[4] == :>= || F.args[4] == :≥ || complements_error(args, "Invalid function bounds")
            func = F.args[3]

            lb_F = esc_nonconstant(F.args[5])
            ub_F = esc_nonconstant(F.args[1])

            Fhaslb = true
            Fhasub = true
        elseif F.args[2] == :<= || F.args[2] == :≤
            # lb <= x <= ub
            func = F.args[3]
            (F.args[4] != :<= && F.args[4] != :≤) &&
                complements_error(args, "Expected <= operator after function expression.")
            lb_F = esc_nonconstant(F.args[1])
            ub_F = esc_nonconstant(F.args[5])

            Fhaslb = true
            Fhasub = true
        else
            complements_error(args, "Use the form lb <= ... <= ub.")
        end

    elseif isexpr(F, :call)
        if F.args[1] == :>= || F.args[1] == :≥
            # x >= lb
            # func = F.args[2]
            # @assert length(F.args) == 3
            # lb_F = esc_nonconstant(F.args[3])
            # ub_F = Inf
            # Fhaslb = true

            # May also be ub >= x
            # if isa(eval(func), Number)
            func = F.args[3]
            @assert length(F.args) == 3
            ub_F = esc_nonconstant(F.args[2])
            lb_F = -Inf
            Fhaslb = false
            Fhasub = true
            # end

        elseif F.args[1] == :<= || F.args[1] == :≤
            # # x <= ub
            # func = F.args[2]
            # @assert length(F.args) == 3
            # ub_F = esc_nonconstant(F.args[3])
            # lb_F = -Inf
            # Fhasub = true

            # May also be lb <= x
            # if isa(eval(func), Number)
            func = F.args[3]
            @assert length(F.args) == 3
            ub_F = Inf
            lb_F = esc_nonconstant(F.args[2])
            Fhaslb = true
            Fhasub = false
            # end

        else
            # Its a comparsion, but not using <= ... <=
            # complements_error(args, "Unexpected syntax $(string(F)).")
            # No bounds provided - free funciable
            # If it isn't, e.g. something odd like f(x), we'll handle later
            func = F
            lb_F = -Inf
            ub_F = Inf
        end

    else
        # No bounds provided - free funciable
        # If it isn't, e.g. something odd like f(x), we'll handle later
        func = F
        lb_F = -Inf
        ub_F = Inf
    end
    ############################### F      #####################################

    number_bounds = xhaslb + xhasub + Fhaslb + Fhasub
    if number_bounds != 2
        complements_error(args, "The total number of bounds on the function and the variable must be exactly two.")
    end

    var_esc = esc(var)
    func_esc = esc(func)

    if isexpr(var, :call)
        complements_error(args, "The second argument must be a variable, not an expression.")
    end

    code = quote end

    # Bound constraints are NOT added in :smooth method.
    if method == :simple || method == :smooth
        if Fhaslb
            expr = Expr(:call, esc(:(>=)), esc(func), :($lb_F)  )
            push!(code.args, Expr(:macrocall, Symbol("@NLconstraint"), m, expr ) )
        end
        if Fhasub
            expr = Expr(:call, esc(:(<=)), esc(func), :($ub_F)  )
            push!(code.args, Expr(:macrocall, Symbol("@NLconstraint"), m, expr ) )
        end
        if xhaslb
            expr = Expr(:call, esc(:(>=)), esc(var), :($lb_x)  )
            push!(code.args, Expr(:macrocall, Symbol("@NLconstraint"), m, expr ) )
        end
        if xhasub
            expr = Expr(:call, esc(:(<=)), esc(var), :($ub_x)  )
            push!(code.args, Expr(:macrocall, Symbol("@NLconstraint"), m, expr ) )
        end
    end

    # There must be a better way of writing the codes below...
    if Fhaslb && Fhasub
        complements_error(args, "Both bounds on the function expression is currently not supported.")

    elseif xhaslb && xhasub

        # Additional variables are defined
        # v defined
        v = copy(var)
        v.args[1] = Symbol(join(["v", gensym(), var.args[1]]))
        v_expr = Expr(:call, esc(:>=), esc(v), 0)
        push!(code.args, Expr(:macrocall, Symbol("@variable"), m, v_expr))

        # w defined
        w = copy(var)
        w.args[1] = Symbol(join(["w", gensym(), var.args[1]]))
        w_expr = Expr(:call, esc(:>=), esc(w), 0)
        push!(code.args, Expr(:macrocall, Symbol("@variable"), m, w_expr))

        # v - w = func
        vwdiff = Expr(:call, esc(:(-)), esc(v), esc(w))
        # expr = Expr(:comparison, vwdiff, esc(:(==)), esc(func))
        expr = Expr(:call, esc(:(==)), vwdiff, esc(func))
        push!(code.args, Expr(:macrocall, Symbol("@constraint"), m, expr))

        # v * (x-lb) = 0
        c1 = Expr(:call, esc(:(-)), esc(var), :($lb_x))
        expr = get_complementarity(c1, esc(v), method)
        push!(code.args, Expr(:macrocall, Symbol("@NLconstraint"), m, expr))

        # w * (ub - x) = 0
        c1 = Expr(:call, esc(:(-)), :($ub_x), esc(var))
        expr = get_complementarity(c1, esc(w), method)
        push!(code.args, Expr(:macrocall, Symbol("@NLconstraint"), m, expr))

    elseif Fhaslb && xhaslb
        c1 = Expr(:call, esc(:(-)), esc(var), :($lb_x))
        c2 = Expr(:call, esc(:(-)), esc(func), :($lb_F))
        expr = get_complementarity(c1, c2, method)
        push!(code.args, Expr(:macrocall, Symbol("@NLconstraint"), m, expr ) )

    elseif Fhaslb && xhasub
        c1 = Expr(:call, esc(:(-)), :($ub_x), esc(var))
        c2 = Expr(:call, esc(:(-)), esc(func), :($lb_F))
        expr = get_complementarity(c1, c2, method)
        push!(code.args, Expr(:macrocall, Symbol("@NLconstraint"), m, expr ) )

    elseif Fhasub && xhaslb
        c1 = Expr(:call, esc(:(-)), esc(var), :($lb_x))
        c2 = Expr(:call, esc(:(-)), :($ub_F), esc(func))
        expr = get_complementarity(c1, c2, method)
        push!(code.args, Expr(:macrocall, Symbol("@NLconstraint"), m, expr ) )

    elseif Fhasub && xhasub
        c1 = Expr(:call, esc(:(-)), :($ub_x), esc(var))
        c2 = Expr(:call, esc(:(-)), :($ub_F), esc(func))
        expr = get_complementarity(c1, c2, method)
        push!(code.args, Expr(:macrocall, Symbol("@NLconstraint"), m, expr ) )

    else
        complements_error(args, "NO VALID CASE")
    end

    return code

end


text:
#
# This file is a part of MolecularGraph.jl
# Licensed under the MIT License http://opensource.org/licenses/MIT
#

@testset "graph.view.reversegraph" begin

@testset "reversegraph" begin
    graph = digraph(5, [(1, 2), (2, 3), (3, 4), (4, 5)])
    rev = reversegraph(graph)

    node = getnode(rev, 4)
    @test isa(node, AbstractNode)

    edge = getedge(rev, 3, 4)
    @test isa(edge, DirectedEdge)
    @test edge.source == 4
    @test edge.target == 3

    nodes = nodesiter(rev)
    (n, state) = iterate(nodes)
    @test in(n[1], [1, 2, 3, 4, 5])
    @test isa(n[2], AbstractNode)

    edges = edgesiter(rev)
    (e, state) = iterate(edges)
    @test in(e[1], [1, 2, 3, 4])
    @test isa(e[2], DirectedEdge)

    succs = outneighbors(rev, 4)
    (s, state) = iterate(succs)
    @test s[1] == 3
    @test s[2] == 3

    preds = inneighbors(rev, 4)
    (p, state) = iterate(preds)
    @test p[1] == 5
    @test p[2] == 4

    @test indegree(rev, 5) == 0
    @test outdegree(rev, 1) == 0

    ntype = nodetype(rev)
    @test ntype <: AbstractNode
    etype = edgetype(rev)
    @test etype <: DirectedEdge

    mapg = digraph(rev)
    @test nodecount(mapg) == 5
    @test edgecount(mapg) == 4
    @test isa(mapg, DiGraph)
end

end # graph.view.reversegraph


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license
module Libdl
# Just re-export Base.Libc.Libdl:
export DL_LOAD_PATH, RTLD_DEEPBIND, RTLD_FIRST, RTLD_GLOBAL, RTLD_LAZY, RTLD_LOCAL,
    RTLD_NODELETE, RTLD_NOLOAD, RTLD_NOW, dlclose, dlopen, dlopen_e, dlsym, dlsym_e,
    dlpath, find_library, dlext, dllist

import Base.Libc.Libdl: DL_LOAD_PATH, RTLD_DEEPBIND, RTLD_FIRST, RTLD_GLOBAL, RTLD_LAZY, RTLD_LOCAL,
                        RTLD_NODELETE, RTLD_NOLOAD, RTLD_NOW, dlclose, dlopen, dlopen_e, dlsym, dlsym_e,
                        dlpath, find_library, dlext, dllist

end # module


text:
# this file is included by wrap.jl and provides several functions for building docstrings

"Build docstring for a function from a Doxygen XML node"
function build_function(node::EzXML.Node)
    io = IOBuffer()

    # code block with function definition
    print(io, "    ", text(node, "name"), "(")
    nspace = position(io)
    params = findall("param", node)
    if isempty(params)
        print(io, ")")
    end
    for param in params
        param !== first(params) && print(io, " "^nspace) # align
        islast = param === last(params)
        print(io, text(param, "type"))
        # declname not always present (with void)
        lastchar = islast ? ")" : ",\n"
        if !isempty(findall("declname", param))
            print(io, " ", text(param, "declname"), lastchar)
        else
            print(io, lastchar)
        end
    end
    println(io, " -> ", text(node, "type"))

    # brief description, always 1 present, sometimes only whitespace
    desc = strip(nodecontent(findfirst("briefdescription", node)))
    if !isempty(desc)
        println(io, "\n", text(node, "briefdescription"))
    end

    # parameters
    params = findall("detaileddescription/para/parameterlist/parameteritem", node)
    if !isempty(params)
        println(io, "\n### Parameters")
        for param in params
            print(io, "* **", text(param, "parameternamelist"), "**: ")
            println(io, text(param, "parameterdescription"))
        end
    end

    # returns
    return_elems = findall("detaileddescription/para/simplesect[@kind='return']", node)
    if !isempty(return_elems)
        println(io, "\n### Returns")
        println(io, text(return_elems[1], "para"))
    end

    String(take!(io))
end

"Build a one line docstring consisting of only the brief description from an XML node"
function brief_description(node::EzXML.Node)
    # brief description, always 1 present, sometimes only whitespace
    strip(nodecontent(findfirst("briefdescription", node)))
end

"Compose a Markdown docstring based on a Doxygen XML element"
function build_docstring(node::EzXML.Node)
    kind = node["kind"]
    if kind == "function"
        build_function(node)
    elseif kind in ("enum", "define", "typedef")
        brief_description(node)
    else
        # "friend" and "variable" kinds remain
        # but we leave them out, not needed
        ""
    end
end

"Return the text of a subelement `el` of `node`"
function text(node::EzXML.Node, el::AbstractString)
    s = findfirst(el, node)
    s === nothing ? "" : strip(nodecontent(s))
end

"Wrap the single- or multi-line docstring in appropriate quotes"
function addquotes(docstr::AbstractString)
    if '\n' in docstr
        string("\"\"\"\n", docstr, "\"\"\"")
    else
        # single line docstring
        repr(rstrip(docstr, '.'))
    end
end

"Get the C name out of a expression"
function cname(ex)
    if @capture(ex, function f_(args__) ccall((a_, b_), xs__) end)
        String(eval(a))
    else
        # TODO make MacroTools.namify work for structs and macros
        if MacroTools.isexpr(ex) && ex.head === :struct
            String(ex.args[2])
        elseif MacroTools.isexpr(ex) && ex.head === :macrocall
            # if the enum has a type specified
            String(ex.args[3].args[1])
        else
            String(namify(ex))
        end
    end
end

"Based on a name, find the best XML node to generate docs from"
function findnode(name::String, doc::EzXML.Document)
    # Names are not unique. We know that kind='friend' (not sure what it is)
    # does not give good docs and is never the only one, so we skip those.
    # First we use XPath to find all nodes with this name and not kind='friend'.
    memberdef = "/doxygen/compounddef/sectiondef/memberdef"
    nofriend = "not(@kind='friend')"  # :-(
    nodes = findall("$memberdef[name='$name' and $nofriend]", doc)

    if length(nodes) == 0
        return nothing
    elseif length(nodes) == 1
        return first(nodes)
    else
        # If we get multiple nodes back, we have to select the best one.
        # Looking at the documentation, sometimes there are two similar docstrings,
        # but one comes from a .cpp file, as seen in location's file attribute,
        # and the other comes from a .h file (.c is the third option).
        # Even though this is a C binding, the .cpp node includes the argument names
        # which makes for an easier to read docstring, since they can be referenced
        # to the argument names in the parameters list.
        # Therefore if .cpp is one of the options, go for that.

        # ExXML uses libxml2 which only supports XPath 1.0, meaning
        # ends-with(@file,'.cpp') is not available, but according to
        # https://stackoverflow.com/a/11857166/2875964 we can rewrite this as
        cpp = "'.cpp' = substring(@file, string-length(@file) - string-length('.cpp') +1)"

        for node in nodes
            cppnode = findfirst("location[$cpp]/..", node)
            if cppnode !== nothing
                return cppnode
            end
        end
        # .cpp not present, just pick the first
        return first(nodes)
    end
end