text:
# *****************************************************************************
# Written by Ritchie Lee, ritchie.lee@sv.cmu.edu
# *****************************************************************************
# Copyright ã 2015, United States Government, as represented by the
# Administrator of the National Aeronautics and Space Administration. All
# rights reserved.  The Reinforcement Learning Encounter Simulator (RLES)
# platform is licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License. You
# may obtain a copy of the License at
# http://www.apache.org/licenses/LICENSE-2.0. Unless required by applicable
# law or agreed to in writing, software distributed under the License is
# distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied. See the License for the specific language
# governing permissions and limitations under the License.
# _____________________________________________________________________________
# Reinforcement Learning Encounter Simulator (RLES) includes the following
# third party software. The SISLES.jl package is licensed under the MIT Expat
# License: Copyright (c) 2014: Youngjun Kim.
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to
# deal in the Software without restriction, including without limitation the
# rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
# sell copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED
# "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT
# NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
# PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
# HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
# *****************************************************************************

module NMACStats

export nmac_stats

import Compat.ASCIIString

using ..SaveHelpers

const NMAC_STATS_ROUND_NDECIMALS = 2

function nmac_stats{T<:AbstractString}(infiles::Vector{T}, txtfile::AbstractString = "nmac_stats.txt")

  stats = Dict{ASCIIString,Any}()

  for file = infiles

    d = trajLoad(file)
    run_type = d["run_type"]

    if !haskey(stats, run_type)
      stats[run_type] = Dict{ASCIIString,Any}()
      stats[run_type]["nmac_count"] = Int64(0)
      stats[run_type]["nmac_encs_rewards"] = (Int64, Float64)[]
      stats[run_type]["total_count"] = Int64(0)
    end

    @show file

    if sv_nmac(d) #nmac occurred

      #increment count
      stats[run_type]["nmac_count"] += 1

      enc_id, enctype = sv_encounter_id(d)

      #add encounter to list
      if enctype != "invalid"
        push!(stats[run_type]["nmac_encs_rewards"],
              (enc_id, sv_reward(d)))
      else
        error("nmac_stats: Invalid encounter id")
      end
    end

    #increment total_count
    stats[run_type]["total_count"] += 1
  end

  #open file
  f = open(txtfile, "w")

  #sort the vectors
  for run_type = keys(stats)

    enc_rewards = stats[run_type]["nmac_encs_rewards"]
    sort!(enc_rewards, by = x -> x[2], rev = true)
    sorted_ids = [tup[1] for tup in enc_rewards]
    sorted_rewards = [tup[2] for tup in enc_rewards]

    if isopen(f)
      println(f, "run type=$run_type")
      println(f, "nmac count=$(stats[run_type]["nmac_count"])")
      println(f, "total count=$(stats[run_type]["total_count"])")
      println(f, "sorted nmac ids=$(sorted_ids)")
      println(f, "sorted nmac rewards=$(round(sorted_rewards, NMAC_STATS_ROUND_NDECIMALS))")
      println(f, "")
    end

  end

  close(f)

  return stats
end

end #module

max_stars_repo_name: sisl/RLESCAS

max_issues_repo_name: sisl/RLESCAS

max_forks_repo_name: sisl/RLESCAS


text:
module Pers

if VERSION>=VersionNumber(0,7)
   using Dates
end

export Per,+,-

struct Per

    id::String
    year::Int64
    month::Int64
    last_day::Date

end

function Per(s::String)

    @assert length(s)==7 "Date String Wrong size, correct format -> yyyy-mm"
    @assert s[5]=='-' "Date String incorrect format, correct format -> yyyy-mm"
    @assert parse(Int64,s[6:7])<=12 "Date String incorrect format, month > 12, correct format -> yyyy-mm"

    year=parse(Int64,s[1:4])
    month=parse(Int64,s[6:7])

    return Per(s,year,month,Dates.lastdayofmonth(Date(year,month,1)))

end

Per(d::Date)=Per(string(Dates.year(d))*"-"*(("0"*string(Dates.month(d)))[end-1:end]),Dates.year(d),Dates.month(d),Dates.lastdayofmonth(d))

import Base.+
import Base.-
import Base.<
import Base.>
import Base.==

Base.:<(p1::Per,p2::Per)=p1.last_day<p2.last_day
Base.:>(p1::Per,p2::Per)=p1.last_day>p2.last_day

==(p1::Per, p2::Per)=p1.last_day==p2.last_day

Base.:+(p1::Per, i::Int)=Per(p1.last_day+Dates.Month(i))
Base.:+(i::Int,p1::Per)=Per(p1.last_day+Dates.Month(i))

Base.:-(p1::Per, i::Int)=Per(p1.last_day-Dates.Month(i))
Base.:-(p1::Per, p2::Per)=(p1.year-p2.year)*12+(p1.month-p2.month)

import Base.convert

convert(::Type{Pers.Per}, str::String) =Per(str)

end # module

ext: jl

max_stars_repo_name: AntonioLoureiro/Pers

max_issues_repo_name: AntonioLoureiro/Pers

max_forks_repo_name: AntonioLoureiro/Pers


text:
module FormInput

using Genie, Stipple, StippleUI, StippleUI.API
import Genie.Renderer.Html: HTMLString, normal_element, template

export textfield, numberfield, textarea, filefield

Genie.Renderer.Html.register_normal_element("q__input", context = @__MODULE__)
Genie.Renderer.Html.register_normal_element("q__file", context = @__MODULE__)

function textfield( label::String = "",
                    fieldname::Union{Symbol,Nothing} = nothing,
                    args...;
                    wrap::Function = StippleUI.DEFAULT_WRAPPER,
                    kwargs...)
  wrap() do
    q__input(args...; attributes([:label => label, :fieldname => fieldname, kwargs...], StippleUI.API.ATTRIBUTES_MAPPINGS)...)
  end
end

function numberfield( label::String = "",
                      fieldname::Union{Symbol,Nothing} = nothing,
                      args...;
                      wrap::Function = StippleUI.DEFAULT_WRAPPER,
                      kwargs...)
  wrap() do
    q__input(args...; attributes([:label => label, :fieldname => fieldname, kwargs...],
              merge(StippleUI.API.ATTRIBUTES_MAPPINGS, Dict("fieldname" => "v-model.number")))...)
  end
end

function textarea(label::String = "",
                  fieldname::Union{Symbol,Nothing} = nothing,
                  args...;
                  wrap::Function = StippleUI.DEFAULT_WRAPPER,
                  kwargs...)

  textfield(label, fieldname, args...; type="textarea", wrap=wrap, kwargs...)
end

function filefield( label::String = "",
                    fieldname::Union{Symbol,Nothing} = nothing,
                    args...;
                    wrap::Function = StippleUI.DEFAULT_WRAPPER,
                    kwargs...)

  wrap() do
    q__file(args...; attributes(
                        [:label => label, :fieldname => fieldname, kwargs...],
                        Dict("fieldname" => "v-model",
                              "itemaligned" => "item-aligned", "stacklabel" => "stack-label", "bottomslots" => "bottom-slots",
                              "hidebottomspace" => "hide-bottom-space", "hidehint" => "hide-hint", "displayvalue" => "display-value",
                              "reactiverules" => "reactive-rules", "lazyrules" => "lazy-rules", "rules" => ":rules",
                              "counterlabel" => "counter-label", "maxfiles" => "max-files", "maxtotalsize" => "max-total-size", "maxfilesize" => "max-file-size",
                              "errormessage" => "error-message", "noerroricon" => "no-error-icon", "clearicon" => "clear-icon",
                              "inputstyle" => "input-style", "inputclass" => "input-class", "bgcolor" => "bg-color",
                              "labelslot" => "label-slot", "labelcolor" => "label-color"
                              )
                      )...)
  end
end

end


text:
## define a model
T = [-0.5 0.5;
    1 -1]
C = [1; -3]

r₁(x) = - (x.>1) + (x.<=1)
r₂(x) = x.*0
R₁(x) = (x.>1).*(2 .- x) + (x.<=1).*x

r = (
    r = function (x)
        [r₁(x) r₂(x)]
    end,
    R = function (x)
        [R₁(x) r₂(x)]
    end
)

bounds = [0 12; -Inf Inf]
model = SFFM.Model( T, C, r, Bounds = bounds)

max_stars_repo_name: angus-lewis/SFFM

max_issues_repo_name: angus-lewis/SFFM

max_forks_repo_name: angus-lewis/SFFM


text:
# %% imports
using Random
using PyCall
gym = pyimport("gym")

include("../../../Modules/Networks/net_arch.jl")
include("../../../Modules/Poisson/PoissonStateSpike.jl")
include("../../../Modules/simviz.jl")

# %% sim functions

function count_spikes(arch::Array{Int64}, spiked::Array{Int64}, layer::Int64)
    count_layer = zeros(arch[layer])
    interval = [sum(arch[1:layer-1])+1, sum(arch[1:layer])]
    spikes_in = filter(x-> interval[1] <= x <= interval[2], spiked)
    spikes_in .-= sum(arch[1:layer-1])
    count_layer[spikes_in] .+= 1
    return count_layer
end

function left_right(count_spikes::Array{Float64}, n_left::Int64, n_right::Int64)
    sl = sum(count_spikes[1:n_left])
    sr = sum(count_spikes[n_left+1:end])
    (sl >= sr) ? action = 0 : action = 1
    return action, sl, sr
end

# %%

function play_episode(net::Network, n_steps::Int64 = 200)
    env = gym.make("CartPole-v1")
    #env.seed(0)
    obs = env.reset()
    rng_action = Random.MersenneTwister(0)
    win_size = 40 # ms
    spike_hist = Array{Int64}[]
    vhist = Float64[]
    tot_reward = 0
    N_rand = 0
    @inbounds for j in 1:n_steps
        count_tot = zeros(net.arch[end])
        spikes = int_spikes(obs) # state2spike function
        @inbounds for msec in 1:win_size
            spiked = step_LIF!(net, spikes, false)
            influence = 2
            if msec >= length(net.arch)-1
                count_tot .+= influence*count_spikes(net.arch, spiked, length(net.arch))
                influence *= 0.95
                push!(spike_hist,spiked)
                push!(vhist,net.v[1])
            end
        end
        n_left = div(net.arch[end],2)
        n_right = net.arch[end] - n_left
        action, sl, sr = left_right(count_tot, n_left, n_right)
        if sl == sr
            action = rand(rng_action, 0:1)
            N_rand += 1
        end
        #push!(score_hist[1], sl)
        #push!(score_hist[2], sr)
        obs_new, reward, done, _ = env.step(action)
        tot_reward += 1
        if done
            break
        end
        obs = obs_new
    end
    println("Rand = ", N_rand, "/", tot_reward)
    env.close()
    env = nothing
    Base.GC.gc()
    return tot_reward,spike_hist,vhist
end

# %% tests
params = YAML.load(open("Julia/CartPole/LIF_base/cfglif4820.yml"))
n = Network([4,8,20], params)
r, shist,vhist = play_episode(n)

# %% viz
spike_histogram(n.arch, shist)
#%%
weight_viz(n,"histogram")
#%%

plot(vhist, legend=false)

max_stars_repo_name: SuReLI/RL-STDP

max_issues_repo_name: SuReLI/RL-STDP

max_forks_repo_name: SuReLI/RL-STDP


text:
using BenchmarkTools
using Evolutionary
using Random
Random.seed!(2);

rosenbrock(x) = (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2

N = 3
cmaes = CMAES(;mu=20, lambda=100)
ga = GA(populationSize=100, ɛ=0.1, selection=rouletteinv, crossover=intermediate(0.25), mutation=domainrange(fill(0.5,N)))
es = ES(initStrategy=IsotropicStrategy(N), recombination=average, srecombination=average, mutation=gaussian, smutation=gaussian, μ=10, ρ=3, λ=100, selection=:plus)
de = DE(populationSize = 100)

suite = BenchmarkGroup()

suite["methods"] = BenchmarkGroup(["ga", "cmaes", "es", "de"])
suite["methods"]["cmaes"] = @benchmarkable Evolutionary.optimize(rosenbrock, randn($N), $cmaes)
suite["methods"]["ga"] = @benchmarkable Evolutionary.optimize(rosenbrock, randn($N), $ga)
suite["methods"]["es"] = @benchmarkable Evolutionary.optimize(rosenbrock, randn($N), $es)
suite["methods"]["de"] = @benchmarkable Evolutionary.optimize(rosenbrock, randn($N), $de)

# parameters
loadparams!(suite, BenchmarkTools.load("bmsetup.json")[1], :evals, :samples);
# tune!(suite);
# BenchmarkTools.save("bmsetup.json", params(suite));
# exit()

# benchmark
results = run(suite, verbose = true)
BenchmarkTools.save("results.json", results)


text:
using Gadfly

set_default_plot_size(21cm, 8cm)

# Issue #678

x = repeat(collect(1:10), inner=[10]) .-0.5
y = repeat(collect(1:10), outer=[10]) .-0.5
palettef1 = Scale.lab_gradient(colorant"green", colorant"white", colorant"red")
palettef2 = Scale.lab_gradient(["blue", "ghostwhite","red"]...)
p1 = plot(x=x,y=y,color=x+y, Geom.rectbin, Scale.ContinuousColorScale(palettef1))
p2 = plot(x=x,y=y,color=x+y, Geom.rectbin, Scale.ContinuousColorScale(palettef2))
hstack(p1, p2)


text:
using Test
using SparseIR
using AssociatedLegendrePolynomials: Plm

@testset "augment.jl" begin
    @testset "MatsubaraConstBasis with stat = $stat" for stat in (fermion, boson)
        beta = 2.0
        value = 1.1
        b = MatsubaraConstBasis(stat, beta; value)
        shift::Int = Dict(fermion => 1, boson => 0)[stat]
        n = 2 .* collect(1:10) .+ shift
        @test b.uhat(n) ≈ fill(value, 1, length(n))
    end

    @testset "LegendreBasis with stat = $stat" for stat in (fermion, boson)
        β = 1.0
        Nl = 10
        cl = sqrt.(2 * collect(0:(Nl - 1)) .+ 1)
        basis = SparseIR.LegendreBasis(stat, β, Nl; cl)
        @test size(basis) == (Nl,)

        τ = Float64[0, 0.1 * β, 0.4 * β, β]
        uval = basis.u(τ)

        ref = Matrix{Float64}(undef, Nl, length(τ))
        for l in 0:(Nl - 1)
            x = 2τ / β .- 1
            ref[l + 1, :] .= cl[l + 1] * (√(2l + 1) / β) * Plm.(l, 0, x)
        end
        @test uval ≈ ref

        sign = stat == fermion ? -1 : 1

        # G(iv) = 1/(iv-pole)
        # G(τ) = -e^{-τ*pole}/(1 + e^{-β*pole}) [F]
        #        = -e^{-τ*pole}/(1 - e^{-β*pole}) [B]
        pole = 1.0
        τ_smpl = TauSampling(basis)
        gτ = -exp.(-τ_smpl.sampling_points * pole) / (1 - sign * exp(-β * pole))
        gl_from_τ = fit(τ_smpl, gτ)

        matsu_smpl = MatsubaraSampling(basis)
        giv = 1 ./ ((im * π / β) * matsu_smpl.sampling_points .- pole)
        gl_from_matsu = fit(matsu_smpl, giv)

        #println(maximum(abs, gl_from_τ-gl_from_matsu))
        #println(maximum(abs, gl_from_τ))
        #println("gl_from_τ", gl_from_τ[1:4])
        #println("gl_from_matsu", gl_from_matsu[1:4])
        @test isapprox(
            gl_from_τ, gl_from_matsu; atol=1e-10 * maximum(abs, gl_from_matsu), rtol=0
        )
    end
end


text:
# this file is inspired from StaticArrays.jl
# https://github.com/JuliaArrays/StaticArrays.jl

import Base: @propagate_inbounds, @_inline_meta, @pure

abstract type TupleVector{N,T} <: AbstractVector{T} end

# Being a member of TupleMatrixLike or TupleVectorLike implies that Val(A)
# returns a static Val instance (none of the dimensions are Dynamic). The converse may not be true.
# These are akin to aliases like StridedArray and in similarly bad taste, but the current approach
# in Base necessitates their existence.
const TupleMatrixLike{n,T} = Union{
    LinearAlgebra.Transpose{T, <:TupleVector{n,T}},
    LinearAlgebra.Adjoint{T, <:TupleVector{n,T}},
    LinearAlgebra.Diagonal{T, <:TupleVector{n,T}},
}
const TupleVectorLike{n,T} = Union{TupleVector{n,T}, TupleMatrixLike{n,T}}

@pure Base.length(::T) where T<:TupleVector{N} where N = N
@pure Base.length(::Type{<:TupleVector{N}}) where N = N
@pure Base.lastindex(::T) where T<:TupleVector{N} where N = N
@pure Base.size(::T) where T<:TupleVector{N} where N = (N,)
@pure Base.size(::Type{<:TupleVector{N}}) where N = (N,)
@pure @inline Base.size(t::T,d::Int) where T<:TupleVector{N} where N = d > 1 ? 1 : length(t)
@pure @inline Base.size(t::Type{<:TupleVector},d::Int) = d > 1 ? 1 : length(t)

include("SOneTo.jl")
include("util.jl")
include("traits.jl")
include("Values.jl")
include("Variables.jl")
include("FixedVector.jl")
include("initializers.jl")
include("convert.jl")
include("abstractvector.jl")
include("indexing.jl")
include("broadcast.jl")
include("mapreduce.jl")
include("arraymath.jl")
include("linalg.jl")


text:
const DISCRETE_INPLACE_DEFAULT = (t,u,du) -> fill!(du,zero(eltype(u)))
const DISCRETE_OUTOFPLACE_DEFAULT = (t,u) -> u

type DiscreteProblem{uType,tType,isinplace,F,C} <: AbstractDiscreteProblem{uType,tType,isinplace}
  f::F
  u0::uType
  tspan::Tuple{tType,tType}
  callback::C
end

function DiscreteProblem(f,u0,tspan; iip = isinplace(f,3),callback = nothing)
  DiscreteProblem{typeof(u0),eltype(tspan),iip,typeof(f),typeof(callback)}(f,u0,tspan,callback)
end

function DiscreteProblem(u0,tspan::Tuple; callback = nothing)
  iip = typeof(u0) <: AbstractArray
  if iip
    f = DISCRETE_INPLACE_DEFAULT
  else
    f = DISCRETE_OUTOFPLACE_DEFAULT
  end
  DiscreteProblem{typeof(u0),eltype(tspan),iip,typeof(f),typeof(callback)}(f,u0,tspan,callback)
end


text:
"""
Case 7:
This case study a three bus system with 2 machine located at bus 2 and 3.
The generator at bus 3 uses the model of a one d- one q- machine, and has a 5-mass shaft and a turbine governor.
The generator at bus 2 uses the model of a one d- one q- machine, and single mass shaft.
The fault disconnects a circuit between buses 1 and 2, doubling its impedance.
"""

##################################################
############### LOAD DATA ########################
##################################################

include(joinpath(TEST_FILES_DIR, "data_tests/test07.jl"))

##################################################
############### SOLVE PROBLEM ####################
##################################################

tspan = (0.0, 20.0);

# Define Fault: Change of YBus
Ybus_change = NetworkSwitch(
    1.0, #change at t = 1.0
    Ybus_fault,
) #New YBus

@testset "Test 07 5-Mass-shaft model ResidualModel" begin
    path = (joinpath(pwd(), "test-07"))
    !isdir(path) && mkdir(path)
    try
        # Define Simulation Problem
        sim = Simulation!(
            ResidualModel,
            threebus_sys, #system,
            path,
            tspan, #time span
            Ybus_change, #Type of Fault
        ) #initial guess

        # Test Initial Condition
        diff = [0.0]
        res = get_init_values_for_comparison(sim)
        for (k, v) in test07_x0_init
            diff[1] += LinearAlgebra.norm(res[k] - v)
        end

        @test (diff[1] < 1e-3)

        # Obtain small signal results for initial conditions
        small_sig = small_signal_analysis(sim)
        eigs = small_sig.eigenvalues
        @test small_sig.stable

        # Test Eigenvalues
        @test LinearAlgebra.norm(eigs - test07_eigvals) < 1e-3

        #Solve problem
        @test execute!(sim, IDA(), dtmax = 0.001) == PSID.SIMULATION_FINALIZED
        results = read_results(sim)

        # Obtain data for angles
        series = get_state_series(results, ("generator-103-1", :δ))
        series2 = get_state_series(results, ("generator-103-1", :δ_hp))
        series3 = get_state_series(results, ("generator-103-1", :δ_ip))
        series4 = get_state_series(results, ("generator-103-1", :δ_ex))
    finally
        @info("removing test files")
        rm(path, force = true, recursive = true)
    end
end

@testset "Test 07 5-Mass-shaft model MassMatrixModel" begin
    path = (joinpath(pwd(), "test-07"))
    !isdir(path) && mkdir(path)
    try
        # Define Simulation Problem
        sim = Simulation!(
            MassMatrixModel,
            threebus_sys, #system,
            path,
            tspan, #time span
            Ybus_change, #Type of Fault
        ) #initial guess

        # Test Initial Condition
        diff = [0.0]
        res = get_init_values_for_comparison(sim)
        for (k, v) in test07_x0_init
            diff[1] += LinearAlgebra.norm(res[k] - v)
        end

        @test (diff[1] < 1e-3)

        # Obtain small signal results for initial conditions
        small_sig = small_signal_analysis(sim)
        eigs = small_sig.eigenvalues
        @test small_sig.stable

        # Test Eigenvalues
        @test LinearAlgebra.norm(eigs - test07_eigvals) < 1e-3

        # Solve problem
        @test execute!(sim, Rodas4(), dtmax = 0.001) == PSID.SIMULATION_FINALIZED
        results = read_results(sim)

        # Obtain data for angles
        series = get_state_series(results, ("generator-103-1", :δ))
        series2 = get_state_series(results, ("generator-103-1", :δ_hp))
        series3 = get_state_series(results, ("generator-103-1", :δ_ip))
        series4 = get_state_series(results, ("generator-103-1", :δ_ex))
    finally
        @info("removing test files")
        rm(path, force = true, recursive = true)
    end
end


text:
"""
This package offers macros that delegate functions over one or more fields of a type;
and macros that delegate operations through fields to return a value of the same type.

**exports**

        @delegate_onefield,                    #     apply functions over field
        @delegate_onefield_twovars,            #          (return type from func)
        @delegate_onefield_threevars,          #          (return type from func)
        @delegate_onefield_astype,             #     and reobtain the same type
        @delegate_onefield_twovars_astype,     #          (return type from arg)
        @delegate_onefield_threevars_astype,   #          (return type from arg)
                                               #
        @delegate_twofields,                   #     apply functions over fields
        @delegate_twofields_twovars,           #          (return type from func)
        @delegate_twofields_astype,            #     and reobtain the same type
        @delegate_twofields_twovars_astype,    #          (return type from args)
                                               #
        @delegate_threefields,                 #     apply functions over fields
        @delegate_threefields_twovars,         #          (return type from func)
        @delegate_threefields_astype,          #     and reobtain the same type
        @delegate_threefields_twovars_astype   #          (return type from args)
"""
module TypedDelegation

export  @delegate_type,                        #
        @delegate_type_astype,                 #
                                               #
        @delegate_onefield,                    #     apply functions over field
        @delegate_onefield_twovars,            #          (return type from func)
        @delegate_onefield_threevars,          #          (return type from func)
        @delegate_onefield_astype,             #     and reobtain the same type
        @delegate_onefield_twovars_astype,     #          (return type from arg)
        @delegate_onefield_threevars_astype,   #          (return type from arg)
                                               #
        @delegate_twofields,                   #     apply functions over fields
        @delegate_twofields_twovars,           #          (return type from func)
        @delegate_twofields_astype,            #     and reobtain the same type
        @delegate_twofields_twovars_astype,    #          (return type from args)
                                               #
        @delegate_threefields,                 #     apply functions over fields
        @delegate_threefields_twovars,         #          (return type from func)
        @delegate_threefields_astype,          #     and reobtain the same type
        @delegate_threefields_twovars_astype   #          (return type from args)

#=
      internal macros, compositable
=#

"""
    @fieldtypes( aType )

works like aType.types
"""
macro fieldtypes( aType )
    esc( :( getfield( $aType, :types ) ) )
end

"""
    @fieldsyms( aType )

works like fieldnames(aType)
"""
macro fieldsyms( aType )
    esc( :( fieldnames( $aType ) ) )
end

"""
    @fieldsyms( aType, nfields )

yields the symbols for the first nfields declared in Type (in order)
"""
macro fieldsyms( aType, nfields )
    return quote
        local allsyms = @fieldsyms( $aType )
        local syms = allvals[1:$nfields]
        syms
    end
end

"""
    @getfield( varOfType, symfieldOfType )

works like getfield(varOfType, symfieldOfType)
"""
macro getfield( varOfType, symfieldOfType)
    esc( :( getfield( ($varOfType), ($symfieldOfType) )) )
end

"""
    @getfields( varOfType )

yields a tuple of field values for varOfType (in order)
"""
macro getfields( varOfType )
    return quote
        local syms = @fieldsyms( typeof( $varOfType ) )
        local vals = Vector( 0 )
        for item in syms
            push!( vals, @getfield( $varOfType, item ) )
        end
        (vals...,)
    end
end

"""
    @fieldvalues( varOfType )

yields the values of the fields in varOfType (in order)
"""
macro fieldvalues( varOfType )
    esc( :( @getfields( $varOfType ) ) )
end

"""
    @fieldvalues( varOfType, nfields )

yields the values of the first nfields declared in Type (in order)
"""
macro fieldvalues( varOfType, nfields )
    return quote
        local allvals = @fieldvalues( $varOfType )
        local vals = allvals[1:$nfields]
        vals
    end
end

#=
       delegation using the type itself
=#

macro delegate_type(sourceType, usedType, targetedFuncs)
    typesname  = esc( :($sourceType) )
    targetname = esc( :($usedType) )
    funcnames  = targetedFuncs.args
    nfuncs = length(funcnames)
    fdefs = Array{Expr}(undef, nfuncs)
    for i in 1:nfuncs
    funcname = esc(funcnames[i])
    fdefs[i] = quote
                 ($funcname)(::Type{$typesname}) =
                   ($funcname)($targetname)
               end
    end
    return Expr(:block, fdefs...)
end

macro delegate_type_astype(sourceType, usedType, targetedFuncs)
    typesname  = esc( :($sourceType) )
    targetname = esc( :($usedType) )
    funcnames  = targetedFuncs.args
    nfuncs = length(funcnames)
    fdefs = Array{Expr}(undef, nfuncs)
    for i in 1:nfuncs
    funcname = esc(funcnames[i])
    fdefs[i] = quote
                 ($funcname)(::Type{$typesname}) =
                   ($typesname)( ($funcname)($targetname) )
               end
    end
    return Expr(:block, fdefs...)
end

#=
         delegation using one field of a type
         when the result type differs from the type of the params
=#

"""
    @delegate_onefield(sourceType, sourcefield, targetedFuncs)

This returns a value of same types as the `targetedFuncs` result types.

    import Base: string, show

    struct MyInt16
      value::Int16
    end

    @delegate_onefield( MyInt16, value, [string, show]);

    three = MyInt16(3);
    seven = MyInt16(7);

    string(three) == "3"           # true
    show(seven)                    # 7
"""
macro delegate_onefield(sourceType, sourcefield, targetedFuncs)
    typesname  = esc( :($sourceType) )
    fieldname  = esc(Expr(:quote, sourcefield))
    funcnames  = targetedFuncs.args
    nfuncs = length(funcnames)
    fdefs = Array{Expr}(undef, nfuncs*2)
    for i in 1:nfuncs
        funcname = esc(funcnames[i])
        fdefs[i] = quote
                     ($funcname)(a::($typesname)) =
                       ($funcname)(getfield(a,($fieldname)))
                   end
        fdefs[i+nfuncs] = quote
                     ($funcname)(a::($typesname), args...) =
                       ($funcname)(getfield(a,($fieldname)), args...)
                   end
    end
    return Expr(:block, fdefs...)
end

"""
    @delegate_onefield_twovars(sourceType, sourcefield, targetedFuncs)

This returns a value of same types as the `targetedOps` result types.

A macro for field delegation over a function{T<:TheType}(arg1::T, arg2::T)

    import Base: (<), (<=)

    struct MyInt16
      value::Int16
    end;

    @delegate_onefield_twovars( MyInt16, value, [ (<), (<=) ] );

    three = MyInt16(3);
    seven = MyInt16(7);

    three <  seven                 # true
    !(seven <= three)              # true
"""
macro delegate_onefield_twovars(sourceType, sourcefield, targetedFuncs)
    typesname  = esc( :($sourceType) )
    fieldname  = esc(Expr(:quote, sourcefield))
    funcnames  = targetedFuncs.args
    nfuncs = length(funcnames)
    fdefs = Array{Expr}(undef, nfuncs*2)
    for i in 1:nfuncs
        funcname = esc(funcnames[i])
        fdefs[i] = quote
                     ($funcname)(a::($typesname), b::($typesname)) =
                       ($funcname)(getfield(a,($fieldname)),
                                   getfield(b,($fieldname)))
                   end
        fdefs[i+nfuncs] = quote
                     ($funcname)(a::($typesname), b::($typesname), args...) =
                       ($funcname)(getfield(a,($fieldname)),
                                   getfield(b,($fieldname)), args...)
                   end
    end
    return Expr(:block, fdefs...)
end

"""
    @delegate_onefield_threevars(sourceType, sourcefield, targetedFuncs)
"""
macro delegate_onefield_threevars(sourceType, sourcefield, targetedFuncs)
    typesname  = esc( :($sourceType) )
    fieldname  = esc(Expr(:quote, sourcefield))
    funcnames  = targetedFuncs.args
    nfuncs = length(funcnames)
    fdefs = Array{Expr}(undef, nfuncs*2)
    for i in 1:nfuncs
        funcname = esc(funcnames[i])
        fdefs[i] = quote
                     ($funcname)(a::($typesname), b::($typesname), c::($typesname)) =
                       ($funcname)(getfield(a,($fieldname)),
                                   getfield(b,($fieldname)),
                                   getfield(c,($fieldname)))
                   end
        fdefs[i+nfuncs] = quote
                     ($funcname)(a::($typesname), b::($typesname), c::($typesname), args...) =
                       ($funcname)(getfield(a,($fieldname)),
                                   getfield(b,($fieldname)),
                                   getfield(c,($fieldname)), args...)
                   end
    end
    return Expr(:block, fdefs...)
end

#=
         delegation using one field of a type
         when the result has the same type as the params
=#

"""
    @delegate_onefield_astype(sourceType, sourcefield, targetedOps)

This returns a value of the same type as the `sourceType` by rewrapping the result.

    import Base: abs, (-)

    struct MyInt16
      value::Int16
    end

    @delegate_onefield_astype( MyInt16, value, [abs, (-)]);

    three = MyInt16(3);
    seven = MyInt16(7);

    abs(three) == three            # true
    -(seven) === MyInt16(-7)       # true
"""
macro delegate_onefield_astype(sourceType, sourcefield, targetedOps)
    typesname  = esc( :($sourceType) )
    fieldname  = esc(Expr(:quote, sourcefield))
    funcnames  = targetedOps.args
    nfuncs = length(funcnames)
    fdefs = Array{Expr}(undef, nfuncs*2)
    for i in 1:nfuncs
        funcname = esc(funcnames[i])
        fdefs[i] = quote
                     ($funcname)(a::($typesname)) =
                       ($typesname)( ($funcname)(getfield(a,($fieldname))) )
                   end
        fdefs[i+nfuncs] = quote
                     ($funcname)(a::($typesname), args...) =
                       ($typesname)( ($funcname)(getfield(a,($fieldname)), args...) )
                   end
    end
    return Expr(:block, fdefs...)
end

"""
    @delegate_onefield_twovars_astype(sourceType, sourcefield, targetedOps)

This returns a value of the same type as the `sourceType` by rewrapping the result.

    import Base: (+), (-), (*)

    struct MyInt16
      value::Int16
    end

    @delegate_onefield_twovars_astype( MyInt16, value, [ (+), (*) ] );

    three = MyInt16(3);
    seven = MyInt16(7);

    import Base: hypot

    struct RightTriangle
      legA::Float64
      legB::Float64
    end;

    @delegate_twofields( RightTriangle, legA, legB, [ hypot, ] );

    myRightTriangle  = RightTriangle( 3.0, 4.0 );
    hypot(myRightTriangle)   #  5.0

    three + seven == MyInt16(3+7)  # true
    three * seven == MyInt16(3*7)  # true
"""
macro delegate_onefield_twovars_astype(sourceType, sourcefield, targetedOps)
    typesname  = esc( :($sourceType) )
    fieldname  = esc(Expr(:quote, sourcefield))
    funcnames  = targetedOps.args
    nfuncs = length(funcnames)
    fdefs = Array{Expr}(undef, nfuncs*2)
    for i in 1:nfuncs
        funcname = esc(funcnames[i])
        fdefs[i] = quote
                     ($funcname)(a::($typesname), b::($typesname)) =
                       ($typesname)( ($funcname)(getfield(a,($fieldname)),
                                                 getfield(b,($fieldname))) )
                   end
        fdefs[i+nfuncs] = quote
                     ($funcname)(a::($typesname), b::($typesname), args...) =
                       ($typesname)( ($funcname)(getfield(a,($fieldname)),
                                                 getfield(b,($fieldname)), args...) )
                   end
    end
    return Expr(:block, fdefs...)
end

"""
    @delegate_onefield_threevars_astype(sourceType, sourcefield, targetedFuncs)
"""
macro delegate_onefield_threevars_astype(sourceType, sourcefield, targetedFuncs)
    typesname  = esc( :($sourceType) )
    fieldname  = esc(Expr(:quote, sourcefield))
    funcnames  = targetedFuncs.args
    nfuncs = length(funcnames)
    fdefs = Array{Expr}(undef, nfuncs*2)
    for i in 1:nfuncs
        funcname = esc(funcnames[i])
        fdefs[i] = quote
                     ($funcname)(a::($typesname), b::($typesname), c::($typesname)) =
                       ($typesname)( ($funcname)(getfield(a,($fieldname)),
                                                 getfield(b,($fieldname)),
                                                 getfield(c,($fieldname))) )
                   end
        fdefs[i+nfuncs] = quote
                     ($funcname)(a::($typesname), b::($typesname), c::($typesname), args...) =
                       ($typesname)( ($funcname)(getfield(a,($fieldname)),
                                                 getfield(b,($fieldname)),
                                                 getfield(c,($fieldname)), args...) )
                   end
    end
    return Expr(:block, fdefs...)
end

#=
        delegation using two fields of a type
        when the result type differs from the type of the params
=#

"""    import Base: hypot

    struct RightTriangle
      legA::Float64
      legB::Float64
    end;

    @delegate_twofields( RightTriangle, legA, legB, [ hypot, ] );

    myRightTriangle  = RightTriangle( 3.0, 4.0 );
    hypot(myRightTriangle)   #  5.0
    @delegate_twofields(sourceType, firstfield, secondfield, targetedFuncs)

This returns a value of same types as the `targetedFuncs` result types.

    import Base: hypot

    struct RightTriangle
      legA::Float64
      legB::Float64
    end;

    @delegate_twofields( RightTriangle, legA, legB, [ hypot, ] );

    myRightTriangle  = RightTriangle( 3.0, 4.0 );
    hypot(myRightTriangle)   #  5.0
"""
macro delegate_twofields(sourceType, firstfield, secondfield, targetedFuncs)
    typesname  = esc( :($sourceType) )
    field1name = esc(Expr(:quote, firstfield))
    field2name = esc(Expr(:quote, secondfield))
    funcnames  = targetedFuncs.args
    nfuncs = length(funcnames)
    fdefs = Array{Expr}(undef, nfuncs*2)
    for i in 1:nfuncs
        funcname = esc(funcnames[i])
        fdefs[i] = quote
                     ($funcname)(a::($typesname)) =
                       ($funcname)(getfield(a, ($field1name)), getfield(a, ($field2name)))
                   end
        fdefs[i+nfuncs] = quote
                     ($funcname)(a::($typesname), args...) =
                       ($funcname)(getfield(a, ($field1name)), getfield(a, ($field2name)), args...)
                   end
    end
    return Expr(:block, fdefs...)
end

#=
        delegation using two fields of a type
        when the result has the same type as the params
=#

"""
    @delegate_twofields_astype(sourceType, firstfield, secondfield, targetedOps)

This returns a value of the same type as the `sourceType` by rewrapping the result.

    function renormalize(a::Float64, b::Float64)
      hi = a + b
      t = hi - a
      lo = (a - (hi - t)) + (b - t)
      return hi, lo
    end

    struct HiLo
      hi::Float64
      lo::Float64
    end;

    @delegate_twofields_astype( HiLo, hi, lo, [ renormalize, ] );

    myHiLo = renormalize( HiLo(12.555555555, 8000.333333333) );
    show(myHiLo)     # HiLo(8012.888888888,4.440892098500626e-14)
"""
macro delegate_twofields_astype(sourceType, firstfield, secondfield, targetedOps)
    typesname  = esc( :($sourceType) )
    field1name = esc(Expr(:quote, firstfield))
    field2name = esc(Expr(:quote, secondfield))
    funcnames  = targetedOps.args
    nfuncs = length(funcnames)
    fdefs = Array{Expr}(undef, nfuncs*2)
    for i in 1:nfuncs
        funcname = esc(funcnames[i])
        fdefs[i] = quote
                     ($funcname)(a::($typesname)) =
                        ($typesname)( ($funcname)(getfield(a, ($field1name)), getfield(a, ($field2name)))... )
                   end
        fdefs[i+nfuncs] = quote
                     ($funcname)(a::($typesname), args...) =
                        ($typesname)( ($funcname)(getfield(a, ($field1name)), getfield(a, ($field2name)), args...)... )
                   end
    end
    return Expr(:block, fdefs...)
end

"""
    @delegate_twofields_twovars(sourceType, firstfield, secondfield, targetedFuncs)

This returns a value of same types as the `targetedFuncs` result types.

A macro for field delegation over a function{T<:TheType}(arg1::T, arg2::T)
  that evalutes two fields of `TheType` from arg1 and also from arg2.

    result = targetedFunc( arg1.firstfield, arg1.secondfield,
                           arg2.firstfield, arg2.secondfield )

    import Base: mean

    type MyInterval
      lo::Float64
      hi::Float64

      function MyInterval(lo::Float64, hi::Float64)
         mn, mx = ifelse( hi<lo, (hi, lo), (lo, hi) )
         new( mn, mx )
      end
    end

    MyInterval(lo::T, hi::T) where T<:Real =
        MyInterval( Float64(lo), Float64(hi) )

    @delegate_twofields_twovars( MyInterval, lo, hi, [ mean, ])

    function mean( x::MyInterval )
        return x.lo * 0.5 + x.hi * 0.5
    end

    function mean( a::MyInterval, b::MyInterval )
        return mean( a ) * 0.5 + mean( b ) * 0.5
    end

    one_three = MyInterval(1.0, 3.0);
    two_four  = MyInterval(2.0, 4.0);

    mean( one_three, two_four ) == 2.5 # true
"""
macro delegate_twofields_twovars(sourceType, firstfield, secondfield, targetedFuncs)
    typesname  = esc( :($sourceType) )
    field1name = esc(Expr(:quote, firstfield))
    field2name = esc(Expr(:quote, secondfield))
    funcnames  = targetedFuncs.args
    nfuncs = length(funcnames)
    fdefs = Array{Expr}(undef, nfuncs*2)
    for i in 1:nfuncs
        funcname = esc(funcnames[i])
        fdefs[i] = quote
                     ($funcname)(a::($typesname), b::($typesname)) =
                         ($funcname)(getfield(a, ($field1name)), getfield(a, ($field2name)),
                                     getfield(b, ($field1name)), getfield(b, ($field2name)))
                   end
        fdefs[i+nfuncs] = quote
                     ($funcname)(a::($typesname), b::($typesname), args...) =
                         ($funcname)(getfield(a, ($field1name)), getfield(a, ($field2name)),
                                     getfield(b, ($field1name)), getfield(b, ($field2name)),
                                     args...)
                   end
    end
    return Expr(:block, fdefs...)
end

"""
    @delegate_twofields_twovars_astype(sourceType, firstfield, secondfield, targetedOps)

This returns a value of the same type as the `sourceType` by rewrapping the result.

A macro for field delegation over a function{T<:TheType}(arg1::T, arg2::T)
  that evalutes two fields of `TheType` from arg1 and also from arg2
  and applies itself over them, obtaining field_values for constructive
  generation of the result; a new realization of TheType( field_values... ).

    TheType( targetedFunc( arg1.firstfield, arg1.secondfield,
                           arg2.firstfield, arg2.secondfield )... )

    import Base: union, intersect

    type MyInterval
      lo::Float64
      hi::Float64

      function MyInterval(lo::Float64, hi::Float64)
         mn, mx = ifelse( hi<lo, (hi, lo), (lo, hi) )
         new( mn, mx )
      end
    end;

    MyInterval(lo::T, hi::T) where T<:Real =
        MyInterval( Float64(lo), Float64(hi) )

    @delegate_twofields_twovars_astype( MyInterval, lo, hi, [ union, ])

    function union( a::T, b::T ) where {T<:MyInterval}
        lo = min( a.lo, b.lo )
        hi = max( a.hi, b.hi )
        return T( lo, hi )
    end

    one_three = MyInterval(1.0, 3.0);
    two_four  = MyInterval(2.0, 4.0);

    union( one_three, two_four ) == MyInterval(1.0, 4.0) # true
"""
macro delegate_twofields_twovars_astype(sourceType, firstfield, secondfield, targetedOps)
    typesname  = esc( :($sourceType) )
    field1name = esc(Expr(:quote, firstfield))
    field2name = esc(Expr(:quote, secondfield))
    funcnames  = targetedOps.args
    nfuncs = length(funcnames)
    fdefs = Array{Expr}(undef, nfuncs*2)
    for i in 1:nfuncs
        funcname = esc(funcnames[i])
        fdefs[i] = quote
                     ($funcname)(a::($typesname), b::($typesname)) =
                        ($typesname)( ($funcname)(getfield(a, ($field1name)), getfield(a, ($field2name)),
                                                  getfield(b, ($field1name)), getfield(b, ($field2name)))... )
                   end
        fdefs[i+nfuncs] = quote
                     ($funcname)(a::($typesname), b::($typesname), args...) =
                        ($typesname)( ($funcname)(getfield(a, ($field1name)), getfield(a, ($field2name)),
                                                  getfield(b, ($field1name)), getfield(b, ($field2name)),
                                                  args...)... )
                   end
    end
    return Expr(:block, fdefs...)
end

#=
        delegation using three fields of a type
        when the result type differs from the type of the params
=#

"""
    @delegate_threefields(sourceType, firstfield, secondfield, targetedFuncs)

This returns a value of same types as the `targetedFuncs` result types.

    import Base: norm

    norm{R<:Real}(xs::Vararg{R,3}) = norm([xs...])

    struct XYZ
      x::Float64
      y::Float64
      z::Float64
    end;

    @delegate_threefields( XYZ, x, y, z, [ norm, ] );

    pointA  = XYZ( 3.0, 4.0, 5.0 );

    norm(pointA)   #  7.0710678+
"""
macro delegate_threefields(sourceType, firstfield, secondfield, thirdfield, targetedFuncs)
    typesname  = esc( :($sourceType) )
    field1name = esc(Expr(:quote, firstfield))
    field2name = esc(Expr(:quote, secondfield))
    field3name = esc(Expr(:quote, thirdfield))
    funcnames  = targetedFuncs.args
    nfuncs = length(funcnames)
    fdefs = Array{Expr}(undef, nfuncs*2)
    for i in 1:nfuncs
        funcname = esc(funcnames[i])
        fdefs[i] = quote
                     ($funcname)(a::($typesname)) =
                       ($funcname)(getfield(a, ($field1name)),
                                   getfield(a, ($field2name)),
                                   getfield(a, ($field3name)))
                   end
        fdefs[i+nfuncs] = quote
                     ($funcname)(a::($typesname), args...) =
                       ($funcname)(getfield(a, ($field1name)),
                                   getfield(a, ($field2name)),
                                   getfield(a, ($field3name)), args...)
                   end
    end
    return Expr(:block, fdefs...)
end

#=
        delegation using three fields of a type
        when the result has the same type as the params
=#

"""
    @delegate_threefields_astype(sourceType, firstfield, secondfield, targetedOps)

This returns a value of the same type as the `sourceType` by rewrapping the result.

    import Base: normalize

    normalize{R<:Real}(xs::Vararg{R,3}) = normalize([xs...])

    struct XYZ
      x::Float64
      y::Float64
      z::Float64
    end;

    @delegate_threefields_astype( XYZ, x, y, z, [ normalize, ] );

    pointA  = XYZ( 3.0, 4.0, 5.0 );

    normalize(pointA)   #  XYZ( 0.424264+, 0.565685+, 0.707107- )
"""
macro delegate_threefields_astype(sourceType, firstfield, secondfield, thirdfield, targetedOps)
    typesname  = esc( :($sourceType) )
    field1name = esc(Expr(:quote, firstfield))
    field2name = esc(Expr(:quote, secondfield))
    field3name = esc(Expr(:quote, thirdfield))
    funcnames  = targetedOps.args
    nfuncs = length(funcnames)
    fdefs = Array{Expr}(undef, nfuncs*2)
    for i in 1:nfuncs
        funcname = esc(funcnames[i])
        fdefs[i] = quote
                     ($funcname)(a::($typesname)) =
                        ($typesname)( ($funcname)(getfield(a, ($field1name)),
                                                  getfield(a, ($field2name)),
                                                  getfield(a, ($field3name)))... )
                   end
        fdefs[i+nfuncs] = quote
                     ($funcname)(a::($typesname), args...) =
                        ($typesname)( ($funcname)(getfield(a, ($field1name)),
                                                  getfield(a, ($field2name)),
                                                  getfield(a, ($field3name)), args...)... )
                   end
    end
    return Expr(:block, fdefs...)
end

"""
    @delegate_threefields_twovars(sourceType, firstfield, secondfield, thirdfield, targetedFuncs)

This returns a value of same types as the `targetedFuncs` result types.

A macro for field delegation over a function{T<:TheType}(arg1::T, arg2::T)
  that evalutes three fields of `TheType` from arg1 and also from arg2.

    result = targetedFunc( arg1.firstfield, arg1.secondfield, arg1.thirdfield,
                           arg2.firstfield, arg2.secondfield, arg2.thirdfield )

    import Base: norm, normalize, cross, sin

    normalize(xs::Vararg{R,3}) where R<:Real  = normalize([xs...])
    cross(xs::Vararg{R,6}) where R<:Real = cross([xs[1:3]...], [xs[4:6]...])

    struct XYZ
      x::Float64
      y::Float64
      z::Float64
    end;

    @delegate_threefields_astype( XYZ, x, y, z, [ normalize, ] );
    @delegate_threefields_twovars( XYZ, x, y, z, [ cross, ] );

    function sin( pointA::XYZ, pointB::XYZ )
        norm( cross( normalize(pointA), normalize(pointB) ) )
    end

    pointA  = XYZ( 3.0, 4.0, 5.0 );
    pointB  = XYZ( 5.0, 4.0, 3.0 );

    sin(pointA, pointB) #  0.391918+
"""
macro delegate_threefields_twovars(sourceType, firstfield, secondfield, thirdfield, targetedFuncs)
    typesname  = esc( :($sourceType) )
    field1name = esc(Expr(:quote, firstfield))
    field2name = esc(Expr(:quote, secondfield))
    field3name = esc(Expr(:quote, thirdfield))
    funcnames  = targetedFuncs.args
    nfuncs = length(funcnames)
    fdefs = Array{Expr}(undef, nfuncs*2)
    for i in 1:nfuncs
        funcname = esc(funcnames[i])
        fdefs[i] = quote
                     ($funcname)(a::($typesname), b::($typesname)) =
                         ($funcname)(getfield(a, ($field1name)), getfield(a, ($field2name)), getfield(a, ($field3name)),
                                     getfield(b, ($field1name)), getfield(b, ($field2name)), getfield(b, ($field2name)))
                   end
        fdefs[i+nfuncs] = quote
                     ($funcname)(a::($typesname), b::($typesname), args...) =
                         ($funcname)(getfield(a, ($field1name)), getfield(a, ($field2name)), getfield(a, ($field3name)),
                                     getfield(b, ($field1name)), getfield(b, ($field2name)), getfield(b, ($field2name)),
                                     args...)
                   end
    end
    return Expr(:block, fdefs...)
end

"""
    @delegate_threefields_twovars_astype(sourceType, firstfield, secondfield, thirdfield, targetedOps)

This returns a value of the same type as the `sourceType` by rewrapping the result.

A macro for field delegation over a function{T<:TheType}(arg1::T, arg2::T)
  that evalutes three fields of `TheType` from arg1 and also from arg2
  and applies itself over them, obtaining field_values for constructive
  generation of the result; a new realization of TheType( field_values... ).

  TheType( targetedFunc( arg1.firstfield, arg1.secondfield, arg1.thirdfield,
                         arg2.firstfield, arg2.secondfield, arg2.thirdfield )... )

    import Base: cross

    cross(xs::Vararg{R,6}) where R<:Real = cross([xs[1:3]...], [xs[4:6]...])

    struct XYZ
      x::Float64
      y::Float64
      z::Float64
    end

    @delegate_threefields_twovars_astype( XYZ, x, y, z, [ cross, ] );

    pointA  = XYZ( 3.0, 4.0, 5.0 );
    pointB  = XYZ( 5.0, 4.0, 3.0 );

    cross(pointA, pointB) #  XYZ(-8.0, 16.0, -8.0)
"""
macro delegate_threefields_twovars_astype(sourceType, firstfield, secondfield,  thirdfield, targetedOps)
    typesname  = esc( :($sourceType) )
    field1name = esc(Expr(:quote, firstfield))
    field2name = esc(Expr(:quote, secondfield))
    field3name = esc(Expr(:quote, thirdfield))
    funcnames  = targetedOps.args
    nfuncs = length(funcnames)
    fdefs = Array{Expr}(undef, nfuncs*2)
    for i in 1:nfuncs
        funcname = esc(funcnames[i])
        fdefs[i] = quote
                     ($funcname)(a::($typesname), b::($typesname)) =
                         ($typesname)( ($funcname)(getfield(a, ($field1name)), getfield(a, ($field2name)), getfield(a, ($field3name)),
                                                   getfield(b, ($field1name)), getfield(b, ($field2name)), getfield(b, ($field3name)))... )
                   end
        fdefs[i+nfuncs] = quote
                     ($funcname)(a::($typesname), b::($typesname), args...) =
                         ($typesname)( ($funcname)(getfield(a, ($field1name)), getfield(a, ($field2name)), getfield(a, ($field3name)),
                                                   getfield(b, ($field1name)), getfield(b, ($field2name)), getfield(b, ($field3name)),
                                                   args...)... )
                   end
    end
    return Expr(:block, fdefs...)
end

#=
    earlier implementations
    description and logic from John Myles White
        <https://gist.github.com/johnmyleswhite/5225361>
    additional macro work from DataStructures.jl
        <https://github.com/JuliaLang/DataStructures.jl/blob/master/src/delegate.jl?
    and delegation with n-ary ops from Toivo Henningsson
        <https://groups.google.com/forum/#!msg/julia-dev/MV7lYRgAcB0/-tS50TreaPoJ?
=#

end # moduleTypedDelegation


text:
function f(a, b=1 ; c=2, d=3)
  [a, b, c, d]
end

@assert f(100) == [100, 1, 2, 3]
@assert f(100, 200) == [100, 200, 2, 3]
@assert f(100, d=10) == [100, 1, 2, 10]
@assert f(100, 200, d=400, c=300) == [100, 200, 300, 400]

max_stars_repo_name: idrougge/ple

max_issues_repo_name: idrougge/ple

max_forks_repo_name: idrougge/ple


text:
#=

Slice sampling utilities

Ignacio Quintero Mächler

t(-_-t)

September 23 2017

=#

"""
    loop_slice_sampler(lhf         ::Function,
                       p           ::Array{Float64,1},
                       fp          ::Array{Float64,1},
                       nnps        ::Array{Int64,1},
                       nps         ::Array{Int64,1},
                       phid        ::Array{Int64,1},
                       mvps        ::Array{Array{Int64,1},1},
                       mvhfs       ::Array{Array{Int64,1},1},
                       w           ::Array{Float64,1},
                       npars       ::Int64,
                       niter       ::Int64,
                       nthin       ::Int64,
                       screen_print::Int64)

Run slice sampling.
"""
function loop_slice_sampler(lhf         ::Function,
                            p           ::Array{Float64,1},
                            fp          ::Array{Float64,1},
                            nnps        ::Array{Int64,1},
                            nps         ::Array{Int64,1},
                            phid        ::Array{Int64,1},
                            mvps        ::Array{Array{Int64,1},1},
                            nngps       ::Array{Array{Bool,1},1},
                            mvhfs       ::Array{Array{Int64,1},1},
                            w           ::Array{Float64,1},
                            npars       ::Int64,
                            niter       ::Int64,
                            nthin       ::Int64,
                            screen_print::Int64)

  # maximum number of multivariate updates
  maxmvu = if iszero(lastindex(mvps)) && iszero(lastindex(mvhfs))
    zero(1)
  elseif iszero(lastindex(mvhfs))
    maximum(map(length,mvps))
  elseif iszero(lastindex(mvps))
    maximum(map(length,mvhfs))
  else
    maximum((maximum(map(length,mvps)), maximum(map(length,mvhfs))))
  end

  nlogs = fld(niter,nthin)
  its   = Array{Float64,1}(undef, nlogs)
  hlog  = Array{Float64,1}(undef, nlogs)
  ps    = Array{Float64,2}(undef, nlogs,npars)
  Lv    = Array{Float64,1}(undef, maxmvu)
  Rv    = Array{Float64,1}(undef, maxmvu)

  lthin, lit = 0, 0

  # preallocate pp
  pp = copy(p)

  # preallocate fpp
  fpp = copy(fp)

  # start iterations
  prog = Progress(niter, screen_print, "running slice-sampler...", 20)

  hc = lhf(p, fp)

  for it in Base.OneTo(niter)

    #=
    univariate updates
    =#
    # nonnegative parameters
    for j in nnps
     S     = hc - randexp()
     L, R  = find_nonneg_int(p, pp, fp, j, S, lhf, w[j], npars)
     p, hc = sample_int(p, pp, fp, j, L, R, S, lhf, npars)
    end

    # real line parameters
    for j in nps
      S     = hc - randexp()
      L, R  = find_real_int(p, pp, fp, j, S, lhf, w[j], npars)
      p, hc = sample_int(p, pp, fp, j, L, R, S, lhf, npars)
    end

    # hidden factors
    for j in phid
     S     = hc - randexp()
     L, R  = find_nonneg_int(p, pp, fp, fpp, j, S, lhf, w[j], npars)
     p, fp, hc = sample_int(p, pp, fp, fpp, j, L, R, S, lhf, npars)
    end

    #=
    multivariate updates
    =#
    # non hidden factors
    for (i,mvp) in enumerate(mvps)
      S = hc - randexp()
      @views find_rect(p, pp, fp, mvp, nngps[i], S, lhf, w, Lv, Rv, npars)
      p, hc  = sample_rect(p, pp, fp, Lv, Rv, mvp, S, lhf, npars)
    end

    # hidden factors
    for mvp in mvhfs
      S = hc - randexp()
      @views find_rect(p, pp, fp, fpp, mvp, S, lhf, w, Lv, Rv, npars)
      p, fp, hc = sample_rect(p, pp, fp, fpp, Lv, Rv, mvp, S, lhf, npars)
    end

    # log samples
    lthin += 1
    if lthin == nthin
      @inbounds begin
        lit += 1
        setindex!(its,  it, lit)
        setindex!(hlog, hc, lit)
        setindex!(ps,   p,  lit, :)
      end
      lthin = 0
    end

    next!(prog)
  end

  return its, hlog, ps
end

"""
    w_sampler(lhf         ::Function,
              p           ::Array{Float64,1},
              fp          ::Array{Float64,1},
              nnps        ::Array{Int64,1},
              nps         ::Array{Int64,1},
              phid        ::Array{Int64,1},
              mvps        ::Array{Array{Int64,1},1},
              mvhfs       ::Array{Array{Int64,1},1},
              npars       ::Int64,
              optimal_w   ::Float64,
              screen_print::Int64,
              nburn       ::Int64,
              ntakew      ::Int64,
              winit       ::Float64)

Run slice sampler for burn-in and to estimate appropriate w's.
"""
function w_sampler(lhf         ::Function,
                   p           ::Array{Float64,1},
                   fp          ::Array{Float64,1},
                   nnps        ::Array{Int64,1},
                   nps         ::Array{Int64,1},
                   phid        ::Array{Int64,1},
                   mvps        ::Array{Array{Int64,1},1},
                   nngps        ::Array{Array{Bool,1},1},
                   mvhfs       ::Array{Array{Int64,1},1},
                   npars       ::Int64,
                   optimal_w   ::Float64,
                   screen_print::Int64,
                   nburn       ::Int64,
                   ntakew      ::Int64,
                   winit       ::Float64)

  if nburn < ntakew
    ntakew = nburn
  end

  # maximum number of multivariate updates
  maxmvu = if iszero(lastindex(mvps)) && iszero(lastindex(mvhfs))
    zero(1)
  elseif iszero(lastindex(mvhfs))
    maximum(map(length,mvps))
  elseif iszero(lastindex(mvps))
    maximum(map(length,mvhfs))
  else
    maximum((maximum(map(length,mvps)), maximum(map(length,mvhfs))))
  end

  Lv = Array{Float64,1}(undef, maxmvu)
  Rv = Array{Float64,1}(undef, maxmvu)

  w  = fill(winit, npars)
  ps = Array{Float64,2}(undef, nburn, npars)

  # posterior
  hc = lhf(p, fp)

  # preallocate pp and fpp
  pp  = copy(p)
  fpp = copy(fp)

  prog = Progress(nburn, screen_print, "estimating optimal widths...", 20)

  for it in Base.OneTo(nburn)

    #=
    univariate updates
    =#
    # nonnegative parameters
    for j in nnps
     S     = hc - randexp()
     L, R  = find_nonneg_int(p, pp, fp, j, S, lhf, w[j], npars)
     p, hc = sample_int(p, pp, fp, j, L, R, S, lhf, npars)
    end

    # real line parameters
    for j in nps
      S     = hc - randexp()
      L, R  = find_real_int(p, pp, fp, j, S, lhf, w[j], npars)
      p, hc = sample_int(p, pp, fp, j, L, R, S, lhf, npars)
    end

    # hidden factors
    for j in phid
     S     = hc - randexp()
     L, R  = find_nonneg_int(p, pp, fp, fpp, j, S, lhf, w[j], npars)
     p, fp, hc = sample_int(p, pp, fp, fpp, j, L, R, S, lhf, npars)
    end

    #=
    multivariate updates
    =#
    # non hidden factors
    for (i, mvp) in enumerate(mvps)
      S = hc - randexp()
      @views find_rect(p, pp, fp, mvp, nngps[i], S, lhf, w, Lv, Rv, npars)
      p, hc  = sample_rect(p, pp, fp, Lv, Rv, mvp, S, lhf, npars)
    end

    # hidden factors
    for mvp in mvhfs
      S = hc - randexp()
      @views find_rect(p, pp, fp, fpp, mvp, S, lhf, w, Lv, Rv, npars)
      p, fp, hc = sample_rect(p, pp, fp, fpp, Lv, Rv, mvp, S, lhf, npars)
    end

    @inbounds setindex!(ps, p, it, :)

    next!(prog)
  end

  sps = nburn-ntakew

  ps = ps[(nburn-ntakew+1):nburn,:]

  w = optimal_w .* (reduce(max, ps, dims=1) .- reduce(min, ps, dims=1))
  w = reshape(w, size(w,2))

  return (p, fp, w)::Tuple{Array{Float64,1},Array{Float64,1},Array{Float64,1}}
end

"""
    find_nonneg_int(p    ::Array{Float64},
                    pp   ::Array{Float64},
                    j    ::Int64,
                    S    ::Float64,
                    postf::Function,
                    w    ::Float64)

Estimate a nonnegative slice interval.
"""
function find_nonneg_int(p    ::Array{Float64,1},
                         pp   ::Array{Float64,1},
                         fp   ::Array{Float64,1},
                         j    ::Int64,
                         S    ::Float64,
                         postf::Function,
                         w    ::Float64,
                         npars::Int64)

  unsafe_copyto!(pp, 1, p, 1, npars)

  L::Float64 = pp[j] - w*rand()
  R::Float64 = L + w

  if L <= 0.0
    L = 1e-30
  end

  # left extreme
  pp[j] = L
  while S < postf(pp, fp)
    L -= w
    if L <= 0.0
      L = 1e-30
      break
    end
    pp[j] = L
  end

  # right extreme
  pp[j] = R
  while S < postf(pp, fp)
    R    += w
    pp[j] = R
  end

  return (L, R)::NTuple{2,Float64}
end

"""
    find_nonneg_int(p    ::Array{Float64,1},
                    pp   ::Array{Float64,1},
                    fp   ::Array{Float64,1},
                    fpp  ::Array{Float64,1},
                    j    ::Int64,
                    S    ::Float64,
                    postf::Function,
                    w    ::Float64)

Estimate a non_negative slice interval for hidden factors.
"""
function find_nonneg_int(p    ::Array{Float64,1},
                         pp   ::Array{Float64,1},
                         fp   ::Array{Float64,1},
                         fpp  ::Array{Float64,1},
                         j    ::Int64,
                         S    ::Float64,
                         postf::Function,
                         w    ::Float64,
                         npars::Int64)

  unsafe_copyto!(pp,  1, p,  1, npars)
  unsafe_copyto!(fpp, 1, fp, 1, npars)

  L = fpp[j] - w*rand()
  R = L + w

  if L <= 0.0
    L = 1e-30
  end

  # left extreme
  fpp[j] = L
  while S < postf(pp, fpp)
    L -= w
    if L <= 0.0
      L = 1e-30
      break
    end
    fpp[j] = L
  end

  # right extreme
  fpp[j] = R
  while S < postf(pp, fpp)
    R    += w
    fpp[j] = R
  end

  return (L, R)::NTuple{2,Float64}
end

"""
    find_real_int(p    ::Array{Float64},
                  pp   ::Array{Float64},
                  fp   ::Array{Float64},
                  j    ::Int64,
                  S    ::Float64,
                  postf::Function,
                  w    ::Float64)

Estimate a nonnegative slice interval.
"""
function find_real_int(p    ::Array{Float64,1},
                       pp   ::Array{Float64,1},
                       fp   ::Array{Float64,1},
                       j    ::Int64,
                       S    ::Float64,
                       postf::Function,
                       wj   ::Float64,
                       npars::Int64)

  unsafe_copyto!(pp, 1, p, 1, npars)

  L = pp[j] - wj*rand()
  R = L + wj

  # left extreme
  pp[j] = L::Float64
  while S < postf(pp, fp)
    L    -= wj::Float64
    pp[j] = L::Float64
  end

  # right extreme
  pp[j] = R::Float64
  while S < postf(pp, fp)
    R    += wj::Float64
    pp[j] = R::Float64
  end

  return (L, R)::NTuple{2,Float64}
end

"""
    sample_int(p    ::Array{Float64,1},
               pp   ::Array{Float64,1},
               fp   ::Array{Float64,1},
               j    ::Int64,
               L    ::Float64,
               R    ::Float64,
               S    ::Float64,
               postf::Function)

Take one sample within the interval of the slice.
"""
function sample_int(p    ::Array{Float64,1},
                    pp   ::Array{Float64,1},
                    fp   ::Array{Float64,1},
                    j    ::Int64,
                    L    ::Float64,
                    R    ::Float64,
                    S    ::Float64,
                    postf::Function,
                    npars::Int64)

  @inbounds begin
    unsafe_copyto!(pp, 1, p, 1, npars)

    while true
      pp[j] = (L + rand()*(R-L))::Float64

      hc = postf(pp, fp)::Float64
      if S < hc
        unsafe_copyto!(p, 1, pp, 1, npars)
        return (p, hc)::Tuple{Array{Float64,1}, Float64}
      end

      if pp[j] < p[j]
        L = pp[j]::Float64
      else
        R = pp[j]::Float64
      end
    end

  end
end

"""
    sample_int(p    ::Array{Float64,1},
               pp   ::Array{Float64,1},
               fp   ::Array{Float64,1},
               fpp  ::Array{Float64,1},
               j    ::Int64,
               L    ::Float64,
               R    ::Float64,
               S    ::Float64,
               postf::Function)

Take one sample within the interval of the slice for hidden factors.
"""
function sample_int(p    ::Array{Float64,1},
                    pp   ::Array{Float64,1},
                    fp   ::Array{Float64,1},
                    fpp  ::Array{Float64,1},
                    j    ::Int64,
                    L    ::Float64,
                    R    ::Float64,
                    S    ::Float64,
                    postf::Function,
                    npars::Int64)

  @inbounds begin
    unsafe_copyto!(pp,  1, p,  1, npars)
    unsafe_copyto!(fpp, 1, fp, 1, npars)

    while true
      fpp[j] = (L + rand()*(R-L))::Float64

      hc = postf(pp, fpp)::Float64
      if S < hc
        unsafe_copyto!(p,  1, pp,  1, npars)
        unsafe_copyto!(fp, 1, fpp, 1, npars)
        return (p, fp, hc)::Tuple{Array{Float64,1}, Array{Float64,1}, Float64}
      end

      if fpp[j] < fp[j]
        L = fpp[j]::Float64
      else
        R = fpp[j]::Float64
      end
    end

  end
end

"""
    find_rect(p    ::Array{Float64,1},
              pp   ::Array{Float64,1},
              fp   ::Array{Float64,1},
              mvp  ::Array{Int64,1},
              nngp  ::Array{Bool,1},
              S    ::Float64,
              postf::Function,
              w    ::Array{Float64,1},
              Lv   ::Array{Float64,1},
              Rv   ::Array{Float64,1},
              npars::Int64)

Returns left and right vertices of an hyper-rectangle.
"""
function find_rect(p    ::Array{Float64,1},
                   pp   ::Array{Float64,1},
                   fp   ::Array{Float64,1},
                   mvp  ::Array{Int64,1},
                   nngp  ::Array{Bool,1},
                   S    ::Float64,
                   postf::Function,
                   w    ::Array{Float64,1},
                   Lv   ::Array{Float64,1},
                   Rv   ::Array{Float64,1},
                   npars::Int64)

  @inbounds begin

    unsafe_copyto!(pp, 1, p, 1, npars)

    # randomly start rectangle
    for (i,j) in enumerate(mvp)
      Lv[i] = pp[j] - w[j]*rand()
      if nngp[i] && Lv[i] <= 0.0
        Lv[i] = 1e-30
      end
      Rv[i] = Lv[i] + w[j]
    end

    # left extremes
    for (i,j) in enumerate(mvp)
      pp[j] = Lv[i]::Float64
      while S < postf(pp, fp)
        Lv[i] -= w[j]::Float64
        if nngp[i] && Lv[i] <= 0.0
          Lv[i] = 1e-30
          break
        end
        pp[j]  = Lv[i]::Float64
      end
    end

    # right extremes
    for (i,j) in enumerate(mvp)
      pp[j] = Rv[i]::Float64
      while S < postf(pp, fp)
        Rv[i] += w[j]::Float64
        pp[j]  = Rv[i]::Float64
      end
    end
  end

  return nothing
end

"""
    find_rect(p    ::Array{Float64,1},
              pp   ::Array{Float64,1},
              fp   ::Array{Float64,1},
              mvp  ::Array{Int64,1},
              nngp  ::Array{Bool,1},
              S    ::Float64,
              postf::Function,
              w    ::Array{Float64,1},
              Lv   ::Array{Float64,1},
              Rv   ::Array{Float64,1},
              npars::Int64)

Returns left and right vertices of an hyper-rectangle for hidden factors.
"""
function find_rect(p    ::Array{Float64,1},
                   pp   ::Array{Float64,1},
                   fp   ::Array{Float64,1},
                   fpp  ::Array{Float64,1},
                   mvp  ::Array{Int64,1},
                   S    ::Float64,
                   postf::Function,
                   w    ::Array{Float64,1},
                   Lv   ::Array{Float64,1},
                   Rv   ::Array{Float64,1},
                   npars::Int64)

  @inbounds begin

    unsafe_copyto!(pp,  1, p, 1, npars)
    unsafe_copyto!(fpp, 1, fp, 1, npars)

    # randomly start rectangle
    for (i,j) in enumerate(mvp)
      Lv[i] = fpp[j] - w[j]*rand()
      if Lv[i] <= 0.0
        Lv[i] = 1e-30
      end
      Rv[i] = Lv[i] + w[j]
    end

    # left extremes
    for (i,j) in enumerate(mvp)
      fpp[j] = Lv[i]::Float64
      while S < postf(pp, fpp)
        Lv[i] -= w[j]::Float64
        if Lv[i] <= 0.0
          Lv[i] = 1e-30
          break
        end
        fpp[j] = Lv[i]::Float64
      end
    end

    # right extremes
    for (i,j) in enumerate(mvp)
      fpp[j] = Rv[i]::Float64
      while S < postf(pp, fpp)
        Rv[i] += w[j]::Float64
        fpp[j]  = Rv[i]::Float64
      end
    end
  end

  return nothing
end

"""
    sample_rect(p    ::Array{Float64,1},
                pp   ::Array{Float64,1},
                fp   ::Array{Float64,1},
                Lv   ::Array{Float64,1},
                Rv   ::Array{Float64,1},
                mvp  ::Array{Float64,1},
                S    ::Float64,
                postf::Function,
                npars::Int64)

Sample from an hyper-rectangle with shrinking.
"""
function sample_rect(p    ::Array{Float64,1},
                     pp   ::Array{Float64,1},
                     fp   ::Array{Float64,1},
                     Lv   ::Array{Float64,1},
                     Rv   ::Array{Float64,1},
                     mvp  ::Array{Int64,1},
                     S    ::Float64,
                     postf::Function,
                     npars::Int64)

  @inbounds begin

    unsafe_copyto!(pp, 1, p, 1, npars)

    while true

      # multivariate uniform sample
      for (i,j) in enumerate(mvp)
        pp[j] = Lv[i] + rand()*(Rv[i] - Lv[i])
      end

      # posterior
      hc = postf(pp, fp)

      if S < hc
        unsafe_copyto!(p, 1, pp, 1, npars)
        return(p, hc)::Tuple{Array{Float64,1}, Float64}
      end

      for (i,j) in enumerate(mvp)
        if pp[j] < p[j]
            Lv[i] = pp[j]::Float64
        else
            Rv[i] = pp[j]::Float64
        end
      end

    end
  end

  return nothing
end

"""
    sample_rect(p    ::Array{Float64,1},
                pp   ::Array{Float64,1},
                fp   ::Array{Float64,1},
                Lv   ::Array{Float64,1},
                Rv   ::Array{Float64,1},
                mvp  ::Array{Float64,1},
                S    ::Float64,
                postf::Function,
                npars::Int64)

Sample from an hyper-rectangle with shrinking for hidden factors.
"""
function sample_rect(p    ::Array{Float64,1},
                     pp   ::Array{Float64,1},
                     fp   ::Array{Float64,1},
                     fpp  ::Array{Float64,1},
                     Lv   ::Array{Float64,1},
                     Rv   ::Array{Float64,1},
                     mvp  ::Array{Int64,1},
                     S    ::Float64,
                     postf::Function,
                     npars::Int64)

  @inbounds begin

    unsafe_copyto!(pp,  1, p,  1, npars)
    unsafe_copyto!(fpp, 1, fp, 1, npars)

    while true

      # multivariate uniform sample
      for (i,j) in enumerate(mvp)
        fpp[j] = Lv[i] + rand()*(Rv[i] - Lv[i])
      end

      # posterior
      hc = postf(pp, fpp)

      if S < hc
        unsafe_copyto!(p, 1, pp, 1, npars)
        unsafe_copyto!(fp, 1, fpp, 1, npars)
        return(p, fp, hc)::Tuple{Array{Float64,1}, Array{Float64,1}, Float64}
      end

      for (i,j) in enumerate(mvp)
        if fpp[j] < fp[j]
            Lv[i] = fpp[j]::Float64
        else
            Rv[i] = fpp[j]::Float64
        end
      end

    end
  end

  return nothing
end

max_stars_repo_name: hmorlon/JPANDA

max_issues_repo_name: hmorlon/JPANDA

max_forks_repo_name: hmorlon/JPANDA


text:
### These were missing from CUDArt:

using CUDArt
using CUDArt: ContiguousArray

import Base: isequal, convert, reshape, resize!, copy!, isempty, fill!, pointer, issparse, deepcopy_internal

issparse(::CudaArray)=false

# For profiling:
copysync!{T}(dst::ContiguousArray{T}, src::ContiguousArray{T}; stream=null_stream)=(copy!(dst,src;stream=stream);gpusync();dst)
fillsync!{T}(X::CudaArray{T}, val; stream=null_stream)=(fill!(X,val;stream=stream);gpusync();X)

isequal(A::CudaArray,B::CudaArray) = (typeof(A)==typeof(B) && size(A)==size(B) && isequal(to_host(A),to_host(B)))

convert{A<:CudaArray,T,N}(::Type{A}, a::Array{T,N})=CudaArray(a)

reshape(a::CudaArray, dims::Dims)=reinterpret(eltype(a), a, dims)
reshape(a::CudaArray, dims::Int...)=reshape(a, dims)

function resize!(a::CudaVector, n::Integer)
    if n < length(a)
        a.dims = (n,)
    elseif n > length(a)
        b = CudaArray(eltype(a), convert(Int,n))
        copysync!(b, 1, a, 1, min(n, length(a)))
        free(a.ptr)
        a.ptr = b.ptr
        a.dims = b.dims
    end
    return a
end

# Generalizing low level copy using linear indexing to/from gpu
# arrays:

typealias BaseArray{T,N} Union{Array{T,N},SubArray{T,N},CudaArray{T,N}}

function copysync!{T<:Real}(dst::BaseArray{T}, di::Integer,
                            src::BaseArray{T}, si::Integer,
                            n::Integer; stream=null_stream)
    eltype(src) <: eltype(dst) || error("$(eltype(dst)) != $(eltype(src))")
    isbits(T) || error("$T not a bits type")
    if si+n-1 > length(src) || di+n-1 > length(dst) || di < 1 || si < 1
        throw(BoundsError())
    end
    esize = sizeof(eltype(src))
    nbytes = n * esize
    dptr = pointer(dst) + (di-1) * esize
    sptr = pointer(src) + (si-1) * esize
    CUDArt.rt.cudaMemcpyAsync(dptr, sptr, nbytes, CUDArt.cudamemcpykind(dst, src), stream)
    gpusync(); return dst
end

isempty(a::CudaArray)=(length(a)==0)

# This one has to be defined like this because of a conflict with the CUDArt version:
#fillsync!(A::AbstractCudaArray,x::Number)=(isempty(A)||cudnnSetTensor(A, x);gpusync();A)
#fillsync!(A::CudaArray,x::Number)=(isempty(A)||cudnnSetTensor(A, x);gpusync();A)

pointer{T}(x::CudaArray{T}, i::Integer) = pointer(x) + (i-1)*sizeof(T)

convert{A<:CudaArray}(::Type{A}, a::Array)=CudaArray(a)
convert{A<:CudaArray}(::Type{A}, a::SparseMatrixCSC)=CudaArray(full(a))
convert{A<:Array}(::Type{A}, a::CudaArray)=to_host(a)
convert{A<:SparseMatrixCSC}(::Type{A}, a::CudaArray)=sparse(to_host(a))

deepcopy_internal(x::CudaArray, s::ObjectIdDict)=(haskey(s,x)||(s[x]=copy(x));s[x])
cpucopy_internal(x::CudaArray, s::ObjectIdDict)=(haskey(s,x)||(s[x]=to_host(x));s[x])
gpucopy_internal(x::CudaArray, s::ObjectIdDict)=deepcopy_internal(x,s)
gpucopy_internal{T<:Number}(x::Array{T}, s::ObjectIdDict)=(haskey(s,x)||(s[x]=CudaArray(x));s[x])

# AbstractArray methods:

Base.getindex{T}(x::CudaArray{T}, i::Integer)=copysync!(T[0], 1, x, i, 1)[1]
Base.setindex!{T}(x::CudaArray{T}, v, i::Integer)=copysync!(x, i, T[convert(T,v)], 1, 1)
Base.endof(x::CudaArray)=length(x)

# Finding memory usage:

_getbytes(x::DataType,d)=sizeof(Int)
_getbytes(x::NTuple,d)=sum(map(y->_getbytes(y,d), x))
_getbytes(x::AbstractCudaArray,d)=(haskey(d,x) ? 0 : (d[x]=1; length(x) * sizeof(eltype(x))))
_getbytes(x::Symbol,d)=sizeof(Int)

function _getbytes(x::DenseArray,d)
    haskey(d,x) && return 0; d[x]=1
    total = sizeof(x)
    if !isbits(eltype(x))
        for i = 1:length(x)
            isassigned(x,i) || continue
            isize = _getbytes(x[i],d)
            # @show (typeof(x), i, isize)
            total += isize
        end
    end
    # @show (typeof(x), total)
    # @show (eltype(x), size(x), total)
    return total
end

function _getbytes(x,d)
    total = sizeof(x)
    isbits(x) && return total
    haskey(d,x) && return 0; d[x]=1
    fieldNames = fieldnames(x)
    if fieldNames != ()
        for fieldName in fieldNames
            isdefined(x, fieldName) || continue
            f = x.(fieldName)
            fieldBytes = _getbytes(f,d)
            fieldSize = (isa(f, Union{AbstractArray,CudaArray}) ? (eltype(f),size(f)) : ())
            # @show (typeof(x), fieldName, fieldSize, fieldBytes)
            total += fieldBytes
        end
    end
    return total
end

getbytes(x)=_getbytes(x, ObjectIdDict())

# Changing the way CudaArray prints:

if false # this takes up extra memory
# if !isdefined(:_CudaArray)
    import Base: size, linearindexing, getindex, writemime, summary
    using Base: with_output_limit, showarray, dims2string
    type _CudaArray{T,N} <: AbstractArray{T,N}; a::CudaArray{T,N}; end
    _CudaArray{T,N}(a::CudaArray{T,N})=_CudaArray{T,N}(a)
    size(a::_CudaArray)=size(a.a)
    linearindexing(::_CudaArray)=Base.LinearFast()
    getindex{T}(a::_CudaArray{T},i::Int)=getindex(a.a,i)
    summary{T,N}(a::_CudaArray{T,N})=string(dims2string(size(a)), " CudaArray{$T,$N}")
    writemime(io::IO, ::MIME"text/plain", v::CudaArray)=with_output_limit(()->showarray(io, _CudaArray(v), header=true, repr=false))
end

# To be able to load/save CudaArrays:
if !isdefined(:_CudaArraySave)
    using JLD
    type _CudaArraySave; a::Array; end
    JLD.writeas(c::CudaArray) = _CudaArraySave(to_host(c))
    JLD.readas(d::_CudaArraySave) = CudaArray(d.a)
end