text:
#=
Copyright (c) 2015, Intel Corporation
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
- Redistributions of source code must retain the above copyright notice,
  this list of conditions and the following disclaimer.
- Redistributions in binary form must reproduce the above copyright notice,
  this list of conditions and the following disclaimer in the documentation
  and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
THE POSSIBILITY OF SUCH DAMAGE.
=#

module LambdaHandling

import ..DebugMsg
DebugMsg.init()

using CompilerTools
using ..Helper
import ..Helper.toLHSVar
using CompilerTools.AstWalker
using Core.Inference: to_tuple_type

import Base.show

export VarDef, LambdaVarInfo, toRHSVar, lookupLHSVarByName, lookupVariableName
export getDesc, setDesc, getType, setType, getReturnType, setReturnType, getVarDef
export isInputParameter, isVarArgParameter, getInputParameters, setInputParameters, getInputParametersAsExpr
export isEscapingVariable, getEscapingVariables, addEscapingVariable, setEscapingVariable, unsetEscapingVariable
export isVariableDefined, isLocalVariable, getLocalVariables, getLocalVariablesNoParam, addLocalVariable, addTempVariable
export getBody, getReturnType, setReturnType, getInputParametersAsLHSVar
export lambdaToLambdaVarInfo, LambdaVarInfoToLambda, lambdaTypeinf, prependStatements
export getRefParams, getArrayParams, updateAssignedDesc, getEscapingVariablesAsLHSVar
export getStaticParameterValue
export mergeLambdaVarInfo, replaceExprWithDict!, countVariables
export ISCAPTURED, ISASSIGNED, ISASSIGNEDBYINNERFUNCTION, ISCONST, ISASSIGNEDONCE

# Possible values of VarDef descriptor that can be OR'ed together.
const ISCAPTURED = 1
const ISASSIGNED = 2
const ISASSIGNEDBYINNERFUNCTION = 4
const ISCONST = 8
const ISASSIGNEDONCE = 16

const digits = Char[x for x in "0123456789" ]

if VERSION >= v"0.5.0-dev+3875"
typealias DescType UInt8
else
typealias DescType Int64
end
getDefaultDesc() = convert(DescType, 0)

const emptyVarName = Symbol("")

"""
Represents a variable's definition, including its (optional) symbol, type, desc, and (optional) id.
"""
type VarDef
  name :: Symbol
  typ  :: Type
  desc :: DescType
  id   :: Int

  function VarDef(n::Symbol, t::Type, d::DescType, i::Int)
    new(n, t, d, i)
  end

  # no id, assuming -1.
  function VarDef(n::Symbol, t::Type, d::DescType)
    new(n, t, d, -1)
  end

  # no id, assuming -1.
  # no desc, assuming 0
  function VarDef(n::Symbol, t::Type)
    new(n, t, convert(DescType,0), -1)
  end

  # no symbol name, meaning GenSym (or SSAValue)
  function VarDef(t::Type, d::DescType, i::Int)
    new(emptyVarName, t, d, i)
  end
end

type LambdaVarInfo
  input_params  :: Array{Symbol,1}          # input parameters always have a name symbol
  vararg_params :: Array{Symbol,1}          # record which parameter is vararg
  escaping_vars :: Array{Symbol,1}          # escaping variables always have a name symbol
  var_defs      :: Array{VarDef}            # all variable definitions
  static_parameter_names :: Array{Any,1}
  static_parameter_values :: Array{Any,1}
  return_type   :: Any
  orig_info     :: Union{LambdaInfo,Void}

  function LambdaVarInfo()
    new(Symbol[], Symbol[], Symbol[], VarDef[], Any[], Any[], Void, nothing)
  end

  function LambdaVarInfo(li::LambdaVarInfo)
    new(copy(li.input_params), copy(li.vararg_params), copy(li.escaping_vars), copy(li.var_defs),
        copy(li.static_parameter_names), copy(li.static_parameter_values),
        li.return_type, li.orig_info == nothing ? nothing : li.orig_info)
  end

if VERSION >= v"0.5.0-dev+3875"
  function LambdaVarInfo(lambda::LambdaInfo, body)
    @dprintln(3, "lambda = ", lambda)
    return_typ = lambda.rettype
    input_params = Symbol[]
    var_defs = VarDef[]
    escaping_vars = Symbol[]
    allnames = Set{Symbol}()
    for i = 1:length(lambda.slotnames)
        name = lambda.slotnames[i]
        typ = lambda.slottypes[i]
        desc = lambda.slotflags[i]
        if in(name, allnames)
            name = gensym(string(name, "@", i))
        elseif startswith(string(name), digits) # numerical names are very confusing, fix them!
            name = gensym(string("@", name))
        end
        push!(allnames, name)
        vd = VarDef(name, typ, desc, i)
        push!(var_defs, vd)
        # skip #self# in parameters
        if i > 1 && i <= lambda.nargs
            push!(input_params, name)
        end
    end

    for i = 1:length(lambda.ssavaluetypes)
        vd = VarDef(emptyVarName, lambda.ssavaluetypes[i], convert(DescType, ISASSIGNED | ISASSIGNEDONCE), i-1)
        push!(var_defs, vd)
    end
    vararg_params = lambda.isva ? Symbol[input_params[end]] : Symbol[]
    static_parameter_names = Any[ x for x in lambda.sparam_syms]
    static_parameter_values = Any[ x for x in lambda.sparam_vals]
    @dprintln(3, "LambdaVarInfo.sparam_vals = ", lambda.sparam_vals)
    x = new(input_params, vararg_params, Symbol[], var_defs, static_parameter_names, static_parameter_values, return_typ, lambda)
    @dprintln(3, "LambdaVarInfo = ", x)
    return x
  end
else
  function LambdaVarInfo(lambda::LambdaInfo, ast)
    li = LambdaVarInfo(ast)
    li.orig_info = lambda
    return li
  end

  function LambdaVarInfo(lambda::Expr)
    @assert (lambda.head == :lambda) "Expect a lambda Expr as argument to LambdaVarInfo(..)"
    assert(length(lambda.args) == 3)

    input_params = Symbol[]
    vararg_params = Symbol[]
    var_defs = VarDef[]
    escaping_vars = Symbol[]

    #parameters
    for i = 1:length(lambda.args[1])
        p = lambda.args[1][i]
        if isa(p, Symbol)
          push!(input_params, p)
        elseif isa(p, SymbolNode)
          push!(input_params, p.name)
        elseif isa(p, Expr) && (p.head == :(::))
          @dprintln(3, "Got typed parameter: ", p)
          assert(isa(p.args[1], Symbol))
          push!(input_params, p.args[1])
          if isa(p.args[2], Expr) && p.args[2].head == :(...)
            push!(vararg_params, p.args[1])
          end
        end
    end

    meta = lambda.args[2]
    assert(length(meta) >= 3)
    #local variables (including parameter and its type)
    locals = meta[1]
    for i = 1:length(locals)
        name = locals[i][1]
        typ  = locals[i][2]
        desc = locals[i][3]
        @assert (isa(name, Symbol)) "Expect local variable name to be Symbol, but got " * string(name)
        @assert (isa(typ, Type)) "Expect local variable type, but got " * string(typ)
        @assert (isa(desc, DescType)) "Expect local variable desc to be " * string(DescType) * " but got " * string(typeof(desc))
        push!(var_defs, VarDef(name, typ, desc))
    end
    #escaping variables
    escapings = meta[2]
    for i = 1:length(escapings)
        name = escapings[i][1]
        typ  = escapings[i][2]
        desc = escapings[i][3]
        @assert (isa(name, Symbol)) "Expect local variable name to be Symbol, but got " * string(name)
        @assert (isa(typ, Type)) "Expect local variable type, but got " * string(typ)
        @assert (isa(desc, DescType)) "Expect local variable desc to be " * string(DescType) * " but got " * string(typeof(desc))
        push!(var_defs, VarDef(name, typ, desc))
        push!(escaping_vars, name)
    end
    #local GenSyms
    gensyms = meta[3]
    if isa(gensyms, Array)
        for i = 1:length(gensyms)
            typ = gensyms[i]
            desc = ISASSIGNED | ISASSIGNEDONCE
            id = i-1
            push!(var_defs, VarDef(typ, desc, id))
        end
    end
    static_param_values = length(meta) > 3 ? Any[ x for x in meta[3]] : Any[]
    static_param_names = Any[ gensym("T") for x in static_param_values]

    @assert (isa(lambda.args[3], Expr) && (lambda.args[3].head == :body)) "Expect a body Expr, but got " * string(lambda[3])
    return_typ = lambda.args[3].typ

    # we do not store original info
    new(input_params, vararg_params, escaping_vars, var_defs, static_param_names, static_param_values, return_typ, nothing)
  end
end

end

symbolToLHSVar(s :: Symbol, linfo :: LambdaVarInfo) = s

if VERSION >= v"0.5.0-dev+3875"

matchVarDef(x :: Symbol, vd :: VarDef) = vd.name == x
matchVarDef(x :: SlotNumber, vd :: VarDef) = vd.name != emptyVarName && vd.id == x.id
matchVarDef(x :: TypedSlot, vd :: VarDef) = vd.name != emptyVarName && vd.id == x.id
matchVarDef(x :: SSAValue, vd :: VarDef) = vd.name == emptyVarName && vd.id == x.id

toRHSVar(x :: Symbol, typ, linfo :: LambdaVarInfo) = toRHSVar(lookupLHSVarByName(x, linfo), linfo)
toRHSVar(x :: SlotNumber, typ, linfo :: LambdaVarInfo) = TypedSlot(x.id, typ)
toRHSVar(x :: TypedSlot, typ, linfo :: LambdaVarInfo) = x
toRHSVar(x :: SSAValue, typ, linfo :: LambdaVarInfo) = x
toRHSVar(vd :: VarDef) = vd.name == emptyVarName ? SSAValue(vd.id) : TypedSlot(vd.id, vd.typ)

toLHSVar(x :: Symbol, linfo :: LambdaVarInfo) = lookupLHSVarByName(x, linfo)
toLHSVar(x, linfo) = toLHSVar(x)
toLHSVar(vd :: VarDef) = vd.name == emptyVarName ? SSAValue(vd.id) : SlotNumber(vd.id)

else

matchVarDef(x :: Symbol, vd :: VarDef) = vd.name == x
matchVarDef(x :: SymbolNode, vd :: VarDef) = vd.name == x.name
matchVarDef(x :: GenSym, vd :: VarDef) = vd.name == emptyVarName && vd.id == x.id

toRHSVar(x :: Symbol, typ, linfo :: LambdaVarInfo) = SymbolNode(x, typ)
toRHSVar(x :: SymbolNode, typ, linfo :: LambdaVarInfo) = x
toRHSVar(x :: GenSym, typ, linfo :: LambdaVarInfo) = x
toRHSVar(vd :: VarDef) = vd.name == emptyVarName ? GenSym(vd.id) : SymbolNode(vd.name, vd.typ)

toLHSVar(x, linfo) = toLHSVar(x)
toLHSVar(vd :: VarDef) = vd.name == emptyVarName ? GenSym(vd.id) : vd.name
end

"""
Consolidate GenSym ids in LambdaVarInfo.
Return a dictionary to be used with replaceExprWithDict!(..)
"""
function consolidateLambdaVarInfo!(li::LambdaVarInfo)
    max_id = -1
    count  = Dict{Int,VarDef}()
    for vd in li.var_defs
        if vd.name == emptyVarName
            if vd.id > max_id
                max_id = vd.id
            end
            count[vd.id] = vd
        end
    end
    dict = Dict{LHSVar,Any}()
    @assert (max_id + 1 >= length(count)) "duplicate entries in " * string(li)
    if max_id < 0 || max_id + 1 == length(count)
        return dict
    end
    # We preserve all ids that are less than count
    for i = 0:length(count)-1
        if !haskey(count, i)
            next_id = length(count)
            while next_id <= max_id
                if haskey(count, next_id) break end
                next_id = next_id + 1
            end
            @assert (next_id <= max_id) "cannot find replacement GenSym for " * string(i)
            vd = count[next_id]
            vd.id = i
            count[i] = vd
            delete!(count, next_id)
            dict[GenSym(next_id)] = GenSym(i)
        end
    end
    @dprintln(3, "consolidateLambdaVarInfo! li=", li)
    @dprintln(3, "consolidateLambdaVarInfo! returns ", dict)
    return dict
end

function toRHSVar(x, linfo)
    return toRHSVar(x, getType(x, linfo), linfo)
end

"""
Pretty print a LambdaVarInfo.
"""
function show(io :: IO, li :: LambdaVarInfo)
  println(io, "Inputs = ")
  for x in li.input_params
    print(io, x)
    if in(x, li.vararg_params)
      print(io, "...")
    end
    print(io, ", ")
  end
  println(io)
  if !isempty(li.static_parameter_names)
    println(io, "Static Parameters = ", zip(li.static_parameter_names, li.static_parameter_values))
  end
  if !isempty(li.escaping_vars)
    println(io, "Escaping Variable Names = ", li.escaping_vars)
  end
  if li.return_type != nothing
    println(io, "Return type = ", li.return_type)
  end
  if !isempty(li.var_defs)
    println(io, "VarDefs")
    for i in li.var_defs
      println(io, "    ", i)
    end
  end
end

"""
Returns true if the given variable is an input parameter.
"""
function isInputParameter(x::LHSVar, li :: LambdaVarInfo)
    name = nothing
    for vd in li.var_defs
        if matchVarDef(x, vd)
            name = vd.name
        end
    end
    return !is(name, nothing) && in(name, li.input_params)
end

"""
Return the (fresh copy of) input parameters as array of Symbols
"""
function getInputParameters(li :: LambdaVarInfo)
  return copy(li.input_params)
end

function getInputParametersAsLHSVar(li :: LambdaVarInfo)
  return [toLHSVar(x, li) for x in li.input_params]
end

"""
Return the (fresh copy of) input parameters as array of Any, so vararg parameters are preserved.
"""
function getInputParametersAsExpr(li :: LambdaVarInfo)
  params = Array(Any, length(li.input_params))
  for i = 1:length(li.input_params)
    p = li.input_params[i]
    if in(p, li.vararg_params)
      params[i] = Expr(:(::), p, Expr(:(...), Any))
    else
      params[i] = p
    end
  end
  return params
end

"""
Set the input parameters to be the given array of Symbols.
The caller has to ensure they are also added to locals in the LambdaVarInfo.
"""
function setInputParameters(params :: Array{Symbol, 1}, li :: LambdaVarInfo)
  li.input_params = copy(params)
end

"""
Returns the type of a variable.
"""
function getType(x::Union{Symbol,RHSVar}, li::LambdaVarInfo)
    for vd in li.var_defs
        if matchVarDef(x, vd)
            return vd.typ
        end
    end
    error("Variable ", x, " is not found in ", li)
    #@dprintln(3, "Variable ", x, " is not found in ", li)
    #return Void
end

getType(x::Number, li::LambdaVarInfo) = typeof(x)

"""
Set the type for a local variable.
"""
function setType(x :: Union{Symbol,RHSVar}, typ :: Type, li :: LambdaVarInfo)
    for vd in li.var_defs
        if matchVarDef(x, vd)
            vd.typ = typ
            return
        end
    end
    error("Variable ", x, " is not found in ", li)
end

"""
Returns the descriptor for a local variable.
"""
function getDesc(x :: Union{Symbol,RHSVar}, li :: LambdaVarInfo)
    for vd in li.var_defs
        if matchVarDef(x, vd)
            return vd.desc
        end
    end
    error("Variable ", x, " is not found in ", li)
end

"""
Returns the descriptor for a local variable.
"""
function getVarDef(x :: Union{Symbol,RHSVar}, li :: LambdaVarInfo)
    for vd in li.var_defs
        if matchVarDef(x, vd)
            return vd
        end
    end
    error("Variable ", x, " is not found in ", li)
end

"""
Set the descriptor for a local variable.
"""
function setDesc(x :: Union{Symbol,RHSVar}, desc, li :: LambdaVarInfo)
    for vd in li.var_defs
        if matchVarDef(x, vd)
            vd.desc = convert(DescType, desc)
            return
        end
    end
    error("Variable ", x, " is not found in ", li)
end

"""
Returns true if the given variable is either a local (including input parameters) or escaping variable.
"""
function isVariableDefined(x::Union{Symbol,RHSVar}, li :: LambdaVarInfo)
    name = nothing
    for vd in li.var_defs
        if matchVarDef(x, vd)
            name = vd.name
        end
    end
    return name != nothing
end

""""
Returns true if the given variable is a local variable (including input parameters, but exclude escaping variables).
"""
function isLocalVariable(x::Union{Symbol,RHSVar}, li :: LambdaVarInfo)
    name = nothing
    for vd in li.var_defs
        if matchVarDef(x, vd)
            name = vd.name
        end
    end
    return name != nothing && !in(name, li.escaping_vars)
end

"""
Returns true if the given variable is an escaping variable.
"""
function isEscapingVariable(x::Union{Symbol,RHSVar}, li :: LambdaVarInfo)
    name = nothing
    for vd in li.var_defs
        if matchVarDef(x, vd)
            name = vd.name
        end
    end
    return name != nothing && in(name, li.escaping_vars)
end

"""
Returns an array of local variables and GenSyms, excluding escaping variables.
"""
function getLocalVariables(li :: LambdaVarInfo)
    locals = LHSVar[]
    for vd in li.var_defs
        if !in(vd.name, li.escaping_vars)
            push!(locals, toLHSVar(vd))
        end
    end
    return locals
end

"""
Returns an array of local variables and GenSyms, excluding escaping variables and parameters.
"""
function getLocalVariablesNoParam(li :: LambdaVarInfo)
    locals = LHSVar[]
    for vd in li.var_defs
        if !in(vd.name, li.escaping_vars) && !in(vd.name, li.input_params)
            push!(locals, toLHSVar(vd))
        end
    end
    return locals
end

"""
Returns the (fresh copy of) escaping variables as an array of Symbols.
"""
getEscapingVariables(li :: LambdaVarInfo) = copy(li.escaping_vars)
getEscapingVariablesAsLHSVar(li :: LambdaVarInfo) = map(x -> toLHSVar(x, li), li.escaping_vars)

"""
Adds a new local variable with the given Symbol "s", type "typ", descriptor "desc" in LambdaVarInfo "li".
Throw error if the variable has already existed.
Return the newly added variable (RHSVar) as result.
"""
function addLocalVariable(s :: Symbol, typ :: Type, desc, li :: LambdaVarInfo, assert_unique = true)
    @dprintln(3, "addLocalVariable s = ", s, " type = ", typ, " desc = ", desc, " li = ", li)
    max_id = 0
    for vd in li.var_defs
        if assert_unique && matchVarDef(s, vd)
            error("Variable ", s, " already exists in ", li)
        end
        if vd.name != emptyVarName && vd.id > max_id
            max_id = vd.id
        end
    end
    vd = VarDef(s, typ, convert(DescType, desc), max_id + 1)
    push!(li.var_defs, vd)
    return toRHSVar(vd)
end

"""
Add the given variable to be an escaping variable. Will error if
the given variable already exists locally, or is a GenSym
(or SSAValue in 0.5) or parameter.  Return the variable itself.
"""
function addEscapingVariable(s :: Symbol, typ :: Type, desc, li :: LambdaVarInfo)
    x = addLocalVariable(s, typ, desc, li)
    setEscapingVariable(x, li)
    return x
end

"""
Set the given variable to be escaping. Will error if the given variable
doesn't already exist, or is a GenSym (or SSAValue in 0.5) or parameter.
Return true if it is not already escaping, or false otherwise.
"""
function setEscapingVariable(x :: Union{Symbol,RHSVar}, li :: LambdaVarInfo)
    for vd in li.var_defs
        if matchVarDef(x, vd)
            @assert (vd.name != emptyVarName) "Variable " * string(x) * " does not have a Symbol name, and cannot be made escaping"
            @assert (!in(vd.name, li.input_params)) "Variable " * string(x) * " is an input parameter, and cannot be made escaping"
            if !in(vd.name, li.escaping_vars)
                push!(li.escaping_vars, vd.name)
                return true
            else
                return false
            end
        end
    end
    error("Variable ", x, " is not found in ", li)
end

"""
Set a local variable to be not escaping. Will error if the given variable
doesn't already exist, or is a GenSym (or SSAValue in 0.5) or parameter.
Return false if it is not already escaping, or true otherwise.
"""
function unsetEscapingVariable(x :: Union{Symbol,RHSVar}, li :: LambdaVarInfo)
    for vd in li.var_defs
        if matchVarDef(x, vd)
            @assert (vd.name != emptyVarName) "Variable " * string(x) * " does not have a Symbol name, and cannot be made escaping"
            @assert (!in(vd.name, li.input_params)) "Variable " * string(x) * " is an input parameter, and cannot be made escaping"
            escapes = Symbol[]
            found = false
            for y in li.escaping_vars
                if y == vd.name
                    found = true
                else
                    push!(escapes, y)
                end
            end
            li.escaping_vars = found ? escapes : li.escaping_vars
            return found
        end
    end
    error("Variable ", x, " is not found in ", li)
end

"""
Add a new temporary (GenSym in 0.4 or SSAValue in 0.5) variable to the LambdaVarInfo.
Returns the newly added variable.
"""
function addTempVariable(typ, li :: LambdaVarInfo)
    max_id = -1
    for vd in li.var_defs
        if vd.name == emptyVarName && vd.id > max_id
            max_id = vd.id
        end
    end
    vd = VarDef(typ, convert(DescType, ISASSIGNEDONCE | ISASSIGNED), max_id + 1)
    push!(li.var_defs, vd)
    return toRHSVar(vd)
end

function findAndFixAliasing!(expr :: ANY, fix_aliases = true, AstWalkFunc = nothing)
  function find_alias(expr :: ANY, object_set, top_level_number :: Int64, is_top_level :: Bool, read :: Bool)
    if isa(expr, Expr)
        oid = object_id(expr)
        if in(oid, object_set)
            @dprintln(2, "findAndFixAliasing found aliased Expr = ", expr)
            if fix_aliases
                return deepcopy(expr)
            else
                return expr
            end
        else
            push!(object_set, oid)
        end
    end
    return CompilerTools.AstWalker.ASTWALK_RECURSE
  end

  if expr == nothing
    return nothing
  end

  expr_set = Set()
  @dprintln(3, "findAndFixAliasing!: ", expr, " AstWalkFunc = ", AstWalkFunc)
  if isa(expr,Array)
    for i = 1:length(expr)
      if AstWalkFunc == nothing
        expr[i] = CompilerTools.AstWalker.AstWalk(expr[i], find_alias, expr_set)
      else
        expr[i] = AstWalkFunc(expr[i], find_alias, expr_set)
      end
    end
  else
    if AstWalkFunc == nothing
      expr = CompilerTools.AstWalker.AstWalk(expr, find_alias, expr_set)
    else
      expr = AstWalkFunc(expr, find_alias, expr_set)
    end
  end
  return expr
end

"""
Replace the symbols in an expression "expr" with those defined in the
dictionary "dict".  Return the result expression, which may share part of the
input expression, and the input "expr" may be modified inplace and shall not be used
after this call.
"""
function replaceExprWithDict!(expr :: ANY, dict :: Dict{LHSVar, Any}, AstWalkFunc = nothing)
  function update_sym(expr :: LambdaInfo, dict, top_level_number :: Int64, is_top_level :: Bool, read :: Bool)
    expr
  end
  function update_sym(expr :: ANY, dict, top_level_number :: Int64, is_top_level :: Bool, read :: Bool)
    if isfunctionhead(expr)
        linfo, body = lambdaToLambdaVarInfo(expr)
        new_dict = Dict{LHSVar,Any}()
        for (k,v) in dict
            if !isLocalVariable(k, linfo) && isEscapingVariable(k, linfo)
                new_v = toLHSVar(v,linfo)
                if isa(new_v, Symbol)
                    if !isLocalVariable(new_v, linfo)
                        @dprintln(3, "replaceExprWithDict! nested lambda replace ", k, " with ", v)
                        new_dict[k] = v
                        typ = getType(k, linfo)
                        desc = getDesc(k, linfo)
                        addEscapingVariable(new_v, typ, desc, linfo)
                    else
                        @dprintln(3, "replaceExprWithDict! nested lambda cannot replace ", k, " with ", v, " which is not Symbol")
                    end
                else
                    push!(new_dict, k, v)
                    @dprintln(3, "replaceExprWithDict! nested lambda replace ", k, " with ", v)
                end
            else
                @dprintln(3, "replaceExprWithDict! nested lambda skip replacing ", k, " because it is not escaping")
            end
        end
        replaceExprWithDict!(body, new_dict, AstWalkFunc)
        return expr
    elseif isa(expr, LHSVar)
      if haskey(dict, expr)
        return deepcopy(toLHSVar(dict[expr]))
      end
    elseif isa(expr, TypedVar)
      lhsVar = toLHSVar(expr)
      if haskey(dict, lhsVar)
        new_val = dict[lhsVar]
        if isa(new_val, TypedVar)
            return new_val
        elseif isa(new_val, LHSRealVar)
            return toTypedVar(new_val, expr.typ)
        else
            return deepcopy(toLHSVar(new_val))
        end
      end
    end
    return CompilerTools.AstWalker.ASTWALK_RECURSE
  end

  if expr == nothing
    return nothing
  end

  findAndFixAliasing!(expr, true, AstWalkFunc)

  @dprintln(3, "replaceExprWithDict!: ", expr, " dict = ", dict, " AstWalkFunc = ", AstWalkFunc)
  if isa(expr,Array)
    for i = 1:length(expr)
      if AstWalkFunc == nothing
        expr[i] = CompilerTools.AstWalker.AstWalk(expr[i], update_sym, dict)
      else
        expr[i] = AstWalkFunc(expr[i], update_sym, dict)
      end
    end
  else
    if AstWalkFunc == nothing
      expr = CompilerTools.AstWalker.AstWalk(expr, update_sym, dict)
    else
      expr = AstWalkFunc(expr, update_sym, dict)
    end
  end
  return expr
end

"""
Merge "inner" LambdaVarInfo into "outer", and "outer" is changed as result.  Note
that the input_params, static_parameter_names, and escaping_defs of "outer" do
not change, other fields are merged. The GenSyms in "inner" will need to adjust
their indices as a result of this merge. We return a dictionary that maps from
old GenSym to new GenSym for "inner", which can be used to adjust the body Expr
of "inner" lambda using "replaceExprWithDict" or "replaceExprWithDict!".

The "extraDesc" can be used to add extra desc flag to all local variables to
the inner lambda.
"""
function mergeLambdaVarInfo(outer :: LambdaVarInfo, inner :: LambdaVarInfo, extra_desc = 0)
    dict = Dict{LHSVar, Any}()
    @dprintln(3,"outer = ", outer)
    @dprintln(3,"inner = ", inner)
    for vd in inner.var_defs
        if vd.name == emptyVarName
            # temporary var, add to outer
            dict[toLHSVar(vd)] = addTempVariable(vd.typ, outer)
        else
            # named var
            # first, all input parameters are skipped
            name = vd.name
            typ  = vd.typ
            desc = vd.desc
            if !in(name, inner.input_params)
                if in(name, inner.escaping_vars)
                    # if escaping, must exist as variable in outer
                    @assert (isLocalVariable(name, outer) || isEscapingVariable(name, outer)) "Variable " * string(vd.name) * " is not found in outer info " * string(outer)
                    # clear ISCAPTURED and ISASSIGNEDBYINNERFUNCTION is not correct since this variable could still be escaping in other context
                    # desc = desc & (~ ISASSIGNEDBYINNERFUNCTION)
                    # desc = desc & (~ ISCAPTURED)
                    # setDesc(name, desc, outer)
                    # must replace escaping variables too, because the actual variable representation would be different in 0.5 between outer and inner
                    dict[toLHSVar(vd)] = toRHSVar(lookupLHSVarByName(name, outer), outer)
                else
                    # fresh name if there is name conflict
                    if isLocalVariable(name, outer) || isEscapingVariable(name, outer)
                        name = gensym(string(name))
                    end
                    dict[toLHSVar(vd)] = addLocalVariable(name, typ, desc | extra_desc, outer)
                end
            end
        end
    end
    return dict
end

function lambdaToLambdaVarInfo(lambda :: LambdaInfo)
    ast = Base.uncompressed_ast(lambda)
    linfo = LambdaVarInfo(lambda, ast)
    body = getBody(ast, linfo.return_type)
    return linfo, body
end

function lambdaToLambdaVarInfo(lambda :: Expr)
    @assert (lambda.head == :lambda) "Expect a :lambda Expr, but got " * string(lambda)
    linfo = LambdaVarInfo(lambda)
    body = getBody(lambda)
    return linfo, body
end

"""
Returns the type of the lambda as stored in LambdaVarInfo "li" and as extracted during lambdaToLambdaVarInfo.
"""
function getReturnType(li :: LambdaVarInfo)
  return li.return_type
end

"""
set the return type of LambdaVarInfo "li".
"""
function setReturnType(ret_typ, li :: LambdaVarInfo)
  @assert (ret_typ != nothing)
  li.return_type = ret_typ
end

if VERSION > v"0.5.0-dev+3260"
function lambdaTypeinf(lambda :: LambdaInfo, typs; optimize = true)
    throw(string("Force inference LambdaInfo in 0.5 not yet supported"))
end
function lambdaTypeinf(lambda :: Function, typs; optimize = true)
  lambdaTypeinf(typeof(lambda), typs, optimize = optimize)
end
else
"""
Force type inference on a LambdaStaticData object.
Return both the inferred AST that is to a "code_typed(Function, (type,...))" call,
and the inferred return type of the input method.
"""
function lambdaTypeinf(lambda :: LambdaStaticData, typs; optimize = true)
  types::Any = to_tuple_type(Tuple{typs...})
  if optimize
    (tree, ty) = Core.Inference.typeinf(lambda, types, Core.svec(), lambda, true, true)
  else
    (tree, ty) = Core.Inference.typeinf_uncached(lambda, types, Core.svec(), optimize = false)
  end
  lambda.ast = tree
  return lambda, ty
end

function lambdaTypeinf(lambda :: Function, typs; optimize = true)
  m = methods(lambda, (typs...))
  assert(length(m) > 0)
  lambdaTypeinf(m[1].func.code, typs, optimize = optimize)
end
end

function lambdaTypeinf(ftyp :: Type, typs; optimize = true)
    typ = Tuple{ftyp, typs...}
    types::Any = to_tuple_type(typ)
    env = Core.Inference.svec()
#    println("lambdaTypeinf typeof(ftyp) = ", typeof(ftyp))
#    println("lambdaTypeinf typeof(ftyp.name) = ", typeof(ftyp.name))
#    println("lambdaTypeinf typeof(ftyp.name.mt) = ", typeof(ftyp.name.mt))
#    println("lambdaTypeinf typeof(ftyp.name.mt.defs) = ", typeof(ftyp.name.mt.defs))
#    println("ftyp = ", ftyp)
#    println("ftyp.name = ", ftyp.name)
#    println("ftyp.name.mt = ", ftyp.name.mt)
#    println("ftyp.name.mt.defs = ", ftyp.name.mt.defs)
if VERSION > v"0.5.0-dev+3260"
#    println("ftyp = ", ftyp, " typs = ", typs, " methods = ", Base.methods(ftyp))
    meth = Base._methods_by_ftype(typ, -1)
    for m in meth
      if m[1] <: Tuple && m[1].parameters[2:end] == typ.parameters[2:end]
#        println("meth in _methods ", meth)
        m = meth[1]
        lambda = Core.Inference.func_for_method_checked(m[3], types)
        (tree, ty) = Core.Inference.typeinf_uncached(lambda, m[1], m[2], optimize = optimize)
#        println("lambdaTypeinf typeof(tree) = ", typeof(tree), " typeof(ty) = ", typeof(ty))
        return tree, ty
      end
    end
    error("Expected one method from call to Base._methods in lambdaTypeinf to match type ", typ, ", but none: ", meth)
else
    lambda = Core.Inference.func_for_method(ftyp.name.mt.defs, typs, env)
#    println("lambdaTypeinf typeof(lambda) = ", typeof(lambda))
    (tree, ty) = Core.Inference.typeinf_uncached(lambda, types, Core.svec(), optimize = optimize)
#    println("lambdaTypeinf typeof(tree) = ", typeof(tree), " typeof(ty) = ", typeof(ty))
    lambda.ast = tree
    return lambda, ty
end
end

if VERSION > v"0.5.0-dev+3260"
function LambdaVarInfoToLambda(li :: LambdaVarInfo, body::Array{Any,1}, AstWalkFunc = nothing)
    @assert (length(li.escaping_vars)==0) "LambdaVarInfo has escaping variables, and cannot be converted back to LambdaInfo " * string(li)
    lambda = li.orig_info
    @assert (lambda != nothing) "li.orig_info is not found " * string(li)
    @dprintln(3, "LambdaVarInfoToLambda, LambdaVarInfo = ", li, " body = ", body)
    nparams = length(li.input_params)
    dict = Dict{LHSVar,Any}()
    vars = Any[nothing for i = 1:length(li.var_defs)]
    ssas = Any[nothing for i = 1:length(li.var_defs)]
    vars_pending = VarDef[]
    ssas_pending = VarDef[]
    # parameters come first
    nvars = 0
    nssas = 0
    for vd in li.var_defs
      dprintln(3, "Processing ", vd)
      if vd.name == Symbol("#self#")
        dprintln(3, "Found self, putting in vars[1]")
        assert(vars[1] == nothing)
        vars[1] = vd
        nvars += 1
      elseif vd.name != emptyVarName
        dprintln(3, "Found non-empty varname")
        nvars += 1
        j = findfirst(li.input_params, vd.name)
        if j > 0
          # parameter
          assert(vars[j + 1] == nothing)
          vars[j+1] = vd
          dprintln(3, "is a param, putting in slot ", j+1)
        else
          dprintln(3, "is NOT a param ", vd.id, " ", nparams+1)
          # non-parameter variable
          if vd.id > nparams+1 && vd.id <= length(vars) && [vd.id] == nothing
            vars[vd.id] = vd
          else # place already taken
            push!(vars_pending, vd)
          end
        end
      else # ssavalues
        dprintln(3, "Found ssavalue")
        nssas += 1
        if ssas[vd.id+1] == nothing
          ssas[vd.id+1] = vd
        else
          push!(ssas_pending, vd)
        end
      end
    end
    lambda.slotnames = Array(Any, nvars)
    lambda.slottypes = Array(Any, nvars)
    lambda.slotflags = Array(DescType, nvars)
    @dprintln(3, "nvars = ", nvars)
    @dprintln(3, "vars = ", vars)
    @dprintln(3, "vars_pending = ", vars_pending)
    @dprintln(3, "ssas = ", ssas)
    @dprintln(3, "ssas_pending = ", ssas_pending)
    j = 1
    for i = 1:nvars
      if vars[i] == nothing
        assert(j <= length(vars_pending))
        vars[i] = vars_pending[j]
        j += 1
      end
      vd = vars[i]
      if vd.id != i
        # id mismatch, must replace
        dict[toLHSVar(vd)] = TypedSlot(i, vd.typ)
        vd.id = i
      end
      li.var_defs[i] = vd
      lambda.slotnames[i] = vd.name
      lambda.slottypes[i] = vd.typ
      lambda.slotflags[i] = vd.desc
    end
    @dprintln(3, "lambda.nargs = ", lambda.nargs, " nparams = ", nparams)
    lambda.nargs = nparams + 1
    # Julia may try to add a Const to this array during an optimization pass and so this can't be Array of Type.
    lambda.ssavaluetypes = Array(Any, nssas)
    j = 1
    for i = 1:nssas
      if ssas[i] == nothing
        assert(j <= length(ssas_pending))
        ssas[i] = ssas_pending[j]
        j += 1
      end
      vd = ssas[i]
      if vd.id != i - 1
        dict[toLHSVar(vd)] = GenSym(i - 1)
        vd.id = i - 1
      end
      li.var_defs[nvars + i] = vd
      lambda.ssavaluetypes[i] = vd.typ
    end
    lambda.rettype = li.return_type
    body = replaceExprWithDict!(body, dict, AstWalkFunc)
    @dprintln(3, "body = ", body)
    #lambda.code =  ccall(:jl_compress_ast, Any, (Any,Any), lambda, body)
    lambda.code =  body
    @dprintln(3, "lambda.slotnames = ", lambda.slotnames)
    @dprintln(3, "lambda.slottypes = ", lambda.slottypes)
    @dprintln(3, "lambda.slotflags = ", lambda.slotflags)
    @dprintln(3, "lambda.ssavaluetypes = ", lambda.ssavaluetypes)
    @dprintln(3, "LambdaVarInfoToLambda, lambda = ", lambda)
    return lambda
end

else # if VERSION

function createMeta(li::LambdaVarInfo)
    locals = Any[]
    escapes = Any[]
    max_id = -1
    count = 0
    for vd in li.var_defs
        if vd.name == emptyVarName && vd.id > max_id
            max_id = vd.id
            count+=1
        end
    end
    gensyms = Array(Any, max_id + 1)
    for vd in li.var_defs
        if vd.name == emptyVarName
            gensyms[vd.id + 1] = vd.typ
        elseif in(vd.name, li.escaping_vars)
            push!(escapes, Any[vd.name, vd.typ, vd.desc])
        else
            push!(locals, Any[vd.name, vd.typ, vd.desc])
        end
    end
    meta = Any[locals, escapes, gensyms]
    if isa(li.static_parameter_names, Array)
        push!(meta, li.static_parameter_names)
    end
    return meta
end

"""
Convert our internal storage format, LambdaVarInfo, back into a lambda expression.
This takes a LambdaVarInfo and a body as input parameters.
This body can be a body expression or you can pass "nothing" if you want but then you will probably need to set the body in args[3] manually by yourself.
"""
function LambdaVarInfoToLambda(li :: LambdaVarInfo, body::Array{Any,1}, AstWalkFunc = nothing)
  expr = Expr(:body)
  expr.args = body
  if li.return_type != nothing
      expr.typ = li.return_type
  end
  return Expr(:lambda, getInputParametersAsExpr(li), createMeta(li), expr)
end

end # if VERSION

function LambdaVarInfoToLambda(li :: LambdaVarInfo, body_expr :: Expr, AstWalkFunc = nothing)
  assert(body_expr.head == :body)
  LambdaVarInfoToLambda(li, body_expr.args, AstWalkFunc)
end

"""
Update the descriptor part of the VarDef dealing with whether the variable is assigned or not in the function.
Takes the LambdaVarInfo and a dictionary that maps symbols names to the number of times they are statically assigned in the function.
"""
function updateAssignedDesc(li :: LambdaVarInfo, symbol_assigns :: Dict{Symbol,Int})
  # For each VarDef
  for vd in li.var_defs
    # If that VarDef's symbol is in the dictionary.
    if vd.name != emptyVarName && haskey(symbol_assigns, vd.name)
      # Get how many times the symbol is assigned to.
      num_assigns = symbol_assigns[vd.name]
      # Remove the parts of the descriptor dealing with the number of assignments.
      vd.desc = vd.desc & (~ (ISASSIGNED | ISASSIGNEDONCE))
      if num_assigns > 1
        # If more than one assignment then OR on ISASSIGNED.
        vd.desc = vd.desc | ISASSIGNED
      elseif num_assigns == 1
        # If exactly one assignment then OR on ISASSIGNED and ISASSIGNEDONCE
        vd.desc = vd.desc | ISASSIGNED | ISASSIGNEDONCE
      end
    end
  end
end

if VERSION > v"0.5.0-dev+3260"
function getBody(lambda :: LambdaInfo)
  return getBody(Base.uncompressed_ast(lambda), lambda.rettype)
end
else
function getBody(lambda :: LambdaStaticData)
  return getBody(Base.uncompressed_ast(lambda))
end
end

function getBody(lambda :: Expr, rettype = nothing)
  assert(lambda.head == :lambda)
  body = lambda.args[3]
  if rettype != nothing
    body.typ = rettype
  end
  return body
end

function getBody(lambda :: Array{Any,1}, rettype)
  body = Expr(:body)
  body.args = lambda
  body.typ = rettype
  return body
end

function prependStatements(body :: Expr, stmts:: Array{Any, 1})
  assert(body.head == :body)
  body.args = [stmts; body.args]
  return body
end

"""
Returns an array of Symbols corresponding to those parameters to the method that are going to be passed by reference.
In short, isbits() types are passed by value and !isbits() types are passed by reference.
"""
function getRefParams(li :: LambdaVarInfo)
  ret = LHSVar[]

  for x in li.input_params
    typ = getType(x, li)
    @dprintln(3,"x = ", x, " type = ", typ)
    if !isbits(typ)
        push!(ret, lookupLHSVarByName(x, li))
    end
  end

  return ret
end

function getArrayParams(li :: LambdaVarInfo)
  ret = LHSVar[]

  input_vars = li.input_params
  var_types  = li.var_defs

  for x in li.input_params
      typ = getType(x, li)
      if isArrayType(typ)
        push!(ret, lookupLHSVarByName(x, li))
      end
  end

  return ret
end

"""
Get the name of a local variable, either as Symbol or GenSym
"""
function lookupVariableName(x::Union{Symbol,RHSVar}, li::LambdaVarInfo)
    for vd in li.var_defs
        if matchVarDef(x, vd)
            return vd.name == emptyVarName ? toLHSVar(vd) : vd.name
        end
    end
    error("Variable ", x, " is not found in ", li)
end

"""
Lookup a Symbol in local variables, and return that variable (as LHSVar).
"""
function lookupLHSVarByName(x::Symbol, li::LambdaVarInfo)
    for vd in li.var_defs
        if matchVarDef(x, vd)
            return toLHSVar(vd)
        end
    end
    error("Variable ", x, " is not found in ", li)
end

"""
Holds symbols and gensyms that are seen in a given AST when using the specified callback to handle non-standard Julia AST types.
"""
type CountVariableState
  used :: Set{LHSVar}

  function CountVariableState()
    new(Set{LHSVar}())
  end
end

"""
Adds symbols and gensyms to their corresponding sets in CountSymbolState when they are seen in the AST.
"""
function count_variables(x::RHSVar,
                       state::CountVariableState,
                       top_level_number,
                       is_top_level,
                       read)
    push!(state.used, toLHSVar(x))
    return CompilerTools.AstWalker.ASTWALK_RECURSE
end

function count_variables(x::ANY,
                       state::CountVariableState,
                       top_level_number,
                       is_top_level,
                       read)
    return CompilerTools.AstWalker.ASTWALK_RECURSE
end

function countVariables(body, AstWalkFunc = nothing)
  css = CountVariableState()
  if AstWalkFunc == nothing
    CompilerTools.AstWalker.AstWalk(body, count_variables, css)
  else
    AstWalkFunc(body, count_variables, css)
  end
  return css.used
end

"""
Eliminates unused symbols from the LambdaVarInfo var_defs.
Takes a LambdaVarInfo to modify, the body to scan using AstWalk and an optional callback to AstWalk for custom AST types.
"""
function eliminateUnusedLocals!(li :: LambdaVarInfo, body, AstWalkFunc = nothing)
  used = countVariables(body, AstWalkFunc)
  @dprintln(3,"used = ", used)
  var_defs = VarDef[]
  escapings = Symbol[]
  for vd in li.var_defs
    if !in(vd.name, li.input_params) && vd.name != Symbol("#self#")
      lhsVar = toLHSVar(vd)
      if !in(lhsVar, used)
        continue
      end
    end
    if in(vd.name, li.escaping_vars)
        push!(escapings, vd.name)
    end
    push!(var_defs, vd)
  end
  li.escaping_vars = escapings
  li.var_defs = var_defs
  dict = consolidateLambdaVarInfo!(li)
  @dprintln(3,"dict = ", dict)
  body = replaceExprWithDict!(body, dict, AstWalkFunc)
  @dprintln(3,"updated body = ", body)
  return body
end

function isVarArgParameter(s::Symbol, linfo::LambdaVarInfo)
    @assert (in(s, linfo.input_params))
    in(s, linfo.vararg_params)
end

function stripCaptureFlag(linfo :: LambdaVarInfo)
  for i = 1:length(linfo.var_defs)
    linfo.var_defs[i].desc &= (~ISCAPTURED)
  end
end

function getStaticParameterValue(i :: Int, linfo :: LambdaVarInfo)
  @assert (i >= 1 && i <= length(linfo.static_parameter_values)) string("getStaticParameterValue out-of-range: ", i, " in ", linfo.static_parameter_values)
  linfo.static_parameter_values[i]
end
end


text:
includet("barnesHut.jl")
includet("visualization.jl")
global secs = 3

# this is a helper function called by other specific simulator functions
function basicSimulation(C::Node)
    global timesInSec = 0.00

    # keep a local dataframe to update
    df = writeStatsToFrame(C)
    # write this stats to a csv file
    CSV.write("data/"*string(timesInSec)*".csv", df)

    frames = secs*24
    # simulate new positions for "frames" number of frames
    for i in 1:frames
        C = makeTreeFromFrame(df)
        creep(C)
        df = writeStatsToFrame(C)
        global timesInSec += frameRate
        CSV.write("data/"*string(timesInSec)*".csv", df)
    end
    plotting()
end

# this function is a helper function that returns minBounds and maxBounds based on the type
function getBoundsBasedOnTypes(type::String, vector::Array{Float64,1}, minX::Float64, maxX::Float64)
    minBounds = SVector{3, Float64}(zeros(Float64,3))
    maxBounds = SVector{3, Float64}(zeros(Float64,3))
    if type == "random"
        minBounds = SVector{3, Float64}(minX,minX,minX)
        maxBounds = SVector{3, Float64}(maxX,maxX,maxX)
    elseif type == "disk"
        normal = vector
        # calcuate the possible max and min value of z for initializing a tree
        if normal[3] == 0 && normal[1] != 0 && normal[2] != 0
            minBounds = SVector{3, Float64}(-maxX,-maxX,minX)
            maxBounds = SVector{3, Float64}(maxX,maxX,maxX)
        else
            minZ = min(-minX*normal[1]-minX*normal[2],-maxX*normal[1]-maxX*normal[2])/normal[3]
            maxZ = max(-minX*normal[1]-minX*normal[2],-maxX*normal[1]-maxX*normal[2])/normal[3]
            if minZ == maxZ
                minZ = -maxX
                maxZ = maxX
            end
            minBounds = SVector{3, Float64}(minX,minX,minZ)
            maxBounds = SVector{3, Float64}(maxX,maxX,maxZ)
        end
    elseif type == "line"
        # calcuate the possible max and min value of y and z for initializing a tree
        dir = vector
        minY = min(minX/dir[1] , maxX/dir[1])*dir[2]
        maxY = max(minX/dir[1] , maxX/dir[1])*dir[2]
        minZ = min(minX/dir[1] , maxX/dir[1])*dir[3]
        maxZ = max(minX/dir[1] , maxX/dir[1])*dir[3]
        minBounds = SVector{3, Float64}(minX,minY,minZ)
        maxBounds = SVector{3, Float64}(maxX,maxY,maxZ)
    end
    return minBounds, maxBounds
end

# this function gets random position based on the type of the galaxy
function getRandomPosition(type::String, num, vector::Array{Float64,1}, minn::Int64,maxx::Int64)
    if type == "random"
        position = [rand(minn:maxx)+rand(Float64), rand(minn:maxx)+rand(Float64) ,rand(minn:maxx)+rand(Float64)]
    elseif type == "disk"
        normal = vector
        x = 0.0
        y = 0.0
        z = 0.0
        if normal[3]==0 && normal[1]!= 0 && normal[2]!= 0
            x = rand(minn:maxx)+rand(Float64)
            y = (0-normal[1]*x)/normal[2]
            z = rand(minn:maxx)+rand(Float64)
        elseif normal[3]==0 && normal[1]== 0 && normal[2]!= 0
            x = rand(minn:maxx)+rand(Float64)
            y = (0-normal[1]*x)/normal[2]
            z = rand(minn:maxx)+rand(Float64)
        elseif normal[2]==0 && normal[1]!= 0 && normal[3]!= 0
            x = rand(minn:maxx)+rand(Float64)
            y = rand(minn:maxx)+rand(Float64)
            z = (0-normal[1]*x)/normal[3]
        elseif normal[1]==0 && normal[2]!= 0 && normal[3]!= 0
            x = rand(minn:maxx)+rand(Float64)
            y = rand(minn:maxx)+rand(Float64)
            z = (0-normal[2]*y)/normal[3]
        else
            x = rand(minn:maxx)+rand(Float64)
            y = rand(minn:maxx)+rand(Float64)
            z = (0-x*normal[1]-y*normal[2])/normal[3]
        end
        position = [x,y,z]
    elseif type == "line"
        dir = vector
        x = 0.0
        y = 0.0
        z = 0.0
        if dir[1]!=0
            x = rand(minn:maxx)+rand(Float64)
            y = dir[2]*x/dir[1]
            z = dir[3]*x/dir[1]
        elseif dir[2]!=0
            y = rand(minn:maxx)+rand(Float64)
            x = dir[1]*y/dir[2]
            z = dir[3]*y/dir[2]
        end
        position = [x,y,z]
    end
    return position
end
# this function is a helper function that adds bodies based on different system
function addBodies(C::Node, type::String, num::Int64, vector::Array{Float64,1}, minn::Int64,maxx::Int64)
    mid = (maxx - minn)/5
    if type == "random"
        # add random bodies
        for i in 1:Int(floor(0.9num))
            mass = exp(randn())*M
            position = getRandomPosition("random", num, [0.0,0.0,0.0], minn,maxx)
            newBody = Body(i, mass, position, SVector{3, Float64}(zeros(Float64,3)), SVector{3, Float64}(zeros(Float64,3)))
            insertBody(C,newBody)
        end
        for i in 1:(num - Int(floor(0.9num)))-1
            mass = (rand(1:12)+rand(Float64))*10M
            position = getRandomPosition("random", num, [0.0,0.0,0.0], Int64(ceil(minn+mid)),Int64(floor(maxx-mid)))
            newBody = Body(i, mass, position, SVector{3, Float64}(zeros(Float64,3)), SVector{3, Float64}(zeros(Float64,3)))
            insertBody(C,newBody)
        end
        mass = (rand(1:10)+rand(Float64))*(maxx-minn)*num*M
        position = rand()/100 .+C.centerOfMass
        newBody = Body(num, mass, position, SVector{3, Float64}(zeros(Float64,3)), SVector{3, Float64}(zeros(Float64,3)))
        insertBody(C,newBody)

    elseif type == "disk"
        for i in 1:Int(floor(0.9num))
            mass = exp(randn())*M
            position = getRandomPosition("disk", num, vector, minn,maxx)
            newBody = Body(i, mass, position, SVector{3, Float64}(zeros(Float64,3)), SVector{3, Float64}(zeros(Float64,3)))
            insertBody(C,newBody)
        end
        for i in 1:(num - Int(floor(0.9num)))-1
            mass = (rand(1:10)+rand(Float64))*10M
            position = getRandomPosition("disk", num, vector, Int64(ceil(minn+mid)),Int64(floor(maxx-mid)))
            newBody = Body(i, mass, position, SVector{3, Float64}(zeros(Float64,3)), SVector{3, Float64}(zeros(Float64,3)))
            insertBody(C,newBody)
        end

        mass = (rand(1:10)+rand(Float64))*(maxx-minn)*num*M

        position = rand()/100 .+C.centerOfMass

        newBody = Body(-1, mass, position, SVector{3, Float64}(zeros(Float64,3)), SVector{3, Float64}(zeros(Float64,3)))

        insertBody(C,newBody)

    elseif type == "line"
        for i in 1:Int(floor(0.9num))
            mass = exp(randn())*M
            position = getRandomPosition("line", num, vector, minn,maxx)
            newBody = Body(i, mass, position, SVector{3, Float64}(zeros(Float64,3)), SVector{3, Float64}(zeros(Float64,3)))
            insertBody(C,newBody)
        end
        for i in 1:(num - Int(floor(0.9num)))-1
            mass = (rand(1:10)+rand(Float64))*10M
            position = getRandomPosition("line", num, vector, Int64(ceil(minn+mid)),Int64(floor(maxx-mid)))
            newBody = Body(i, mass, position, SVector{3, Float64}(zeros(Float64,3)), SVector{3, Float64}(zeros(Float64,3)))
            insertBody(C,newBody)
        end
        mass = (rand(1:10)+rand(Float64))*100M
        position = getRandomPosition("line", num, vector, Int64(ceil(minn+2mid)),Int64(floor(maxx-2mid)))
        newBody = Body(num, mass, position, SVector{3, Float64}(zeros(Float64,3)), SVector{3, Float64}(zeros(Float64,3)))
        insertBody(C,newBody)
    end
end

# this is a random galaxy simulator similating a space of 1000 pc * 1000 pc *1000 pc
# optional argument: num -> number of bodies in the system
function randomGalaxy(num::Int64=100)
    minBounds, maxBounds = getBoundsBasedOnTypes("random",[0.0,0.0,0.0], -1.0, 1001.0)
    # making an empty tree
    C= Node(
           false, Array{Node,1}(undef, 8),0,
           Body(), 0.0, SVector{3, Float64}(zeros(3)),
           minBounds,maxBounds
           )
    # add random bodies
    addBodies(C,"random",num, [0.0,0.0,0.0],1,1000)
    # assign initial velocity for all the body
    caluculateInitialSpeed(C)

    # tuning for different cases
    if num== 2
        global timeScalar = nothing
        global timeScalar = 60*60*60*24*365*7e7
        global strengthOfInteraction = nothing
        global strengthOfInteraction = 3
    else
        global timeScalar = nothing
        global timeScalar = 60*60*60*24*365*7e14
        global strengthOfInteraction = nothing
        global strengthOfInteraction = 2e33
    end
    basicSimulation(C)

end

# this is a disk-shaped galaxy simulator similating a space of 1000 pc * 1000 pc *1000 pc
# optional arguments: normal-> the normal vector of the plane;
#       num -> number of bodies in the system
function diskGalaxy(normal::Array{Float64,1} = [1.0,1.0,0.0], num::Int64=100)
    if length(normal)!=3
        error("Normal vector must have 3 components")
    end

    # invert the sign to positive
    if normal[1] < 0 && normal[2] <0 && normal[3] <0
        normal = -normal
    end

    minBounds, maxBounds = getBoundsBasedOnTypes("disk",normal, -1.0, 1001.0)
    # making an empty tree
    C= Node(
           false, Array{Node,1}(undef, 8),0,
           Body(), 0.0, SVector{3, Float64}(zeros(3)),
           minBounds,maxBounds
           )
    # add bodies in a plane
    addBodies(C,"disk",num,normal,1,1000)
    # assign initial velocity for all the body
    changeRotatingPlane(normal)
    #caluculateInitialSpeed(C)
    caluculateInitialSpeed(C)
    if num== 2
        global timeScalar = nothing
        global timeScalar = 60*60*60*24*365*7e7
        global strengthOfInteraction = nothing
        global strengthOfInteraction = 3
    else
        global timeScalar = nothing
        global timeScalar = 60*60*60*24*365*7e14
        global strengthOfInteraction = nothing
        global strengthOfInteraction = 2e33
        # or 4e33
    end
    # addEnvironment(C,"disk",10*num,normal,1,1000)
    basicSimulation(C)
end

# this is a line-shaped galaxy simulator similating a space of 1000 pc * 1000 pc *1000 pc
# optional arguments: dir-> the direction of the line;
#                       num -> number of bodies on the line
function lineGalaxy(dir::Array{Float64,1}=[1.0,1.0,0.0], num::Int64=100)
    if length(dir)!=3
        error("Direction vector must have 3 components")
    end

    # invert the sign to positive
    if dir[1] < 0 && dir[2] <0 && dir[3] <0
        dir = -dir
    end

    minBounds, maxBounds = getBoundsBasedOnTypes("line",dir, -1.0, 1001.0)
    # making an empty tree
    C= Node(
           false, Array{Node,1}(undef, 8),0,
           Body(), 0.0, SVector{3, Float64}(zeros(3)),
           minBounds,maxBounds
           )
    # add bodies in a line
    addBodies(C,"line",num,dir,1,1000)
    # assign initial velocity for all the body
    caluculateInitialSpeed(C)
    if num== 2
        global timeScalar = nothing
        global timeScalar = 60*60*60*24*365*7e7
        global strengthOfInteraction = nothing
        global strengthOfInteraction = 3
    else
        global timeScalar = nothing
        global timeScalar = 60*60*60*24*365*5e5
        global strengthOfInteraction = nothing
        global strengthOfInteraction = 0.1
    end
    basicSimulation(C)
end

# this is a system of two galaxies simulator similating a space of 1000 pc * 1000 pc *1000 pc
# optional arguments: num1 -> number of bodies in the first galaxy,
#                     num2 -> number of bodies in the second galaxy,
#                     velocity1 -> the velocity of the first galaxy,
#                     velocity2 -> the velocity of the second galaxy,
#                     type1 -> the type of the first galaxy,
#                     type2 -> the type of the second galaxy,
#                     vector1 -> the normal vector or direction vector of the first galaxy,
#                     vector2 -> the normal vector or direction vector of the second galaxy]
function twoCollapseGalaxy(num1::Int64=100,num2::Int64=100, velocity1::SVector{3,Float64}=SVector{3,Float64}(0.0,0.0,0.0),
    velocity2::SVector{3,Float64}=SVector{3,Float64}(0.0,0.0,0.0), type1::String="disk",type2::String="disk",
    vector1::Array{Float64,1}=[1.0,1.0,0.0], vector2::Array{Float64,1}=[1.0,1.0,0.0])
    if length(velocity2)!=3 || length(velocity1)!=3
        error("Velocity vector(s) must have 3 components")
    elseif length(vector1) != 3 || length(vector2) != 3
        error("direcion/normal vector(s) must have 3 components")
    elseif type1 != "disk" && type1 != "random" && type1 != "line"
        error("Type1 not supported")
    elseif type2 != "disk" && type2 != "random" && type2 != "line"
        error("Type2 not supported")
    end

    global timeScalar = nothing
    global timeScalar = 60*60*60*24*365*1e14
    global strengthOfInteraction = nothing
    global strengthOfInteraction = 2e33

    changeRotatingPlane(vector1)
    #initialize a new empty tree
    minBounds1, maxBounds1 = getBoundsBasedOnTypes(type1,vector1, -10.0, 510.0)
    minBounds2, maxBounds2 = getBoundsBasedOnTypes(type2,vector2, 490.0, 1010.0)

    netMinBounds = SVector{3, Float64}(min(minBounds1[1],minBounds2[1]),
                                       min(minBounds1[2],minBounds2[2]),
                                       min(minBounds1[3],minBounds2[3]))

    netMaxBounds = SVector{3, Float64}(max(maxBounds1[1],maxBounds2[1]),
                                       max(maxBounds1[2],maxBounds2[2]),
                                       max(maxBounds1[3],maxBounds2[3]))

    # C1 is the first galaxy
    C1= Node(
           false, Array{Node,1}(undef, 8),0,
           Body(), 0.0, SVector{3, Float64}(zeros(3)),
           minBounds1, maxBounds1
           )

    # C2 is the first galaxy
    C2= Node(
          false, Array{Node,1}(undef, 8),0,
          Body(), 0.0, SVector{3, Float64}(zeros(3)),
          minBounds2, maxBounds2
          )

    # add bodies to the first tree
    if type1 == "disk"
        addBodies(C1,"disk",num1, vector1,0,500)
    elseif type1 == "random"
        addBodies(C1,"random",num1, [0.0,0.0,0.0],0,500)
    elseif type1 == "line"
        addBodies(C1,"line",num1, vector1,0,500)
    end

    if type2 == "disk"
        addBodies(C2,"disk",num2,vector2,500,1000)
    elseif type2 == "random"
        addBodies(C2,"random",num2, [0.0,0.0,0.0],500,1000)
    elseif type2 == "line"
        addBodies(C2,"line",num1, vector2,500,1000)
    end

    # calculate the velocity of the second galaxy
    C = mergeTwoNodes(C1, C2)

    velocity1 = sqrt(G*C.mass/dis(C.centerOfMass,C1.centerOfMass)/distanceRate) .*getDirectionOfVelocity(C.centerOfMass,C1.centerOfMass)
    velocity1 = SVector{3, Float64}(0.0,0.0,0.0)

    caluculateInitialSpeedPlus(C1,velocity1)

    velocity2 = sqrt(G*C.mass/dis(C.centerOfMass,C2.centerOfMass)/distanceRate) .*getDirectionOfVelocity(C.centerOfMass,C2.centerOfMass)
    velocity2 = SVector{3, Float64}(0.0,0.0,0.0)

    caluculateInitialSpeedPlus(C2,velocity2)
    C = mergeTwoNodes(C1, C2)

    basicSimulation(C)
end

function jupiterSring()
    minBounds, maxBounds = getBoundsBasedOnTypes("random",[0.0,0.0,0.0], 499.999999995, 500.000000005)
    # making an empty tree
    C= Node(
           false, Array{Node,1}(undef, 8),0,
           Body(), 0.0, SVector{3, Float64}(zeros(3)),
           minBounds,maxBounds
           )
    # add 100 small random bodies
    for i in 1:100
        mass = rand()*M*5e-17
        position = [500+rand()/250000000,500+rand()/250000000,500+rand()/250000000 ]
        newBody = Body(i, mass, position, SVector{3, Float64}(zeros(Float64,3)), SVector{3, Float64}(zeros(Float64,3)))
        insertBody(C,newBody)
    end

    # add jupitar
    newBody = Body(0, 0.001M, [500.0,500.0,500.0], SVector{3, Float64}(zeros(Float64,3)), SVector{3, Float64}(zeros(Float64,3)))
    insertBody(C,newBody)
    # assign initial velocity for all the body
    caluculateInitialSpeed(C)

    # tuning for different cases

    global timeScalar = nothing
    global timeScalar = 60*60*60
    global strengthOfInteraction = nothing
    global strengthOfInteraction = 1

    basicSimulation(C)

end

# this function add 10X random stars
function addEnvironment(C::Node, type::String, num::Int64, vector::Array{Float64,1}, minn::Int64,maxx::Int64)
    for i in 1:num
        mass = rand()*M
        position = getRandomPosition("disk", num, vector, minn,maxx)
        newBody = Body(i, mass, position, SVector{3, Float64}(zeros(Float64,3)), SVector{3, Float64}(zeros(Float64,3)))
        insertBody(C,newBody)
    end
    printSpeed(C)
    caluculateInitialSpeedForEnvironment(C, vector)

end
# this function changes theta (s/d)
function changeTheta(theta::Float64)
    global THETA = nothing
    global THETA = theta
end

# this function changes how long the simulation is
function changeDuration(factor::Float64)
    oldSecs =  secs
    global secs = nothing
    global secs = oldSecs*factor
    global secsToPlot = nothing
    global secsToPlot = oldSecs*factor

end

# this function changes how long the simulation is
function changeDuration(factor::Int64)
    oldSecs =  secs
    global secs = nothing
    global secs = oldSecs*factor
    global secsToPlot = nothing
    global secsToPlot = oldSecs*factor

end

function getDuration()
    return secsToPlot
end

function changeSpeed(factor::Float64)
    oldSpeed = timeScalar
    global timeScalar = nothing
    global timeScalar = oldSpeed*factor
end

function changeRotatingPlane(vector::Array{Float64,1})
    if length(vector) != 3
        error("vector must be 3-dimensional")
    end
    global rotatingPlane = nothing
    global rotatingPlane = vector
end

function changeEpsilon(e::Float64)
    global epsilon = nothing
    global epsilon = e
end

function checkEpsilon()
    return epsilon
end

function printSpeed(C::Node)
    if C.hasChildren && C.numOfChild != 1
        for i in 1:8
            printSpeed(C.children[i])
        end
    elseif C.hasChildren && C.numOfChild == 1
        println(C.body.velocity)
    end
end

ext: jl

max_stars_repo_name: KeKeJin/Barnes-Hut-Galaxy-Simulator

max_issues_repo_name: KeKeJin/Barnes-Hut-Galaxy-Simulator

max_forks_repo_name: KeKeJin/Barnes-Hut-Galaxy-Simulator


text:
using PyPlot, ForneyLab
import Distributions

function generateNoisyMeasurements( z_start::Vector{Float64},
                                    u::Vector{Float64},
                                    A::Matrix{Float64},
                                    b::Vector{Float64},
                                    Σz::Matrix{Float64},
                                    Σx::Matrix{Float64})
    # Simulate linear state-space model
    # z[t] = A*z[t-1] + b*u[t] + N(0,Σz)
    # x[t] = N(z[t],Σx)

    # Return noisy observations [x[1],...,x[n]]

    n = length(u)
    d = length(z_start)

    # Sanity checks
    (n>0) || error("u should contain at least one value")
    (d>0) || error("z_start has an invalid dimensionality")
    (size(A) == (d,d)) || error("Transition matrix A does not have correct dimensions")
    (length(b) == d) || error("b does not have the correct dimensionality")
    (size(Σz) == (d,d)) || error("Covariance matrix Σz does not have correct dimensions")
    (size(Σx) == (d,d)) || error("Covariance matrix Σx does not have correct dimensions")

    x = Vector[] # x will be a list of n vectors
    z = copy(z_start)
    chol_Σz = cholesky(Σz)
    chol_Σx = cholesky(Σx)
    proc_noise = convert(Matrix,chol_Σz) * randn(d)
    for i=1:n
        proc_noise = convert(Matrix,chol_Σz) * randn(d)
        z = A*z + b*u[i] + proc_noise
        obs_noise = convert(Matrix,chol_Σx) * randn(length(z))
        push!(x, z + obs_noise)
    end

    return x
end

function plotCartPrediction(prediction::ProbabilityDistribution{Multivariate, GaussianMeanVariance},
                            measurement::ProbabilityDistribution{Multivariate, GaussianMeanVariance},
                            corr_prediction::ProbabilityDistribution{Multivariate, GaussianMeanVariance})
    # Make a fancy plot of the Kalman-filtered cart position
    p = Distributions.Normal(mean(prediction)[1], var(prediction)[1])
    m = Distributions.Normal(mean(measurement)[1], var(measurement)[1])
    c = Distributions.Normal(mean(corr_prediction)[1], var(corr_prediction)[1])

    plot_range = range(mean(c)-8*sqrt(var(p)), mean(c)+8*sqrt(var(p)), length=300)
    PyPlot.figure(figsize=(15,5))
    bg_img = imread("figures/cart-bg.png")
    height = floor((plot_range[end] - plot_range[1])/3)
    imshow(bg_img, zorder=0, extent=[plot_range[1], plot_range[end], 0., height])
    plot(plot_range, height*Distributions.pdf.(p, plot_range), "r-")
    plot(plot_range, height*Distributions.pdf.(m, plot_range), "b-")
    plot(plot_range, height*Distributions.pdf.(c, plot_range), "g-")
    legend(["Prediction "*L"p(z[n]|z[n-1],u[n])",
            "Noisy measurement "*L"p(z[n]|x[n])",
            "Corrected prediction "*L"p(z[n]|z[n-1],u[n],x[n])"], prop=Dict("size"=>14), loc=1)
    fill_between(plot_range, 0, height*Distributions.pdf.(p, plot_range), color="r", alpha=0.1)
    fill_between(plot_range, 0, height*Distributions.pdf.(m, plot_range), color="b", alpha=0.1)
    fill_between(plot_range, 0, height*Distributions.pdf.(c, plot_range), color="g", alpha=0.1)
    xlim([plot_range[1],plot_range[end]]); ylim([0.,height])
    ax=gca()
    ax.yaxis.set_visible(false)
    xlabel("Position")
end

function plotCartPrediction2(predictive_mean, predictive_cov, mean, cov,measurement_position,measurement_cov)
    # Make a fancy plot of the Kalman-filtered cart position
    p = Distributions.Normal(predictive_mean, predictive_cov)
    m = Distributions.Normal(measurement_position, measurement_cov)
    c = Distributions.Normal(mean, cov)
    plot_range = range(mean-8*sqrt(predictive_cov), mean+8*sqrt(predictive_cov), length=300)
    PyPlot.figure(figsize=(15,5))
    bg_img = imread("figures/cart-bg.png")
    height = floor((plot_range[end] - plot_range[1])/3)
    imshow(bg_img, zorder=0, extent=[plot_range[1], plot_range[end], 0., height])
    plot(plot_range, height*Distributions.pdf.(p, plot_range), "r-")
    plot(plot_range, height*Distributions.pdf.(m, plot_range), "b-")
    plot(plot_range, height*Distributions.pdf.(c, plot_range), "g-")
    legend(["Prediction "*L"p(z[n]|z[n-1],u[n])",
            "Noisy measurement "*L"p(z[n]|x[n])",
            "Corrected prediction "*L"p(z[n]|z[n-1],u[n],x[n])"], prop=Dict("size"=>14), loc=1)
    fill_between(plot_range, 0, height*Distributions.pdf.(p, plot_range), color="r", alpha=0.1)
    fill_between(plot_range, 0, height*Distributions.pdf.(m, plot_range), color="b", alpha=0.1)
    fill_between(plot_range, 0, height*Distributions.pdf.(c, plot_range), color="g", alpha=0.1)
    xlim([plot_range[1],plot_range[end]]); ylim([0.,height])
    ax=gca()
    ax.yaxis.set_visible(false)
    xlabel("Position")
end

max_stars_repo_name: Yikeru/BMLIP

max_issues_repo_name: Yikeru/BMLIP

max_forks_repo_name: Yikeru/BMLIP


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule pprof_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("pprof")
JLLWrappers.@generate_main_file("pprof", UUID("cf2c5f97-e748-59fa-a03f-dda3c62118cb"))
end  # module pprof_jll


text:
struct FileFactor{T<:Unsigned} <: AbstractFactor{T}
    name::String
    length::Int64
    levels::AbstractVector{<:AbstractString}
    datapath::String
end

Base.length(factor::FileFactor{T}) where {T<:Unsigned} = factor.length

function slice(factor::FileFactor{T}, fromobs::Integer, toobs::Integer, slicelength::Integer) where {T<:Unsigned}
    slicelength = verifyslicelength(fromobs, toobs, slicelength)
    sizeT = sizeof(T)
    iostream = open(factor.datapath)
    seek(iostream, sizeT * (fromobs - 1))
    buffer = Vector{T}(slicelength)
    map(Seq(Vector{T}, (buffer, iostream, fromobs, toobs, slicelength), nextdatachunk), SubArray{T,1,Array{T,1},Tuple{UnitRange{Int64}},true}) do slice
        view(slice, 1:length(slice))
    end

    # len = toobs - fromobs + 1
    # slicelength = verifyslicelength(fromobs, toobs, slicelength)
    # data = Vector{UInt8}(sizeof(T) * len)
    # open(factor.datapath) do f
    #     seek(f, sizeof(T) * (fromobs - 1))
    #     readbytes!(f, data)
    # end
    # factordata = reinterpret(T, data)
    # slice(factordata, 1, len, slicelength)
end

function Base.map(factor::FileFactor, dataframe::AbstractDataFrame)
    dataframe[getname(factor)]
end

max_stars_repo_name: vikashdwd/JuML

max_issues_repo_name: vikashdwd/JuML

max_forks_repo_name: vikashdwd/JuML


text:
# Functions in modules

export
    CuFunction

const CuFunction_t = Ptr{Cvoid}

"""
    CuFunction(mod::CuModule, name::String)

Acquires a function handle from a named function in a module.
"""
struct CuFunction
    handle::CuFunction_t
    mod::CuModule

    "Get a handle to a kernel function in a CUDA module."
    function CuFunction(mod::CuModule, name::String)
        handle_ref = Ref{CuFunction_t}()
        @apicall(:cuModuleGetFunction, (Ptr{CuFunction_t}, CuModule_t, Ptr{Cchar}),
                                       handle_ref, mod, name)
        new(handle_ref[], mod)
    end
end

Base.unsafe_convert(::Type{CuFunction_t}, fun::CuFunction) = fun.handle

Base.:(==)(a::CuFunction, b::CuFunction) = a.handle == b.handle
Base.hash(fun::CuFunction, h::UInt) = hash(mod.handle, h)


text:
using GibbsTypePriors
using Test, Nemo

include("test_common_functions.jl")
include("test_Cnk.jl")
include("test_Cnk_variable_prec.jl")
include("test_Vnk.jl")
include("test_Vnk_variable_prec.jl")
include("test_Pkn.jl")
include("test_Pkn_variable_prec.jl")
include("test_Expect_Kn.jl")

max_stars_repo_name: konkam/GibbsTypePriors

max_issues_repo_name: konkam/GibbsTypePriors

max_forks_repo_name: konkam/GibbsTypePriors


text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

export drawrectlattice, drawhexneighbors, drawhexregion, drawhexrect, drawhexrectcolors
export hex3cluster, hex3RGB, hexRGBW, hex12RGB
export drawhex3RGB, drawhex12RGB

drawrectlattice() = Vis.drawcells(Repeat.RectangularBasis(),50.0,SMatrix{2,4,Int64}(0,0,0,1,1,0,1,1))

"""draw the 2 ring neighbors of the hex cell at coordinates (0,0)"""
drawhexneighbors() =  Vis.drawcells(Repeat.HexBasis1(),50,Repeat.neighbors(Repeat.HexBasis1,(0,0),2))

"""draw hex cell at coordinates (0,0) and the 1 and 2 ring neighbors"""
drawhexregion() = Vis.drawcells(Repeat.HexBasis1(),50,Repeat.region(Repeat.HexBasis1,(0,0),2))

"""draw hex cells that fit within a rectangular box centered at coordinates (0,0). Use fill color yellow."""
function drawhexrect()
    cells = Repeat.hexcellsinbox(2,2)
    Vis.drawcells(Repeat.HexBasis1(),50,cells,color = repeat(["yellow"],length(cells)))
end

"""draw hex cells that fit within a rectangular box centered at coordinates (0,0). Use random fill colors selected for maximum distinguishability."""
function drawhexrectcolors()
    cells = Repeat.hexcellsinbox(4,4)
    Vis.drawcells(Repeat.HexBasis1(),30,cells)
end

""" Create a LatticeCluser with three elements at (0,0),(-1,0),(-1,1) coordinates in the HexBasis1 lattice"""
function hex3cluster()
    clusterelts = SVector((0,0),(-1,0),(-1,1))
    eltlattice = HexBasis1()
    clusterbasis = LatticeBasis(( -1,2),(2,-1))
    return LatticeCluster(clusterbasis,eltlattice,clusterelts)
end

""" Create a ClusterWithProperties with three types of elements, R,G,B """
function hex3RGB()
    clusterelements = SVector((0,0),(-1,0),(-1,1))
    colors = [colorant"red",colorant"green",colorant"blue"]
    names = ["R","G","B"]
    eltlattice = HexBasis1()
    clusterbasis = LatticeBasis(( -1,2),(2,-1))
    lattice = LatticeCluster(clusterbasis,eltlattice,clusterelements)
    properties =  DataFrame(Color = colors, Name = names)
    return ClusterWithProperties(lattice,properties)
end

""" Create a ClusterWithProperties with four types of elements, R,G,B,W """
function hexRGBW()
    clusterelements = SVector((0,0),(-1,0),(-1,1),(0,-1))
    colors = [colorant"red",colorant"green",colorant"blue",colorant"white"]
    names = ["R","G","B","W"]
    eltlattice = HexBasis1()
    clusterbasis = LatticeBasis((0,2),(2,-2))
    lattice = LatticeCluster(clusterbasis,eltlattice,clusterelements)
    properties =  DataFrame(Color = colors, Name = names)
    return ClusterWithProperties(lattice,properties)
end

function hex12RGB()
    clusterelements = SVector(
        (-1,1),(0,1),
        (-1,0),(0,0),(1,0),
        (-1,-1),(0,-1),(1,-1),(2,-1),
        (0,-2),(1,-2),(2,-2)
    )
    red = colorant"red"
    grn = colorant"green"
    blu = colorant"blue"
    colors = [
        grn,blu,
        blu,red,grn,
        red,grn,blu,red,
        blu,red,grn
        ]
    names = [
        "G3","B0",
        "B2","R1","G2",
        "R0","G1","B1","R3",
        "B3","R2","G0"
    ]
    eltlattice = HexBasis3()
    clusterbasis = LatticeBasis((2,2),(-3,2))
    lattice = LatticeCluster(clusterbasis,eltlattice,clusterelements)
    properties =  DataFrame(Color = colors, Name = names)
    return ClusterWithProperties(lattice,properties)
end

""" draw 3 repeats of hex3RGB cluster """
drawhex3RGB() = Vis.draw(hex3RGB(),[0 1 0; 0 0 1])

""" draw 3 repeats of hex12RGB cluster """
drawhex12RGB() = Vis.draw(hex12RGB(),[0 1 0; 0 0 1])


text:
# --------------------------------------------------------------------------
# ACE.jl: Julia implementation of the Atomic Cluster Expansion
# Copyright (c) 2019 Christoph Ortner <christophortner0@gmail.com>
# All rights reserved.
# --------------------------------------------------------------------------

@testset "Cylindrical Coordinates" begin

@info("Testset Cylindrical Coordinates")
##

using StaticArrays, Test
using LinearAlgebra
using JuLIP.Testing: print_tf

using ACE
using ACE.Bonds: CylindricalCoordinateSystem, cylindrical,
                  cartesian, CylindricalCoordinates

##
for ntest = 1:7
   R = 1.0 .+ rand(SVector{3, Float64})
   R̂ = R/norm(R)
   o = R/2
   C = CylindricalCoordinateSystem(R, o)
   for mtest = 1:7
      r = rand(SVector{3, Float64}) .- 0.5
      c = cylindrical(C, r)
      r1 = cartesian(C, c)
      print_tf(@test r ≈ r1)
      cref = CylindricalCoordinates(c.cosθ, c.sinθ, c.r, - c.z)
      rref = r - 2 * dot(r-o, R̂) * R̂
      print_tf(@test rref ≈ cartesian(C, cref))
   end
end

##

end


text:
module DirectionFinding

using LinearAlgebra, UnicodePlots

export Antenna,
       AntennaArray,
       Position, PositionS,
       SignalSource,
       Simulation,
       origin,
       halfsine_train,
       lineararray, circulararray, cornerarray, randomarray,
       cbf,
       MUSIC,
       findpeaks,
       beamwidth

const C = 299.79e6  # speed of light

abstract type AbstractPosition end

struct Position <: AbstractPosition
    x :: Float64
    y :: Float64
    z :: Float64
end

Position() = Position(0, 0, 0)
Position(x, y) = Position(x, y, 0)
origin() = Position()

struct PositionS <: AbstractPosition
    r :: Float64
    θ :: Float64
    ϕ :: Float64
end

function tocartesian(p::PositionS)
    r , θ, ϕ = p.r, p.θ, p.ϕ
    return Position.(r.*cos.(ϕ)*sin.(θ), r.*sin.(ϕ)*sin.(θ), r.*cos.(θ))
end

function tospherical(p::Position)
    x, y, z = p.x, p.y, p.z
    r= sqrt.(x.^2 + y.^2 + z.^2)
    return PositionS.(r, acos.(z/r), atan.(y, x))
end

struct Antenna{R}
    position :: Position
    gain     :: R
end

Antenna(a::Position) = Antenna(a, isotropic)

Antenna(x::T, y::T, gain = isotropic) where {T<:Number} = Antenna(Position(x, y), gain)

function Base.getproperty(a::Antenna, s::Symbol)
    s === :x && return a.position.x
    s === :y && return a.position.y
    s === :z && return a.position.z
    getfield(a, s)
end

isotropic(θ, ϕ) = 1.0

struct AntennaArray
    array    :: Vector{Antenna}
    ref      :: Position
end

AntennaArray(av) = AntennaArray(av, origin())

Base.length(a::AntennaArray) = length(a.array)
Base.getindex(a::AntennaArray, args... ; kwargs...) = getindex(a.array, args..., kwargs...)

""" lineararray(M, d, ref ; gain, σ²) :: AntennaArray

Returns a planar, linear array of `M` antennas spaced `d` meters with reference
point `ref` (defaults to the origin). Each antenna has gain `gain` (defaults to
isotropic). σ² specifies the variance of a (Gaussian) random variation in antenna
position."""
function lineararray(M, d, ref = origin(); gain = isotropic, σ² = 0)
    s() = sqrt(σ²)*randn()
    return AntennaArray([Antenna(0+s(), dd+s(), gain)
                         for dd in range((-d/2)*(M-1), (d/2)*(M-1), length=M)], ref)
end

""" circulararray(M, r, ref ; gain, σ²) ::AntennaArray

Returns a planar, circular array with `M` antennas on a circle of radius `r`, centered
on reference point `ref`."""
function circulararray(M, r, ref = origin(); gain = isotropic, σ² = 0)
    s() = sqrt(σ²)*randn()
    angles = range(0, step = 2π/M, length=M)
    return AntennaArray([Antenna(x+s(), y+s(), gain)
                         for (x,y) = zip(r.*cos.(angles), r.*sin.(angles))], ref)
end

""" cornerarray(d, ref ; gain, σ² = 0)

Returns a planar array made up of four antennas in the corner of a square of side `d`
centered on reference point `ref`. """
function cornerarray(d, ref = origin(); gain = isotropic, σ² = 0)
    return(circulararray(4, d/sqrt(2), ref, gain=gain, σ²=σ²))
end

""" randomarray(M, σ², ref ; gain)

Returns a planar array with `M` antennas on (Gaussian) random positions around reference
point `ref`. The variance of the positions is `σ²`."""
function randomarray(M, σ² = 1, ref = origin() ; gain = isotropic)
    s() = sqrt(σ²)*randn()
    return AntennaArray([Antenna(ref.x+s(), ref.y+s(), gain) for _=1:M])
end

function Base.show(io::IO, a::AntennaArray)
    M = length(a.array)
    xcoord = [ant.x for ant in a.array]
    ycoord = [ant.y for ant in a.array]
    zcoord = [ant.z for ant in a.array]
    println(io, "An antenna array with $M elements.")
    for i in 1:M
        println(io, "   $i: x = $(xcoord[i]), y = $(ycoord[i]), z = $(zcoord[i])")
    end
    println("Reference point: x = $(a.ref.x), y = $(a.ref.y), z = $(a.ref.z)")
    # plot
    plt = scatterplot(xcoord, ycoord, marker = '📡')
    print(io, scatterplot!(plt, [a.ref.x], [a.ref.y], marker = 'R', color=:red))
end

"""Half-sine pulse at rate Rp with random amplitude."""
function halfsine_train(Rp)
    Tp = 1/Rp
    p(t) = sin(π*Rp*mod(t,Tp))
    a() = rand([-1,1])
    t -> Complex(a()*p(t), a()*p(t))
end

abstract type Source end

""" SignalSource(θ, sbb)

Returns a signal source with the following fields:

    `θ   :: Float64`  : Angle between a line orthogonal to the array's
                        reference point and the source.
    `sbb :: Function` : function that returns the transmitted lowpass
                        complex signal evaluated at snapshot instants `t`.  """
struct SignalSource{T, F} <: Source
    θ      :: T
    sbb    :: F
end

""" SignalSource(t)

Return the lowpass complex signal produced by the source `s` at time `t`.  """
(s::SignalSource)(t) = s.sbb(t)

"""Simulation

Fields:
    array :: an AbstractArray
    fc :: simulation carrier frequency
    N₀ :: noise spectral density in each receiver
    sources :: a tuple of the sources in the simulation
    phaseshifts : vector of functions `f[i](θ)` = phase shift of antenna `i`
                  relative to array reference point when source is at angle `θ`.  """
struct Simulation{T, S, M}
    array       :: T
    fc          :: Float64
    N₀          :: Float64
    sources     :: S
    phaseshifts :: Vector{M}  # array manifold
end

function Simulation(a, fc, N₀, s::Source...)
    λ = C/fc
    phaseshifts = [θ -> exp(1im*phaseshift(antenna, a.ref, θ, λ)) for antenna in a.array]
    Simulation(a, fc, N₀, s, phaseshifts)
end

""" projection(a, r, θ)

Returns the projection of antenna `a` on the wavefront of a source
at angle `θ` when it touches reference point `r`."""
function projection(a::Antenna, r::Position, θ)
    ϕ = θ + π/2
    if ϕ ≈ π/2 || ϕ ≈ 3π/2
        proj = Position(a.x, r.y)
    else
        m = tan(ϕ); n = -1.0/m;
    	x = (m*r.x-r.y-n*a.x+a.y)/(m-n)
    	y = m*(x-r.x)+r.y
		proj = Position(x, y)
    end
end

""" phaseshift(a, r, θ, λ)

Returns the phase shift seen by `a` relative to `r` when the
source is at angle `θ` and with wavelength `λ`.  """
function phaseshift(a::Antenna, r::Position, θ, λ)
    θ < 0 && (θ += 2π)
    p = projection(a, r, θ)
    # determine distance between antenna and projection
    d = sqrt((a.x-p.x)^2 + (a.y-p.y)^2)
    ps = d*2π/λ
    # Determine if `a` is in front of or behind `r` relative to the wavefront.
    # If in front, the phase shift is positive.
    if θ ≈ 0
		a.x > r.x && return ps
        return -ps
	elseif θ ≈ π
        a.x < r.x && return ps
        return -ps
	else
		m = tan(θ+π/2);
        th = m*(a.x-r.x)+r.y
		if θ <= π
            a.y < th && return -ps
            return ps
		else
            a.y > th && return -ps
            return ps
		end
	end
end

""" snapshot(sim, t)

Return a snapshot of the signal received by the array in simulation `sim` at time `t`" """
function snapshot(sim::Simulation, t)
    y = zeros(ComplexF64, length(sim.array))  # initialize received array
    for (source, manifold) in zip(sim.sources, sim.phaseshifts)
        v = randn(ComplexF64, length(sim.array))
        # received signal with unit-variance noise
        manifold = [f(source.θ) for f in sim.phaseshifts]
        y .+= source(t).*manifold + sqrt(sim.N₀/2)*randn(ComplexF64, length(manifold))
    end
    y
end

""" cbf(sim, fs, K)

Return the spectral power estimate of simulation `sim` using classical beamforming.
`fs` is the sampling frequency and `K` is the number of snapshots.  """
function cbf(sim::Simulation, fs, K ; ϕ = range(-π/2, π/2, length=100))
    P = length(sim.sources)
    #TODO: support more than one source
    P == 1 || error("Simulation must contain only one source")
    # snapshot instants
    T = range(0, length=K-1, step=1/fs)  # snapshot times
    max = 0.0
    doa = 0.0
    for θ in ϕ
        pbf = 0.0
        manifold = [f(θ) for f in sim.phaseshifts]'
        for t in T
            pbf += abs2(manifold*snapshot(sim, t))
        end
        if pbf > max
            doa = θ
            max = pbf
        end
    end
    doa
end

""" MUSIC(sim, fs, K)

Run the MUSIC algorithm on the simulation `sim`, sampling `K` times
at frequency `fs`.

Returns the MUSIC pseudospectrum of the simulation.  """
function MUSIC(sim::Simulation, fs, K)
    P = length(sim.sources)
    M = length(sim.array)
    # snapshot instants
    T = range(0, length=K-1, step=1/fs)  # snapshot times
    # generate received snapshots and calculate covariance matrix
    Ry = zeros(ComplexF64, M, M)  # initialize matrix
    for t in T
        y = snapshot(sim, t)
        Ry .+= y*y'
    end
    F = svd(Ry, full=true)  # in descending order
    Un = F.U[:,P+1:M]
    θ -> 1.0./sum(abs2.([f(θ) for f in sim.phaseshifts]'*Un))
end

""" findpeaks(mps, r)
Find peaks in function `mps`, defined on the range `r`.  """
function findpeaks(mps, r ; n = 1)
    ps = mps.(r)  # evaluate pseudospectrum
    peaks = [(height = 0.0, angle = 0.0)]  # largest peaks seen
    for i in 2:length(r)-1
        if (ps[i] > ps[i-1]) && (ps[i] > ps[i+1])
            for j in 1:n
                if ps[i] > peaks[j].height
                    x = (height=ps[i], angle=r[i])
                    if length(peaks) < n
                        push!(peaks,x)
                    else
                        peaks[j] = x
                    end
                    break
                end
            end
            sort!(peaks, by = x->x.height)
        end
    end
    return sort([x.angle for x in peaks])
end

""" beamwidth(S::Simulation ; c = 300e6)

Return the beamwidth (in radians) of the array in simulation `S`. As currently
implemented, the calculation is accurate only for linear arrays."""
beamwidth(S::Simulation; c = 300e6) = 0.89*(c/S.fc)/abs(maximum(S.array.d)-minimum(S.array.d))

beamwidth(a::AntennaArray, f ; c = 300e6) = 0.89*(c/f)/abs(maximum(a.d)-minimum(a.d))

end


text:
include(joinpath("..", "quadruped.jl"))

# test kinematics
q = rand(nq)

@assert norm(kinematics_1(model, q, body = :torso, mode = :ee) - [q[1] + model.l_torso * sin(q[3]); q[2] - model.l_torso * cos(q[3])]) ≈ 0.0
@assert norm(kinematics_1(model, q, body = :torso, mode = :com) - [q[1] + model.d_torso * sin(q[3]); q[2] - model.d_torso * cos(q[3])]) ≈ 0.0
@assert norm(kinematics_1(model, q, body = :thigh_1, mode = :ee) - [q[1] + model.l_thigh1 * sin(q[4]); q[2] - model.l_thigh1 * cos(q[4])]) ≈ 0.0
@assert norm(kinematics_1(model, q, body = :thigh_1, mode = :com) - [q[1] + model.d_thigh1 * sin(q[4]); q[2] - model.d_thigh1 * cos(q[4])]) ≈ 0.0
@assert norm(kinematics_1(model, q, body = :thigh_2, mode = :ee) - [q[1] + model.l_thigh2 * sin(q[6]); q[2] - model.l_thigh2 * cos(q[6])]) ≈ 0.0
@assert norm(kinematics_1(model, q, body = :thigh_2, mode = :com) - [q[1] + model.d_thigh2 * sin(q[6]); q[2] - model.d_thigh2 * cos(q[6])]) ≈ 0.0

@assert norm(kinematics_2(model, q, body = :calf_1, mode = :ee) - [q[1] + model.l_thigh1 * sin(q[4]) + model.l_calf1 * sin(q[5]); q[2] - model.l_thigh1 * cos(q[4]) - model.l_calf1 * cos(q[5])]) ≈ 0.0
@assert norm(kinematics_2(model, q, body = :calf_1, mode = :com) - [q[1] + model.l_thigh1 * sin(q[4]) + model.d_calf1 * sin(q[5]); q[2] - model.l_thigh1 * cos(q[4]) - model.d_calf1 * cos(q[5])]) ≈ 0.0
@assert norm(kinematics_2(model, q, body = :calf_2, mode = :ee) - [q[1] + model.l_thigh2 * sin(q[6]) + model.l_calf2 * sin(q[7]); q[2] - model.l_thigh2 * cos(q[6]) - model.l_calf2 * cos(q[7])]) ≈ 0.0
@assert norm(kinematics_2(model, q, body = :calf_2, mode = :com) - [q[1] + model.l_thigh2 * sin(q[6]) + model.d_calf2 * sin(q[7]); q[2] - model.l_thigh2 * cos(q[6]) - model.d_calf2 * cos(q[7])]) ≈ 0.0
@assert norm(kinematics_2(model, q, body = :thigh_3, mode = :ee) - [q[1] + model.l_torso * sin(q[3]) + model.l_thigh3 * sin(q[8]); q[2] - model.l_torso * cos(q[3]) - model.l_thigh3 * cos(q[8])]) ≈ 0.0
@assert norm(kinematics_2(model, q, body = :thigh_3, mode = :com) - [q[1] + model.l_torso * sin(q[3]) + model.d_thigh3 * sin(q[8]); q[2] - model.l_torso * cos(q[3]) - model.d_thigh3 * cos(q[8])]) ≈ 0.0
@assert norm(kinematics_2(model, q, body = :thigh_4, mode = :ee) - [q[1] + model.l_torso * sin(q[3]) + model.l_thigh4 * sin(q[10]); q[2] - model.l_torso * cos(q[3]) - model.l_thigh4 * cos(q[10])]) ≈ 0.0
@assert norm(kinematics_2(model, q, body = :thigh_4, mode = :com) - [q[1] + model.l_torso * sin(q[3]) + model.d_thigh4 * sin(q[10]); q[2] - model.l_torso * cos(q[3]) - model.d_thigh4 * cos(q[10])]) ≈ 0.0

@assert norm(kinematics_3(model, q, body = :calf_3, mode = :ee) - [q[1] + model.l_torso * sin(q[3]) + model.l_thigh3 * sin(q[8]) + model.l_calf3 * sin(q[9]); q[2] - model.l_torso * cos(q[3]) - model.l_thigh3 * cos(q[8]) - model.l_calf3 * cos(q[9])]) ≈ 0.0
@assert norm(kinematics_3(model, q, body = :calf_3, mode = :com) - [q[1] + model.l_torso * sin(q[3]) + model.l_thigh3 * sin(q[8]) + model.d_calf3 * sin(q[9]); q[2] - model.l_torso * cos(q[3]) - model.l_thigh3 * cos(q[8]) - model.d_calf3 * cos(q[9])]) ≈ 0.0
@assert norm(kinematics_3(model, q, body = :calf_4, mode = :ee) - [q[1] + model.l_torso * sin(q[3]) + model.l_thigh4 * sin(q[10]) + model.l_calf4 * sin(q[11]); q[2] - model.l_torso * cos(q[3]) - model.l_thigh4 * cos(q[10]) - model.l_calf4 * cos(q[11])]) ≈ 0.0
@assert norm(kinematics_3(model, q, body = :calf_4, mode = :com) - [q[1] + model.l_torso * sin(q[3]) + model.l_thigh4 * sin(q[10]) + model.d_calf4 * sin(q[11]); q[2] - model.l_torso * cos(q[3]) - model.l_thigh4 * cos(q[10]) - model.d_calf4 * cos(q[11])]) ≈ 0.0

k1(z) = kinematics_1(model, z, body = :torso, mode = :ee)
@assert norm(ForwardDiff.jacobian(k1,q) - jacobian_1(model, q, body = :torso, mode = :ee)) ≈ 0.0
k1(z) = kinematics_1(model, z, body = :torso, mode = :com)
@assert norm(ForwardDiff.jacobian(k1,q) - jacobian_1(model, q, body = :torso, mode = :com)) ≈ 0.0
k1(z) = kinematics_1(model, z, body = :thigh_1, mode = :ee)
@assert norm(ForwardDiff.jacobian(k1,q) - jacobian_1(model, q, body = :thigh_1, mode = :ee)) ≈ 0.0
k1(z) = kinematics_1(model, z, body = :thigh_1, mode = :com)
@assert norm(ForwardDiff.jacobian(k1,q) - jacobian_1(model, q, body = :thigh_1, mode = :com)) ≈ 0.0
k1(z) = kinematics_1(model, z, body = :thigh_2, mode = :ee)
@assert norm(ForwardDiff.jacobian(k1,q) - jacobian_1(model, q, body = :thigh_2, mode = :ee)) ≈ 0.0
k1(z) = kinematics_1(model, z, body = :thigh_2, mode = :com)
@assert norm(ForwardDiff.jacobian(k1,q) - jacobian_1(model, q, body = :thigh_2, mode = :com)) ≈ 0.0

k2(z) = kinematics_2(model, z, body = :calf_1, mode = :ee)
@assert norm(ForwardDiff.jacobian(k2,q) - jacobian_2(model, q, body = :calf_1, mode = :ee)) ≈ 0.0
k2(z) = kinematics_2(model, z, body = :calf_1, mode = :com)
@assert norm(ForwardDiff.jacobian(k2,q) - jacobian_2(model, q, body = :calf_1, mode = :com)) ≈ 0.0
k2(z) = kinematics_2(model, z, body = :calf_2, mode = :ee)
@assert norm(ForwardDiff.jacobian(k2,q) - jacobian_2(model, q, body = :calf_2, mode = :ee)) ≈ 0.0
k2(z) = kinematics_2(model, z, body = :calf_2, mode = :com)
@assert norm(ForwardDiff.jacobian(k2,q) - jacobian_2(model, q, body = :calf_2, mode = :com)) ≈ 0.0
k2(z) = kinematics_2(model, z, body = :thigh_3, mode = :ee)
@assert norm(ForwardDiff.jacobian(k2,q) - jacobian_2(model, q, body = :thigh_3, mode = :ee)) ≈ 0.0
k2(z) = kinematics_2(model, z, body = :thigh_3, mode = :com)
@assert norm(ForwardDiff.jacobian(k2,q) - jacobian_2(model, q, body = :thigh_3, mode = :com)) ≈ 0.0
k2(z) = kinematics_2(model, z, body = :thigh_4, mode = :ee)
@assert norm(ForwardDiff.jacobian(k2,q) - jacobian_2(model, q, body = :thigh_4, mode = :ee)) ≈ 0.0
k2(z) = kinematics_2(model, z, body = :thigh_4, mode = :com)
@assert norm(ForwardDiff.jacobian(k2,q) - jacobian_2(model, q, body = :thigh_4, mode = :com)) ≈ 0.0

k3(z) = kinematics_3(model, z, body = :calf_3, mode = :ee)
@assert norm(ForwardDiff.jacobian(k3, q) - jacobian_3(model, q, body = :calf_3, mode = :ee)) ≈ 0.0
k3(z) = kinematics_3(model, z, body = :calf_3, mode = :com)
@assert norm(ForwardDiff.jacobian(k3, q) - jacobian_3(model, q, body = :calf_3, mode = :com)) ≈ 0.0
k3(z) = kinematics_3(model, z, body = :calf_4, mode = :ee)
@assert norm(ForwardDiff.jacobian(k3, q) - jacobian_3(model, q, body = :calf_4, mode = :ee)) ≈ 0.0
k3(z) = kinematics_3(model, z, body = :calf_4, mode = :com)
@assert norm(ForwardDiff.jacobian(k3, q) - jacobian_3(model, q, body = :calf_4, mode = :com)) ≈ 0.0

q̇ = rand(nq)
eigen(M_func(model, q))

tmp_q(z) = _dLdq̇(model, z, q̇)
tmp_q̇(z) = _dLdq̇(model, q, z)
norm(ForwardDiff.jacobian(tmp_q̇,q̇) - M_func(model, q))


text:
"""
    TanhDerivative(x)

Derivative of the hyperbolic tangent activation function.

    Usage:
```julia-repl
TanhDerivative(x)
```
Arguments:
- `x` : Input value.
"""
function TanhDerivative(x::Float64)
    x₁ = 1 - tanh2(x)
    return x₁
end


text:
try
    using ProgressMeter
catch
    import Pkg
    Pkg.add("ProgressMeter")
    using ProgressMeter
end

struct RideShare
    passengers_from::Vector{NTuple{2,Int}}
    passengers_to::Vector{NTuple{2,Int}}
    time_windows::Vector{NTuple{2,Int}}
    drivers::Vector{NTuple{2,Int}}
    speed::Float64
    gas::Float64
    rate::Float64
end

function get_initial_solution(rs::RideShare, prob=0.5)::Vector{Vector{Int}}
    driver_idxs = eachindex(rs.drivers)
    solution = [Int[] for x in driver_idxs]
    for x in eachindex(rs.passengers_from)
        if rand() > prob
            push!(solution[rand(driver_idxs)], x)
        else
            pushfirst!(solution[rand(driver_idxs)], x)
        end
    end
    solution
end

function iterate_solution(solution::Vector{Vector{Int}})
    new_solution = deepcopy(solution)
    possible_drivers = findall(x -> x > 0,  length.(solution))
    s_i = new_solution[rand(possible_drivers)]
    if rand(Bool) # swap
        s_j = new_solution[rand(possible_drivers)]
        pos_i = rand(eachindex(s_i))
        pos_j = rand(eachindex(s_j))
        s_i[pos_i], s_j[pos_j] = s_j[pos_j], s_i[pos_i]
    else # insert-delete
        s_j = rand(new_solution)
        pos_i = rand(eachindex(s_i));
        t = s_i[pos_i]
        deleteat!(s_i, pos_i)
        pushfirst!(s_j, t)
    end
    new_solution
end

dist_L1(a, b) = abs(b[1] - a[1]) + abs(b[2] - a[2])

function rideshare_optimality(rs::RideShare, solution::Vector{Vector{Int}})
    revenue = 0.0
    trips_completed = 0
    for x in eachindex(rs.drivers)
        t = 0.0
        driver_position = rs.drivers[x]
        for y in solution[x]
            curr_pass_from = rs.passengers_from[y]
            curr_pass_to = rs.passengers_to[y]
            current_time_window = rs.time_windows[y]
            revenue -= dist_L1(driver_position, curr_pass_from) * rs.gas
            t += dist_L1(driver_position, curr_pass_from) / rs.speed
            if t > current_time_window[2]
                driver_position = curr_pass_from
                continue
            end
            t = max(t, current_time_window[1])
            revenue += dist_L1(curr_pass_from, curr_pass_to) * rs.rate
            t += dist_L1(curr_pass_from, curr_pass_to) / rs.speed
            driver_position = curr_pass_to
            trips_completed += 1
        end
    end
    if revenue > 0.0
        revenue * (trips_completed / length(rs.passengers_to))
    else
        revenue * (1 - trips_completed / length(rs.passengers_to))
    end
end

rs_samples=[
    RideShare(
        [(15, 22),(23,13),(25,23),(6,30),(13,28),(24,26),(4,2),
         (30,23),(28,24),(28,18),(20,2),(28,7),(14,22),(20,21),
         (27,10),(8,7),(24,19),(24,6),(7,7),(7,22)], # passengers_from
        [(17,17),(22,10),(24,26),(3,27),(14,23),(22,28),(5,7),
         (27,22), (31,28),(23,23),(23,-3),(23,9),(19,18),(19,21),
         (30,6),(7,2),(29,23),(20,6),(10,10),(6,18)], # passengers_to
        [(91,101),(44,54),(22,32),(68,78),(108,118),(182,192),(19,29),
         (36,46),(183,193), (48,58),(209,219),(145,155),(141,151),(59,69),
         (236,246),(69,79),(177,187),(15,25),(145,155),(110,120)], # time_windows
        [(18,22),(9,20),(6,29),(30,28),(26,8)], #drivers
        1.0, #speed
        0.5, #gas
        2.5 # rate
        ),
    RideShare(
        [(4,30),(4,25),(21,1),(19,4),(3,14),(19,24),(29,17),
         (11,23),(17,29),(24,22)], # passengers_from
        [(11,36),(11,20),(22,6),(29,5),(-7,7),(20,31),(26,22),
         (10,23),(20,35),(24,28)], # passengers_to
        [(152,177),(254,279),(50,75),(255,280),(296,321),(76,101),
         (228,253),(134,159),(82,107),(201,226)], # time_windows
        [(27,1),(7,14),(20,11)], # drivers
        1.0, #speed
        0.5, #gas
        2.5 # rate
        )
    ]


text:
using Test

include("solver.jl")

@testset "(220, 284) are amicable numbers" begin

  @test divisors(220) == [1, 2, 4, 5, 10, 11, 20, 22, 44, 55, 110]
  @test divisors(284) == [1, 2, 4, 71, 142]

  @test is_amicable(220, 284)
end

@testset "amicable exception" begin
  @test_throws ArgumentError is_amicable(3.5, -10)
end

@testset "Σ(amicable) ≤ 10_000" begin
  s, _v = find_amicable(10_000)

  @test sum(s) == 31626
end

max_stars_repo_name: pascal-p/Euler

max_issues_repo_name: pascal-p/Euler

max_forks_repo_name: pascal-p/Euler


text:
using ParaDwell
using CSSUtil, Interact, Blink, CSV, DataFrames, PlotlyJS, StatsBase, JSON
using Dates

# Constructor script
env = ParaDwell.loadPaths()

try
    close(ui.mainWindow)
catch
end

ParaDwell.OS_TileSummries(env)
project = ParaDwell.Project()
project.dat = Dict()
project.tileRegister = CSV.read(joinpath(env.paths[:projects],".OS_TileRegister.csv"), DataFrame)

include("UI/MainUI_Launch.jl")
include("Util/DataLinks.jl")
include("Util/JsonHandling.jl")
include("Util/ExportData.jl")
include("Util/GeomOperations.jl")
include("DataMethods/DataLinks.jl")
include("DataMethods/LinkEpcData.jl")
include("DataMethods/LinkOsBhaData.jl")
include("DataMethods/GeoRef.jl")
include("ParametricModel/ProcessGeom.jl")
include("ParametricModel/LinkMultiDwellingPlans.jl")
include("ParametricModel/ParametricModel.jl")
include("ParametricModel/ExtrapolateEpc.jl")

ui = LaunchMainUI(env);
include("UI/ClickTile.jl")
ClickTile(env,project,ui)
project = SelectTile(env,project,ui,"HY40NE")
project = SelectTile(env,project,ui,"HY40NW")
project = SelectTile(env,project,ui,"HY41SW")
project = SelectTile(env,project,ui,"HY41SE")
# project = SelectTile(env,project,ui,"NZ16SE")

@time project = LoadStockData(env,project)

@time project = buildarchetypes(project)

@time fillblockconfig.((project.dat["master"],),2:4)

@time abovebelowadj(project.dat["master"])

@time gen_archcode(project.dat["master"],1.)

# exportstockdata(project)


text:
export  PlaneStress

@doc raw"""
    PlaneStress

Creates a plane stress element

- `H`: Linear elasticity matrix, $3\times3$
- `E`: Young's modulus
- `ν`: Poisson's ratio
- `ρ`: density
- `σ0`: stress at the **last** time step
- `σ0_`: (for internal use), stress to be updated in `commitHistory`
- `ε0`: strain at the **last** time step
- `ε0_`: (for internal use), strain to be updated in `commitHistory`

# Example
```julia
prop = Dict("name"=> "PlaneStress", "rho"=> 0.0876584, "E"=>0.07180760098, "nu"=>0.4)
mat = PlaneStress(prop)
```
"""
mutable struct PlaneStress
    H::Array{Float64}
    E::Float64
    ν::Float64
    ρ::Float64
    σ0::Array{Float64} # stress at last time step
    σ0_::Array{Float64} # σ0 to be updated in `commitHistory`
    ε0::Array{Float64}
    ε0_::Array{Float64}
end

"""
    PlaneStress(prop::Dict{String, Any})

`prop` should contain at least the following three fields: `E`, `nu`, `rho`
"""
function PlaneStress(prop::Dict{String, Any})
    E = prop["E"]; ν = prop["nu"]; ρ = prop["rho"]
    H = zeros(3,3)

    H[1,1] = E/(1. -ν*ν)
    H[1,2] = H[1,1]*ν
    H[2,1] = H[1,2]
    H[2,2] = H[1,1]
    H[3,3] = E/(2.0*(1.0+ν))
    σ0 = zeros(3); σ0_ = zeros(3); ε0 = zeros(3); ε0_ = zeros(3)
    PlaneStress(H, E, ν, ρ, σ0, σ0_,ε0,ε0_)
end

function getStress(self::PlaneStress, strain::Array{Float64}, Dstrain::Array{Float64}, Δt::Float64 = 0.0)
    sigma = self.H * strain
    self.σ0_ = copy(sigma)
    self.ε0_ = copy(strain)
    return sigma, self.H
end

function getTangent(self::PlaneStress)
    self.H
end

function commitHistory(self::PlaneStress)
    self.σ0 = self.σ0_
    self.ε0 = self.ε0_
end


text:
module StructLH

using ArgCheck, DocStringExtensions

export @common_fields
export retrieve_property, has_property, retrieve_child_property
export merge_object_arrays!, reduce_object_vector
export ApplyFctResult, apply_fct_to_object, obj_name, obj_type, fct_value, children
export describe, describe_object, show_description
export NodeInfo, struct2dict, dict2struct!
# Reductions
export reduce_scalar_vector, reduce_array_vector

include("helpers.jl");
include("retrieve.jl")
include("reduce.jl");
include("struct2dict.jl")
include("traverse_object.jl")

"""
    merge_object_arrays!

Merge all arrays (and vectors) from object `oSource` into the corresponding arrays
in another object `oTg` (at given index values `idxV`).
If target does not have corresponding field: behavior is governed by `skipMissingFields`.

# Arguments
- idxV
    `Vector{Integer}` or other iterable with integer results. Indexes the first dimension of each object to be copied.
    The assignment is: `oSource.x[i1,:] => oTg.x[idxV[i1],:]`.
"""
function merge_object_arrays!(oSource, oTg, idxV,
    skipMissingFields :: Bool; dbg :: Bool = false)

    nameV = propertynames(oSource);
    for name in nameV
        xSrc = getproperty(oSource, name);
        if isa(xSrc,  Array)
            # Does target have this field?
            if hasproperty(oTg, name)
                xTg = getproperty(oTg, name);
                if dbg
                    @assert size(xSrc, 1) == length(idxV)
                    @assert size(xSrc)[2:end] == size(xTg)[2:end] "Size mismatch: $(size(xSrc)) vs $(size(xTg))"
                end
                # The n-dim array code also works for Vectors, but is less efficient.
                if isa(xSrc, Vector)
                    xTg[idxV] .= xSrc;
                elseif isa(xSrc, Matrix)
                    xTg[idxV, :] .= xSrc;
                else
                    # For multidimensional arrays (we don't know the dimensions!)
                    # we need to loop over "rows". This is expensive.
                    for (i1, idx) in enumerate(idxV)
                        # This selects target "row" `idx`
                        tgView = selectdim(xTg, 1, idx);
                        # Copy source "row" `i1` into target row (in place, hence [:])
                        tgView[:] = selectdim(xSrc, 1, i1);
                    end
                end
            elseif !skipMissingFields
                error("Missing field $name in target object")
            end
        end
    end
    return nothing
end

end # module