text:
using Revise, Test, ForwardDiff, Parameters, Setfield, Plots, LinearAlgebra
using BifurcationKit, Test
const BK = BifurcationKit

norminf(x) = norm(x, Inf)
####################################################################################################
function COm(u, p)
	@unpack q1,q2,q3,q4,q5,q6,k = p
	x, y, s = u
	z = 1-x-y-s
	[
		2q1 * z^2 - 2q5 * x^2 - q3 * x * y,
		q2 * z - q6 * y - q3 * x * y,
		q4 * (z - k * s)
	]
end
jet = BK.getJet(COm; matrixfree=false)

par_com = (q1 = 2.5, q2 = 2.0, q3 = 10., q4 = 0.0675, q5 = 1., q6 = 0.1, k = 0.4)

z0 = [0.001137, 0.891483, 0.062345]

opts_br = ContinuationPar(pMin = 0.5, pMax = 2.0, ds = 0.002, dsmax = 0.01, nInversion = 6, detectBifurcation = 3, maxBisectionSteps = 25, nev = 3, maxSteps = 20000)
	@set! opts_br.newtonOptions.verbose = true
	br, = @time continuation(jet[1], jet[2], z0, par_com, (@lens _.q2), opts_br;
		recordFromSolution = (x, p) -> (x = x[1], y = x[2], s = x[3]),
		plot = false, verbosity = 3, normC = norminf,
	bothside = true)
	show(br)

plot(br, plotfold=false, markersize=4, legend=:topright, ylims=(0,0.16))
####################################################################################################
@set! opts_br.newtonOptions.verbose = true
@set! opts_br.newtonOptions.maxIter = 10
opts_br = @set opts_br.newtonOptions.tol = 1e-12

sn, = newton(jet[1:2]..., br, 2; options = opts_br.newtonOptions, bdlinsolver = MatrixBLS())

hp, = newton(jet[1:2]..., br, 1; options = NewtonPar( opts_br.newtonOptions; maxIter = 10),startWithEigen=true, d2F = jet[3])

BK.hopfNormalForm(jet..., hp, opts_br.newtonOptions.linsolver)
hpnf = computeNormalForm(jet..., br, 1)

sn_codim2, = continuationFold(jet[1:2]..., br, 2, (@lens _.k), ContinuationPar(opts_br, pMax = 3.2, pMin = 0., detectBifurcation = 0, dsmin=1e-5, ds = -0.001, dsmax = 0.05, nInversion = 6, detectEvent = 2, detectFold = false) ; plot = true,
	verbosity = 3,
	normC = norminf,
	updateMinAugEveryStep = 1,
	startWithEigen = true,
	recordFromSolution = (u,p; kw...) -> (x = u.u[1] ),
	bothside=true,
	bdlinsolver = MatrixBLS()
	)

using Test
@test sn_codim2.specialpoint[1].printsol.k 	≈ 0.971397 rtol = 1e-4
@test sn_codim2.specialpoint[1].printsol.q2 ≈ 1.417628 rtol = 1e-4
@test sn_codim2.specialpoint[3].printsol.k 	≈ 0.722339 rtol = 1e-4
@test sn_codim2.specialpoint[3].printsol.q2 ≈ 1.161199 rtol = 1e-4

plot(sn_codim2)#, real.(sn_codim2.BT), ylims = (-1,1), xlims=(0,2))

plot(sn_codim2, vars=(:q2, :x), branchlabel = "Fold", plotstability = false);plot!(br,xlims=(0.8,1.8))

hp_codim2, = continuation(jet[1:2]..., br, 1, (@lens _.k), ContinuationPar(opts_br, pMin = 0., pMax = 2.8, detectBifurcation = 0, ds = -0.0001, dsmax = 0.08, dsmin = 1e-4, nInversion = 6, detectEvent = 2, detectLoop = true, maxSteps = 50, detectFold=false) ; plot = true,
	verbosity = 3,
	normC = norminf,
	tangentAlgo = BorderedPred(),
	detectCodim2Bifurcation = 2,
	updateMinAugEveryStep = 1,
	startWithEigen = true,
	recordFromSolution = (u,p; kw...) -> (x = u.u[1] ),
	d2F = jet[3], d3F = jet[4],
	bothside = true,
	bdlinsolver = MatrixBLS())

@test hp_codim2.branch[5].l1 |> real 		≈ 33.15920 rtol = 1e-1
@test hp_codim2.specialpoint[1].printsol.k 	≈ 0.305879 rtol = 1e-3
@test hp_codim2.specialpoint[1].printsol.q2 ≈ 0.924255 rtol = 1e-3
@test hp_codim2.specialpoint[2].printsol.k 	≈ 0.235550 rtol = 1e-4
@test hp_codim2.specialpoint[2].printsol.q2 ≈ 0.896099 rtol = 1e-4

plot(sn_codim2, vars=(:q2, :x), branchlabel = "Fold", plotcirclesbif = true)
	plot!(hp_codim2, vars=(:q2, :x), branchlabel = "Hopf",plotcirclesbif = true)
	plot!(br,xlims=(0.6,1.5))

plot(sn_codim2, vars=(:k, :q2), branchlabel = "Fold")
	plot!(hp_codim2, vars=(:k, :q2), branchlabel = "Hopf",)

plot(hp_codim2, vars=(:q2, :x), branchlabel = "Hopf")
####################################################################################################


text:
# julia --project=. plotscript.jl
#
# This will take a long time to start so won't
# recommend doing this.
using Plots
gr

# x = 1:10
# y = rand(10)
# p = plot(x,y)
# gui()
# display(p)
# println("Hello World!")
# readline()

n = 50 # points
x_vec = fill(0.0, n)
y_vec = fill(0.0, n)
xmin = -10.0
xmax = 10.0
x_range = LinRange(xmin, xmax, n)
# Sigmoid function
# Generalized logistic function
function fx(x, w, b, alpha)
    #y = 1/(1 + exp(-w*x + b))
    y = (1+ exp(-w*x)+b)^(-alpha)
    #y = exp(-w*x)
    return y
end

"""
    glf(x, A, K, B, mu, Q, C, M)

Generalized logistic function
"""
function glf(x, A, K, B, mu, Q, C, M)
    Y = A + (K-A)/(C + Q*exp(-B*(x-M)))^(1/mu)
end
# w = 1
# b = 20# 0 - 10
# alp = 0.3
A = 0.01
K = 1.0
Q = 1.0 #0.5
C = 1.0

# B = 1.8# 2 # 5 # 0 to 10 # top angle
# mu = 1.8 # Close rate  # bottom angle
# M = 5.8 # shift left and right
B = 1.8 # 2 # 5 # 0 to 10 # top angle
mu = 1.7 # Close rate  # bottom angle
M = 5.8 # shift left and right
for (ind,x) in enumerate(x_range)
    y_vec[ind] = glf(x, A, K, B, mu, Q, C, M)
end

# scale to 14.7 to 9200 psi
pmin = 14.7
pmax = 9200
pdist = pmax - pmin
dist = xmax - xmin
scale = pdist/dist
for (ind,x) in enumerate(x_range)
    x_vec[ind] = (x+10)*scale + pmin
end

pt = plot(x_vec,y_vec, legend=(0.1,0.9))

# B = 2 # 2 # 5 # 0 to 10 # top angle
# mu = 2 # Close rate  # bottom angle
# M = 5 # shift left and right
# for (ind,x) in enumerate(x_range)
#     y_vec[ind] = glf(x, A, K, B, mu, Q, C, M)
# end

# scale to 14.7 to 9200 psi
# pmin = 14.7
# pmax = 9200
# pdist = pmax - pmin
# dist = xmax - xmin
# scale = pdist/dist
# for (ind,x) in enumerate(x_range)
#     x_vec[ind] = (x+10)*scale + pmin
# end
#
# plot!(pt, x_vec, y_vec)

ps = [  14.7,   3000,    5000,    5316,    5800,   6284,   6500,       6700, 7167,  7601, 7900, 8200, 9204]
ks = [0.01, 0.01,    0.01,    0.015,    0.04,      0.08,     0.16,   0.30,   0.65,  0.85, 0.95, 1, 1]

plot!(pt, ps, ks)

ext: jl


text:
# License for this file: MIT (expat)
# Copyright 2017-2020, DLR Institute of System Dynamics and Control

"""
    module Modia3D.PathPlanning

Module for path planning and designing reference paths for Modia3D

# Main developers

Andrea Neumayr and Martin Otter, [DLR - Institute of System Dynamics and Control](https://www.dlr.de/sr/en)
"""
module PathPlanning

export PTP_path, pathEndTime, getPosition!, getPosition, getIndex, plotPath
export ReferencePath, scheduleReferenceMotion, calculateRobotMovement, getRefPathPosition, getRefPathInitPosition
export ptpJointSpace
export ptpJointSpace!

export getPathComponent

import Base
import Modia3D
import Modia3D.Basics
import Modia3D.Composition
import Modia3D.Shapes

using ModiaLang

using Unitful
using StaticArrays
using OrderedCollections

include("pathPlanning.jl")
include("referencePathInternal.jl")
include("referencePathUser.jl")

end


text:
using BinaryProvider # requires BinaryProvider 0.3.0 or later
using Libdl: dlext

# Parse some basic command-line arguments
const verbose = "--verbose" in ARGS
const prefix = Prefix(get([a for a in ARGS if a != "--verbose"], 1, joinpath(@__DIR__, "usr")))
const LIBDIR = Sys.iswindows() ? Sys.BINDIR : joinpath(Sys.BINDIR, Base.LIBDIR)

const products = [
    LibraryProduct(LIBDIR, "librelic_gmp_pbc_bls381", :librelic_gmp_pbc_bls381),
    LibraryProduct(LIBDIR, "librelic_gmp_pbc_bn254", :librelic_gmp_pbc_bn254),
]

# Download binaries from hosted location
const version = "0.4.0-613-g613795f7-3"
const bin_prefix = "https://github.com/blenessy/RelicToolkitBuilder/releases/download/$version"

# Listing of files generated by BinaryBuilder:
const download_info = Dict(
    Linux(:i686, :glibc) => ("$bin_prefix/RelicToolkit.v$version.i686-linux-gnu.tar.gz", "e3019742022ba470fc8666fed87e510c55de2277932cb4c64dafab1915ff03e6"),
    Windows(:i686) => ("$bin_prefix/RelicToolkit.v$version.i686-w64-mingw32.tar.gz", "f60e0993059febdf7cf6b23121345e13a71aa5ffba0a391373c467aca94aac56"),
    MacOS(:x86_64) => ("$bin_prefix/RelicToolkit.v$version.x86_64-apple-darwin14.tar.gz", "60e2b5b7aacacad61d1fd18f476733fecb1749b23afb9ddafb5c8bd949e7681e"),
    Linux(:x86_64, :glibc) => ("$bin_prefix/RelicToolkit.v$version.x86_64-linux-gnu.tar.gz", "cd097d550e9a13ab4968f6f7b2d61a9bef291687d66833f248b1d9f3532e646d"),
    Windows(:x86_64) => ("$bin_prefix/RelicToolkit.v$version.x86_64-w64-mingw32.tar.gz", "83f17928a0d3fc7370d489a6787c0a77ec0f5d84fee73aa058d7e1e5cfa6c91b"),
)

# First, check to see if we're all satisfied
if any(!satisfied(p; verbose=verbose) for p in products)
    try
        # Download and install binaries
        url, tarball_hash = choose_download(download_info)
        install(url, tarball_hash; prefix=prefix, force=true, verbose=verbose)
        # The libs link towards libgmp-10.dll - this is an easy workaround
        if Sys.iswindows()
            cp(joinpath(LIBDIR, "libgmp.dll"), joinpath(Sys.BINDIR, "libgmp-10.dll"), force=true)
        end
        # Workaround: Copy libs to Sys.BINDIR to make the Windows libs work
        # tl;dr: there seems to bug when native libs depend on other native libs (libgmp-10.dll in our case)
        for p in products
            for n in p.libnames
                cp(joinpath(libdir(prefix), "$n.$dlext"), joinpath(LIBDIR, "$n.$dlext"), force=true)
            end
        end
    catch e
        if typeof(e) <: ArgumentError
            error("Your platform $(Sys.MACHINE) is not supported by this package!")
        else
            rethrow(e)
        end
    end
    # Finally, write out a deps.jl file
    write_deps_file(joinpath(@__DIR__, "deps.jl"), products, verbose=verbose)
end


text:
using IntUtils

function main()
    n = 1
    while true
        if num_factors(sum(1:n)) > 500
            break
        end
        n += 1
    end
    @printf("The first triangle number with over 500 factors is %d\n", sum(1:n))
end

main()


text:
function get_method_instance(f, typs) ::MethodInstance
  world = ccall(:jl_get_world_counter, UInt, ())
  tt = typs isa Type ? Tuple{typeof(f), typs.parameters...} : Tuple{typeof(f), typs...}
  results = Base._methods_by_ftype(tt, -1, world)
  @assert length(results) == 1 "get_method_instance should return one method, instead returned $(length(results)) methods: $results"
  (_, _, meth) = results[1]
  # TODO not totally sure what jl_match_method is needed for - I think it's just extracting type parameters like `where {T}`
  (ti, env) = ccall(:jl_match_method, Any, (Any, Any), tt, meth.sig)::SimpleVector
  meth = Base.func_for_method_checked(meth, tt)
  linfo = ccall(:jl_specializations_get_linfo, Ref{MethodInstance}, (Any, Any, Any, UInt), meth, tt, env, world)
end

function get_code_info(method_instance::MethodInstance; optimize=true) ::Tuple{CodeInfo, Type}
  world = ccall(:jl_get_world_counter, UInt, ())
  # TODO inlining=false would make analysis easier to follow, but it seems to break specialization on function types
  params = Core.Inference.InferenceParams(world)
  optimize = true
  cache = false # not sure if cached copies use the same params
  (_, code_info, return_typ) = Core.Inference.typeinf_code(method_instance, optimize, cache, params)
  (code_info, return_typ)
end

"Does this look like error reporting code ie not worth looking inside?"
function is_error_path(expr)
  expr == :throw ||
  expr == :throw_boundserror ||
  expr == :error ||
  expr == :assert ||
  (expr isa QuoteNode && is_error_path(expr.value)) ||
  (expr isa Expr && expr.head == :(.) && is_error_path(expr.args[2])) ||
  (expr isa GlobalRef && is_error_path(expr.name)) ||
  (expr isa MethodInstance && is_error_path(expr.def.name))
end

"Is it pointless to look inside this expression?"
function should_ignore(expr::Expr)
  is_error_path(expr.head) ||
  (expr.head == :call && is_error_path(expr.args[1])) ||
  (expr.head == :invoke && is_error_path(expr.args[1]))
end

"Return all function calls in the method whose argument types can be determined statically"
function get_child_calls(method_instance::MethodInstance)
  code_info, return_typ = get_code_info(method_instance, optimize=true)
  calls = Set{MethodInstance}()

  function walk_expr(expr)
    if isa(expr, MethodInstance)
      push!(calls, expr)
    elseif isa(expr, Expr)
      if !should_ignore(expr)
        foreach(walk_expr, expr.args)
      end
    end
  end
  foreach(walk_expr, code_info.code)

  calls
end

"A node in the call graph"
struct CallNode
  call::MethodInstance
  parent_calls::Set{MethodInstance}
  child_calls::Set{MethodInstance}
end

"Return as much of the call graph of `method_instance` as can be determined statically"
function get_call_graph(method_instance::MethodInstance, max_calls=1000::Int64) ::Vector{CallNode}
  all = Dict{MethodInstance, CallNode}()
  ordered = Vector{MethodInstance}()
  unexplored = Set{Tuple{Union{Void, MethodInstance}, MethodInstance}}(((nothing, method_instance),))
  for _ in 1:max_calls
    if isempty(unexplored)
      return [all[call] for call in ordered]
    end
    (parent, call) = pop!(unexplored)
    child_calls= get_child_calls(call)
    parent_calls = parent == nothing ? Set() : Set((parent,))
    all[call] = CallNode(call, parent_calls, child_calls)
    push!(ordered, call)
    for child_call in child_calls
      if !haskey(all, child_call)
        push!(unexplored, (call, child_call))
      else
        push!(all[child_call].parent_calls, call)
      end
    end
  end
  error("get_call_graph reached $max_calls calls and gave up")
end

@generated function show_structure(x)
  quote
    @show x
    $([:(isdefined(x, $(Expr(:quote, fieldname))) ? (@show x.$fieldname) : nothing) for fieldname in fieldnames(x)]...)
    x
  end
end

function trace_static(filter::Function, warn::Function, f::Function, typs)
  for call_node in get_call_graph(get_method_instance(f, typs))
    if filter(call_node.call)
      analyse((a...) -> warn(Warning(a...)), StaticCall(call_node.call))
    end
  end
end

warntrace_static(filter::Function, f::Function, typs) = trace_static(filter, warning_printer(), f, typs)
warntrace_static(f::Function, typs) = warntrace_static((_) -> true, f, typs)

@eval begin
  macro trace_static(ex0)
    Base.gen_call_with_extracted_types($(Expr(:quote, :warntrace_static)), ex0)
  end

  macro trace_static(filter, ex0)
    expr = Base.gen_call_with_extracted_types($(Expr(:quote, :warntrace_static)), ex0)
    insert!(expr.args, 2, esc(filter))
    expr
  end
end

function warnings_static(f)
  warnings = Warning[]
  trace_static((_) -> true, w -> push!(warnings, w), f, ())
  return warnings
end


text:
module BSONify

export as_type
export TypeValueMismatch
export FieldsMismatch
export TypeIsNotAMapping

include("tobson.jl")
include("frombson.jl")

end # module

max_stars_repo_name: OliverColeman/BSONify

max_forks_repo_name: OliverColeman/BSONify


text:
"""
    za = agg(z,nas,nat)

Aggregate field `z` to an array `za` of size `(nas,nas,nat)`
"""
function agg(zi,nas,nat)
        nss=size(zi);
        if(length(nss)>=3)
           ns=nss[1]; nt=nss[3];
        else
           ns=nss[1]; nt=1;
        end
        sdim=div(ns,nas);
        tdim=div(nt,nat);
        ir=1:sdim; jr=ir; kr=1:tdim;
        rs=0:sdim:(ns-sdim);
        rt=0:tdim:(nt-tdim);
	z=zeros(Float64,nas,nas,nat);
        n=zeros(Int64,nas,nas,nat);
        for i=1:sdim
           for j=1:sdim
              for k=1:tdim
		 n=n+.~isnan.(zi[i.+rs,j.+rs,k.+rt]);
		 z=z+zi[i.+rs,j.+rs,k.+rt].*.~isnan.(zi[i.+rs,j.+rs,k.+rt]) ;
              end
           end
        end
        z=z./n;
end


text:
for (op,dotop) in ((:(+),:(.+)), (:(-),:(.-)), (:(*),:(.*)), (:(/),:(./)))
  @eval begin
      function ($op){P}(x::ArbFloat{P}, y::Array{ArbFloat{P},1})
          return ($dotop)(x, y)
      end
      function ($op){P}(x::Array{ArbFloat{P},1}, y::ArbFloat{P})
          return ($dotop)(x, y)
      end
      function ($op){P}(x::Array{ArbFloat{P},1}, y::Array{ArbFloat{P},1})
          return ($dotop)(x, y)
      end
  end
end


text:
using Dates, DelimitedFiles, Plots, Random, ScoreDrivenModels

# Define dates and load historical Affluent Natural Energy data
dates = collect(Date(1961):Month(1):Date(2000, 12))
y = vec(readdlm("../test/data/ane_northeastern.csv"))
y_train = y[1:400]
y_test = y[401:460]

# Set RNG seed to guarantee consistent results
Random.seed!(123)

# Specify GAS model: a lognormal model with time-varying μ, constant σ, and lags 1 and 12
gas = Model([1, 2, 11, 12], [1, 2, 11, 12], LogNormal, 0.0; time_varying_params=[1])

# Obtain initial parameters to start the GAS recursion
initial_params = dynamic_initial_params(y_train, gas)

# Fit specified model to historical data using initial parameters
f = ScoreDrivenModels.fit!(gas, y_train; initial_params=initial_params)

# Print estimation statistics
fit_stats(f)

# Simulate 1000 future scenarios and obtain the 5% and 95% quantiles in each time period
forec = forecast(y_train, gas, 60; S=1000, initial_params=initial_params)

# Plot results
plot(dates[401:460], forec.observation_scenarios, color="grey", width=0.05, label="", ylims=(0, 70))
plot!(dates[360:460], y[360:460], label="ANE", color="black", xlabel="Months", ylabel="GWmed", legend=:topright)
plot!(dates[401:460], forec.observation_quantiles, label=["Quantiles" "" ""], color="red", line=:dash)


text:
using Devectorize
function cost(expanded_target_batch, output_layer_state, batchsize)
@devec CE = -sum(expanded_target_batch .* log(output_layer_state + exp(-30))) ./ batchsize
return CE
end

max_stars_repo_name: kirnap/Machine-Learning-in-Julia

max_issues_repo_name: kirnap/Machine-Learning-in-Julia

max_forks_repo_name: kirnap/Machine-Learning-in-Julia


text:
#=
Given a list of integers, write a function that returns the largest sum of non-adjacent numbers. Numbers can be 0 or negative.

For example, [2, 4, 6, 2, 5] should return 13, since we pick 2, 6, and 5. [5, 1, 1, 5] should return 10, since we pick 5 and 5.

Follow-up: Can you do this in O(N) time and constant space?
=#

using Test
include("Solutions/problem9_sum_notadjacentnumbers.jl")
@test sum_notadjacent([2, 4, 6, 2, 5]) == 13
@test sum_notadjacent([5, 1, 1, 5]) == 10
@test sum_notadjacent([-3, -4, -7, 1, 0, 1]) == 2

max_stars_repo_name: DominiqueCaron/daily-coding-problem

max_issues_repo_name: DominiqueCaron/daily-coding-problem

max_forks_repo_name: DominiqueCaron/daily-coding-problem


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "sed"
version = v"4.8.1"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://ftp.gnu.org/gnu/sed/sed-4.8.tar.xz", "f79b0cfea71b37a8eeec8490db6c5f7ae7719c35587f21edb0617f370eeff633")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/sed-*
if [[ "${target}" == *-mingw* ]]; then
    # Fix error
    #    sed/sed-compile.o: In function `sprintf':
    #    /opt/x86_64-w64-mingw32/x86_64-w64-mingw32/sys-root/include/stdio.h:366: undefined reference to `__chk_fail'
    # See https://github.com/msys2/MINGW-packages/issues/5868#issuecomment-544107564
    export LIBS="-lssp"
fi

./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target}
make -j${nproc} SUBDIRS="po ."
make install SUBDIRS="po ."
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms(; experimental=true)

# The products that we will ensure are always built
products = [
    ExecutableProduct("sed", :sed)
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
    Dependency("Libiconv_jll"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
module DisjunctiveProgramming

using JuMP, IntervalArithmetic, Symbolics

export add_disjunction
export @disjunction

include("reformulate.jl")
include("disjunction.jl")
include("macro.jl")

end # module


text:
"""
    applyrule(data::AbstractSimData, rule::Rule{R,W}, state, index::Tuple{Int,Int}) -> cell value(s)

Apply a rule to the cell state and return values to write to the grid(s).

This is called in `maprule!` methods during the simulation,
not by the user. Custom `Rule` implementations must define this method.

## Arguments

- `data` : [`AbstractSimData`](@ref)
- `rule` : [`Rule`](@ref)
- `state`: the value(s) of the current cell
- `index`: a (row, column) tuple of Int for the current cell coordinates

Returns the value(s) to be written to the current cell(s) of
the grids specified by the `W` type parameter.
"""
function applyrule end

"""
    applyrule!(data::AbstractSimData, rule::{R,W}, state, index::Tuple{Int,Int}) -> Nothing

Apply a rule to the cell state and manually write to the grid data array.
Used in all rules inheriting from [`SetCellRule`](@ref).

This is called in internal `maprule!` methods during the simulation, not by
the user. Custom [`SetCellRule`](@ref) implementations must define this method.

Only grids specified with the `W` type parameter will be writable from `data`.

## Arguments

- `data` : [`AbstractSimData`](@ref)
- `rule` : [`Rule`](@ref)
- `state`: the value(s) of the current cell
- `index`: a (row, column) tuple of Int for the current cell coordinates - `t`: the current time step
"""
function applyrule! end

"""
    modifyrule(rule::Rule, data::AbstractSimData) -> Rule

Precalculates rule fields at each timestep. Define this method if a [`Rule`](@ref)
has fields that need to be updated over time.

`Rule`s are immutable (it's faster and works on GPU), so `modifyrule` is
expected to return a new rule object with changes applied to it. Setfield.jl or
Acessors.jl may help with updating the immutable struct.

The default behaviour is to return the existing rule without change. Updated rules
are discarded after use, and the `rule` argument is always the original object passed in.

# Example

We define a rule with a parameter that is the total sum of the grids current,
and update it for each time-step using `modifyrule`.

This could be used to simulate top-down control e.g. a market mechanism in a
geographic model that includes agricultural economics.

```jldoctest
using DynamicGrids, Setfield
struct MySummedRule{R,W,T} <: CellRule{R,W}
    gridsum::T
end
function modifyrule(rule::MySummedRule{R,W}, data::AbstractSimData) where {R,W}
    Setfield.@set rule.gridsum = sum(data[R])
end

# output
modifyrule (generic function with 1 method)
```
"""
function modifyrule end

"""
    neighbors(x::Union{Neighborhood,NeighborhoodRule}}) -> iterable

Returns an indexable iterator for all cells in the neighborhood,
either a `Tuple` of values or a range.

Custom `Neighborhood`s must define this method.
"""
function neighbors end

"""
    neighborhood(x::Union{NeighborhoodRule,SetNeighborhoodRule}}) -> Neighborhood

Returns a rules neighborhood object
"""
function neighborhood end

"""
    kernel(hood::AbstractKernelNeighborhood) => iterable

Returns the kernel object, an array or iterable matching the length
of the neighborhood.
"""
function kernel end

"""
    kernelproduct(rule::NeighborhoodRule})
    kernelproduct(hood::AbstractKernelNeighborhood)
    kernelproduct(hood::Neighborhood, kernel)

Returns the vector dot product of the neighborhood and the kernel,
although differing from `dot` in that the dot product is not take for
vector members of the neighborhood - they are treated as scalars.
"""
function kernelproduct end

"""
    offsets(x::Union{Neighborhood,NeighborhoodRule}}) -> iterable

Returns an indexable iterable over all cells, containing `Tuple`s of
the index offset from the central cell.

Custom `Neighborhood`s must define this method.
"""
function offsets end

"""
    positions(x::Union{Neighborhood,NeighborhoodRule}}, cellindex::Tuple) -> iterable

Returns an indexable iterable, over all cells as `Tuple`s of each
index in the main array. Useful in [`SetNeighborhoodRule`](@ref) for
setting neighborhood values, or for getting values in an Aux array.
"""
function positions end

"""
    add!(data::WritableGridData, x, I...)

Add the value `x` to a grid cell.

## Example useage

```jldoctest
using DynamicGrids
rule = SetCell{:a}() do data, a, cellindex
    dest, is_inbounds = inbounds(data, (jump .+ cellindex)...)

    # Update spotted cell if it's on the grid
    is_inbounds && add!(data[:a], state, dest...)
end

# output
SetCell{:a,:a} :
    f = var"#1#2"
```
"""
function add! end

"""
    sub!(data::WritableGridData, x, I...)

Subtract the value `x` from a grid cell. See `add!` for example usage.
"""
function sub! end

"""
    min!(data::WritableGridData, x, I...)

Set a gride cell to the minimum of `x` and the current value. See `add!` for example usage.
"""
function min! end

"""
    max!(data::WritableGridData, x, I...)

Set a gride cell to the maximum of `x` and the current value. See `add!` for example usage.
"""
function max! end

"""
    and!(data::WritableGridData, x, I...)
    and!(A::AbstractArray, x, I...)

Set the grid cell `c` to `c & x`. See `add!` for example usage.
"""
function and! end

"""
    or!(data::WritableGridData, x, I...)
    or!(A::AbstractArray, x, I...)

Set the grid cell `c` to `c | x`. See `add!` for example usage.
"""
function or! end

"""
    xor!(data::WritableGridData, x, I...)
    xor!(A::AbstractArray, x, I...)

Set the grid cell `c` to `xor(c, x)`. See `add!` for example usage.
"""
function xor! end

"""
    inbounds(data::AbstractSimData, I::Tuple) -> Tuple{NTuple{2,Int}, Bool}
    inbounds(data::AbstractSimData, I...) -> Tuple{NTuple{2,Int}, Bool}

Check grid boundaries for a coordinate before writing in [`SetCellRule`](@ref).

Returns a `Tuple` containing a coordinates `Tuple` and a `Bool` - `true`
if the cell is inside the grid bounds, `false` if not.

[`BoundaryCondition`](@ref) of type [`Remove`](@ref) returns the coordinate and `false`
to skip coordinates that boundary outside of the grid.

[`Wrap`](@ref) returns a tuple with the current position or it's
wrapped equivalent, and `true` as it is allways in-bounds.
"""
function inbounds end

"""
    isinbounds(data, I::Tuple) -> Bool
    isinbounds(data, I...) -> Bool

Check that a coordinate is within the grid, usually in [`SetCellRule`](@ref).

Unlike [`inbounds`](@ref), [`BoundaryCondition`](@ref) status is ignored.
"""
function isinbounds end

"""
    radius(rule, [key]) -> Int

Return the radius of a rule or ruleset if it has one, otherwise zero.
"""
function radius end

"""
    init(obj) -> Union{AbstractArray,NamedTUple}

Retrieve the mask from an [`Output`](@ref), [`Extent`](@ref) or [`AbstractSimData`](@ref) object.
"""
function init end

"""
    mask(obj) -> AbstractArray

Retrieve the mask from an [`Output`](@ref), [`Extent`](@ref) or [`AbstractSimData`](@ref) object.
"""
function mask end

"""
    aux(obj, [key])

Retrieve auxilary data `NamedTuple` from an [`Output`](@ref),
[`Extent`](@ref) or [`AbstractSimData`](@ref) object.

Given `key` specific data will be returned. `key` should be a
`Val{:symbol}` for type stability and zero-cost access inside rules.
`Symbol` will also work, but may be slow.
"""
function aux end

"""
    tspan(obj) -> AbstractRange

Retrieve the time-span `AbstractRange` from an [`Output`](@ref),
[`Extent`](@ref) or [`AbstractSimData`](@ref) object.
"""
function tspan end

"""
    timestep(obj)

Retrieve the timestep size from an [`Output`](@ref),
[`Extent`](@ref), [`Ruleset`](@ref) or [`AbstractSimData`](@ref) object.

This will be in whatever type/units you specify in `tspan`.
"""
function timestep end

"""
    currentframe(simdata::AbstractSimData) -> Int

Retrieve the current simulation frame a [`AbstractSimData`](@ref) object.
"""
function currentframe end

"""
    currenttime(simdata::AbstractSimData)

Retrieve the current simulation time from a [`AbstractSimData`](@ref) object.

This will be in whatever type/units you specify in `tspan`.
"""
function currenttime end

"""
    currenttimestep(simdata::AbstractSimData)

Retrieve the current timestep from a [`AbstractSimData`](@ref) object.

This may be different from the `timestep`. If the timestep is `Month`,
`currenttimestep` will return `Seconds` for the length of the specific month.
"""
function currenttimestep end


text:
import LinearAlgebra
using LinearAlgebra
using SparseArrays
using DelimitedFiles

using jInv
### matrix free S (slicing) matrix
export normalize

function normalize(xyz::Array{Float64},eps::Float64 = 1e-5)
	n = div(length(xyz),3);
	x = xyz[1:n];
	y = xyz[(n+1):2*n];
	z = xyz[(2*n+1):3*n];
	norms = sqrt.(x.^2 .+ y.^2 .+ z.^2 .+ eps);
	inv_n = 1.0./norms;
	inv_np3 = inv_n.^3;
	nx = x./norms;
	ny = y./norms;
	nz = z./norms;
	xyz_norm = [nx;ny;nz];
	#d(x/sqrt(x^2+y^2+z^2))_dx = 1/sqrt(x^2+y^2+z^2) - x^2/(x^2+y^2+z^2)^{3/2}
	#d(x/sqrt(x^2+y^2+z^2))_dy =  - x*y/(x^2+y^2+z^2)^{3/2}
	dxdx = sparse(Diagonal(inv_n .- (x.^2).*inv_np3));
	dxdy = sparse(Diagonal(      .- (x.*y).*inv_np3));
	dxdz = sparse(Diagonal(      .- (x.*z).*inv_np3));
	dydx = dxdy;
	dydy = sparse(Diagonal(inv_n .- (y.^2).*inv_np3));
	dydz = sparse(Diagonal( 	 .- (y.*z).*inv_np3));
	dzdx = dxdz;
	dzdy = dydz;
	dzdz = sparse(Diagonal(inv_n .- (z.^2).*inv_np3));
	J = [dxdx dxdy dxdz ; dydx dydy dzdy ; dzdx dzdy dzdz];
	return xyz_norm, J
end

export getData
function getData(m::Array,pFor::PointCloudParam,doClear::Bool=false)
d = 0;
Minv = pFor.Mesh;
n = Minv.n;
npc = pFor.npcAll;
samplingBinning = pFor.samplingBinning;
ind_array = pFor.P;
n_points = 0;
for k=1:length(ind_array)
	n_points+=length(ind_array[k]);
end
d = zeros(Float32,3*n_points);
Jacobian = convert(SparseMatrixCSC{Float64,Int32},spzeros(length(d),length(m)));

if pFor.method == MATFree
	Af   = getFaceAverageMatrix(Minv)
	A1 = Af[:,1:Int(size(Af,2)/3)];
	A2 = Af[:,Int(size(Af,2)/3)+1:2*Int(size(Af,2)/3)];
	A3 = Af[:,2*Int(size(Af,2)/3)+1:Int(size(Af,2))];
	Af = blockdiag(A1,A2,A3);

	D1 = kron(sparse(1.0I, Minv.n[3], Minv.n[3]),kron(sparse(1.0I, Minv.n[2], Minv.n[2]),ddx(Minv.n[1])))
	D2 = kron(sparse(1.0I, Minv.n[3], Minv.n[3]),kron(ddx(Minv.n[2]),sparse(1.0I, Minv.n[1], Minv.n[1])))
	D3 = kron(ddx(Minv.n[3]),kron(sparse(1.0I, Minv.n[2], Minv.n[2]),sparse(1.0I, Minv.n[1], Minv.n[1])))
	D = [D1 D2 D3];
	global count = 1;
	for i=1:npc
		ind = ind_array[i];
		nx = A1*(D1'*m);
		ny = A2*(D2'*m);
		nz = A3*(D3'*m);
		nx = nx[ind];
		ny = ny[ind];
		nz = nz[ind];
		#normals = Af*D*m;
		#normals,J = normalize([nx;ny;nz]);
		normals = [nx;ny;nz];
		d[count:(count + 3*length(ind) - 1)] = normals;
		I2 = collect(1:length(ind));
		J2 = ind;
		V2 = ones(size(ind));
		P = sparse(I2,J2,V2,length(ind),length(m));
		P = blockdiag(P,P,P);
		#J = J*P*Af*D';
		J = P*Af*D';
		Jacobian[count:(count + 3*length(ind) - 1),:] = J;
		count = count + 3*length(ind);
	end
	pFor.Jacobian = Jacobian;
	return d,pFor;
end

if pFor.method == RBFBasedSimple || pFor.method == RBF10BasedSimple
	if pFor.method == RBF10BasedSimple
		numParamOfRBF = 10;
	else
		numParamOfRBF = 5;
	end
        println("numparam:",numParamOfRBF);
	println("m size:",size(m));
	println("m:",m)
	theta_phi = pFor.theta_phi_rad;
	b = pFor.b;
	println("theta phi:",theta_phi); println("b:",b)
	mrot,Jrot = rotateAndMoveRBFsimple(m,pFor.Mesh,theta_phi,b;computeJacobian = 1,numParamOfRBF=numParamOfRBF);
	println("mrot size:",size(mrot));
	println("mrot:",mrot)
	d,JacT = getPCDataRBF(pFor, mrot,theta_phi,b,numParamOfRBF);
	# multiply Jacobian with Jrot, and then take the transpose.
	pFor.Jacobian = Jrot'*JacT;
end

if pFor.method == RBFBased || pFor.method == RBF10Based || pFor.method == RBF5Based
	if pFor.method == RBF10Based
		numParamOfRBF = 10;
	else
		numParamOfRBF = 5;
	end
	#nRBF 				= div(length(m)- 5*pFor.npcAll,numParamOfRBF) ;
	#(m1,theta_phi,b) 	= splitRBFparamAndRotationsTranslations(m,nRBF,pFor.npcAll,numParamOfRBF);
	#theta_phi 			= theta_phi[pFor.workerSubIdxs,:];
	#b 					= b[pFor.workerSubIdxs,:];
	#mrot,Jrot 			= rotateAndMoveRBF(m1,Minv,theta_phi,b;computeJacobian = 1,numParamOfRBF=numParamOfRBF);
	theta_phi = 0; b =0;
	nRBF = div(length(m),numParamOfRBF) ;
	println("size of m:",size(m));
	d,JacT 				= getPCDataRBF(pFor, m,theta_phi,b,numParamOfRBF);
	# multiply Jacobian with Jrot, and then take the transpose.

	JacobianT = convert(SparseMatrixCSC{Float64,Int32},spzeros(length(m),prod(size(d))));
	IpIdxs = getIpIdxs(pFor.workerSubIdxs,nRBF,pFor.npcAll,numParamOfRBF);
	#JacobianT[IpIdxs,:] = Jrot'*JacT;
	#JacobianT[IpIdxs,:] = JacT;
	pFor.Jacobian = JacT';
end
return d,pFor;
end

function getPCDataRBF(pFor, m, theta_phi,b,numParamOfRBF = 5)
Mesh = pFor.Mesh;
n = pFor.Mesh.n;
h = pFor.Mesh.h;
ndips = pFor.npcAll;
volCell = prod(Mesh.h);
npc = ndips;
eps = 1e-4/Mesh.h[1];
Af   = getFaceAverageMatrix(Mesh)
A1 = Af[:,1:Int(size(Af,2)/3)];
A2 = Af[:,Int(size(Af,2)/3)+1:2*Int(size(Af,2)/3)];
A3 = Af[:,2*Int(size(Af,2)/3)+1:Int(size(Af,2))];
Af = blockdiag(A1,A2,A3);

D1 = kron(sparse(1.0I, Mesh.n[3], Mesh.n[3]),kron(sparse(1.0I, Mesh.n[2], Mesh.n[2]),ddx(Mesh.n[1])))
D2 = kron(sparse(1.0I, Mesh.n[3], Mesh.n[3]),kron(ddx(Mesh.n[2]),sparse(1.0I, Mesh.n[1], Mesh.n[1])))
D3 = kron(ddx(Mesh.n[3]),kron(sparse(1.0I, Mesh.n[2], Mesh.n[2]),sparse(1.0I, Mesh.n[1], Mesh.n[1])))
D = [D1 D2 D3];

#P is a sparse matrix of size k \times n, where k = #points on point cloud , n = volume size of the mesh
Jacobians = Array{SparseMatrixCSC{Float64,Int32}}(undef, npc);
count = 1;
ind_array = pFor.P;

n_points = 0;
for k=1:length(ind_array)
	n_points+=length(ind_array[k]);
end
d = zeros(Float32,3*n_points);

for i=1:npc
	sigmaH = getDefaultHeavySide();
	u,I1,J1,V1 = ParamLevelSetModelFunc(Mesh,m[:,i];computeJacobian = 1,sigma = sigmaH,bf = 1,numParamOfRBF = numParamOfRBF);
	J1 = sparse(I1,J1,V1,prod(n),length(m[:,i]));
	#u,I1,J1,V1 = ParamLevelSetModelFunc(Mesh,m;computeJacobian = 1,sigma = sigmaH,bf = 1,numParamOfRBF = numParamOfRBF);
	#J1 = sparse(I1,J1,V1,prod(n),length(m));
	ind = ind_array[i];
	nx = A1*(D1'*u);
	ny = A2*(D2'*u);
	nz = A3*(D3'*u);
	nx = nx[ind];
	ny = ny[ind];
	nz = nz[ind];
	#normals = Af*D*m;
	#normals,J = normalize([nx;ny;nz]);
	normals = [nx;ny;nz];
	d[count:(count + 3*length(ind) - 1)] = normals;
	I2 = collect(1:length(ind));
	J2 = ind;
	V2 = ones(size(ind));
	P = sparse(I2,J2,V2,length(ind),length(u));
	P = blockdiag(P,P,P);
	#J = J*P*Af*D';
	J = P*Af*D';
	#Jacobian[count:(count + 3*length(ind) - 1),:] = J1'*J';
	count = count + 3*length(ind);
	Jacobians[i] = J1'*J';

end
JacT = blockdiag(Jacobians...);
return d,JacT
end

max_stars_repo_name: BGUCompSci/ShapeReconstructionPaLS

max_issues_repo_name: BGUCompSci/ShapeReconstructionPaLS

max_forks_repo_name: BGUCompSci/ShapeReconstructionPaLS


text:
# julia script to test PUMI interface

using PumiInterface
include("funcs2.jl")
apf.declareNames();  # declare global variable names
# apf.initilize mesh
dmg_name = "cube.dmg"
smb_name = "tet-mesh-1.smb"
#dmg_name = "apf.reorder_a.dmg"
#smb_name = "apf.reorder_a.smb"
#dmg_name = ".null"
#smb_name = ".smb"
downward_counts_tmp, num_entities_tmp, m2_ptr, mshape_ptr = apf.init(dmg_name, smb_name)
m_ptr = apf.getMeshPtr()
mshape_ptr = apf.getMeshShapePtr()

vertN_ptr = getVertNumbering()
edgeN_ptr = getEdgeNumbering()
faceN_ptr = getFaceNumbering()
elN_ptr = getElNumbering()

global const downward_counts = downward_counts_tmp
global const num_entities = num_entities_tmp

println("downward_counts = ", downward_counts)
println("num_entities = ", num_entities)

checkVars();

resetVertIt();
resetEdgeIt();
resetFaceIt();
resetElIt();

checkNums()

resetVertIt();
resetEdgeIt();
resetFaceIt();
resetElIt();

coords = Array{Float64}(3, 1)   # pass an array 3 by n (3 coordinates each for n points)
(m,n) = size(coords)
apf.getVertCoords(coords, m, n)

coords = zeros(3,2)
(m,n) = size(coords)
apf.getEdgeCoords(coords, m, n)

# get face coordinates
numP = downward_counts[3, 1]
coords = zeros(3,numP)
(m,n) = size(coords)
apf.getFaceCoords(coords, m, n)
println("from julia, face coords = ", coords)
print("\n")

# get element coordinates
numP = downward_counts[4,1]
coords = zeros(3, numP)
(m,n) = size(coords)
apf.getElCoords(coords, m,n)
println("from julia, element coords = ", coords)

resetElIt()

print("\n")

# get all vertex coordinates
vertCoords = getAllVertCoords()
println("verts coords = ", vertCoords)

resetVertIt()
resetEdgeIt()
resetFaceIt()
resetElIt()

i = getVertNumber()
println(" first vertex number = " , i)
incrementVertItn(3);
i = getVertNumber();
println(" 4th vertex number = ", i )

i = getEdgeNumber()
println(" first edge number = " , i)
incrementEdgeItn(3);
i = getEdgeNumber();
println(" 4th Edge number = ", i )

i = getFaceNumber()
println(" first face number = " , i)
incrementFaceItn(3);
i = getFaceNumber();
println(" 4th Face number = ", i )

i = getElNumber()
println(" first element number = " , i)
incrementElItn(3);
i = getElNumber();
println(" 4th element number = ", i )

entity = getVert()
println("typeof(entity) in main = ", typeof(entity) )
i = getVertNumber2(entity)
println("vertex number = ", i)

i = apf.getMeshDimension(m_ptr);
println("mesh dimension from julia = ", i )

i = apf.countNodesOn(mshape_ptr, 2);
println("number of nodes on type 2 = ", i)

numbering_ptr = apf.createNumberingJ(m_ptr, "numberingJ1", mshape_ptr, 1)
println("created numbering")

apf.numberJ(numbering_ptr, entity, 0, 0, 5)
i = apf.getNumberJ(numbering_ptr, entity, 0, 0)
println("number of entity = ", i)

i = apf.getNumberJ(vertN_ptr, entity, 0, 0)
println("global number of vertex entity = ", i)

entity = getEdge()
i = apf.getNumberJ(edgeN_ptr, entity, 0, 0)
println("global number of edge entity = ", i)

entity = getFace()
i = apf.getNumberJ(faceN_ptr, entity, 0, 0)
println("global number of face entity = ", i)

entity = getEl()
i = apf.getNumberJ(elN_ptr, entity, 0, 0)
println("global number of element entity = ", i)

i = apf.getMeshDimension(m2_ptr)
println("mesh2 dimension = ", i)

downwards, npoints = apf.getDownward(m_ptr, entity, 0)
println("npoints = ", npoints)
j = 1
for entity2 in downwards
  num = apf.getNumberJ(vertN_ptr, entity2, 0, 0);
  println("entity $j on current element has number", num)
  j += 1
end

entity_type = apf.getType(m_ptr, entity)
println("type of current entity = ", entity_type)
i = apf.countNodesOn(mshape_ptr, entity_type)
println(" number of nodes on current entity = ", i)

apf.printNumberingName(faceN_ptr)

resetVertIt()
incrementIt(3)
entity = getEl()
i = apf.getNumberJ(elN_ptr, entity, 0, 0)
println("global number of vertex entity = ", i)

tag_ptr = apf.createDoubleTag( m_ptr, "tag1", 2)
apf.setDoubleTag( m_ptr, entity, tag_ptr, [3.5, 2.5])
data_ret = zeros(1,2)
data_ret[1] = 2.0
apf.getDoubleTag( m_ptr, entity, tag_ptr, data_ret)
println("data_ret = ", data_ret)

entity = getVert()
n = apf.countAdjacent(m_ptr, entity, 3)
#println("vertex has ", i, " 3d regions ")

adj = apf.getAdjacent(n)

(adj, n) = apf.getAdjacentFull(m_ptr, entity, 3)
println("vertex has ", n, " 3d regions ")
println("size(adj) = ", size(adj))
for k in adj
  n = apf.getNumberJ(elN_ptr, k, 0, 0)
  m = apf.countDownwards(m_ptr, k)
  println("element number = ", n)
  println("downward adjacency counts = ", m)
end

apf.hasNodesIn(mshape_ptr, 0)
apf.hasNodesIn(mshape_ptr, 1)

i = apf.countNodesOn(mshape_ptr, apf.getType(m_ptr, entity))
println("number of nodes on vertex = ", i)

for i = 0:4
  j = apf.countNodesOn(mshape_ptr, i)
  println("entity of type ", i, "has ", j, " nodes")
end

eshape_ptr = apf.getEntityShape(mshape_ptr, 3)
i = apf.countNodes(eshape_ptr)
println("number of nodes on entity of type 3 = ", i)

i = apf.countAllNodes(mshape_ptr, 3)
println("number of nodes on entity of type 3 = ", i)

apf.printEdgeVertNumbers(edgeN_ptr, vertN_ptr)

apf.writeVtkFiles("vtkTest", m_ptr)

entity = getFace()
etype = apf.getType(m_ptr, entity)
eshape_ptr = apf.getEntityShape(mshape_ptr, etype)

mel_ptr = apf.createMeshElement(m_ptr, entity)
i = apf.countIntPoints(mel_ptr, 5)
println("face requires ", i, " points for 5th order accurate integration")
numN = apf.countNodes(eshape_ptr)

for k=1:i
  coords = apf.getIntPoint(mel_ptr, 5, k)
  weight = apf.getIntWeight(mel_ptr, 5,k)
  sh_vals = apf.getValues2(eshape_ptr, coords)
  shd_vals = apf.getLocalGradients2(eshape_ptr, coords)
  jac = apf.getJacobian(mel_ptr, coords)
  println("point ", k, " has coordinates ", coords, " and weight ", weight)
  println("shape function values = ", sh_vals)
  println("shape function derivative values = ", shd_vals)
  println("jacobian = \n", jac)
  print("\n")
end

resetVertIt()
entity = getVert()
coords = zeros(3,1)
apf.getVertCoords(entity, coords, 3, 1)
println("coords = ", coords)

resetEdgeIt()
entity = getEdge()
coords = zeros(3,2)
apf.getEdgeCoords(entity, coords, 3, 2)
println("coords = ", coords)

resetFaceIt()
entity = getFace()
coords = zeros(3,3)
apf.getFaceCoords(entity, coords, 3, 4)


text:
module _PermEn
export PermEn
using Combinatorics: permutations
using Statistics: std, var, mean
    """
        Perm, Pnorm, cPE = PermEn(Sig)

    Returns the permuation entropy estimates `Perm`, the normalised
    permutation entropy `Pnorm` and the conditional permutation entropy `cPE`
    for `m` = [1,2] estimated from the data sequence `Sig`  using the default
    parameters: embedding dimension = 2, time delay = 1, logarithm = base 2,
    normalisation = w.r.t #symbols (`m`-1)
    Note: using the standard PermEn estimation, `Perm` = 0 when `m` = 1.
    Note: It is recommeneded that signal length > 5m! (see [8] and Amigo et al., Europhys. Lett. 83:60005, 2008)

        Perm, Pnorm, cPE = PermEn(Sig, m)

    Returns the permutation entropy estimates `Perm` estimated from the data
    sequence `Sig` using the specified embedding dimensions = [1,...,`m`]
    with other default parameters as listed above.

        Perm, Pnorm, cPE = PermEn(Sig::AbstractArray{T,1} where T<:Real; m::Int=2, tau::Int=1, Typex::String="none", tpx::Union{Real,Nothing}=nothing, Logx::Real=2, Norm::Bool=false)

    Returns the permutation entropy estimates `Perm` for dimensions = [1,...,`m`]
    estimated from the data sequence `Sig` using the specified 'keyword' arguments:

    # Arguments:
    `m`     - Embedding Dimension, an integer > 1 \n
    `tau`   - Time Delay, a positive integer\n
    `Logx`  - Logarithm base, a positive scalar (enter 0 for natural log) \n
    `Norm`  - Normalisation of PermEn value:\n
              false -  normalises w.r.t log(# of permutation symbols [m-1]) - default
              true  -  normalises w.r.t log(# of all possible permutations [m!])
              * Note: Normalised permutation entropy is undefined for m = 1.
              ** Note: When Typex = 'uniquant' and Norm = true, normalisation
              is calculated w.r.t. log(tpx^m)\n
    `Typex`  - Permutation entropy variation, one of the following:
              {`"none", "uniquant", "finegrain", "modified", "ampaware",
              "weighted", "edge"}
              See the EntropyHub guide for more info on PermEn variations.    \n
    `tpx`   - Tuning parameter for associated permutation entropy variation.\n
              [uniquant]  'tpx' is the L parameter, an integer > 1 (default = 4).
              [finegrain] 'tpx' is the alpha parameter, a positive scalar (default = 1)
              [ampaware]  'tpx' is the A parameter, a value in range [0 1] (default = 0.5)
              [edge]      'tpx' is the r sensitivity parameter, a scalar > 0 (default = 1)
              See the EntropyHub guide for more info on PermEn variations.    \n

    # See also `XPermEn`, `MSEn`, `XMSEn`, `SampEn`, `ApEn`, `CondEn`

    # References:
        [1] Christoph Bandt and Bernd Pompe,
                "Permutation entropy: A natural complexity measure for time series."
                Physical Review Letters,
                88.17 (2002): 174102.

        [2] Xiao-Feng Liu, and Wang Yue,
                "Fine-grained permutation entropy as a measure of natural
                complexity for time series."
                Chinese Physics B
                18.7 (2009): 2690.

        [3] Chunhua Bian, et al.,
                "Modified permutation-entropy analysis of heartbeat dynamics."
                Physical Review E
                85.2 (2012) : 021906

        [4] Bilal Fadlallah, et al.,
                "Weighted-permutation entropy: A complexity measure for time
                series incorporating amplitude information."
                Physical Review E
                87.2 (2013): 022911.

        [5] Hamed Azami and Javier Escudero,
                "Amplitude-aware permutation entropy: Illustration in spike
                detection and signal segmentation."
                Computer methods and programs in biomedicine,
                128 (2016): 40-51.

        [6] Zhiqiang Huo, et al.,
                "Edge Permutation Entropy: An Improved Entropy Measure for
                Time-Series Analysis,"
                45th Annual Conference of the IEEE Industrial Electronics Soc,
                (2019), 5998-6003

        [8] Zhe Chen, et al.
                "Improved permutation entropy for measuring complexity of time
                series under noisy condition."
                Complexity
                1403829 (2019).

        [9] Maik Riedl, Andreas Müller, and Niels Wessel,
                "Practical considerations of permutation entropy."
                The European Physical Journal Special Topics
                222.2 (2013): 249-262.

    """
    function PermEn(Sig::AbstractArray{T,1} where T<:Real; m::Int=2, tau::Int=1,
        Typex::String="none", tpx::Union{Real,Nothing}=nothing, Logx::Real=2, Norm::Bool=false)

    Logx == 0  ? Logx = exp(1) : nothing

    N = size(Sig)[1]
    (N>10) ? nothing :   error("Sig:   must be a numeric vector")
    (m > 1) ? nothing :  error("m:     must be an integer > 1")
    (tau > 0) ? nothing :  error("tau:   must be an integer > 0")
    (Logx>0) ? nothing : error("Logx:     must be a positive number > 0")
    (lowercase(Typex) in ["none","uniquant","finegrain","modified","ampaware","weighted","edge"]) ?
     nothing : error("Typex:    must be one of the following strings -
        'uniquant','finegrain','modified','ampaware','weighted','edge'")
    (isnothing(tpx) || tpx>0) ? nothing : error("tpx:   the value of tpx relates to 'Type'.
        See the EntropyHub guide for further info on the 'tpx' value.")

    Sx = zeros(N,m)
    Perm = zeros(m)
    Pnorm = zeros(m)
    for k = 1:m
        Nx = N-(k-1)*tau
        Sx[1:Nx,k] = Sig[1+(k-1)*tau:N]
        Temp = sortind(Sx[1:Nx,1:k])
        Px = collect(permutations(collect(1:k)))
        Counter = zeros(length(Px))

        if lowercase(Typex) == "uniquant"
            Temp = sort(Sx[1:Nx,1:k], dims=2)
            S = zeros(size(Temp))
            if isnothing(tpx)
                tpx = 4;
            elseif tpx <= 1 || typeof(tpx)!=Int
                error("When Typex = 'Uniquant', L parameter (tpx) must be an integer > 1")
            end
            delta = (maximum(Sig)-minimum(Sig))/tpx
            S[:,1] = map(x -> searchsortedfirst(
                minimum(Sig):delta:maximum(Sig),x), Temp[:,1]) .- 1
            #S[findall(S[:,1].==0),1] .+= 1
            S[S[:,1].==0,1] .+= 1
            if k > 1
                S[:,2:k] = S[:,1] .+ floor.((Temp[:,2:k] .- Temp[:,1])/delta)
            end
            Px = unique(S,dims=1)
            Counter = zeros(size(Px,1))
            for n = 1:size(Px,1)
                Counter[n] = sum(all(S .- transpose(Px[n,:]).==0,dims=2))
            end
            Counter = Counter[Counter.!=0]
            Ppi = Counter/sum(Counter)
            Norm = true  # Might need to change back to 3

        elseif lowercase(Typex) == "finegrain"
            if k > 1
                if isnothing(tpx)
                    tpx = 1;
                elseif tpx <= 0
                    error("When Typex = 'finegrain', Alpha parameter (tpx) must be greater than 0")
                end
                q =  floor.(maximum(abs.(diff(Sx[1:Nx,1:k],dims=2)),dims=2)./
                                (tpx*std(abs.(diff(Sig)),corrected=false)))
                Temp = hcat(Temp, q)
                Px = unique(Temp,dims=1)
                Counter = zeros(size(Px,1))
                for n = 1:size(Px,1)
                    Counter[n] = sum(all(Temp .- transpose(Px[n,:]).== 0,dims=2))
                end
                Counter = Counter[Counter.!=0]
                Ppi = Counter./sum(Counter)
                #clear q n qt
            else
                Ppi = 1
            end

        elseif lowercase(Typex) == "modified"
            Tx = (diff(sort(Sx[1:Nx,1:k],dims=2),dims=2).==0)
            for km = 1:k-1
                Temp[Tx[:,km],km+1] = Temp[Tx[:,km],km];
            end
            Px = unique(Temp,dims=1)
            Counter = zeros(size(Px,1))
            for n = 1:size(Px,1)
                Counter[n] = sum(all(Temp .- transpose(Px[n,:]) .==0,dims=2))
            end
            Counter = Counter[Counter.!=0]
            Ppi = Counter/sum(Counter)
            #clear Tx km

        elseif lowercase(Typex) == "weighted"
            if k > 1
                Wj = var(Sx[1:Nx,1:k],corrected=false,dims=2)
                for n = 1:size(Px,1)
                    Counter[n] = sum(Wj[all(Temp .- transpose(Px[n]) .==0,dims=2)])
                end
                Counter = Counter[Counter.!=0]
                Ppi = Counter/sum(Wj)
                #clear Wj n
            else
                Ppi = 1;
            end

        elseif lowercase(Typex) == "ampaware"
            if k > 1
                if isnothing(tpx)
                    tpx = 0.5;
                elseif tpx<0 || tpx>1
                    error("When Typex = 'ampaware', the A parameter must be in the range [0 1]")
                end
                AA = sum(abs.(Sx[1:Nx,1:k]),dims=2)
                AB = sum(abs.(diff(Sx[1:Nx,1:k],dims=2)),dims=2)
                Ax = (tpx*AA/k) + ((1-tpx)*AB/(k-1));
                for n = 1:size(Px,1)
                    Counter[n] = sum(Ax[all(Temp.-transpose(Px[n]).==0,dims=2)])
                end
                Counter = Counter[Counter.!=0]
                Ppi = Counter/sum(Ax);
            else
                Ppi = 1;
            end
            #clear AA AB Ax

        elseif lowercase(Typex) == "edge"
            if isnothing(tpx)
                tpx = 1;
            elseif tpx <=0
                error("When Typex = 'Edge', the r sensitivity parameter (tpx) must be > 0")
            end
            if k > 1
                for n = 1:size(Px,1)
                    Sy = Sx[1:Nx,1:k]
                    Tx = diff(Sy[all(Temp .- transpose(Px[n]) .==0,dims=2)[:],:],dims=2)
                    Counter[n] = sum(mean(hypot.(Tx,1),dims=2).^tpx)
                end
                Counter = Counter[Counter.!=0]
                Ppi = Counter/sum(Counter)
            else Ppi = 1;
            end

        else
            for n = 1:size(Px,1)
                Counter[n] = sum(all(Temp .- transpose(Px[n]).==0,dims=2));
                #sum(all(Temp .- transpose(Px[n]).==0,dims=2));
            end
            Counter = Counter[Counter.!=0]
            Ppi = Counter/sum(Counter)
        end

        if round(sum(Ppi),digits=3) != 1
            @warn ("Potential error with probability calculation")
        end

        Perm[k] = -sum(Ppi.*(log.(Logx, Ppi)));
        if Norm
            if Norm && lowercase(Typex)=="uniquant"
                Pnorm[k] = Perm[k]/(log(Logx, tpx^k));
            else
                Pnorm[k] = Perm[k]/(log(Logx, factorial(k)));
            end
        else
            Pnorm[k] = Perm[k]/(k-1);
        end
    end
    cPE = diff(Perm);

    return Perm, Pnorm, cPE
    end

    function sortind(X)
        Y = zeros(Int, size(X))
        for k = 1:length(X[:,1])
            Y[k,:] = sortperm(X[k,:])
        end
        return Y
    end

end

"""
Copyright 2021 Matthew W. Flood, EntropyHub

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

    http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.

For Terms of Use see https://github.com/MattWillFlood/EntropyHub
"""


text:
module MarkovChainUtil

using ..Resampling
using ..Reproduce
import StatsBase

export
    env_settings!,
    get_experience,
    get_experience_interpolation,
    policy_settings,
    max_is_dict,
    get_policies

const policy_settings = Dict(
    "easy"=>([0.5, 0.5], [0.1, 0.9]),
    "hard"=>([0.9, 0.1], [0.1, 0.9]),
    "hardest"=>([0.99, 0.01], [0.01, 0.99]))

const max_is_dict = Dict(
    "easy"=>1.8,
    "hard"=>9.0,
    "hardest"=>99.0)

function env_settings!(s::ArgParseSettings)
    @add_arg_table s begin
        "--chainsize"
        arg_type=Int64
        default=10
        "--policy"
        help = "The policy settings"
        default="easy"
    end
end

function get_policies(policy)
    return policy_settings[policy][1], policy_settings[policy][2]
end

function get_policies(parsed::Dict)
    return get_policies(parsed["policy"])
end

function get_experience(buffer_size, μ, π; rng=Random.GLOBAL_RNG, chain_size=10, state_mod=identity)

    env = MarkovChain(chain_size)
    _, prev_state = start!(env; rng=rng)
    ϕ = state_mod(prev_state)

    er_names = [:s_t, :s_tp1, :a_t, :r, :γ_tp1, :terminal, :mu_t, :pi_t, :ρ]
    er_types = [typeof(ϕ), typeof(ϕ), Int64, Float64, Float64, Bool, Float64, Float64, Float64]
    ER = ExperienceReplay(buffer_size, er_types; column_names=er_names)
    WER = WeightedExperienceReplay(buffer_size, er_types; column_names=er_names)

    for i = 1:buffer_size

        action = StatsBase.sample(rng, get_actions(env), StatsBase.Weights(μ))
        _, state, reward, terminal = step!(env, action)

        experience = (state_mod(prev_state), state_mod(state), action, reward, 0.9, terminal,
                      μ[action], π[action], π[action]/μ[action])
        add!(ER, experience)
        add!(WER, experience, π[action]/μ[action])

        prev_state = copy(state)
        if terminal
            _, prev_state = start!(env; rng=rng)
        end
    end

    return ER, WER
end

function get_experience_interpolation(buffer_size, μ, π, sample_policy;
                                      rng=Random.GLOBAL_RNG, chain_size=10,
                                      state_mod=identity)

    env = MarkovChain(chain_size)
    _, prev_state = start!(env; rng=rng)
    ϕ = state_mod(prev_state)

    er_names = [:s_t, :s_tp1, :a_t, :r, :γ_tp1, :terminal, :mu_t, :pi_t, :ρ]
    er_types = [typeof(ϕ), typeof(ϕ), Int64, Float64, Float64, Bool, Float64, Float64, Float64]
    WER = WeightedExperienceReplay(buffer_size, er_types; column_names=er_names)

    for i = 1:buffer_size

        action = StatsBase.sample(rng, get_actions(env), StatsBase.Weights(μ))
        _, state, reward, terminal = step!(env, action)

        experience = (state_mod(prev_state), state_mod(state), action, reward, 0.9, terminal,
                      μ[action], π[action], π[action]/sample_policy[action])
        add!(WER, experience, sample_policy[action]/μ[action])

        prev_state = copy(state)
        if terminal
            _, prev_state = start!(env; rng=rng)
        end
    end

    return WER
end

end


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

using Markdown
using InteractiveUtils

# ╔═╡ aab4db99-8eaf-4b70-add7-d470e6b6b6e2
using Plots, LaTeXStrings, DifferentialEquations, Distributions

# ╔═╡ 460194e8-add2-11eb-0f92-27b8fffb5046
md"# Chapter 1 Figures. Statistics"

# ╔═╡ de27faec-4e2b-4520-9d3f-1ac15344bab1
random_count() = floor(Int, n_min + rand() * n_max)

# ╔═╡ a8fa0a9a-4fb7-4380-87dd-2534687409c4
begin
	affect1!(integrator) = integrator.u += 1
	affect2!(integrator) = integrator.u -= 1
	prod1Jump(A) = ConstantRateJump((u, p, t) -> A, affect1!)
	deg1Jump = ConstantRateJump((u, p, t) -> u, affect2!)
	u01 = 0
	t1 = (0.0, 15.0)
	dprob1 = DiscreteProblem(u01, t1)
	jprob1 = JumpProblem(dprob1, Direct(), prod1Jump(100), deg1Jump)
	continous(A) = t -> -A*ℯ^(-t) + A
end

# ╔═╡ 343e56a2-4837-49a6-aa69-5a285d362c03
sol1 = solve(jprob1, SSAStepper())

# ╔═╡ 238a13fa-6d94-448d-a236-fdde5b78ad98
begin
	plot(sol1, label="Stochastic", legend=:topleft)
	p1 = plot!(continous(100), sol1.t, label="Continuous", xlabel=L"t", ylabel=L"n", ylimit=(0,145))
end

# ╔═╡ 67a22053-294a-47d1-9793-2188398de91c
begin
	last_ones = []
	num_samples=10000
	for _ ∈ 1:num_samples
		push!(last_ones, solve(jprob1, SSAStepper()).u[end])
	end
	last_ones
end

# ╔═╡ 67f8f4d0-0e6d-4857-9ac6-ea3b9f1eca3f
begin
	# poisson(mean)
	histogram(last_ones, label="Distribution of final values", color=:orange, orientation=:h)
	p2 = plot!(rand(Poisson(100), num_samples), seriestype=:stephist, label="Sampled Poisson Distribution", color=:blue, orientation=:h, ylimit=(0,145), ticks=nothing, yaxis=false, xaxis=false)
end

# ╔═╡ 6c9d9245-317a-47ca-83c4-b00523eaf3fd
pf = plot(p1, p2)

# ╔═╡ 2be2984e-1315-41a5-b7c4-63e79fad18f2
savefig(pf, "ch1_stochastic.pdf")

# ╔═╡ Cell order:
# ╠═aab4db99-8eaf-4b70-add7-d470e6b6b6e2
# ╟─460194e8-add2-11eb-0f92-27b8fffb5046
# ╠═de27faec-4e2b-4520-9d3f-1ac15344bab1
# ╠═a8fa0a9a-4fb7-4380-87dd-2534687409c4
# ╠═343e56a2-4837-49a6-aa69-5a285d362c03
# ╠═238a13fa-6d94-448d-a236-fdde5b78ad98
# ╠═67a22053-294a-47d1-9793-2188398de91c
# ╠═67f8f4d0-0e6d-4857-9ac6-ea3b9f1eca3f
# ╠═6c9d9245-317a-47ca-83c4-b00523eaf3fd
# ╠═2be2984e-1315-41a5-b7c4-63e79fad18f2


text:
# water to soil routine
# snow pack
# water content in snow
# meltwater
# refreezing meltwater
# boolean (snow at start of timestep?)

# Q_box0, Q_box1, Q_box2 - Q from box1 och box2
# SUZ, SLZ - storage in box1 och box2
# Q_gen - Q before MAXBAS

# Packages

using Distributions
using PyPlot
using CSV

# Function for computing ordinates of unit hydrograph

function compute_hbv_ord(maxbas, tstep)

  maxbas = maxbas / tstep

  triang = Distributions.TriangularDist(0, maxbas)

  triang_cdf = Distributions.cdf(triang, 0:ceil(Int64, maxbas + 2))

  hbv_ord = diff(triang_cdf)

  return(hbv_ord)

end

# Compute monthly potential evapotranspiration

function epot_monthly(datevec)

    # Monthly values of potential evapotranspiration (mm/day)

    epot_month = [0.05, 0.14, 0.46, 1.5, 3.01, 4.15, 3.66, 2.72, 1.42, 0.43, 0.03, 0.0]

    # Assign montly values to days

    epot = zeros(Float64, length(datevec))

    for i in eachindex(datevec)

        imonth = Dates.month(datevec[i])

        epot[i] = epot_month[imonth]

    end

    return(epot)

end

# Input data

data = CSV.read("C:\\Work\\Studies\\vann\\hbv_light\\ptq.txt",
                types = [DateTime, Float64, Float64, Float64],
                delim = "\t", dateformat = "yyyymmdd", header = true)

date = Array(data[:date])
Prec = Array(data[:Prec])
temp = Array(data[:Temp])
Qobs = Array(data[:Qsim])

Prec = transpose(Prec)
temp = transpose(temp)

pET = epot_monthly(date)

# Settings

UseOldSUZ = true

# Elevation and vegetation zones

lake = 0.0
nelev = size(Prec, 1)
nveg = 1
area = ones(nveg, nelev)
area = area / sum(area)

# Parameter values

PERC = 0.7
K0 = 0.2
K1 = 0.08
K2 = 0.03
UZL = 20.0

par_TT = [-1.0]
par_CFMAX = [5.0]
par_SFCF = [0.8]
par_CFR = [0.05]
par_CWH = [0.1]
par_FC = [250.0]
par_LP = [0.7]
par_BETA = [3.0]

par_MAXBAS = 2.5*24.0
tstep = 24.0

maxTS = size(Prec,2)

# Initilize state variables

state_SP = zeros(nveg, nelev)
state_WC = zeros(nveg, nelev)
state_SM = zeros(nveg, nelev)

# Temporary

runoff_sim = zeros(size(Prec,2))
snow_sim = zeros(size(Prec,2))
SM_sim = zeros(size(Prec,2))
SUZ_sim = zeros(size(Prec,2))
SLZ_sim = zeros(size(Prec,2))
aET_sim = zeros(size(Prec,2))
recharge_sim = zeros(size(Prec,2))
Q_com = zeros(size(Prec,2))

# Compute maxbas weights

maxbas_w = compute_hbv_ord(par_MAXBAS, tstep)

# Initilize state variables

SUZ = 0.0
SLZ = min(PERC / K2, Qobs[1] / K2)
Q_box0 = 0.0

for t = 1:size(Prec, 2)

    pot_E = pET[t] * (1 - lake)

    till_Qsum = 0.0
    avg_act_E = 0.0

    for elevzone = 1:nelev

        P_zone = Prec[elevzone, t]
        T_zone = temp[elevzone, t]

        for vegzone = 1:nveg

            if area[vegzone, elevzone] > 0.0

                SP = state_SP[vegzone, elevzone]
                WC = state_WC[vegzone, elevzone]
                SM = state_SM[vegzone, elevzone]

                tt = par_TT[vegzone]
                CFMAX = par_CFMAX[vegzone]
                SFCF = par_SFCF[vegzone]
                CFR = par_CFR[vegzone]
                CWH = par_CWH[vegzone]
                FC = par_FC[vegzone]
                LP = par_LP[vegzone]
                BETA = par_BETA[vegzone]

                snow = (SP > 0)

                #snow routine

                insoil = 0.0
                if SP > 0.0

                    if P_zone > 0.0
                        if T_zone > tt
                            WC = WC + P_zone
                        else
                            SP = SP + P_zone * SFCF
                        end
                    end # if P_zone

                    if T_zone > tt
                        melt = CFMAX * (T_zone - tt)
                        if melt > SP
                            insoil = SP + WC
                            WC = 0.0
                            SP = 0.0
                        else
                            SP = SP - melt
                            WC = WC + melt
                            if WC >= CWH * SP
                                insoil = WC - CWH * SP
                                WC = CWH * SP
                            end
                        end
                    else
                        refrez = CFR * CFMAX * (tt - T_zone)
                        if refrez > WC
                            refrez = WC
                        end
                        SP = SP + refrez
                        WC = WC - refrez
                    end # if T_zone
                else
                    if T_zone > tt
                        insoil = P_zone
                    else
                        SP = P_zone * SFCF
                    end # if T_zone
                end # if SP

                #soil routine

                till_Q = 0.0
                old_SM = SM
                if insoil > 0.0
                    if insoil < 1.0
                        y = insoil
                    else
                        m = floor(insoil)   # IS THIS CORRECT
                        y = insoil - m
                        for i in 1:m
                            dQdP = (SM / FC) ^ BETA
                            if dQdP > 1.0
                                dQdP = 1.0
                            end
                            SM = SM + 1.0 - dQdP
                            till_Q = till_Q + dQdP
                        end
                    end
                    dQdP = (SM / FC) ^ BETA
                    if dQdP > 1.0
                        dQdP = 1.0
                    end
                    SM = SM + (1 - dQdP) * y
                    till_Q = till_Q + dQdP * y
                end # if insoil

                mean_SM = (SM + old_SM) / 2.0
                if mean_SM < (LP * FC)
                    act_E = pot_E * mean_SM / (LP * FC)
                else
                    act_E = pot_E
                end
                if snow
                    act_E = 0.0
                end
                SM = SM - act_E
                if SM < 0.0
                    SM = 0.0
                end

                avg_act_E = avg_act_E + act_E * area[vegzone, elevzone]
                till_Qsum = till_Qsum + till_Q * area[vegzone, elevzone]
                state_SP[vegzone, elevzone] = SP
                state_WC[vegzone, elevzone] = WC
                state_SM[vegzone, elevzone] = SM

            end # if area

        end # veg loop

    end # elev loop

    till_box1 = till_Qsum

    # generation of runoff
    SUZ = SUZ + till_box1
    if ( SUZ - PERC / (1 - lake) ) < 0.0
        SLZ = SLZ + SUZ * (1.0 - lake)
        SUZ = 0.0
    else
        SLZ = SLZ + PERC
        SUZ = SUZ - PERC / (1.0 - lake)
    end

    tt = mean(par_TT)   # CHECK WITH JAN WHAT THIS MEANS
    SFCF = mean(SFCF)   # CHECK WITH JAN WHAT THIS MEANS

    if mean(temp[:,t]) > tt   # CHECK WITH JAN WHAT THIS MEANS
        SLZ = max(SLZ - pET[t] * lake, 0.0)
        avg_act_E = avg_act_E + min(SLZ, pET[t] * lake)
    end
    if mean(temp[:,t]) <= tt   # CHECK WITH JAN WHAT THIS MEANS
        SLZ = SLZ + SFCF * mean(Prec[:,t]) * lake   # CHECK WITH JAN WHAT THIS MEANS
    else
        SLZ = SLZ + mean(Prec[:,t]) * lake   # CHECK WITH JAN WHAT THIS MEANS
    end
    if UseOldSUZ
        Q_box1 = K1 * SUZ
        if SUZ < UZL
            Q_box0 = 0.0
        else
            Q_box0 = K0 * (SUZ - UZL)
        end
    else
        Q_box1 = min(K1 * SUZ ^ (1 + UZL), SUZ)
    end
    Q_box2 = K2 * SLZ
    SUZ = SUZ - Q_box1 - Q_box0
    SLZ = SLZ - Q_box2
    Q_gen = Q_box1 + Q_box2 + Q_box0

    # transformation of runoff

    for i = 0:(length(maxbas_w)-1)
        if t + i <= maxTS
            Q_com[t + i] = Q_com[t + i] + Q_gen * maxbas_w[i+1]
        end
    end

    avg_SP = 0.0
    avg_SM = 0.0
    avg_WC = 0.0
    for elevzone = 1:nelev
        for vegzone = 1:nveg
            avg_SP = avg_SP + state_SP[vegzone, elevzone] * area[vegzone, elevzone]
            avg_WC = avg_WC + state_WC[vegzone, elevzone] * area[vegzone, elevzone]
            avg_SM = avg_SM + state_SM[vegzone, elevzone] * area[vegzone, elevzone]
        end
    end

    runoff_sim[t] = Q_gen
    snow_sim[t] = (avg_SP + avg_WC) / (1 - lake)
    SM_sim[t] = avg_SM / (1 - lake)
    SUZ_sim[t] = SUZ
    SLZ_sim[t] = SLZ
    aET_sim[t] = avg_act_E
    recharge_sim[t] = till_Qsum

end  # t loop

plot(date, Q_com, label = "Julia code")
plot(date, Qobs, label = "Original code")
ylabel("Runoff (mm/day)")

@show runoff_sim

#=plot(snow_sim)
plot(SM_sim)
plot(SUZ_sim)
plot(SLZ_sim)
plot(aET_sim)
plot(recharge_sim)=#

max_stars_repo_name: jmgnve/Vann

max_forks_repo_name: jmgnve/Vann


text:
module DiscreteChoice

using LinearAlgebra

import ForwardDiff

export optimize,
       armijo,
       tcg,
       BFGS!,
       inv_BFGS,
       SR1!

include("optimizers/gradientdescent.jl")
include("optimizers/newton.jl")
include("optimizers/quasinewton.jl")
include("optimizers/trustregion.jl")
include("utils/hessian.jl")

end # module


text:
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
using Plasm

VOID = [[]], [[1]]

# example 1

V = [[0,0] [1,0] [2,0] [0,1] [1,1] [2,1] [0,2] [1,2] [2,2]];
FV = [[1,2,4],[2,3,5],[3,5,6],[4,5,7],[5,7,8],[6,8,9]];
pattern = repeat([1,2,-3],outer=4);
model = (V,FV);
W,FW = Lar.extrudeSimplicial(model, pattern);
Plasm.view(W,FW)

# example 2

model = Lar.extrudeSimplicial( VOID, ones(10) )
Plasm.view(model)
model = Lar.extrudeSimplicial( model, ones(10) )
Plasm.view(model)
model = Lar.extrudeSimplicial( model, ones(10) )
Plasm.view(model)

# example 3

model = Lar.extrudeSimplicial( VOID, repeat([1,-1],outer=10) )
Plasm.view(model)
model = Lar.extrudeSimplicial( model, repeat([1,-1],outer=10) )
Plasm.view(model)

# example 4

grid_2d = Lar.simplexGrid([3,3])
Plasm.view(grid_2d)
grid_3d = Lar.simplexGrid([2,3,4])
Plasm.view(grid_3d)
V,CV = Lar.simplexGrid([1,1,1])
Plasm.view(V,CV)

# example 5

SK2 = Lar.simplexFacets(CV)
Plasm.view(V, SK2)
SK1 = Lar.simplexFacets(SK2)
Plasm.view(V, SK1)


text:
module Plotting

using PyPlot
using ClimateMARGO.Models
using ClimateMARGO.Diagnostics

export
    plot_controls, plot_emissions, plot_concentrations,
    plot_temperatures, plot_benefits, plot_damages,
    plot_state

include("line_plots.jl")

end


text:
"""
    accumulated_edges(adj)

Return a vector which acts as a mapping table. The index is the vertex index,
value is accumulated numbers of edge (current vertex not included).
"""
accumulated_edges(adj::AbstractVector{<:AbstractVector{<:Integer}}) = [0, cumsum(map(length, adj))...]

Zygote.@nograd accumulated_edges

Zygote.@nograd function generate_cluster(M::AbstractArray{T,N}, accu_edge) where {T,N}
    num_V = length(accu_edge) - 1
    num_E = accu_edge[end]
    cluster = similar(M, Int, num_E)
    @inbounds for i = 1:num_V
        j = accu_edge[i]
        k = accu_edge[i+1]
        cluster[j+1:k] .= i
    end
    cluster
end

"""
    edge_index_table(adj[, directed])

Generate a mapping from vertex pair (i, j) to edge index. The edge indecies are determined by
the sorted vertex indecies.
"""
function edge_index_table(adj::AbstractVector{<:AbstractVector{<:Integer}}, directed::Bool=is_directed(adj))
    table = Dict{Tuple{UInt32,UInt32},UInt64}()
    e = one(UInt64)
    if directed
        for (i, js) = enumerate(adj)
            js = sort(js)
            for j = js
                table[(i, j)] = e
                e += one(UInt64)
            end
        end
    else
        for (i, js) = enumerate(adj)
            js = sort(js)
            js = js[i .≤ js]
            for j = js
                table[(i, j)] = e
                table[(j, i)] = e
                e += one(UInt64)
            end
        end
    end
    table
end

function edge_index_table(vpair::AbstractVector{<:Tuple})
    table = Dict{Tuple{UInt32,UInt32},UInt64}()
    for (i, p) = enumerate(vpair)
        table[p] = i
    end
    table
end

edge_index_table(fg::FeaturedGraph) = edge_index_table(fg.graph, fg.directed)

Zygote.@nograd edge_index_table

### TODO move these to GraphSignals ######
import GraphSignals: FeaturedGraph

function FeaturedGraph(fg::FeaturedGraph;
                        nf=node_feature(fg),
                        ef=edge_feature(fg),
                        gf=global_feature(fg))

    return FeaturedGraph(graph(fg); nf, ef, gf)
end

function check_num_nodes(fg::FeaturedGraph, x::AbstractArray)
    @assert nv(fg) == size(x, ndims(x))
end

function check_num_edges(fg::FeaturedGraph, e::AbstractArray)
    @assert ne(fg) == size(e, ndims(e))
end


text:
function rollout!(solver::iLQRSolver{T,Q,n}, α) where {T,Q,n}
    Z = solver.Z; Z̄ = solver.Z̄
    K = solver.K; d = solver.d;

    Z̄[1].z = [solver.x0; control(Z[1])]

    temp = 0.0
	δx = solver.S[end].q
	δu = solver.S[end].r

    for k = 1:solver.N-1
        δx .= RobotDynamics.state_diff(solver.model, state(Z̄[k]), state(Z[k]))
		δu .= d[k] .* α
		mul!(δu, K[k], δx, 1.0, 1.0)
        ū = control(Z[k]) + δu
        RobotDynamics.set_control!(Z̄[k], ū)

        # Z̄[k].z = [state(Z̄[k]); control(Z[k]) + δu]
        Z̄[k+1].z = [RobotDynamics.discrete_dynamics(Q, solver.model, Z̄[k]);
            control(Z[k+1])]

        max_x = norm(state(Z̄[k+1]),Inf)
        if max_x > solver.opts.max_state_value || isnan(max_x)
            solver.stats.status = STATE_LIMIT
            return false
        end
        max_u = norm(control(Z̄[k+1]),Inf)
        if max_u > solver.opts.max_control_value || isnan(max_u)
            solver.stats.status = CONTROL_LIMIT
            return false
        end
    end
    solver.stats.status = UNSOLVED
    return true
end

"Simulate the forward the dynamics open-loop"
function rollout!(solver::iLQRSolver{<:Any,Q}) where Q
    rollout!(Q, solver.model, solver.Z, SVector(solver.x0))
    for k in eachindex(solver.Z)
        solver.Z̄[k].t = solver.Z[k].t
    end
end


text:
include("Dataset.jl")
import .Dataset:
  get_MINST,
  TrainTest,
  split_X_by_batches,
  test_viz_idx,
  min_distance_digits

include("Model.jl")
import .Model:
  encoder,
  decoder,
  split_encoder_result,
  create_vae

include("Utils.jl")
import .Utils:
  gen_images,
  img,
  safe_img_save

using Test
using Flux
using Flux.Tracker: TrackedReal
using Printf
using Images

function experiment1()
  # CONFIG
  # TODO: break this out, after complete spec is determined
  n_traning_examples = 60000
  batch_size = 128
  n_latent = 50
  max_epoch_idx = 20
  output_dir = joinpath(pwd(), "output/1/")
  if !isdir(output_dir) # make output directory, (git ignored)
    mkdir(output_dir)
  end

  # DATASET
  dataset = get_MINST()
  X_batches = split_X_by_batches(dataset.train_x, batch_size)
  viz_idx = min_distance_digits()
  [@info "viz_idx: " idx for idx in viz_idx]

  x_view = dataset.test_x[:,:,:,viz_idx]
  # save test images (ignored by git)
  x_imgs_test = hcat(img.([x_view[:,:,:,i] for i = 1:size(x_view,4)])...)
  safe_img_save(joinpath(output_dir, "test.png"), x_imgs_test)

  # MODEL
  ps, loss_fn, f, g = create_vae(n_latent, batch_size)
  opt = ADAM()

  for epoch_idx = range(1, stop = max_epoch_idx)
    for (data_idx, dataset_batch) in enumerate(X_batches)
      start = time()
      Flux.train!(loss_fn, ps, [dataset_batch], opt)
      loss_test_set = loss_fn(dataset.test_x[:, :, :, 1:batch_size])
      @info @sprintf("[%d][%d] test loss : %.4f", epoch_idx, data_idx, loss_test_set)
      @info @sprintf("[%d][%d] elapsed %d(s)", epoch_idx, data_idx, time() - start)
      if data_idx % 5 == 1
        epoch_img = joinpath(output_dir, @sprintf("img_epoch=%d_batch=%d.png", epoch_idx, data_idx))
        gen_images(epoch_img, g, f, x_view)
      end
    end
  end
end

experiment1()

max_stars_repo_name: adamwespiser/variational-autoencoders

max_issues_repo_name: adamwespiser/variational-autoencoders

max_forks_repo_name: adamwespiser/variational-autoencoders


text:
"""
    postprocess_lonlat(sol,𝑃::NamedTuple; id=missing, 𝑇=missing)

Copy `sol` to a `DataFrame` & map position to lon,lat coordinates
using "exchanged" 𝑃.XC, 𝑃.YC via `add_lonlat!`
"""
function postprocess_lonlat(sol,𝑃::NamedTuple; id=missing, 𝑇=missing)
    ismissing(id) ? id=collect(1:size(sol,2)) : nothing
    ismissing(𝑇) ? 𝑇=𝑃.𝑇 : nothing

    nd=length(size(sol))
    nt=size(sol,nd)
    nf=size(sol,nd-1)
    id=id*ones(1,size(sol,nd))
    if (size(sol,1)>1)&&(nd>2)
        x=sol[1,:,:]
        y=sol[2,:,:]
        fIndex=sol[end,:,:]
    elseif (nd>2)
        x=[sol[1,i,j][1] for i in 1:nf, j in 1:nt]
        y=[sol[1,i,j][2] for i in 1:nf, j in 1:nt]
        fIndex=[sol[1,i,j][end] for i in 1:nf, j in 1:nt]
    else
        x=sol[1,:]
        y=sol[2,:]
        fIndex=sol[end,:]
        nf=1
    end

    𝑃.XC.grid.nFaces==1 ? fIndex=ones(size(x)) : nothing

    t=[ceil(i/nf)-1 for i in 1:nt*nf]
    t=𝑇[1] .+ (𝑇[2]-𝑇[1])/t[end].*t

    df = DataFrame(ID=Int.(id[:]), x=x[:], y=y[:], fid=Int.(fIndex[:]), t=t[:])
    add_lonlat!(df,𝑃.XC,𝑃.YC)
    return df
end

"""
    add_lonlat!(df::DataFrame,XC,YC)

Add lon & lat to dataframe using "exchanged" XC, YC
"""
function add_lonlat!(df::DataFrame,XC,YC)
    x=df[!,:x];
    y=df[!,:y];
    f=Int.(df[!,:fid]);
    dx,dy=(x - floor.(x) .+ 0.5,y - floor.(y) .+ 0.5);
    i_c = Int32.(floor.(x)) .+ 1;
    j_c = Int32.(floor.(y)) .+ 1;

    lon=zeros(length(x),4)
    [lon[k,:]=XC[f[k]][i_c[k]:i_c[k]+1,j_c[k]:j_c[k]+1][:] for k in 1:length(i_c)]
    lat=zeros(length(x),4)
    [lat[k,:]=YC[f[k]][i_c[k]:i_c[k]+1,j_c[k]:j_c[k]+1][:] for k in 1:length(i_c)]

    k=findall(vec(maximum(lon,dims=2)-minimum(lon,dims=2)) .> 180.0)
    tmp=view(lon,k,:)
    tmp[findall(tmp.<0.0)]=tmp[findall(tmp.<0.0)] .+ 360.0

    df.lon=(1.0 .-dx).*(1.0 .-dy).*lon[:,1]+dx.*(1.0 .-dy).*lon[:,2] +
         (1.0 .-dx).*dy.*lon[:,3]+dx.*dy.*lon[:,4]
    df.lat=(1.0 .-dx).*(1.0 .-dy).*lat[:,1]+dx.*(1.0 .-dy).*lat[:,2] +
         (1.0 .-dx).*dy.*lat[:,3]+dx.*dy.*lat[:,4]

    return df
end

"""
    postprocess_xy()

Copy `sol` to a `DataFrame` & map position to x,y coordinates,
and define time axis for a simple doubly periodic domain
"""
function postprocess_xy(sol,𝑃::NamedTuple; id=missing, 𝑇=missing)
    ismissing(id) ? id=collect(1:size(sol,2)) : nothing
    ismissing(𝑇) ? 𝑇=𝑃.𝑇 : nothing

    nf=size(sol,2)
    nt=size(sol,3)
    nx,ny=𝑃.ioSize[1:2]
    nd=length(size(sol))

    id=id*ones(1,size(sol,nd))
    if (size(sol,1)>1)&&(nd>2)
        x=mod.(sol[1,:,:],Ref(nx))
        y=mod.(sol[2,:,:],Ref(ny))
    elseif (nd>2)
        x=[mod(sol[1,i,j][1],nx) for i in 1:nf, j in 1:nt]
        y=[mod(sol[1,i,j][2],ny) for i in 1:nf, j in 1:nt]
    else
        x=mod.(sol[1,:],Ref(nx))
        y=mod.(sol[2,:],Ref(ny))
    end
    t=[ceil(i/nf)-1 for i in 1:nt*nf]
    #size(𝑃.XC,1)>1 ? fIndex=sol[3,:,:] : fIndex=fill(1.0,size(x))
    t=𝑇[1] .+ (𝑇[2]-𝑇[1])/t[end].*t

    return DataFrame(ID=Int.(id[:]), t=t[:], x=x[:], y=y[:])
end

"""
    set_up_𝑃(k::Int,t::Float64,Γ::Dict,pth::String)

Define the `𝑃` _parameter_ tuple given grid variables `Γ`, vertical level
choice `k`, time `t` in `seconds`, and velocity fields obtained from
files in `pth`.

The two climatological months (`m0`,`m1`) that bracket time `t` are
read to memory (e.g. months 12 & 1 then 1 & 2 and so on).

_Note: the initial implementation approximates every month duration
to 365 days / 12 months for simplicity._
"""
function set_up_𝑃(k::Int,t::Float64,Γ::Dict,pth::String)
    XC=exchange(Γ["XC"]) #add 1 lon point at each edge
    YC=exchange(Γ["YC"]) #add 1 lat point at each edge
    iDXC=1. ./Γ["DXC"]
    iDYC=1. ./Γ["DYC"]
    γ=Γ["XC"].grid
    mon=86400.0*365.0/12.0

    𝑃 = (u0=MeshArray(γ,Float64), u1=MeshArray(γ,Float64),
         v0=MeshArray(γ,Float64), v1=MeshArray(γ,Float64),
         𝑇=[-mon/2,mon/2], 🔄 = update_𝑃!, pth=pth,
         XC=XC, YC=YC, iDXC=iDXC, iDYC=iDYC,
         k=k, msk=Γ["hFacC"][:, k])

    tmp = dict_to_nt(IndividualDisplacements.NeighborTileIndices_cs(Γ))
    𝑃 = merge(𝑃 , tmp)

    𝑃.🔄(𝑃,0.0)
    return 𝑃
end

"""
    update_𝑃!(𝑃::NamedTuple,t::Float64)

Update input data (velocity arrays) and time period (array) inside 𝑃 (𝑃.u0[:], etc, and 𝑃.𝑇[:])
based on the chosen time `t` (in `seconds`).

_Note: for now, it is assumed that (1) input 𝑃.𝑇 is used to infer `dt` between consecutive velocity fields,
(2) periodicity of 12 monthly records, (3) vertical 𝑃.k is selected -- but this could easily be generalized._
"""
function update_𝑃!(𝑃::NamedTuple,t::Float64)
    dt=𝑃.𝑇[2]-𝑃.𝑇[1]

    m0=Int(floor((t+dt/2.0)/dt))
    m1=m0+1
    t0=m0*dt-dt/2.0
    t1=m1*dt-dt/2.0

    m0=mod(m0,12)
    m0==0 ? m0=12 : nothing
    m1=mod(m1,12)
    m1==0 ? m1=12 : nothing

    (U,V)=read_velocities(𝑃.u0.grid,m0,𝑃.pth)
    u0=U[:,𝑃.k]; v0=V[:,𝑃.k]
    u0[findall(isnan.(u0))]=0.0; v0[findall(isnan.(v0))]=0.0 #mask with 0s rather than NaNs
    u0=u0.*𝑃.iDXC; v0=v0.*𝑃.iDYC; #normalize to grid units
    (u0,v0)=exchange(u0,v0,1) #add 1 point at each edge for u and v

    (U,V)=read_velocities(𝑃.u0.grid,m1,𝑃.pth)
    u1=U[:,𝑃.k]; v1=V[:,𝑃.k]
    u1[findall(isnan.(u1))]=0.0; v1[findall(isnan.(v1))]=0.0 #mask with 0s rather than NaNs
    u1=u1.*𝑃.iDXC; v1=v1.*𝑃.iDYC; #normalize to grid units
    (u1,v1)=exchange(u1,v1,1) #add 1 point at each edge for u and v

    𝑃.u0[:]=u0[:]
    𝑃.u1[:]=u1[:]
    𝑃.v0[:]=v0[:]
    𝑃.v1[:]=v1[:]
    𝑃.𝑇[:]=[t0,t1]

end

"""
    initialize_gridded(𝑃::NamedTuple,n_subset::Int=1)

Define initial condition (u0,du) as a subset of grid points
"""
function initialize_gridded(𝑃::NamedTuple,n_subset::Int=1)
    msk=𝑃.msk
    uInitS = Array{Float64,2}(undef, 3, prod(msk.grid.ioSize))

    kk = 0
    for fIndex = 1:length(msk.fSize)
        nx, ny = msk.fSize[fIndex]
        ii1 = 0.5:1.0:nx
        ii2 = 0.5:1.0:ny
        n1 = length(ii1)
        n2 = length(ii2)
        for i1 in eachindex(ii1)
            for i2 in eachindex(ii2)
                if msk[fIndex][Int(round(i1+0.5)),Int(round(i2+0.5))]
                    kk += 1
                    let kk = kk
                        uInitS[1, kk] = ii1[i1]
                        uInitS[2, kk] = ii2[i2]
                        uInitS[3, kk] = fIndex
                    end
                end
            end
        end
    end

    uInitS=uInitS[:,1:kk]
    du=fill(0.0,size(uInitS));

    uInitS=uInitS[:,1:n_subset:end]
    du=du[:,1:n_subset:end]
    return uInitS,du
end

"""
    randn_lonlat(nn=1,seed=missing)

Randomly distributed longitude, latitude positions on the sphere.
"""
function randn_lonlat(nn=1;seed=missing)
    !ismissing(seed) ? rng = MersenneTwister(1234) : rng = MersenneTwister()
    tmp = randn(rng, Float64, (nn, 3))
    tmpn = tmp ./ sqrt.(sum(tmp.^2, dims=2))
    lon = rad2deg.(atan.(tmpn[:,2], tmpn[:,1]))
    lat = 90.0 .- rad2deg.(acos.(tmpn[:,3]))
    return lon, lat
end

"""
    initialize_lonlat(Γ,lon,lat ; msk=missing)

Define initial condition (u0,du) in grid coordinates (Γ) from longitude
& latitude vectors (lon,lat) optionally with a land mask (msk).
"""
function initialize_lonlat(Γ::Dict,lon::Array{Float64,1},lat::Array{Float64,1};msk=missing)
    (f,i,j,w,j_f,j_x,j_y)=InterpolationFactors(Γ,lon,lat)
    ii=findall( ((!isnan).(j_x)).&(j_f.!==0) )
    if !ismissing(msk)
        jj=[msk[Int(j_f[i]),1][ Int(round(j_x[i] .+ 0.5)), Int(round(j_y[i] .+ 0.5)) ] for i in ii]
        ii=ii[findall(jj.>0.0)]
    end
    u0=Array(transpose([j_x[ii] j_y[ii] j_f[ii]]))
    du=similar(u0)
    return u0,du
end

initialize_lonlat(Γ::Dict,lon::Float64,lat::Float64;msk=missing) = initialize_lonlat(Γ,[lon],[lat];msk=msk)

"""
    interp_to_lonlat

Use MeshArrays.Interpolate() to interpolate to e.g. a regular grid (e.g. maps for plotting purposes).

```
using MeshArrays, IndividualDisplacements

lon=[i for i=19.5:1.0:379.5, j=-78.5:1.0:78.5]
lat=[j for i=19.5:1.0:379.5, j=-78.5:1.0:78.5]
(f,i,j,w,_,_,_)=InterpolationFactors(Γ,vec(lon),vec(lat))
IntFac=(lon=lon,lat=lat,f=f,i=i,j=j,w=w)

D=Γ["Depth"]
tmp1=interp_to_lonlat(D,Γ,lon,lat)
tmp2=interp_to_lonlat(D,IntFac)
```
"""
function interp_to_lonlat(X::MeshArray,Γ::Dict,lon,lat)
    (f,i,j,w,_,_,_)=InterpolationFactors(Γ,vec(lon),vec(lat))
    return reshape(Interpolate(X,f,i,j,w),size(lon))
end

function interp_to_lonlat(X::MeshArray,IntFac::NamedTuple)
    @unpack f,i,j,w,lon,lat = IntFac
    return reshape(Interpolate(X,f,i,j,w),size(lon))
end

"""
    interp_to_xy(df::DataFrame,Zin)

Interpolate "exchanged" / "hallo-included" Zin to df[!,:x], df[!,:y] on df[!,:fid]
"""
function interp_to_xy(df::DataFrame,Zin)
    x=df[!,:x];
    y=df[!,:y];
    f=Int.(df[!,:fid]);
    dx,dy=(x - floor.(x) .+ 0.5,y - floor.(y) .+ 0.5);
    i_c = Int32.(floor.(x)) .+ 1;
    j_c = Int32.(floor.(y)) .+ 1;

    Z=zeros(length(x),4)
    [Z[k,:]=Zin[f[k]][i_c[k]:i_c[k]+1,j_c[k]:j_c[k]+1][:] for k in 1:length(i_c)]

    return (1.0 .-dx).*(1.0 .-dy).*Z[:,1]+dx.*(1.0 .-dy).*Z[:,2] +
           (1.0 .-dx).*dy.*Z[:,3]+dx.*dy.*Z[:,4]
end


text:
function mahalanobisSquaredBetweenPairs(pairs::Matrix; covmatrix=nothing)
    n, p = size(pairs)
    newmat = zeros(Float64, n, n)
    if covmatrix === nothing
        covmatrix = cov(pairs)
    end
    try
        invm = inv(covmatrix)
        for i in 1:n
            @inbounds for j in i:n
                newmat[i, j] = ((pairs[i,:] .- pairs[j,:])' * invm * (pairs[i,:] .- pairs[j,:]))
                newmat[j, i] = newmat[i,j]
            end
        end
        return newmat
    catch e
        @warn e
        if det(covmatrix) == 0
            @warn "singular covariance matrix, mahalanobis distances can not be calculated"
        end
        return zeros(Float64, (n, n))
    end
end

"""

    asm2000(setting)

Perform the Setan, Halim and Mohd (2000) algorithm for the given regression setting.

# Arguments
- `setting::RegressionSetting`: RegressionSetting object with a formula and dataset.

# Examples
```julia-repl
julia> reg0001 = createRegressionSetting(@formula(calls ~ year), phones);
julia> asm2000(reg0001)
Dict{Any,Any} with 1 entry:
  "outliers" => [15, 16, 17, 18, 19, 20]
```

# References
Setan, Halim, and Mohd Nor Mohamad. "Identifying multiple outliers in
linear regression: Robust fit and clustering approach." (2000).
"""
function asm2000(setting::RegressionSetting)
    X = designMatrix(setting)
    y = responseVector(setting)
    return asm2000(X, y)
end

function asm2000(X::Array{Float64,2}, y::Array{Float64,1})
    n, p = size(X)
    h = floor((n + p - 1) / 2)
    ltsreg = lts(X, y)

    betas = ltsreg["betas"]
    hsubset = ltsreg["hsubset"]

    predicteds = [sum(X[i,:] .* betas) for i in 1:n]
    resids = y .- predicteds
    stdres = standardize(ZScoreTransform, resids, dims=1)
    stdfit = standardize(ZScoreTransform, predicteds, dims=1)
    pairs = hcat(stdfit, stdres)

    pairs = hcat(resids, predicteds)

    covmatrix = cov(pairs[hsubset, :])
    mahdist = mahalanobisSquaredBetweenPairs(pairs, covmatrix=covmatrix)

    outlierset = Array{Int,1}()

    hcl = hclust(mahdist, linkage=:single)
    majonacrit = majona(hcl)
    clustermappings = cutree(hcl, h=majonacrit)
    uniquemappings = unique(clustermappings)
    for clustid in uniquemappings
        cnt = count(x -> x == clustid, clustermappings)
        if cnt >= h
            outlierset = filter(i -> clustermappings[i] != clustid, 1:n)
        end
    end

    result = Dict()
    result["outliers"] = outlierset
    return result
end

max_stars_repo_name: akadal/LinRegOutliers

max_issues_repo_name: akadal/LinRegOutliers

max_forks_repo_name: akadal/LinRegOutliers


text:
"""
    loaddata(data_type, file_name, f_s, f_center, f_gnss, t_length;
             start_data_idx=1, site_lla=missing, data_start_time=missing)

Loads data from `scN` file and loads into `GNSSData` type struct.

`site_lla` is `[latitude, longitude, height]` in (degrees, degrees, meters).

`data_start_time` is `[year, month, day, hour, minute, seconds]`, where "ss" is
decimal seconds.

Required Arguments:

- `data_type`: one of the following for `Complex{IntN}` datasets
 	* `:sc4`
	* `:sc8`
	* `:sc16`
	* `:sc32`
	* `:sc64`
- `file_name::String`: data file name
- `f_s::Float64`: sampling rate of data in Hz
- `f_center::Float64`: center frequency of the receiver in Hz
- `f_gnss::Float64`: center frequency of the GNSS signal in Hz
- `t_length::Float64`: the length of data in seconds to be loaded

Optional Arguments:

- `skip_to::Float64`: where the first data sample will be loaded in the data
                      `(default = 0. sec)`
	* units are in seconds
- `site_lla`: receiver LLA location in the format `(lat, lon, height)`
	* `lat`, latitude, is in degrees
	* `lon`, longitude, is in degrees
	* `height` is in meters
	* `(default = missing)`
- `data_start_time`: data UTC start time `(year, month, day, hour, minute, seconds)`
	* `seconds` can be `Int` or `Float64`, all others are `Int`
	* `(default = missing)`
- `start_t`: the start time in seconds `(default = 0 seconds)`

Returns:

- `GNSSData` struct
"""
function loaddata(data_type, file_name, f_s, f_center, f_gnss, t_length;
                  skip_to=0., site_loc_lla=missing,
                  data_start_time=missing, start_t=0)
	# Compute number of samples to extract
	sample_num = Int(f_s * t_length)
	# Calculate f_if
	f_if = f_gnss - f_center
	# Determine nADC
	if data_type == :sc4
		data_type = Val(:sc4)
		nADC = 4
		data = Array{Complex{Int8}}(undef, sample_num)
	elseif data_type == :sc8
		data_type = Val(:sc8)
		nADC = 8
		data = Array{Complex{Int8}}(undef, sample_num)
	elseif data_type == :sc16
		data_type = Val(:sc16)
		nADC = 16
		data = Array{Complex{Int16}}(undef, sample_num)
	elseif data_type == :sc32
		data_type = Val(:sc32)
		nADC = 32
		data = Array{Complex{Int32}}(undef, sample_num)
	elseif data_type == :sc64
		data_type = Val(:sc64)
		nADC = 64
		data = Array{Complex{Int64}}(undef, sample_num)
	else
		error("Invalid data type: $(typeof(data_type))")
	end
	# Read data
	start_data_idx = floor(Int, f_s*skip_to)
	data, end_idx, dtype = readdatafile!(data, data_type, file_name,
	                                     sample_num, nADC, start_data_idx)
	# Calculate total data length in seconds
	total_data_length = (end_idx)/f_s
	# Generate time vector
	# t = calctvector(sample_num, f_s)
	return GNSSData(file_name, f_s, f_if, t_length, start_data_idx,
	                data, String(dtype), data_start_time,
	                site_loc_lla, sample_num, total_data_length, nADC, start_t)
end

"""
	reloaddata!(gnss_data::GNSSData, start_data_idx, sample_num)

Reloads portion of data up to the length of the data array inside `gnss_data`.

Required Arguments:

- `gnss_data::GNSSData`: structure that holds loaded data from data files
- `start_data_idx::Int`: where the first data sample will be loaded in the data
- `sample_num::Int`: number of samples after `start_data_idx` to load into
                     `gnss_data`
    * must not be greater than `gnss_data.sample_num`

Modifies in Place and Returns:

- `gnss_data::GNSSData`
"""
function reloaddata!(gnss_data::GNSSData, skip_to)
	f_s = gnss_data.f_s
	nADC = gnss_data.nADC
	start_data_idx = floor(Int, f_s*skip_to)
	sample_num=gnss_data.sample_num
	file_name = gnss_data.file_name
	data_type = gnss_data.data_type
	data = gnss_data.data
	if sample_num <= gnss_data.sample_num
        f = open(file_name, "r")
		byte_num = Int(2*nADC/8)
		seek(f, byte_num*start_data_idx)
		read!(f, data)
		close(f)
	else
		error("`sample_num` is greater than data sample size!")
	end
	return gnss_data
end

"""
	readdatafile!(data, data_type, file_name, sample_num, start_idx=1)

Loads `scN` data files. First N-bit number is real, second is imaginary.

Required Arguments:

- `data::Vector{Complex{IntN}}`: the complex data vector to store the new data
- `data_type::Val{:scN}`: only accepts `Val(:scN)` to indicate data is N-bit
                          complex
- `file_name::String`: data file name
- `sample_num::Int`: `[NOT USED]` number of samples after `start_data_idx` to
                     load

Optional Arguments:

- `start_idx::Int`: where the first data sample will be loaded in the data
                    `(default = 0)`

Returns:

- `data::Vector{Complex{IntN}}`: elements are replaced with new data
- `end_idx::Int`: the last sample index loaded plus one
- `:scN::Symbol`: symbol which describes the type of data
"""
function readdatafile!(data, data_type, file_name, sample_num,
                       nADC, start_idx=0)
	# Open file
	f = open(file_name, "r")
	# Go to start location
	byte_num = 2*Int(nADC/8)
	seek(f, byte_num*start_idx)
	# Load Complex{IntN} values directly from file
	read!(f, data)
	# Get the index value for the end of the file
	end_idx = Int(position(seekend(f))/byte_num)
    close(f)
    return (data, end_idx, gnsstypes[data_type])
end

"""
	readdatafile!(data, data_type::Val{:sc4}, file_name, sample_num,
				  start_idx=1)

Loads `sc4` data files. For each UInt8 number, the LSB is real and MSB is
imaginary.

**NOTE:** this method is slower than the version used for 8-bit complex data.

Required Arguments:

- `data::Vector{Complex{Int}}`: the complex data vector to store the new data
- `data_type::Val{:sc4}`: only accepts `Val(:sc4)` to indicate data is 4-bit
                          complex
- `file_name::String`: data file name
- `sample_num::Int`: `[NOT USED]` number of samples after `start_data_idx` to
                     load

Optional Arguments:

- `start_idx::Int`: where the first data sample will be loaded in the data
	                `(default = 0)`

Returns:

- `data::Vector{Complex{Int}}`: elements are replaced with new data
- `end_idx::Int`: the last sample index loaded plus one
- `:sc4::Symbol`: symbol which describes the type of data
"""
function readdatafile!(data, data_type::Val{:sc4}, file_name, sample_num,
                       nADC, start_idx=0)
	# Open file
	f = open(file_name, "r")
	# Go to start location
	seek(f, start_idx)
	# Load UInt8 values and convert to Complex{Int8}
	data[:] = bytetocomplex.(read(f, sample_num))
	# Get the index value for the end of the file
	end_idx = position(seekend(f))
    close(f)
    return (data, end_idx, :sc4)
end

"""
    bytetocomplex(byte::UInt8)

Converts a `UInt8` value to `Complex{Int8}`. Assumes that the LSB and MSB are the
real and imaginary parts, respectively.

Required Arguments:

- `byte::UInt8`: a single byte to converted to `Complex{Int8}`

Returns:

- `Complex{Int8}`: complex number converted from `byte`
"""
function bytetocomplex(byte::UInt8)
	# Get unsigned I and Q values
	Iuns = byte & 0x0f
	Quns = byte >> 0x04
	# Get I and Q signs
	Isign = (Iuns&0x08) >> 0x03
	Qsign = (Quns&0x08) >> 0x03
	# Get I and Q magnitudes
	Imag = xor(Isign*0x07, Iuns&0x07)
	Qmag = xor(Qsign*0x07, Quns&0x07)
	# Compute the complex number
	return (-1)^(Isign)*Imag + (-1)^(Qsign)*Qmag*1im
end

"""
	FileInfo

`Struct` that holds data file information based off its file name.

Fields:

- `f_s::Float64`: sampling rate of signal in Hz
- `f_if::Float64`: signal IF frequency in Hz
- `f_center::Float64`: receiver center frequency in Hz
- `sig_freq::Float64`: navigation band frequency in Hz
                       (e.g. for L1, `sig_freq = L1_freq`)
- `sigtype::T1`: `[DEPRICATED]`
- `data_type::T2`: the data type of the data file which is either set to
                   `Val(:sc4)` or `Val(:sc8)`
- `timestamp::T3`: six element `Tuple` with format
                   `(year, month, day, hour, minute, second)` where everything
				   other than `second`, which can be either `Int` or `Float64`,
				   is an Int
- `timestamp_JD::Float64`: the Julia date of the timestamp
- `file_name::String`: file name of the data file
"""
struct FileInfo{T1,T2,T3}
	f_s::Float64
	f_if::Float64
	f_center::Float64
	sig_freq::Float64
	sigtype::T1
	data_type::T2
	timestamp::T3
	timestamp_JD::Float64
	file_name::String
end

"""
	data_info_from_name(file_name)

Determine the IF and sampling frequency of the data from the file name. Looks
for the following strings inside the file name:

- `g1b1`: 25 Msps, centered at 1.57 GHz which may contain the followin signals
	* GPS L1
	* Galileo E1bc
	* Beidou B1i
- `g2r2`: `[NOTE SUPPORTED]` 25 Msps, centered at 1.2375 GHz which may contain
          the followin signals
	* GPS L2
	* GLONASS L2
- `g5`: 25 Msps, centered at 1.17645 GHz which may contain the followin signals
	* GPS L5
	* Galileo E5a

Also looks at the file name extension to determine the type of data stored in
the file where

- `.sc4`: 4-bit complex data where for each UInt8 number in the data file, the
  LSB is real and MSB is imaginary
- `.sc8`: 8-bit complex data where for each pair of Int8 numbers, the first
  value is real and the second is imaginary

Required Arguments:

- `file_name::String`: the data file name to parse

Returns:

- `FileInfo` struct
"""
function data_info_from_name(file_name)
	# Get file type (N-bit complex)
	if occursin("sc64", file_name)
		data_type = Val(:sc64)
	elseif occursin("sc32", file_name)
		data_type = Val(:sc32)
	elseif occursin("sc16", file_name)
		data_type = Val(:sc16)
	elseif occursin("sc8", file_name)
		data_type = Val(:sc8)
	elseif occursin("sc4", file_name)
		data_type = Val(:sc4)
	else
		error("File type not supported. Only sc4, sc8, sc16, sc32 and sc64 files are supported.")
	end
	# Determine sampling and IF frequency and frequency center
	f_s, f_if, f_center, sig_freq, sigtype = get_signal_type(file_name)
	timestamp, timestamp_JD = find_and_get_timestamp(file_name)
	return FileInfo(f_s, f_if, f_center, sig_freq, sigtype, data_type,
	                timestamp, timestamp_JD, file_name)
end

"""
	signal_to_sc_data_file(signal::GNSSSignal, sig_freq, file_location,
                           file_prefix; operation="replace")
"""
function signal_to_sc_data_file(signal::GNSSSignal, sig_freq, file_location,
	                            file_prefix; operation="replace")
    if sig_freq == L1_freq
		freq_band = "GPSL1"
	elseif sig_freq == L5_freq
		freq_band = "GPSL5"
	else
		freq_band = "$(round(sig_freq*1e-6, digits=8))"
	end
	f_center = "$(round((sig_freq - signal.f_if)*1e-6, digits=8))"  # MHz
	f_s = "$(round(signal.f_s*1e-6, digits=8))"  # MHz
	if signal.nADC == 8
		var_type = Complex{Int8}
	elseif signal.nADC == 16
		var_type = Complex{Int16}
	elseif signal.nADC == 32
		var_type = Complex{Int32}
	elseif  signal.nADC == 64
		var_type = Complex{Int64}
	else
		error("Support for 8-, 16-, 32-, and 64-bit signals only.")
	end
	file_extension = ".sc$(signal.nADC)"
	file_name = string(file_location, file_prefix, "_", freq_band, "_",
	                   f_center, "_", f_s, file_extension)
	if operation == "replace"
		file = open(file_name, "w")
	elseif operation == "append"
		file = open(file_name, "r+")
		seekend(file)
	else
		error("Unknown value for argument `operation`.")
	end
	write(file, var_type.(signal.data))
	close(file)
	return file_name
end