text:
using BinaryBuilder, Pkg

name = "MKL"
version = v"2021.1.1"

sources = [
    ArchiveSource("https://anaconda.org/intel/mkl/2021.1.1/download/linux-64/mkl-2021.1.1-intel_52.tar.bz2",
                  "bfb0fd056576cad99ae1d9c69ada2745420da9f9cf052551d5b91f797538bda2"; unpack_target = "mkl-x86_64-linux-gnu"),
    ArchiveSource("https://anaconda.org/intel/mkl/2021.1.1/download/linux-32/mkl-2021.1.1-intel_52.tar.bz2",
                  "7b6f55a30886154bd96d4b4c6b7428494a59397b87779b58e5b3de00250343f9"; unpack_target = "mkl-i686-linux-gnu"),
    ArchiveSource("https://anaconda.org/intel/mkl/2021.1.1/download/osx-64/mkl-2021.1.1-intel_50.tar.bz2",
                  "819fb8875909d4d024e2a936c54b561aebd1e3aebe58fc605c70aa1ad9a66b70"; unpack_target = "mkl-x86_64-apple-darwin14"),
    ArchiveSource("https://anaconda.org/intel/mkl/2021.1.1/download/win-32/mkl-2021.1.1-intel_52.tar.bz2",
                  "dba6a12a481407ec55fba9895b68afacb15f044905dcb5e185db341b688e6177"; unpack_target = "mkl-i686-w64-mingw32"),
    ArchiveSource("https://anaconda.org/intel/mkl/2021.1.1/download/win-64/mkl-2021.1.1-intel_52.tar.bz2",
                  "4024391b8a45836d5a7ee92405b7767874b3c3bbf2f490349fda042db3b60dfd"; unpack_target = "mkl-x86_64-w64-mingw32"),
]

# Bash recipe for building across all platforms
script = raw"""
cd ${WORKSPACE}/srcdir/mkl-${target}
if [[ ${target} == *-mingw* ]]; then
    cp -r Library/bin/* ${libdir}
else
    cp -r lib/* ${libdir}
fi

install_license info/licenses/*.txt
"""

platforms = [
    Platform("x86_64", "linux"; libc="glibc"),
    Platform("i686", "linux"; libc="glibc"),
    Platform("x86_64", "macos"),
    Platform("i686", "windows"),
    Platform("x86_64", "windows"),
]

# The products that we will ensure are always built
products = [
    LibraryProduct(["libmkl_core", "mkl_core", "mkl_core.1"], :libmkl_core),
    LibraryProduct(["libmkl_rt", "mkl_rt", "mkl_rt.1"], :libmkl_rt),
]

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

non_reg_ARGS = filter(arg -> arg != "--register", ARGS)
include("../../fancy_toys.jl")
no_autofix_platforms = [Platform("i686", "windows"), Platform("x86_64", "windows"), Platform("x86_64", "macos")]
autofix_platforms = [Platform("x86_64", "linux"), Platform("i686", "linux")]
if any(should_build_platform.(triplet.(no_autofix_platforms)))
    # Need to disable autofix: updating linkage of libmkl_intel_thread.dylib on
    # macOS causes runtime issues:
    # https://github.com/JuliaPackaging/Yggdrasil/issues/915.
    build_tarballs(non_reg_ARGS, name, version, sources, script, no_autofix_platforms, products, dependencies; lazy_artifacts = true, autofix = false)
end
if any(should_build_platform.(triplet.(autofix_platforms)))
    # ... but we need to run autofix on Linux, because here libmkl_rt doesn't
    # have a soname, so we can't ccall it without specifying the path:
    # https://github.com/JuliaSparse/Pardiso.jl/issues/69
    build_tarballs(ARGS, name, version, sources, script, autofix_platforms, products, dependencies; lazy_artifacts = true)
end

ext: jl

max_stars_repo_name: waralex/Yggdrasil

max_issues_repo_name: waralex/Yggdrasil

max_forks_repo_name: waralex/Yggdrasil


text:
using FuseApi
using FuseApi.ExampleFs

ExampleFs.fg()


text:
#========================================================================

	Steady state solution for an Aiyagari model with a diffusion process

	Translated julia code based on code from Ben Moll's Website:
		http://www.princeton.edu/moll/HACTproject.htm

		For Julia 1.0.0
========================================================================#
using LinearAlgebra, SparseArrays, Plots

# Model parameters
γ = 2 # CRRA utility
α = 1/3 # parameter for cobb-douglas production function
δ = 0.05 # captal depreciation
ρ = 0.05 # discount rate
A_prod =0.3 # term in production function

# Parameters for productivity and the poisson process q
z_1 = .2
z_2 = 2*z_1
z = [z_1 z_2]
λ_1 = 1 # prob of transitioning from state 1 to 2
λ_2 = 1  # state 2 to 1
λ = [λ_1 λ_2]
z_ave = (z_1*λ_2 + z_2*λ_1)/(λ_1 + λ_2) # Average value of z

H = 500 # Number of elements in grid space for a
a_min = 0.0
a_max = 20.0

# Grid spaces for variables
a = LinRange(a_min, a_max, H)
a = convert(Array,a)
da = (a_max-a_min)/(H-1)

aa = [a a]
zz = ones(H,1)*z #check later

# Paramters for the simulation
maxit = 100
Ir = 100
crit = 10^(-6)
critS = 1e-5
Δ = 1000 #step size for HJB

# Spaces for the finite difference terms
Vaf, Vab, c= [zeros(H,2) for i in 1:3]
V_r = zeros(H,2,H)

# Create r grid space
r_min = -0.0499
r_max = 0.049
r_grid = LinRange(r_min, r_max, Ir)

# Set up an initial guess
r=r_grid[1]
KD_1 = (α*A_prod/(r+δ))^(1/(1-α))*z_ave
w =(1-α)*A_prod*(KD_1/z_ave)^α

v0 = zeros(H,2)
v0[:,1]=(w*z[1] .+ max(r,0.01).*a).^(1-γ)/(1-γ)/ρ
v0[:,2]=(w*z[2] .+ max(r,0.01).*a).^(1-γ)/(1-γ)/ρ

# Spaces for other terms
w_r =[]
KD= zeros(1,Ir)
S=[]
adot=[]
V_n=[]
g_r=[]

# Construct the matrix to summarize the evolution of z over time
# this matrix will not change over the finite differences and is based on the poisson process

Aswitch = [-sparse(I,H,H)*λ[1] sparse(I,H,H)*λ[2] ; sparse(I,H,H)*λ[1] -sparse(I,H,H)*λ[2]]

# The main loop, this makes sure markets clear

for ir in 1:Ir
	global r =r_grid[ir]
	global v0 = v0

	println("Main Loop iteration $(ir)")

	KD[ir]=(α*A_prod/(r+δ))^(1/(1-α))*z_ave
	global w = (1-α)*A_prod*KD[ir].^α * z_ave^(-α)

	push!(w_r,w)

	if ir >1
		global v0 = V_r[:,:,ir-1]
	end

	global v=v0

# Inner loop, value function iteration

	for n in 1:maxit
		V=v
		push!(V_n, V)
		# Forward difference
		Vaf[1:H-1,:] = (V[2:H,:]-V[1:H-1,:])/da
		Vaf[H,:] = (w*z .+ r.*a_max).^(-γ)

		# Backward difference
		Vab[2:H,:] = (V[2:H,:]-V[1:H-1,:])/da
		Vab[1,:] = (w*z .+ r.*a_min).^(-γ)

		# caculate consumption and savings
		# First the forward difference case
		cf = Vaf.^(-1/γ)
		sf = w*zz + r.*aa -cf

		# consumption and savings for the backward difference
		cb = Vab.^(-1/γ)
		sb = w*zz + r.*aa -cb

		#Consumption and V' at the steady state
		c0 = w*zz + r.*aa
		Va0 = c0.^(-γ)

		# Setup to the Upwind scheme
		If = sf .> 0
		Ib = sb .< 0
		I0 = (1 .- If .- Ib)

		#Now implement the upwind scheme
		Va_Upwind = Vaf.*If + Vab.*Ib + Va0.*I0

		global c= Va_Upwind.^(-1/γ)
		u = c.^(1-γ)/(1-γ)

		X = -min.(sb,0)/da
		Y = -max.(sf,0)/da + min.(sb,0)/da
		Z = max.(sf,0)/da

		A1 = spdiagm(0=>Y[:,1], -1 => X[2:H,1], 1=> Z[1:H-1,1])
		A2 = spdiagm(0=>Y[:,2], -1 => X[2:H,2], 1=> Z[1:H-1,2])
		AA = [A1 spzeros(H,H); spzeros(H,H) A2]

		global A= AA + Aswitch

		B = (1/Δ + ρ)*sparse(I,2*H,2*H)-A

		u_stacked = [u[:,1]; u[:,2]]
		V_stacked = [V[:,1]; V[:,2]]

		b = u_stacked + V_stacked/Δ

		V_stacked = B\b

		V = [V_stacked[1:H] V_stacked[H+1:2*H]]

		V_change = V-v

		global v = V

		dist = findmax(V_change)[1]

		if n>1 && dist < crit
			println("Value Function Converge Iteration = $(n)")
			break
		end
	end
	# Now solve the Fokker-Planck equation
	AT = A'
	b = zeros(2*H,1)

	# Dirty fix to avoid matrix singularities
	i_fix=1
	b[i_fix] =.1
	row = [zeros(1,i_fix-1) 1 zeros(1,H*2-i_fix)]
	AT[i_fix,:] = row

	# Now solve the linear system
	gg = AT\b
	g_sum = gg'*ones(2*H,1)*da
	gg = gg./g_sum

	global g = [gg[1:H] gg[H+1:2*H]]

	push!(g_r, g)
	push!(adot, w*zz+r.*aa-c)
	V_r[:,:,ir] = v

	push!(S, g[:,1]'*a*da + g[:,2]'*a*da)
end

# Graphs
S_max = findmax(S)[1]
amin1=a_min-0.02
aaa=LinRange(amin1,S_max,Ir)
rrr=LinRange(-0.06,0.06,Ir)
KD_new = (α*A_prod./max.(rrr.+δ,0)).^(1/(1-α)).*z_ave

plot(S,r_grid, xlims=(amin1,0.6),ylims=(-0.06,0.06), color=:blue,
	label="\$S(r)\$", legend=:bottomright, xlabel="K",ylabel="r",grid=false)
plot!(KD_new,rrr, label="\$F_{k}(K,1)-\\delta\$")
plot!(zeros(Ir,1).+a_min,rrr,label="\$a=\\underline{a}\$",color=:yellow,line=:dash)
plot!(aaa,ones(Ir,1).*ρ,label="\$r=\\rho\$",color=:purple,line=:dash)
plot!(aaa,ones(Ir,1).*(-δ),label="\$r=-\\delta\$",color=:green, line=:dash)
png("Asset_supply")


text:
using BitBasis
using VAN:get_logp, unpack_gradient

export free_energy, free_energy_local, energy, entropy

function get_input_reg(nbits::Int, samples::AbstractArray)
    get_input_reg(nbits, Vector(packbits(samples)))
end

function get_input_reg(nbits::Int, samples::Vector{<:Integer})
    nbatch = length(samples)
    config = zeros(ComplexF64, nbatch, 1<<nbits)
    for i = 1:nbatch
        config[i,Int(samples[i])+1] = 1
    end
    # transpose the register for better performance
    return ArrayReg(transpose(config))
end

function free_energy(β::Real, H::AbstractBlock{N}, sampler::AbstractSampler, circuit::AbstractBlock{N}, samples) where N
    mean(free_energy_local(β, H, sampler, circuit, samples))
end

function energy(H::AbstractBlock{N}, sampler, circuit::AbstractBlock{N}, samples) where N
    reg = get_input_reg(N, samples)
    mean(real.(expect(H, reg|>circuit)))
end

function energy2(H::AbstractBlock{N}, sampler::AbstractSampler, circuit::AbstractBlock{N}, samples) where N
    reg = get_input_reg(N, samples)
    mean(real.(expect(H*H, reg|>circuit)))
end

function entropy(H::AbstractBlock{N}, sampler::AbstractSampler, circuit::AbstractBlock{N}, samples) where N
    -mean(get_logp(sampler, samples))
end

function purity(H::AbstractBlock{N}, sampler::AbstractSampler, circuit::AbstractBlock{N}, samples) where N
    mean(exp.(get_logp(sampler, samples)))
end

function spectra(H::AbstractBlock{N}, sampler::AbstractSampler, circuit::AbstractBlock{N}, samples) where N
    reg = get_input_reg(N, samples)
    real.(expect(H, reg|>circuit))
end

function free_energy_local(β::Real, H::AbstractBlock{N}, sampler, circuit::AbstractBlock{N}, samples) where N
   reg = get_input_reg(N, samples)
   return transpose(real.(expect(H, reg|>circuit))) + get_logp(sampler, samples)/β
end

# E(x) = <x|U'HU|x>
# d([E(x)]{x∼p} - TS)/dθ = [d(E(x))/dθ]{x∼p}
function grad_θ(H, circuit::AbstractBlock{N}, samples) where N
    reg = get_input_reg(N, samples)
    _, paramsδ = expect'(H, reg=>circuit)
    return paramsδ /size(samples)[end] # divide by batch size since there is a mean in loss
end

function loss_reinforce(β, H, sampler, circuit::AbstractBlock{N}, samples) where N
    f = free_energy_local(β, H, sampler, circuit, samples)
    logp = get_logp(sampler, samples)
    b = mean(f)
    return mean(logp .* (f .- b))
end

function grad_sampler(β, H, sampler, circuit::AbstractBlock{N}, samples) where N
    g = gradient(loss_reinforce, β, H, sampler, circuit, samples)[3]
    unpack_gradient(sampler, g)
end


text:
module KitMSA

export	remove_gapped_cols,
		remove_gapped_sequences,
		remove_close_seqs,
		compute_energy_single_sequence,
		extract_params
		fasta2matrix,
		num2lettr,
		lettr2num,
		vec2string,
		string2vec

using DelimitedFiles
using FastaIO
using StatsBase
using Random

include("clean_MSA.jl")
include("read_write_MSA.jl")
include("MSA_analysis.jl")
include("BM_utils")

end

max_stars_repo_name: sabrinacotogno/KitMSA

max_issues_repo_name: sabrinacotogno/KitMSA

max_forks_repo_name: sabrinacotogno/KitMSA


text:
__precompile__()

module DevTools

using Blink, Media, Lazy, Requires

export BlinkDisplay, pin, top, docs, profiler

include("display/BlinkDisplay.jl")
using .BlinkDisplay
include("profile/profile.jl")
include("codemirror.jl")
include("collab.jl")

profiler() = ProfileView.fetch()

end # module


text:
int main() {
    int x;

    x = 5;
    while (x > 0) {
        printInt(x);
        x--;
    }

    printInt(x);

    return 0;
}

max_stars_repo_name: fredefox/javalette

max_issues_repo_name: fredefox/javalette

max_forks_repo_name: fredefox/javalette


text:
### math.jl
#
# Copyright (C) 2016 Mosè Giordano.
#
# Maintainer: Mosè Giordano <mose AT gnu DOT org>
# Keywords: uncertainty, error propagation, physics
#
# This file is a part of Measurements.jl.
#
# License is MIT "Expat".
#
### Commentary:
#
# This file contains definition of mathematical functions that support
# Measurement objects.
#
# Note: some functions defined here (like all degree-related and reciprocal
# trigonometric functions, fld, cld, hypot, cbrt, abs, mod) are redundant in the
# sense that you would get the correct result also without their definitions,
# but having them defined here avoids some calculations and slightly improves
# performance.  Likewise, multiple methods are provided for functions taking two
# (or more) arguments because when only one argument is of Measurement type we
# can use the simple `result' function for one derivative that is faster than
# the generic method.
#
### Code:

export @uncertain

# This function is to be used by methods of mathematical operations to produce a
# `Measurement' object in output.  Arguments are:
#   * val: the nominal result of operation G(a)
#   * der: the derivative ∂G/∂a of G with respect to the variable a
#   * a: the only argument of G
# In this simple case of unary function, we don't have the problem of correlated
# variables (thus making this method much faster than the next one), so we can
# calculate the uncertainty of G(a) as
#   σ_G = |σ_a·∂G/∂a|
# The list of derivatives with respect to each measurement is updated with
#   ∂G/∂a · previous_derivatives
function result{T<:AbstractFloat}(val::Real, der::Real, a::Measurement{T})
    val, der = promote(val, der)
    newder = similar(a.der)
    @inbounds for tag in keys(a.der)
        if tag[2] != 0.0 # Skip values with 0 uncertainty
            newder = Derivatives(newder, tag=>der*a.der[tag])
        end
    end
    # If uncertainty of "a" is null, the uncertainty of result is null as well,
    # even if the derivative is NaN or infinite.  In any other case, use
    # σ_G = |σ_a·∂G/∂a|.
    σ = (a.err == 0.0) ? 0.0 : abs(der*a.err)
    # The tag is NaN because we don't care about tags of derived quantities, we
    # are only interested in independent ones.
    Measurement(val,  σ, NaN, newder)
end

# This function is similar to the previous one, but applies to mathematical
# operations with more than one argument, so the formula to propagate
# uncertainty is more complicated because we have to take into account
# correlation between arguments.  The arguments are the same as above, but `der'
# and `a' are tuples of the same length (`der' has the derivatives of G with
# respect to the corresponding variable in `a').
#
# Suppose we have a function G = G(a1, a2) of two arguments.  a1 and a2 are
# correlated, because they come from some mathematical operations on really
# independent variables x, y, z, say a1 = a1(x, y), a2 = a2(x, z).  The
# uncertainty on G(a1, a2) is calculated as follows:
#   σ_G = sqrt((σ_x·∂G/∂x)^2 + (σ_y·∂G/∂y)^2 + (σ_z·∂G/∂z)^2)
# where ∂G/∂x is the partial derivative of G with respect to x, and so on.  We
# can expand the previous formula to:
#   σ_G = sqrt((σ_x·(∂G/∂a1·∂a1/∂x + ∂G/∂a2·∂a2/∂x))^2 + (σ_y·∂G/∂a1·∂a1/∂y)^2 +
#               + (σ_z·∂G/∂a2·∂a2/∂z)^2)
function result(val::Real, der::Tuple{Vararg{Real}},
                a::Tuple{Vararg{Measurement}})
    @assert length(der) == length(a)
    a = promote(a...)
    T = typeof(a[1].val)
    newder = similar(a[1].der)
    err::T = zero(T)
    # Iterate over all independent variables.  We first iterate over all
    # variables listed in `a' in order to get all independent variables upon
    # which those variables depend, then we get the `tag' of each independent
    # variable, skipping variables that have been already taken into account.
    @inbounds for y in a
        for tag in keys(y.der)
            if tag ∉ keys(newder) # Skip independent variables already considered
                σ_x = tag[2]
                if σ_x != 0.0 # Skip values with 0 uncertainty
                    ∂G_∂x::T = 0.0
                    # Iteratate over all the arguments of the function
                    for (i, x) in enumerate(a)
                        # Calculate the derivative of G with respect to the
                        # current independent variable.  In the case of the x
                        # independent variable of the example above, we should
                        # get   ∂G/∂x = ∂G/∂a1·∂a1/∂x + ∂G/∂a2·∂a2/∂x
                        ∂a_∂x = derivative(x, tag) # ∂a_i/∂x
                        if ∂a_∂x != 0.0 # Skip values with 0 partial derivative
                            # der[i] = ∂G/∂a_i
                            ∂G_∂x = ∂G_∂x + der[i]*∂a_∂x
                        end
                    end
                    newder = Derivatives(newder, tag=>∂G_∂x)
                    # Add (σ_x·∂G/∂x)^2 to the total uncertainty (squared)
                    err = err + abs2(σ_x*∂G_∂x)
                end
            end
        end
    end
    return Measurement(T(val), sqrt(err), NaN, newder)
end

# "result" function for complex-valued functions (like "besselh").  This takes
# the same argument as the first implementation of "result", but with complex
# "val" and "der".
function result(val::Complex, der::Complex, a::Measurement)
    return complex(result(real(val), real(der), a), result(imag(val), imag(der), a))
end

### @uncertain macro.
"""
    @uncertain f(value ± stddev, ...)

A macro to calculate \$f(value) ± uncertainty\$, with \$uncertainty\$ derived
from \$stddev\$ according to rules of linear error propagation theory.

Function \$f\$ can accept any number of real arguments, the type of the
arguments provided must be `Measurement`.
"""
macro uncertain(expr::Expr)
    f = esc(expr.args[1]) # Function name
    n = length(expr.args) - 1
    if n == 1
        a = esc(expr.args[2]) # Argument, of Measurement type
        return :( result($f($a.val), Calculus.derivative($f, $a.val), $a) )
    else
        a = expr.args[2:end] # Arguments, as an array of expressions
        args = :([])  # Build up array of arguments
        [push!(args.args, :($(esc(a[i])))) for i=1:n] # Fill the array
        argsval =:([])  # Build up the array of values of arguments
        [push!(argsval.args, :($(args.args[i]).val)) for i=1:n] # Fill the array
        return :( result($f($argsval...),
                         (Calculus.gradient(x -> $f(x...), $argsval)...),
                         ($args...)) )
    end
end

### Elementary arithmetic operations:
import Base: +, -, *, /, div, inv, fld, cld

# Addition: +
+(a::Measurement) = a
+(a::Measurement, b::Measurement) = result(a.val + b.val, (1.0, 1.0), (a, b))
+(a::Real, b::Measurement) = result(a + b.val, 1.0, b)
+(a::Measurement, b::Bool) = result(a.val +b, 1.0, a)
+(a::Measurement, b::Real) = result(a.val + b, 1.0, a)

# Subtraction: -
-(a::Measurement) = result(-a.val, -1.0, a)
-(a::Measurement, b::Measurement) = result(a.val - b.val, (1.0, -1.0), (a, b))
-(a::Real, b::Measurement) = result(a - b.val, -1.0, b)
-(a::Measurement, b::Real) = result(a.val - b, 1.0, a)

# Multiplication: *
function *(a::Measurement, b::Measurement)
    aval = a.val
    bval = b.val
    return result(aval*bval, (bval, aval), (a, b))
end
*(a::Bool, b::Measurement) = result(a*b.val, a, b)
*(a::Real, b::Measurement) = result(a*b.val, a, b)
*(a::Measurement, b::Bool) = result(a.val*b, b, a)
*(a::Measurement, b::Real) = result(a.val*b, b, a)

# muladd and fma
import Base: muladd, fma

for f in (:fma, :muladd)
    @eval begin
        # All three arguments are Measurement
        function ($f)(a::Measurement, b::Measurement, c::Measurement)
            x = a.val
            y = b.val
            z = c.val
            return result(($f)(x, y, z), (y, x, one(z)), (a, b, c))
        end

        # First argument is always Measurement
        function ($f)(a::Measurement, b::Measurement, c::Real)
            x = a.val
            y = b.val
            return result(($f)(x, y, c), (y, x), (a, b))
        end

        function ($f)(a::Measurement, b::Real, c::Measurement)
            x = a.val
            z = c.val
            return result(($f)(x, b, z), (b, one(z)), (a, c))
        end

        ($f)(a::Measurement, b::Real, c::Real) =
            result(($f)(a.val, b, c), b, a)

        # Secon argument is always Measurement
        function ($f)(a::Real, b::Measurement, c::Measurement)
            y = b.val
            z = c.val
            return result(($f)(a, y, z), (a, one(z)), (b, c))
        end

        ($f)(a::Real, b::Measurement, c::Real) =
            result(($f)(a, b.val, c), a, b)

        # Third argument is Measurement
        function ($f)(a::Real, b::Real, c::Measurement)
            z = c.val
            return result(($f)(a, b, z), one(z), c)
        end
    end
end

# Division: /, div, fld, cld
function /(a::Measurement, b::Measurement)
    aval = a.val
    oneoverbval = inv(b.val)
    return result(aval*oneoverbval, (oneoverbval, -aval*abs2(oneoverbval)),
                  (a, b))
end
/(a::Real, b::Measurement) = result(a/b.val, -a/abs2(b.val), b)
/(a::Measurement, b::Real) = result(a.val/b, 1/b, a)

# 0.0 as partial derivative for both arguments of "div", "fld", "cld" should be
# correct for most cases.  This has been tested against "@uncertain" macro.
div(a::Measurement, b::Measurement) = result(div(a.val, b.val), (0.0, 0.0), (a, b))
div(a::Measurement, b::Real) = result(div(a.val, b), 0.0, a)
div(a::Real, b::Measurement) = result(div(a, b.val), 0.0, b)

fld(a::Measurement, b::Measurement) = result(fld(a.val, b.val), (0.0, 0.0), (a, b))
fld(a::Measurement, b::Real) = result(fld(a.val, b), 0.0, a)
fld(a::Real, b::Measurement) = result(fld(a, b.val), 0.0, b)

cld(a::Measurement, b::Measurement) = result(cld(a.val, b.val), (0.0, 0.0), (a, b))
cld(a::Measurement, b::Real) = result(cld(a.val, b), 0.0, a)
cld(a::Real, b::Measurement) = result(cld(a, b.val), 0.0, b)

# Inverse: inv
function inv(a::Measurement)
    inverse = inv(a.val)
    return result(inverse, -abs2(inverse), a)
end

# signbit
import Base: signbit

signbit(a::Measurement) = signbit(a.val)

# Power: ^
import Base: ^, exp2

function ^(a::Measurement, b::Measurement)
    aval = a.val
    bval = b.val
    pow = aval^bval
    return result(pow, (aval^(bval - 1.0)*bval, pow*log(aval)), (a, b))
end

function ^{T<:Integer}(a::Measurement, b::T)
    aval = a.val
    return result(aval^b, aval^(b-1)*b, a)
end

function ^{T<:Rational}(a::Measurement,  b::T)
    if isinteger(b)
        return a^trunc(Integer, b)
    else
        aval = a.val
        return result(aval^b, b*aval^(b - 1.0), a)
    end
end

function ^{T<:Real}(a::Measurement,  b::T)
    if isinteger(float(b))
        return a^trunc(Integer, b)
    else
        aval = a.val
        return result(aval^b, b*aval^(b - 1.0), a)
    end
end

^(::Irrational{:e}, b::Measurement) = exp(b)

function ^{T<:Real}(a::T,  b::Measurement)
    res = a^b.val
    return result(res, res*log(a), b)
end

function exp2{T<:AbstractFloat}(a::Measurement{T})
    pow = exp2(a.val)
    return result(pow, pow*log(T(2)), a)
end

### Trigonometric functions

# deg2rad, rad2deg
import Base: deg2rad, rad2deg

deg2rad(a::Measurement) = result(deg2rad(a.val), oftype(a.val, pi)/180, a)
rad2deg(a::Measurement) = result(rad2deg(a.val), 180/oftype(a.val, pi), a)

# Cosine: cos, cosd, cosh
import Base: cos, cosd, cosh

function cos(a::Measurement)
    aval = a.val
    result(cos(aval), -sin(aval), a)
end

function cosd(a::Measurement)
    aval = a.val
    return result(cosd(aval), -deg2rad(sind(aval)), a)
end

function cosh(a::Measurement)
    aval = a.val
    result(cosh(aval), sinh(aval), a)
end

# Sine: sin, sind, sinh
import Base: sin, sind, sinh

function sin(a::Measurement)
    aval = a.val
    result(sin(aval), cos(aval), a)
end

function sind(a::Measurement)
    aval = a.val
    return result(sind(aval), deg2rad(cosd(aval)), a)
end

function sinh(a::Measurement)
    aval = a.val
    result(sinh(aval), cosh(aval), a)
end

# Tangent: tan, tand, tanh
import Base: tan, tand, tanh

function tan(a::Measurement)
    aval = a.val
    return result(tan(aval), abs2(sec(aval)), a)
end

function tand(a::Measurement)
    aval = a.val
    return result(tand(aval), deg2rad(abs2(secd(aval))), a)
end

function tanh(a::Measurement)
    aval = a.val
    return result(tanh(aval), abs2(sech(aval)), a)
end

# Inverse trig functions: acos, acosd, acosh, asin, asind, asinh, atan, atand,
# atan2, atanh
import Base: acos, acosd, acosh, asin, asind, asinh, atan, atand, atan2, atanh

function acos(a::Measurement)
    aval = a.val
    return result(acos(aval), -inv(sqrt(1.0 - abs2(aval))), a)
end

function acosd(a::Measurement)
    aval = a.val
    return result(acosd(aval), -rad2deg(inv(sqrt(1.0 - abs2(aval)))), a)
end

function acosh(a::Measurement)
    aval = a.val
    return result(acosh(aval), inv(sqrt(abs2(aval) - 1.0)), a)
end

function asin(a::Measurement)
    aval = a.val
    return result(asin(aval), inv(sqrt(1.0 - abs2(aval))), a)
end

function asind(a::Measurement)
    aval = a.val
    return result(asind(aval), rad2deg(inv(sqrt(1.0 - abs2(aval)))), a)
end

function asinh(a::Measurement)
    aval = a.val
    return result(asinh(aval), inv(hypot(aval, 1.0)), a)
end

function atan(a::Measurement)
    aval = a.val
    return result(atan(aval), inv(abs2(aval) + 1.0), a)
end

function atand(a::Measurement)
    aval = a.val
    return result(atand(aval), rad2deg(inv(abs2(aval) + 1.0)), a)
end

function atanh(a::Measurement)
    aval = a.val
    return result(atanh(aval), inv(1.0 - abs2(aval)), a)
end

function atan2(a::Measurement, b::Measurement)
    aval = a.val
    bval = b.val
    invdenom = inv(abs2(aval) + abs2(bval))
    return result(atan2(aval, bval),
                  (bval*invdenom, -aval*invdenom),
                  (a, b))
end

function atan2(a::Measurement, b::Real)
    x = a.val
    return result(atan2(x, b), -b/(abs2(x) + abs2(b)), a)
end

function atan2(a::Real, b::Measurement)
    y = b.val
    return result(atan2(a, y), -a/(abs2(a) + abs2(y)), b)
end

# Reciprocal trig functions: csc, cscd, csch, sec, secd, sech, cot, cotd, coth
import Base: csc, cscd, csch, sec, secd, sech, cot, cotd, coth

function csc(a::Measurement)
    aval = a.val
    val = csc(aval)
    return result(val, -val*cot(aval), a)
end

function cscd(a::Measurement)
    aval = a.val
    val = cscd(aval)
    return result(val, -deg2rad(val*cotd(aval)), a)
end

function csch(a::Measurement)
    aval = a.val
    val = csch(aval)
    return result(val, -val*coth(aval), a)
end

function sec(a::Measurement)
    aval = a.val
    val = sec(aval)
    return result(val, val*tan(aval), a)
end

function secd(a::Measurement)
    aval = a.val
    val = secd(aval)
    return result(val, deg2rad(val*tand(aval)), a)
end

function sech(a::Measurement)
    aval = a.val
    val = sech(aval)
    return result(val, val*tanh(aval), a)
end

function cot(a::Measurement)
    aval = a.val
    return result(cot(aval), -abs2(csc(aval)), a)
end

function cotd(a::Measurement)
    aval = a.val
    return result(cotd(aval), -deg2rad(abs2(cscd(aval))), a)
end

function coth(a::Measurement)
    aval = a.val
    return result(coth(aval), -abs2(csch(aval)), a)
end

### Exponential-related

# Exponentials: exp, expm1, exp10, frexp, ldexp
import Base: exp, expm1, exp10, frexp, ldexp

function exp(a::Measurement)
    val = exp(a.val)
    return result(val, val, a)
end

function expm1(a::Measurement)
    aval = a.val
    return result(expm1(aval), exp(aval), a)
end

function exp10{T<:AbstractFloat}(a::Measurement{T})
    val = exp10(a.val)
    return result(val, log(T(10))*val, a)
end

function frexp(a::Measurement)
    x, y = frexp(a.val)
    return (result(x, inv(exp2(y)), a), y)
end

ldexp(a::Measurement, e::Integer) = result(ldexp(a.val, e), ldexp(1.0, e), a)

# Logarithms
import Base: log, log2, log10, log1p

function log(a::Measurement, b::Measurement)
    aval = a.val
    bval = b.val
    val = log(aval, bval)
    loga = log(aval)
    return result(val, (-val/(aval*loga), inv(loga*bval)), (a, b))
end

function log(a::Measurement) # Special case
    aval = a.val
    return result(log(aval), inv(aval), a)
end

function log2{T<:AbstractFloat}(a::Measurement{T}) # Special case
    x = a.val
    return result(log2(x), inv(log(T(2))*x), a)
end

function log10{T<:AbstractFloat}(a::Measurement{T}) # Special case
    aval = a.val
    return result(log10(aval), inv(log(T(10))*aval), a)
end

function log1p(a::Measurement) # Special case
    aval = a.val
    return result(log1p(aval), inv(aval + one(aval)), a)
end

log(::Irrational{:e}, a::Measurement) = log(a)

function log(a::Real, b::Measurement)
    bval = b.val
    return result(log(a, bval), inv(log(a)*bval), b)
end

function log(a::Measurement, b::Real)
    aval = a.val
    res = log(aval, b)
    return result(res, -res/(aval*log(aval)), a)
end

# Hypotenuse: hypot
import Base: hypot

function hypot(a::Measurement, b::Measurement)
    aval = a.val
    bval = b.val
    val = hypot(aval, bval)
    invval = inv(val)
    return result(val,
                  (aval*invval, bval*invval),
                  (a, b))
end

function hypot(a::Real, b::Measurement)
    bval = b.val
    res = hypot(a, bval)
    return result(res, bval*inv(res), b)
end

function hypot(a::Measurement, b::Real)
    aval = a.val
    res = hypot(aval, b)
    return result(res, aval*inv(res), a)
end

# Square root: sqrt
import Base: sqrt

function sqrt(a::Measurement)
    val = sqrt(a.val)
    return result(val, 0.5*inv(val), a)
end

# Cube root: cbrt
import Base: cbrt

function cbrt(a::Measurement)
    aval = a.val
    val = cbrt(aval)
    return result(val, val*inv(3.0*aval), a)
end

### Absolute value, sign and the likes

# Absolute value
import Base: abs, abs2

function abs(a::Measurement)
    aval = a.val
    return result(abs(aval), copysign(1, aval), a)
end

function abs2(a::Measurement)
    x = a.val
    return result(abs2(x), 2*x, a)
end

# Sign: sign, copysign, flipsign
import Base: sign, copysign, flipsign

sign(a::Measurement) = result(sign(a.val), 0.0, a)

function copysign(a::Measurement, b::Measurement)
    aval = a.val
    bval = b.val
    result(copysign(aval, bval),
           (copysign(1, aval)/copysign(1, bval), 0.0),
           (a, b))
end

copysign(a::Measurement, b::Real) = copysign(a, measurement(b))
copysign(a::Signed, b::Measurement) = copysign(measurement(a), b)
copysign(a::Rational, b::Measurement) = copysign(measurement(a), b)
copysign(a::Float32, b::Measurement) = copysign(measurement(a), b)
copysign(a::Float64, b::Measurement) = copysign(measurement(a), b)
copysign(a::Real, b::Measurement) = copysign(measurement(a), b)

function flipsign(a::Measurement, b::Measurement)
    flip = flipsign(a.val, b.val)
    return result(flip, (copysign(1.0, flip), 0.0), (a, b))
end

flipsign(a::Measurement, b::Real) = flipsign(a, measurement(b))
flipsign(a::Signed, b::Measurement) = flipsign(measurement(a), b)
flipsign(a::Float32, b::Measurement) = flipsign(measurement(a), b)
flipsign(a::Float64, b::Measurement) = flipsign(measurement(a), b)
flipsign(a::Real, b::Measurement) = flipsign(measurement(a), b)

### Special functions

# Error function: erf, erfinv, erfc, erfcinv, erfcx, erfi, dawson
import Base: erf, erfinv, erfc, erfcinv, erfcx, erfi, dawson

function erf{T<:AbstractFloat}(a::Measurement{T})
    aval = a.val
    return result(erf(aval), 2*exp(-abs2(aval))/sqrt(T(pi)), a)
end

function erfinv{T<:AbstractFloat}(a::Measurement{T})
    res = erfinv(a.val)
    # For the derivative, see http://mathworld.wolfram.com/InverseErf.html
    return result(res, 0.5*sqrt(T(pi))*exp(abs2(res)), a)
end

function erfc{T<:AbstractFloat}(a::Measurement{T})
    aval = a.val
    return result(erfc(aval), -2*exp(-abs2(aval))/sqrt(T(pi)), a)
end

function erfcinv{T<:AbstractFloat}(a::Measurement{T})
    res = erfcinv(a.val)
    # For the derivative, see http://mathworld.wolfram.com/InverseErfc.html
    return result(res, -0.5*sqrt(T(pi))*exp(abs2(res)), a)
end

function erfcx{T<:AbstractFloat}(a::Measurement{T})
    aval = a.val
    res = erfcx(aval)
    return result(res, 2*(aval*res - inv(sqrt(T(pi)))), a)
end

function erfi{T<:AbstractFloat}(a::Measurement{T})
    aval = a.val
    return result(erfi(aval), 2*exp(abs2(aval))/sqrt(T(pi)), a)
end

function dawson{T<:AbstractFloat}(a::Measurement{T})
    aval = a.val
    res = dawson(aval)
    return result(res, 1.0 - 2.0*aval*res, a)
end

# Factorial and gamma
import Base: factorial, gamma, lgamma, digamma, invdigamma, trigamma, polygamma

function factorial(a::Measurement)
    aval = a.val
    fact = factorial(aval)
    return result(fact, fact*digamma(aval + one(aval)), a)
end

function gamma(a::Measurement)
    aval = a.val
    Γ = gamma(aval)
    return result(Γ, Γ*digamma(aval), a)
end

function lgamma(a::Measurement)
    aval = a.val
    return result(lgamma(aval), digamma(aval), a)
end

function digamma(a::Measurement)
    aval = a.val
    return result(digamma(aval), trigamma(aval), a)
end

function invdigamma(a::Measurement)
    aval = a.val
    res = invdigamma(aval)
    return result(res, inv(trigamma(res)), a)
end

function trigamma(a::Measurement)
    aval = a.val
    return result(trigamma(aval), polygamma(2, aval), a)
end

function polygamma(n::Integer, a::Measurement)
    aval = a.val
    return result(polygamma(n, aval), polygamma(n + 1, aval), a)
end

# Beta function: beta, lbeta
import Base: beta, lbeta

function beta(a::Measurement, b::Measurement)
    aval = a.val
    bval = b.val
    res = beta(aval, bval)
    return result(res,
                  (res*(digamma(aval) - digamma(aval + bval)),
                   res*(digamma(bval) - digamma(aval + bval))),
                  (a, b))
end

function beta(a::Measurement, b::Real)
    aval = a.val
    res = beta(aval, b)
    return result(res, res*(digamma(aval) - digamma(aval + b)), a)
end

beta(a::Real, b::Measurement) = beta(b, a)

function lbeta(a::Measurement, b::Measurement)
    aval = a.val
    bval = b.val
    return result(lbeta(aval, bval),
                  (digamma(aval) - digamma(aval + bval),
                   digamma(bval) - digamma(aval + bval)),
                  (a, b))
end

function lbeta(a::Measurement, b::Real)
    aval = a.val
    return result(lbeta(aval, b), digamma(aval) - digamma(aval + b), a)
end

lbeta(a::Real, b::Measurement) = lbeta(b, a)

# Airy functions
import Base: airy

function airy(k::Integer, a::Measurement)
    aval = a.val
    if k == 0 || k == 2
        return result(airy(k, aval), airy(k + 1, aval), a)
    else
        # Use Airy equation: y'' - xy = 0 => y'' = xy
        return result(airy(k, aval), aval*airy(k - 1, aval), a)
    end
end

# Bessel functions
import Base: besselj0, besselj1, besselj, bessely0, bessely1, bessely, besselh,
besseli, besselix, besselk, besselkx

function besselj0(a::Measurement)
    x = a.val
    return result(besselj0(x), -besselj1(x), a)
end

function besselj1(a::Measurement)
    x = a.val
    return result(besselj1(x), 0.5*(besselj0(x) - besselj(2, x)), a)
end

# XXX: this is necessary to fix a method ambiguity in Julia 0.4.  Remove this
# definition when that version will not be supported anymore
function besselj(nu::Integer, a::Measurement)
    x = a.val
    return result(besselj(nu, x), 0.5*(besselj(nu - 1, x) - besselj(nu + 1, x)), a)
end

# XXX: this is necessary to fix a method ambiguity in Julia 0.4.  Remove this
# definition when that version will not be supported anymore
function besselj(nu::AbstractFloat, a::Measurement)
    x = a.val
    return result(besselj(nu, x), 0.5*(besselj(nu - 1, x) - besselj(nu + 1, x)), a)
end

# XXX: I don't know a closed form expression for the derivative with respect to
# first argument of J_n.  Arguably, there will be more cases where the
# measurement is the second argument, than the first one.  In any case, you can
# use "@uncertain" macro when both arguments are of Measurement type.
function besselj(nu::Real, a::Measurement)
    x = a.val
    return result(besselj(nu, x), 0.5*(besselj(nu - 1, x) - besselj(nu + 1, x)), a)
end

function bessely0(a::Measurement)
    x = a.val
    return result(bessely0(x), -bessely1(x), a)
end

function bessely1(a::Measurement)
    x = a.val
    return result(bessely1(x), 0.5*(bessely0(x) - bessely(2, x)), a)
end

# XXX: this is necessary to fix a method ambiguity in Julia 0.4.  Remove this
# definition when that version will not be supported anymore
function bessely(nu::Integer, a::Measurement)
    x = a.val
    return result(bessely(nu, x), 0.5*(bessely(nu - 1, x) - bessely(nu + 1, x)), a)
end

# XXX: I don't know a closed form expression for the derivative with respect to
# first argument of y_n, see comments about "besselj".
function bessely(nu::Real, a::Measurement)
    x = a.val
    return result(bessely(nu, x), 0.5*(bessely(nu - 1, x) - bessely(nu + 1, x)), a)
end

function besselh(nu::Real, k::Integer, a::Measurement)
    x = a.val
    return result(besselh(nu, k, x),
                  0.5*(besselh(nu - 1, k, x) - besselh(nu + 1, k, x)),
                  a)
end

function besseli(nu::Real, a::Measurement)
    x = a.val
    return result(besseli(nu, x), 0.5*(besseli(nu - 1, x) + besseli(nu + 1, x)), a)
end

function besselix(nu::Real, a::Measurement)
    x = a.val
    return result(besselix(nu, x),
                  0.5*(besseli(nu - 1, x) + besseli(nu + 1, x))*exp(-abs(x)) -
                  besseli(nu, x)*sign(x)*exp(-abs(x)),
                  a)
end

function besselk(nu::Real, a::Measurement)
    x = a.val
    return result(besselk(nu, x), -0.5*(besselk(nu - 1, x) + besselk(nu + 1, x)), a)
end

function besselkx(nu::Real, a::Measurement)
    x = a.val
    return result(besselkx(nu, x),
                  -0.5*(besselk(nu - 1, x) + besselk(nu + 1, x))*exp(x) +
                  besselk(nu, x)*exp(x),
                  a)
end

### Modulo

import Base: mod, rem, mod2pi

# Use definition of "mod" function:
# http://docs.julialang.org/en/stable/manual/mathematical-operations/#division-functions
mod(a::Measurement, b::Measurement) = a - fld(a, b)*b
mod(a::Measurement, b::Real) = mod(a, measurement(b))
mod(a::Real, b::Measurement) = mod(measurement(a), b)

# Use definition of "rem" function:
# http://docs.julialang.org/en/stable/manual/mathematical-operations/#division-functions
rem(a::Measurement, b::Measurement) = a - div(a, b)*b
rem(a::Measurement, b::Real) = rem(a, measurement(b))
rem(a::Real, b::Measurement) = rem(measurement(a), b)

mod2pi(a::Measurement) = result(mod2pi(a.val), 1, a)

### Machine precision

import Base: eps, nextfloat, maxintfloat, typemax

eps{T<:AbstractFloat}(::Type{Measurement{T}}) = eps(T)
eps{T<:AbstractFloat}(a::Measurement{T}) = eps(a.val)

nextfloat(a::Measurement) = nextfloat(a.val)

maxintfloat{T<:AbstractFloat}(::Type{Measurement{T}}) = maxintfloat(T)

typemax{T<:AbstractFloat}(::Type{Measurement{T}}) = typemax(T)

### Rounding
import Base: round, floor, ceil, trunc

round(a::Measurement) = round(a.val)
round{T<:Integer}(::Type{T}, a::Measurement) = round(T, a.val)
floor(a::Measurement) = floor(a.val)
floor{T<:Integer}(::Type{T}, a::Measurement) = floor(T, a.val)
ceil(a::Measurement) = ceil(a.val)
ceil{T<:Integer}(::Type{T}, a::Measurement) = ceil(Integer, a.val)
trunc(a::Measurement) = trunc(a.val)
trunc{T<:Integer}(::Type{T}, a::Measurement) = trunc(T, a.val)

# Widening
import Base: widen

widen{T<:AbstractFloat}(::Type{Measurement{T}}) = Measurement{widen(T)}

# To big float
import Base: big

big{T<:AbstractFloat}(x::Measurement{T}) = convert(Measurement{BigFloat}, x)

# Sum and prod
import Base: sum, prod

# This definition is not strictly needed, because `sum' works out-of-the-box
# with Measurement type, but this makes the function linear instead of quadratic
# in the number of arguments, but `result' is quadratic in the number of
# arguments, so in the end the function goes from cubic to quadratic.  Still not
# ideal, but this is an improvement.
sum{T<:Measurement}(a::AbstractArray{T}) =
    result(sum(value(a)), (ones(length(a))...), (a...))

# Same as above.  I'm not particularly proud of how the derivatives are
# computed, but something like this is needed in order to avoid errors with null
# nominal values: you may think to x ./ prod(x), but that would fail if one or
# more elements are zero.
function prod{T<:Measurement}(a::AbstractArray{T})
    x = value(a)
    return result(prod(x),
                  ntuple(i -> prod(deleteat!(copy(x), i)), length(x)),
                  (a...))
end


text:
# the code was tested on v0.4.2
# reinstall HetaSimulator
# ] add https://github.com/hetalang/HetaSimulator.jl

using HetaSimulator, Plots

p = load_platform(".")

# load scenarios
scn_df = read_scenarios("./scenarios.csv")
add_scenarios!(p, scn_df)

### Load measurements

# see table measurements.csv describing data with unknown variance
# see more detales in https://hetalang.github.io/HetaSimulator.jl/dev/table-formats/measurement/
measurements_df = read_measurements("./measurements.csv")
add_measurements!(p, measurements_df)

# display platform content
p

### simulate all
res = sim(p)

# plot all default
plot(res)
# plotd = plot(res)
# savefig(plotd, "fig.png")

# plot C1, C2 in log scale
plot(res, vars=[:C1,:C2], yscale=:log10, ylim=(1e-3, 1e3))
# plotd = plot(res, vars=[:C1,:C2], yscale=:log10, ylim=(1e-3, 1e3))
# savefig(plotd, "fig.png")

### Fitting

# fitted parameters
to_fit = [
    :kabs => 8.0,
    :Q => 4.0,
    :kel => 2.2,
    :sigma1 => 0.1,
    :sigma2 => 0.1,
    :sigma3 => 0.1,
]
res_optim = fit(p, to_fit) # default fitting

# optimal parameters
optim(res_optim)

# check fitting quality
res = sim(p, parameters_upd = optim(res_optim))
plot(res, yscale=:log10, ylims=(1e-3,1e2))
# plotd = plot(res_optim, vars=[:C1,:C2], yscale=:log10, ylims=(1e-3,1e2))
# savefig(plotd, "fig.png")

### Fitting with parameters table

# read parameters from table
params_df = read_parameters("./parameters.csv")
res_optim = fit(p, params_df)

### additional optimization-specific options

# see also https://nlopt.readthedocs.io/en/latest/NLopt_Algorithms/
res_optim = fit(
    p,
    params_df,
    fit_alg = :LN_SBPLX,
    ftol_abs = 1e-5,
    ftol_rel = 0,
    maxeval = 10^6
)
optim(res_optim)

#= see API docs
fit_alg : fitting algorithm. Default is :LN_NELDERMEAD
ftol_abs : absolute tolerance on function value. Default is 0.0
ftol_rel : relative tolerance on function value. Default is 1e-4
xtol_rel : relative tolerance on optimization parameters. Default is 0.0
xtol_rel : absolute tolerance on optimization parameters. Default is 0.0
maxeval : maximum number of function evaluations. Default is 1e4
maxtime : maximum optimization time (in seconds). Default is 0
=#


text:
using Gadfly
using Colors

include("pf_SQMC.jl")
include("sobol.jl")

function GordonKitagawaUpdate!(x,t, u, xtemp, i)
  b₁, b₂, b₃, b₄, σ = .5, 25, 8, 1.2, sqrt(1);
  xtemp[1,i] = b₁*x[1]+b₂*x[1]/(t+x[1]^2)+b₃*cos(b₄*t)+σ*erfinv(u[1]*2-1);
  return
end

function GordonKitagawaOut(x,t)
  σ = sqrt(.1);
  a = 20
  x[1]^2/a+σ*randn()[1]
end

function GordonKitagawaOutNoNoise(x, t)
  a = 20
  x[1].^2./a
end

function GordonKitagawapxy(yhat,y,t)
  σ = sqrt(.1);
  #w = 1/sqrt(2π)*e^(-(yhat[1]-y[1])^2/2)
  w = -(yhat[1]-y[1])^2/(2*σ^2)
end

function generateRealSequence(f!, g, x0, T, N = length(x0))
  x = Array{Float64,2}(length(x0), T)
  f!(x0, 1, rand(N),x,1)
  y1 = g(x[:,1],1)
  y = Array{Float64,2}(length(y1), T)
  y[:,1] = y1
  for t = 2:T
    f!(x[t-1], t, rand(N),x,t)
    y[:,t] = g(x[t], t)
  end
  x, y
end

function estMean(x,w,t)
  sum(x[:].*exp(w[:]))
end

function runTest(N, method, debug)
  T = 150

  ##Gordon Kitagawa
  f! = GordonKitagawaUpdate!#(x,t,u,xtemp,i) -> GordonKitagawaUpdate!(x[1],t,u,xtemp,i)
  g = GordonKitagawaOutNoNoise#(x,t) -> GordonKitagawaOutNoNoise(x[1])
  gn = GordonKitagawaOut#(x,t) -> GordonKitagawaOut(x[1])
  ginv = GordonKitagawapxy#(yhat, y,t) -> GordonKitagawapxy(yhat,y)
  rands = rand
  ##LTI System
  #σ = .5
  #f! = (x,t,u,xtemp,i) -> xtemp[1,i] = .8x[1]+4*erfinv(u[1]*2-1)
  #g = (x,t) -> 2*x[1]
  #gn = (x,t) -> 2*x[1] + σ*randn()[1]
  #ginv = (yhat, y, t) -> -(yhat[1]-y[1])^2/(2*σ^2)

  ##Estimator
  #est = (x,w,t) -> x[findmax(w)[2]]
  est = (x,w,t) -> sum(x[:].*exp(w[:]))

  x0 = .5
  x, y = generateRealSequence(f!, gn, x0, T)
  xhat = method(f!, g, ginv, rands, x0, y, est, N, debug=debug, xreal=x)
  x, xhat
end

function plotPoints(x, w, y, N, a, τ, t, xreal, xhat)
  c = w[:]-minimum(w)+1
  ##Use for GordonK
  #p = plot(layer(x=collect(1:N), y=x[:], Geom.point, color=c),
  #        layer(x=[1,N],y=ones(2).*sqrt(20*max(y[:,t],0)),Geom.line, Theme(default_color=color(colorant"red"))),
  #        layer(x=[1,N],y=-ones(2).*sqrt(20*max(y[:,t],0)),Geom.line, Theme(default_color=color(colorant"red"))),
  #        layer(x=[1,N],y=ones(2).*xreal[:,t],Geom.line, Theme(default_color=color(colorant"blue"),line_width=2px)),
  #        layer(x=[1,N],y=ones(2).*xhat[:,t],Geom.line, Theme(default_color=color(colorant"black"),line_width=4px)),
  #        Guide.XLabel("Particle "*string(t)), Guide.YLabel("Estimate"), Coord.Cartesian(ymin=-15,ymax=15))
  ##Use for LTI
  p = plot(layer(x=collect(1:N), y=x[:], Geom.point, color=c),
  layer(x=[1,N],y=ones(2)*1/2.*y[:,t],Geom.line, Theme(default_color=color(colorant"red"))),
  layer(x=[1,N],y=ones(2).*xreal[:,t],Geom.line, Theme(default_color=color(colorant"blue"),line_width=2px)),
  layer(x=[1,N],y=ones(2).*xhat[:,t],Geom.line, Theme(default_color=color(colorant"black"),line_width=4px)),
  Guide.XLabel("Particle "*string(t)), Guide.YLabel("Estimate"), Coord.Cartesian(ymin=-10,ymax=10))

  display(p)
  print("here")
  readline(STDIN)
end

function rms(x)
  sqrt(1/length(x)*sum(x.^2))
end

function testSQMC()
  debug = false
  Ns = 2.^(2:7)
  M = 50
  RMS = Array{Float64,2}(length(Ns),M)
  largeRMS = Array{Float64}(length(Ns),2)
  rmsMean = Array{Float64,2}(length(Ns),2)
  rmsVariance = Array{Float64,2}(length(Ns),2)
  for (methodidx,method) in enumerate([SQMC,pf])
    xhat, xreal = 0, 0
    @time for (i, N) in enumerate(Ns)
      rmslocal = Array{Float64,1}(M)
      for j = 1:M
        xreal, xhat = runTest(N,method,debug)
        RMS[i,j] = rms(xreal-xhat)
      end
    end
    rmsMean[:,methodidx] = mean(RMS,2)
    rmsVariance[:,methodidx] = std(RMS,2)./sqrt(M)
    for (i, N) in enumerate(Ns)
      largeRMS[i,methodidx] = length(find(RMS[i,:].>rmsMean[i,methodidx]+2*rmsVariance[i,methodidx]))
    end
  end
  rmsMean, rmsVariance, Ns, largeRMS, RMS
end

function testPlot(rmsMean, rmsVariance, Ns)
  p = plot(
  layer(x=Ns,y=rmsMean[:,1],Geom.line,Theme(default_color=color(colorant"red"))),
  layer(x=Ns,y=rmsMean[:,2],Geom.line,Theme(default_color=color(colorant"blue"))),
  layer(x=Ns,y=rmsMean[:,1]+rmsVariance[:,1]*2,Geom.line,Theme(default_color=color(colorant"red"))),
  layer(x=Ns,y=rmsMean[:,1]-rmsVariance[:,1]*2,Geom.line,Theme(default_color=color(colorant"red"))),
  layer(x=Ns,y=rmsMean[:,2]+rmsVariance[:,2]*2,Geom.line,Theme(default_color=color(colorant"blue"))),
  layer(x=Ns,y=rmsMean[:,2]-rmsVariance[:,2]*2,Geom.line,Theme(default_color=color(colorant"blue"))),
  )
end

function testPlot(N)
  for j = 1:N
    M = 1000
    vals = Array{Float64,2}(M,2)
    ds = linspace(rand()/M,1,M)[1:end]
    for i = 1:M
      o =  d2xy(ds[i])
      vals[i,:] = [o[1], o[2]]
    end
    plot(vals[:,1], vals[:,2])
  end
end

function testPlot2(N)
  for j = 1:N
    M = 1000
    vals = Array{Float64,2}(M,2)
    ds = sort(rand(M))
    for i = 1:M
      o =  d2xy(ds[i])
      vals[i,:] = [o[1], o[2]]
    end
    plot(vals[:,1], vals[:,2])
  end
end


text:
module Tutorial3
using FEniCS

T = 2.0            # final time
num_steps = 10     # number of time steps
dt = T / num_steps # time step size
alpha = 3          # parameter alpha
beta = 1.2
nx = ny = 8
mesh = UnitSquareMesh(nx, ny)
V = FunctionSpace(mesh, "P", 1)
u_D = Expression("1 + x[0]*x[0] + alpha*x[1]*x[1] + beta*t",
degree=2, alpha=alpha, beta=beta, t=0)

bc = DirichletBC(V, u_D, "on_boundary")
u_n = interpolate(u_D, V)
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(beta - 2 - 2*alpha)
F = u*v*dx + dt*dot(grad(u), grad(v))*dx - (u_n + dt*f)*v*dx

a, L = lhs(F),rhs(F)

u=FeFunction(V)
global t = 0

for n = 0:(num_steps-1)
    global t=t+dt
    u_D.pyobject.t=t
    lvsolve(a,L,u,bc)
    u_e= interpolate(u_D,V)
    vv = get_array(u_e)
    ww = get_array(u)
    error = maximum(abs.(vv-ww))
    #@printf "t = %f: error = %10g" t error
    assign(u_n,u)
end
#colorbar()
end#module


text:
"""
Abstract type for Device Formulations (a.k.a Models)

# Example

import PowerSimulations
const PSI = PowerSimulations
struct MyCustomDeviceFormulation <: PSI.AbstractDeviceFormulation
"""
abstract type AbstractDeviceFormulation end

########################### Thermal Generation Formulations ################################
abstract type AbstractThermalFormulation <: AbstractDeviceFormulation end
abstract type AbstractThermalDispatchFormulation <: AbstractThermalFormulation end
abstract type AbstractThermalUnitCommitment <: AbstractThermalFormulation end

abstract type AbstractStandardUnitCommitment <: AbstractThermalUnitCommitment end
abstract type AbstractCompactUnitCommitment <: AbstractThermalUnitCommitment end

struct ThermalBasicUnitCommitment <: AbstractStandardUnitCommitment end
struct ThermalStandardUnitCommitment <: AbstractStandardUnitCommitment end
struct ThermalBasicDispatch <: AbstractThermalDispatchFormulation end
struct ThermalStandardDispatch <: AbstractThermalDispatchFormulation end
struct ThermalDispatchNoMin <: AbstractThermalDispatchFormulation end

struct ThermalMultiStartUnitCommitment <: AbstractCompactUnitCommitment end
struct ThermalCompactUnitCommitment <: AbstractCompactUnitCommitment end
struct ThermalBasicCompactUnitCommitment <: AbstractCompactUnitCommitment end
struct ThermalCompactDispatch <: AbstractThermalDispatchFormulation end

############################# Electric Load Formulations ###################################
abstract type AbstractLoadFormulation <: AbstractDeviceFormulation end
abstract type AbstractControllablePowerLoadFormulation <: AbstractLoadFormulation end
struct StaticPowerLoad <: AbstractLoadFormulation end
struct InterruptiblePowerLoad <: AbstractControllablePowerLoadFormulation end
struct DispatchablePowerLoad <: AbstractControllablePowerLoadFormulation end

########################### Hybrid Generation Formulations ################################
abstract type AbstractHybridFormulation <: AbstractDeviceFormulation end
abstract type AbstractStandardHybridFormulation <: AbstractHybridFormulation end
struct BasicHybridDispatch <: AbstractHybridFormulation end
struct StandardHybridDispatch <: AbstractStandardHybridFormulation end

############################ Hydro Generation Formulations #################################
abstract type AbstractHydroFormulation <: AbstractDeviceFormulation end
abstract type AbstractHydroDispatchFormulation <: AbstractHydroFormulation end
abstract type AbstractHydroUnitCommitment <: AbstractHydroFormulation end
abstract type AbstractHydroReservoirFormulation <: AbstractHydroDispatchFormulation end
struct HydroDispatchRunOfRiver <: AbstractHydroDispatchFormulation end
struct HydroDispatchReservoirBudget <: AbstractHydroReservoirFormulation end
struct HydroDispatchReservoirStorage <: AbstractHydroReservoirFormulation end
struct HydroDispatchPumpedStorage <: AbstractHydroReservoirFormulation end
struct HydroCommitmentRunOfRiver <: AbstractHydroUnitCommitment end
struct HydroCommitmentReservoirBudget <: AbstractHydroUnitCommitment end
struct HydroCommitmentReservoirStorage <: AbstractHydroUnitCommitment end

############################ Regulation Device Formulations ################################
abstract type AbstractRegulationFormulation <: AbstractDeviceFormulation end
struct ReserveLimitedRegulation <: AbstractRegulationFormulation end
struct DeviceLimitedRegulation <: AbstractRegulationFormulation end

########################### Renewable Generation Formulations ##############################
abstract type AbstractRenewableFormulation <: AbstractDeviceFormulation end
abstract type AbstractRenewableDispatchFormulation <: AbstractRenewableFormulation end
struct RenewableFullDispatch <: AbstractRenewableDispatchFormulation end
struct RenewableConstantPowerFactor <: AbstractRenewableDispatchFormulation end

############################ Storage Generation Formulations ###############################
abstract type AbstractStorageFormulation <: AbstractDeviceFormulation end
abstract type AbstractEnergyManagement <: AbstractStorageFormulation end
struct BookKeeping <: AbstractStorageFormulation end
struct BatteryAncillaryServices <: AbstractStorageFormulation end
struct EnergyTarget <: AbstractEnergyManagement end

"""
Abstract type for Branch Formulations (a.k.a Models)

# Example
import PowerSimulations
const PSI = PowerSimulations
struct MyCustomBranchFormulation <: PSI.AbstractDeviceFormulation
"""
# Generic Branch Models
abstract type AbstractBranchFormulation <: AbstractDeviceFormulation end

############################### AC Branch Formulations #####################################
"""
Branch type to add unbounded flow variables and use flow constraints
"""
struct StaticBranch <: AbstractBranchFormulation end
"""
Branch type to add bounded flow variables and use flow constraints
"""
struct StaticBranchBounds <: AbstractBranchFormulation end
"""
Branch type to avoid flow constraints
"""
struct StaticBranchUnbounded <: AbstractBranchFormulation end

############################### DC Branch Formulations #####################################
abstract type AbstractDCLineFormulation <: AbstractBranchFormulation end
struct HVDCUnbounded <: AbstractDCLineFormulation end
struct HVDCLossless <: AbstractDCLineFormulation end
struct HVDCDispatch <: AbstractDCLineFormulation end
# Not Implemented
# struct VoltageSourceDC <: AbstractDCLineFormulation end

############################## Network Model Formulations ##################################
# These formulations are taken directly from PowerModels

abstract type AbstractPTDFModel <: PM.AbstractDCPModel end
struct StandardPTDFModel <: AbstractPTDFModel end
struct PTDFPowerModel <: AbstractPTDFModel end

struct CopperPlatePowerModel <: PM.AbstractActivePowerModel end
struct AreaBalancePowerModel <: PM.AbstractActivePowerModel end

#================================================
    # exact non-convex models
    ACPPowerModel, ACRPowerModel, ACTPowerModel

    # linear approximations
    DCPPowerModel, NFAPowerModel

    # quadratic approximations
    DCPLLPowerModel, LPACCPowerModel

    # quadratic relaxations
    SOCWRPowerModel, SOCWRConicPowerModel,
    SOCBFPowerModel, SOCBFConicPowerModel,
    QCRMPowerModel, QCLSPowerModel,

    # sdp relaxations
    SDPWRMPowerModel, SparseSDPWRMPowerModel
================================================#

##### Exact Non-Convex Models #####
import PowerModels: ACPPowerModel

import PowerModels: ACRPowerModel

import PowerModels: ACTPowerModel

##### Linear Approximations #####
import PowerModels: DCPPowerModel

import PowerModels: NFAPowerModel

##### Quadratic Approximations #####
import PowerModels: DCPLLPowerModel

import PowerModels: LPACCPowerModel

##### Quadratic Relaxations #####
import PowerModels: SOCWRPowerModel

import PowerModels: SOCWRConicPowerModel

import PowerModels: QCRMPowerModel

import PowerModels: QCLSPowerModel

abstract type AbstractServiceFormulation end

abstract type AbstractAGCFormulation <: AbstractServiceFormulation end

struct PIDSmoothACE <: AbstractAGCFormulation end

abstract type AbstractReservesFormulation <: AbstractServiceFormulation end

struct GroupReserve <: AbstractReservesFormulation end
struct RangeReserve <: AbstractReservesFormulation end
struct StepwiseCostReserve <: AbstractReservesFormulation end
struct RampReserve <: AbstractReservesFormulation end
struct NonSpinningReserve <: AbstractReservesFormulation end


text:
using LoopVectorization
using ExplainableAI
using Test
using ReferenceTests

@testset "ExplainableAI.jl" begin
    @testset "Utilities" begin
        @info "Running tests on utilities..."
        include("test_utils.jl")
    end
    @testset "Neuron selection" begin
        @info "Running tests on neuron selection..."
        include("test_neuron_selection.jl")
    end
    @testset "Input augmentation" begin
        @info "Running tests on input augmentation..."
        include("test_input_augmentation.jl")
    end
    @testset "Heatmaps" begin
        @info "Running tests on heatmaps..."
        include("test_heatmaps.jl")
    end
    @testset "Canonize" begin
        @info "Running tests on model canonization..."
        include("test_canonize.jl")
    end
    @testset "LRP model checks" begin
        @info "Running tests on LRP model checks..."
        include("test_checks.jl")
    end
    @testset "LRP rules" begin
        @info "Running tests on LRP rules..."
        include("test_rules.jl")
    end
    @testset "Batches" begin
        @info "Running analyzer tests on batches..."
        include("test_batches.jl")
    end
    @testset "VGG11" begin
        @info "Running analyzer tests on VGG11..."
        include("test_vgg11.jl")
    end
end


text:
module FxRobot
using Downloads: download
using Dates
using ZipFile: Reader
using DataFrames
using CSV
using TimeZones
using Chain
import Arrow

export
    unzip_file,ts2dt,dt2ts, # utils.jl
    csv_to_dataframe, # df.jl
    arrow_read,arrow_write, # arrow.jl
    download_monthly_binance_data, # binance.jl
    download_daily_binance_data,
    trans_df

include("utils.jl")
include("./feed/binance.jl")
include("./feed/df.jl")
include("./feed/arrow.jl")
end

max_stars_repo_name: mrchaos/fxrobot

max_issues_repo_name: mrchaos/fxrobot

max_forks_repo_name: mrchaos/fxrobot


text:
using DistCtrl4DistMan
using LinearAlgebra
using Test

# @testset "DistCtrl4DistMan.jl" begin
# end

@testset "ObjectAgentModule.jl" begin
    include("test_data.jl")

    @testset "Test - MAG" begin# Used types
        F = Float64;
        U = UInt64;

        # Load test data
        Gxy_ref, MAG_test_currents, MAG_fi_ref, MAG_F_ref = test_data_MAG();

        n = U(8);
        dx = F(25e-3);
        aa = DistCtrl4DistMan.ActuatorArray(n, n, dx, dx, :coil);
        ball_pos = ( (n-1)*dx/2 - dx/2, (n-1)*dx/2 + 1.3*dx );

        # calcMAGForce() and fi() test values
        fi_test = DistCtrl4DistMan.ObjectAgentModule.mag_fi(MAG_test_currents);
        F_MAG = DistCtrl4DistMan.calcMAGForce(aa, ball_pos, fi_test);

        aL, a_used = DistCtrl4DistMan.genActList(aa, (ball_pos[1], ball_pos[2], F(0)), F(2.5*dx), F(3.5*dx));
        oa_mag = DistCtrl4DistMan.ObjectAgent_MAG("Agent_MAG", ball_pos, F_MAG, aa, aL, a_used, F(1));

        @testset "Test the function generating the Gxy matrix" begin
            @test DistCtrl4DistMan.ObjectAgentModule.genGxy(aa, ball_pos) ≈ Gxy_ref
        end
        @testset "Test the function fi()" begin
            @test fi_test ≈ MAG_fi_ref
        end
        @testset "Test the function calcMagForce()" begin
            @test all(F_MAG .≈ MAG_F_ref)
        end
    end

    @testset "Test - DEP" begin
        # Used types
        F = Float64;
        U = UInt64;

        # Load test data
        Gamma_ref, Lambda_x, Lambda_y, Lambda_z, DEP_test_phases, DEP_F_ref = test_data_DEP();
        n = U(5);
        dx = F(100e-6);
        aa = DistCtrl4DistMan.ActuatorArray(n, n, dx, dx/2, :square);
        center_pos = ( (n-1)*dx/2, (n-1)*dx/2, F(100e-6));

        F_DEP = DistCtrl4DistMan.calcDEPForce(aa, center_pos, DEP_test_phases);
        aL, a_used = DistCtrl4DistMan.genActList(aa, center_pos, F(350e-6), F(350e-6));
        oa_dep = DistCtrl4DistMan.ObjectAgent_DEP("Agent1", center_pos, F_DEP, aa, aL, a_used);

        @testset "Test the function generating the Gamma and Lambda_a marices" begin
            Gamma_test, Lambda_a = DistCtrl4DistMan.ObjectAgentModule.genC_DEP(aa, center_pos);
            @test Gamma_ref ≈ Gamma_test
            @test Lambda_a[1] ≈ Lambda_x
            @test Lambda_a[2] ≈ Lambda_y
            @test Lambda_a[3] ≈ Lambda_z
        end
        @testset "Test the function calculating the DEP force" begin
            @test all(F_DEP .≈ DEP_F_ref)
        end
        @testset "Test the cost function" begin # The cost function (oa_dep,fvk) should be zero since we set oa_dep.vk_r and oa_dep.vk_i to values which should result in the desired force
            # Calculate the cost function (oa_dep.fvk)
            oa_dep.vk_r = cos.(DEP_test_phases'[:]);
            oa_dep.vk_i = sin.(DEP_test_phases'[:]);
            DistCtrl4DistMan.costFun!(oa_dep);

            @test all( abs.(oa_dep.fvk) .< [1,1,1].*1e-6 )
        end
    end

    @testset "Test - ACU" begin
        # Used types
        F = Float64;
        U = UInt64;

        # Load test data
        C1ref, C2ref, th1, th2, objVal1_ref, objVal2_ref, D1ref, D2ref, phases_test, p1_ref, p2_ref = test_data_ACU();

        n = U(8);
        dx = F(10.0e-3);
        aa = DistCtrl4DistMan.ActuatorArray(n, n, dx);

        # two agents at positions close to the center
        z0 = F(-50e-3);
        maxDist1 = 3*dx; #2
        maxDist2 = 3*dx; #3
        oa1_pos = ((n-1)*aa.dx/2-dx, (n-1)*aa.dx/2, z0);
        oa2_pos = ((n-1)*aa.dx/2+dx, (n-1)*aa.dx/2+1.6f0*dx, z0);
        Pdes = F(1000); # desired presure (same for both agents)

        # Generate list of used actuators for each object agent
        aL1, a_used1 = DistCtrl4DistMan.genActList(aa, oa1_pos, maxDist1, maxDist2);
        aL2, a_used2 = DistCtrl4DistMan.genActList(aa, oa2_pos, maxDist1, maxDist2);

        oa1 = DistCtrl4DistMan.ObjectAgent_ACU("Agent1", oa1_pos, Pdes, aa, aL1, a_used1);
        oa2 = DistCtrl4DistMan.ObjectAgent_ACU("Agent2", oa2_pos, Pdes, aa, aL2, a_used2);

        @testset "Test the function generating the C matrix" begin
            @test isapprox(oa1.zvec*oa1.zvec' + 1im*oa1.zvec*[0 1; -1 0]*oa1.zvec', C1ref, atol=1e-6)
            @test isapprox(oa2.zvec*oa2.zvec' + 1im*oa2.zvec*[0 1; -1 0]*oa2.zvec', C2ref, atol=1e-6)
        end
        @testset "Test the cost function" begin
            oa1.vk_r = cos.(th1);
            oa1.vk_i = sin.(th1);
            oa2.vk_r = cos.(th2);
            oa2.vk_i = sin.(th2);

            # Compute the values of the cost function (oa.fvk)
            DistCtrl4DistMan.costFun!(oa1)
            DistCtrl4DistMan.costFun!(oa2)

            @test isapprox(oa1.fvk[1],  objVal1_ref, atol=1e-5)
            @test isapprox(oa2.fvk[1],  objVal2_ref, atol=1e-5)
        end
        @testset "Test the cost function calculating the jacobian" begin
            # Compute the jacobians
            DistCtrl4DistMan.ObjectAgentModule.jac!(oa1);
            DistCtrl4DistMan.ObjectAgentModule.jac!(oa2);

            @test isapprox(oa1.J,  D1ref', atol=1e-6)
            @test isapprox(oa2.J,  D2ref', atol=1e-6)
        end
        @testset "Test the function calculating the pressure" begin
            @test isapprox(DistCtrl4DistMan.calcPressure(aa, oa1_pos, convert.(F, phases_test)), F(p1_ref), atol=1e-4)
            @test isapprox(DistCtrl4DistMan.calcPressure(aa, oa2_pos, convert.(F, phases_test)), F(p2_ref), atol=1e-4)
        end
        @testset "Test the function finding the shared actuators among agents (function from ActuatorArray module)" begin
            DistCtrl4DistMan.resolveNeighbrRelations!([oa1, oa2])
            @test oa1.actList[oa1.neighbors[1][2]] == oa2.actList[oa1.neighbors[1][3]]
            @test oa2.actList[oa2.neighbors[1][2]] == oa1.actList[oa2.neighbors[1][3]]
        end
    end
end

@testset "ActuatorArray tests" begin
    # Used types
    F = Float64;
    U = UInt64;

    # Load test data
    aL1_ref = Array{Tuple{UInt8,UInt8},1}([(1, 3), (1, 4), (1, 5), (1, 6), (2, 2), (2, 3), (2, 4), (2, 5), (2, 6), (2,7), (3, 2), (3, 3), (3, 4), (3, 5), (3, 6), (3, 7), (4, 2), (4, 3), (4, 4), (4, 5), (4, 6), (4, 7), (5, 2), (5, 3),(5, 4), (5, 5), (5, 6), (5, 7), (6, 3), (6, 4), (6, 5), (6, 6)]);
    aL2_ref = Array{Tuple{UInt8,UInt8},1}([(3, 5), (3, 6), (3, 7), (4, 4), (4, 5), (4, 6), (4, 7), (4, 8), (5, 4), (5, 5), (5, 6), (5, 7), (5, 8), (6, 4), (6, 5), (6, 6), (6, 7), (6, 8), (7, 4), (7, 5), (7, 6), (7, 7), (7, 8), (8, 5), (8, 6), (8, 7)]);

    n = U(8);
    dx = F(10.0e-3);
    aa = DistCtrl4DistMan.ActuatorArray(n, n, dx);

    # two agents at positions close to the center
    z0 = F(-50e-3);
    maxDist1 = 3*dx; #2
    maxDist2 = 3*dx; #3
    oa1_pos = ((n-1)*aa.dx/2-dx, (n-1)*aa.dx/2, z0);
    oa2_pos = ((n-1)*aa.dx/2+dx, (n-1)*aa.dx/2+1.6f0*dx, z0);

    # Generate list of used actuators for each object agent
    aL1, a_used1 = DistCtrl4DistMan.genActList(aa, oa1_pos, maxDist1, maxDist2);
    aL2, a_used2 = DistCtrl4DistMan.genActList(aa, oa2_pos, maxDist1, maxDist2);

    @testset "Test the function generating the actuator list" begin
        @test aL1 == aL1_ref
        @test aL2 == aL2_ref
    end
end

@testset "LinSolvers.jl package tests" begin
    n = 100;
    m = 5;

    @testset "Test backward and forward substitution for solution of R^T*Rx=b" begin
        for i=1:20
            R = Matrix(UpperTriangular(rand(n,n) + I));
            A = R'*R;
            x_test = rand(n);
            b_test = A*x_test;
            b = copy(b_test)

            DistCtrl4DistMan.LinSolvers.forward_substitution!(R, b)
            DistCtrl4DistMan.LinSolvers.backward_substitution!(R, b)
            x_sol = b;

            @test  A*x_sol ≈ b_test
        end
    end
    @testset "Test solving of (D^T*D + 1/λ*I)x=D^T*b by kernel_solve!()" begin
        x_sol = zeros(n);
        lsd = DistCtrl4DistMan.LinSystemData(n, m, x_sol);
        λ = 1.0e3;
        for i in 1:20
            Dt = rand(n,m);
            b  = rand(m);
            DistCtrl4DistMan.kernel_solve!(lsd, Dt, λ, b);

            @test  (Dt*Dt' + 1/λ*I)*x_sol ≈ Dt*b
        end
    end
    @testset "Test solving of (D^T*D + 1/λ*I)x=D^T*b by kernel_solve()" begin
        λ = 1.0e3;
        for i in 1:20
            Dt = rand(n,m);
            b  = rand(m);
            x_sol = DistCtrl4DistMan.kernel_solve(Dt, λ, b);

            @test  (Dt*Dt' + 1/λ*I)*x_sol ≈ Dt*b
        end
    end
    @testset "Test solving of (D^T*D + 1/λ*I)x=D^T*b by kernel_solve_prealloc!()" begin
        λ = 1.0e3;
        x_sol = zeros(n);
        A_prealloc = zeros(m+n,m);
        b_prealloc = zeros(m);
        R_prealloc = zeros(m,m);

        for i in 1:20
            Dt = rand(n,m);
            b  = rand(m);
            DistCtrl4DistMan.LinSolvers.kernel_solve_prealloc!(Dt, λ, b, x_sol, A_prealloc, b_prealloc, R_prealloc);

            @test  (Dt*Dt' + 1/λ*I)*x_sol ≈ Dt*b
        end
    end

    @testset "Test qr_r!() function" begin
        for i in 1:20
            A = rand(n,m);
            R  = zeros(m,m);

            # Compute the upper triangular R matrix from QR factorization by qr_r!()
            # Matrix A is modified in qr_r!(), that is why we use copy of it as an argument
            DistCtrl4DistMan.LinSolvers.qr_r!(copy(A), R);

            # Compute the cholesky factor of AᵀA which has to match to R
            C = cholesky(A'*A);

            @test  C.U ≈ R
        end
    end

end

@testset "Test - ADMM" begin
    F = Float64;
    U = UInt64;

    n = U(8);
    dx = F(25e-3);
    aa = DistCtrl4DistMan.ActuatorArray(n, n, dx, dx, :coil);

    # Set the maximum distance so that both agents use all the coils
    maxdist = n*dx;

    for i in 1:20
        # Generate random positions of the agents
        ball_pos1 = ( dx*(n-1)*rand(),  dx*(n-1)*rand());
        ball_pos2 = ( dx*(n-1)*rand(),  dx*(n-1)*rand());

        # Generate the desired forces so that one can actually generate them
        random_ctrls = rand(F, n, n);
        F_des1 = DistCtrl4DistMan.calcMAGForce(aa, ball_pos1, random_ctrls)
        F_des2 = DistCtrl4DistMan.calcMAGForce(aa, ball_pos2, random_ctrls)

        # Initialize the agents
        aL1, a_used1 = DistCtrl4DistMan.genActList(aa, (ball_pos1[1], ball_pos1[2], F(0)), maxdist, maxdist);
        aL2, a_used2 = DistCtrl4DistMan.genActList(aa, (ball_pos2[1], ball_pos2[2], F(0)), maxdist, maxdist);
        oa1 = DistCtrl4DistMan.ObjectAgent_MAG("Agent 1", ball_pos1, F_des1, aa, aL1, a_used1, F(1));
        oa2 = DistCtrl4DistMan.ObjectAgent_MAG("Agent 2", ball_pos2, F_des2, aa, aL2, a_used2, F(1));

        agents = [oa1, oa2];
        DistCtrl4DistMan.admm(agents,
            λ = 1.0, ρ = 1.0,
            log = false,
            maxiter = 1000,
            method = :freedir);

        F_dev1 = tuple((oa1.Gxy'*oa1.zk*oa1.Fdes_sc)...)
        F_dev2 = tuple((oa2.Gxy'*oa2.zk*oa2.Fdes_sc)...)

        @test all(isapprox.(F_dev1, F_des1, atol=1e-4))
        @test all(isapprox.(F_dev2, F_des2, atol=1e-4))
    end
end


text:
function _container_spec(m::M, ax...) where M <: JuMP.AbstractModel
    return JuMP.Containers.DenseAxisArray{JuMP.variable_type(m)}(undef, ax...)
end

function add_variable(ps_m::CanonicalModel,
                      devices::Array{T,1},
                      time_range::UnitRange{Int64},
                      var_name::Symbol,
                      binary::Bool) where {T <: PSY.Device}

    ps_m.variables[var_name] = _container_spec(ps_m.JuMPmodel, [d.name for d in devices], time_range)

   for t in time_range, d in devices
       ps_m.variables[var_name][d.name,t] = JuMP.@variable(ps_m.JuMPmodel, base_name="$(var_name)_{$(d.name),$(t)}", start = 0.0, binary=binary)
   end

   return

end

function add_variable(ps_m::CanonicalModel,
                      devices::Array{T,1},
                      time_range::UnitRange{Int64},
                      var_name::Symbol,
                      binary::Bool,
                      expression::Symbol) where {T <: PSY.Device}

    ps_m.variables[var_name] = _container_spec(ps_m.JuMPmodel, [d.name for d in devices], time_range)

   for t in time_range, d in devices
       ps_m.variables[var_name][d.name,t] = JuMP.@variable(ps_m.JuMPmodel, base_name="$(var_name)_{$(d.name),$(t)}", start = 0.0, binary=binary)
       _add_to_expression!(ps_m.expressions[expression], d.bus.number, t, ps_m.variables[var_name][d.name,t])
   end

   return

end

function add_variable(ps_m::CanonicalModel,
                      devices::Array{T,1},
                      time_range::UnitRange{Int64},
                      var_name::Symbol,
                      binary::Bool,
                      expression::Symbol,
                      sign::Int64) where {T <: PSY.Device}

    ps_m.variables[var_name] = _container_spec(ps_m.JuMPmodel, [d.name for d in devices], time_range)

   for t in time_range, d in devices
       ps_m.variables[var_name][d.name,t] = JuMP.@variable(ps_m.JuMPmodel, base_name="$(var_name)_{$(d.name),$(t)}", start = 0.0, binary=binary)
       _add_to_expression!(ps_m.expressions[expression], d.bus.number, t, ps_m.variables[var_name][d.name,t], sign)
   end

   return

end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiExtensionsV1beta1NetworkPolicyIngressRule <: SwaggerModel
    from::Any # spec type: Union{ Nothing, Vector{IoK8sApiExtensionsV1beta1NetworkPolicyPeer} } # spec name: from
    ports::Any # spec type: Union{ Nothing, Vector{IoK8sApiExtensionsV1beta1NetworkPolicyPort} } # spec name: ports

    function IoK8sApiExtensionsV1beta1NetworkPolicyIngressRule(;from=nothing, ports=nothing)
        o = new()
        validate_property(IoK8sApiExtensionsV1beta1NetworkPolicyIngressRule, Symbol("from"), from)
        setfield!(o, Symbol("from"), from)
        validate_property(IoK8sApiExtensionsV1beta1NetworkPolicyIngressRule, Symbol("ports"), ports)
        setfield!(o, Symbol("ports"), ports)
        o
    end
end # type IoK8sApiExtensionsV1beta1NetworkPolicyIngressRule

const _property_map_IoK8sApiExtensionsV1beta1NetworkPolicyIngressRule = Dict{Symbol,Symbol}(Symbol("from")=>Symbol("from"), Symbol("ports")=>Symbol("ports"))
const _property_types_IoK8sApiExtensionsV1beta1NetworkPolicyIngressRule = Dict{Symbol,String}(Symbol("from")=>"Vector{IoK8sApiExtensionsV1beta1NetworkPolicyPeer}", Symbol("ports")=>"Vector{IoK8sApiExtensionsV1beta1NetworkPolicyPort}")
Base.propertynames(::Type{ IoK8sApiExtensionsV1beta1NetworkPolicyIngressRule }) = collect(keys(_property_map_IoK8sApiExtensionsV1beta1NetworkPolicyIngressRule))
Swagger.property_type(::Type{ IoK8sApiExtensionsV1beta1NetworkPolicyIngressRule }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiExtensionsV1beta1NetworkPolicyIngressRule[name]))}
Swagger.field_name(::Type{ IoK8sApiExtensionsV1beta1NetworkPolicyIngressRule }, property_name::Symbol) =  _property_map_IoK8sApiExtensionsV1beta1NetworkPolicyIngressRule[property_name]

function check_required(o::IoK8sApiExtensionsV1beta1NetworkPolicyIngressRule)
    true
end

function validate_property(::Type{ IoK8sApiExtensionsV1beta1NetworkPolicyIngressRule }, name::Symbol, val)
end


text:
# module VoxelPlots.Goxel

struct Mesh
end

struct Proc
end

function run(proc::Proc, mesh::Mesh)
end

function save(mesh::Mesh, path::String)
end

# module VoxelPlots.Goxel


text:
"Parameters for right-hand side function of explicit Runge-Kutta methods."
struct ParametersERK{DT, TT, D, S, ET <: NamedTuple} <: Parameters{DT,TT}
    equs::ET
    tab::Tableau{TT}
    Δt::TT

    function ParametersERK{DT,D}(equs::ET, tab::Tableau{TT}, Δt::TT) where {DT, TT, D, ET <: NamedTuple}
        new{DT, TT, D, tab.s, ET}(equs, tab, Δt)
    end
end

"Explicit Runge-Kutta integrator cache."
struct IntegratorCacheERK{DT,D,S} <: ODEIntegratorCache{DT,D}
    Q::Vector{Vector{DT}}
    V::Vector{Vector{DT}}

    function IntegratorCacheERK{DT,D,S}() where {DT,D,S}
        Q = create_internal_stage_vector(DT, D, S)
        V = create_internal_stage_vector(DT, D, S)
        new(Q, V)
    end
end

function IntegratorCache{ST}(params::ParametersERK{DT,TT,D,S}; kwargs...) where {ST,DT,TT,D,S}
    IntegratorCacheERK{ST,D,S}(; kwargs...)
end

@inline CacheType(ST, params::ParametersERK{DT,TT,D,S}) where {DT,TT,D,S} = IntegratorCacheERK{ST,D,S}

@doc raw"""
Explicit Runge-Kutta integrator solving the system
```math
\begin{aligned}
V_{n,i} &= v (Q_{n,i}, P_{n,i}) , &
Q_{n,i} &= q_{n} + h \sum \limits_{j=1}^{s} a_{ij} \, V_{n,j} , &
q_{n+1} &= q_{n} + h \sum \limits_{i=1}^{s} b_{i} \, V_{n,i} .
\end{aligned}
```
"""
struct IntegratorERK{DT, TT, D, S, ET} <: AbstractIntegratorRK{DT,TT}
    params::ParametersERK{DT,TT,D,S,ET}
    caches::CacheDict{ParametersERK{DT,TT,D,S,ET}}

    function IntegratorERK(params::ParametersERK{DT,TT,D,S,ET}, caches) where {DT,TT,D,S,ET}
        new{DT, TT, D, S, ET}(params, caches)
    end

    function IntegratorERK{DT,D}(equations::ET, tableau::Tableau{TT}, Δt::TT) where {DT, TT, D, ET <: NamedTuple}
        # get number of stages
        S = tableau.s

        # check if tableau is explicit
        @assert isexplicit(tableau)

        # create params
        params = ParametersERK{DT,D}(equations, tableau, Δt)

        # create cache dict
        caches = CacheDict(params)

        # create integrator
        IntegratorERK(params, caches)
    end

    function IntegratorERK{DT,D}(v::Function, tableau::Tableau{TT}, Δt::TT; kwargs...) where {DT,TT,D}
        IntegratorERK{DT,D}(NamedTuple{(:v,)}((v,)), tableau, Δt; kwargs...)
    end

    function IntegratorERK{DT,D}(v::Function, h::Function, tableau::Tableau{TT}, Δt::TT; kwargs...) where {DT,TT,D}
        IntegratorERK{DT,D}(NamedTuple{(:v,:h)}((v,h)), tableau, Δt; kwargs...)
    end

    function IntegratorERK(equation::ODE{DT}, tableau::Tableau{TT}, Δt::TT; kwargs...) where {DT,TT}
        IntegratorERK{DT, ndims(equation)}(get_functions(equation), tableau, Δt; kwargs...)
    end
end

@inline Base.ndims(::IntegratorERK{DT,TT,D,S}) where {DT,TT,D,S} = D

function integrate_step!(int::IntegratorERK{DT,TT}, sol::AtomicSolutionODE{DT,TT},
                         cache::IntegratorCacheERK{DT}=int.caches[DT]) where {DT,TT}
    # temporary variables
    local tᵢ::TT
    local yᵢ::DT

    # reset atomic solution
    reset!(sol, timestep(int))

    # compute internal stages
    for i in eachstage(int)
        for k in eachindex(cache.Q[i], cache.V[i])
            yᵢ = 0
            for j in 1:i-1
                yᵢ += tableau(int).a[i,j] * cache.V[j][k]
            end
            cache.Q[i][k] = sol.q̄[k] + timestep(int) * yᵢ
        end
        tᵢ = sol.t̄ + timestep(int) * tableau(int).c[i]
        equations(int)[:v](tᵢ, cache.Q[i], cache.V[i])
    end

    # compute final update
    update_solution!(sol.q, sol.q̃, cache.V, tableau(int).b, tableau(int).b̂, timestep(int))
end


text:
abstract type AbstractFixPoint end

"""
  Solvers.fixpoint(f,x0::Real,method=SimpleFixPoint())
does a fixpoint iteration convergence on a series of real numbers.
f is a function that should
the following strategies:
 SimpleFixPoint(): simple fixpoint, evaluates until criteria are met
 AitkenFixPoint(): uses Aitken's delta-squared process to accelerate convergence of the series
"""

struct SimpleFixPoint{T<:Real} <: AbstractFixPoint
    dampingfactor::T
end

SimpleFixPoint(;dampingfactor=1.0) = SimpleFixPoint(dampingfactor)

struct AitkenFixPoint <: AbstractFixPoint end

function fixpoint(f,x0::Real,
    method::AbstractFixPoint = SimpleFixPoint();
    atol=zero(nested_eltype(x0)),
    rtol=8eps(one(nested_eltype(x0))),
    max_iters=100,
    norm = norm,
    antinan=false)

    x0,atol,rtol = promote(x0,atol,rtol)
    return _fixpoint(f,x0,method,atol,rtol,max_iters,norm,antinan)
end

function _fixpoint(f,
        x0::T,
        method::SimpleFixPoint = SimpleFixPoint(),
        atol::T = zero(T),
        rtol::T =8*eps(T),
        max_iters=100,
        norm=norm,
        antinan=false) where {T}

    nan = (0*atol)/(0*atol)

    xi = f(x0)
    itercount = 1 #one iteration is already done
    #if fixed point is already the initial point
    xi == x0 && return xi
    isnan(xi) && return nan

    xold = x0
    while itercount < max_iters
    #("iter = $itercount")

    xi = f(xi)
    #@show xi
    xi == xold && return xi
    if isnan(xi) && abs(xi) == T(Inf)
        if antinan
            return xold
        else
            return nan
        end
    end

    Δx = abs(xi-xold)
    #@show Δx/xi,xi
    if abs(Δx) < max(atol,abs(xi)*rtol)
        return xi
    end

    itercount +=1
    xold = xi
end
    if antinan
        return xold
    else
        return nan
    end

end


text:
#=
Function overloads and custom types for interacting with Mapper2.

Specifically:

    - Defines attributes in the architectural model and how that relates to
        mappability.

    - Custom placement TwoChannel edge that contains a "cost" field for
        annotating important links. Also provides a custom "channel_cost" method
        for this new type of link.

    - Custom routing channel with a cost field.
=#

################################################################################
# Custom architectural definitions for Mapper2
################################################################################

# Custom RuleSet for KiloCore style mappings
abstract type AbstractKC <: RuleSet end

struct KC <: AbstractKC end
struct Asap2 <: AbstractKC end

# Shortcuts for convenience
const TN = TaskgraphNode
const TE = TaskgraphEdge

# Two tasks are considered equivalent if they have the exact same attribute
# metadata. NOTE: These are just strings rather than vectors of strings like
# with components.
function Mapper2.isequivalent(::AbstractKC, a::TN, b::TN)
    return a.metadata[typekey()] == b.metadata[typekey()]
end

# A Component is mappable if it has a right "typekey" in the metadata and
# the corresponding value for that key has length greater than 1. In general,
# the value should be a vector of strings, so this last criteria takes care
# of the instance where a core gets assigned an empty metadata string for some reason.
function Mapper2.ismappable(::AbstractKC, c::Component)
    return haskey(c.metadata, typekey()) && length(c.metadata[typekey()]) > 0
end

# For a task to be able to map to a component, the attribute of the task must
# be in the list of attributes supplied by the component.
function Mapper2.canmap(::AbstractKC, t::TN, c::Component)
    ismappable(c) || return false
    return in(t.metadata[typekey()], c.metadata[typekey()])
end

#####
##### Extensions for routing
#####
function Mapper2.is_source_port(::AbstractKC, p::Port, e::TE)
    port_link_class = p.metadata["link_class"]
    edge_link_class = e.metadata["link_class"]

    return port_link_class == edge_link_class
end

function Mapper2.is_sink_port(::AbstractKC, p::Port, e::TE)
    port_link_class = p.metadata["link_class"]
    edge_link_class = e.metadata["link_class"]

    # Check if this is a circuit_link. If so, preserve the destination index.
    if e.metadata["preserve_dest"] && port_link_class == edge_link_class
        dest_index = get(e.metadata, "dest_index", nothing)
        dest_index == nothing && return true
        return e.metadata["dest_index"] == p.metadata["index"]
    end

    return port_link_class == edge_link_class
end

# Don't perform the "preserve_dest" check
Mapper2.is_sink_port(a::Asap2, p::Port, e::TE) = Mapper2.is_source_port(a,p,e)

function Mapper2.needsrouting(::AbstractKC, edge::TaskgraphEdge)
    return edge.metadata["route_link"]
end

################################################################################
# Placement
################################################################################

# All edges for KC types will be "CostChannel". Even if the profiled_links option
# is not being used, links between memories and memory_processors will still
# be given a higher weight to encourage them to be mapped together.
#
# The main difference here is the "cost" field, which will be multiplied by
# the distance between nodes to emphasize some connections over others.
struct CostChannel <: Mapper2.SA.TwoChannel
    source::Int64
    sink::Int64
    cost::Float64
end

# Extend the "channel_cost" function for KC types.
Base.@propagate_inbounds function Mapper2.SA.channel_cost(sa::SAStruct, channel::CostChannel)
    src = sa.nodes[channel.source]
    dst = sa.nodes[channel.sink]
    distance = getdistance(sa.distance, src, dst)

    return channel.cost * distance
end

# # Constructor for CostChannel. Extracts the "cost" field from the metadata
# # of each taskgraph channel type.
function Mapper2.SA.buildchannel(::AbstractKC, channel, sources, sinks)
    # Quick verification that no fanout is happening. This should never happen for normal
    # KiloCore mappings.
    @assert length(sources) == 1
    @assert length(sinks) == 1

    # Since all the source and sink vectors are of length 1, we can get the
    # source and sink simply by taking the first element.
    source = first(sources)
    sink = first(sinks)
    cost = channel.metadata["cost"]
    return CostChannel(source, sink, cost)
end

########################
# Extensions for Asap2 #
########################

build_aux(map::Map{D,Asap2}) where {D} =
    MutableBinaryMaxHeap(zeros(UInt8, num_edges(map.taskgraph)))

# Asap2 has the quirk where longer circuit switched links run slower
#
# The idea is to keep a record of all the links in the mapping in a binary maxheap, and
# attach this maxheap to the auxiliary slot in the SAStruct. We then apply a global penalty
# on the maximum link length, promoting a lower global link length.
Base.@propagate_inbounds function SA.unsafe_assign(
            sa :: SAStruct{Asap2},
            node :: SA.SANode,
            new_location
        )

    SA.assign(node, new_location)
    # Iterate through input and output channels, update their location in the
    # edge heap.
    for edge in node.outchannels
        # Get the destination node.
        @inbounds dest = sa.nodes[sa.channels[edge].sink]

        # The distance between the two nodes.
        distance = SA.getdistance(sa.distance, new_location, SA.location(dest))

        # Update the heap entry for this edge with the new distance
        update!(sa.aux, edge, distance)
    end
    for edge in node.inchannels
        @inbounds src = sa.nodes[sa.channels[edge].source]
        # The distance between the two nodes.
        distance = SA.getdistance(sa.distance, SA.location(src), new_location)
        update!(sa.aux, edge, distance)
    end
    return nothing
end

# Apply a high global penalty for the longest link length
Mapper2.SA.aux_cost(sa_struct::SAStruct{Asap2}) = 512.0 * first(sa_struct.aux)

################################################################################
# Routing
################################################################################

# Custom Channels
struct CostRoutingChannel <: RoutingChannel
    start_vertices   ::Vector{PortVertices}
    stop_vertices    ::Vector{PortVertices}
    cost             ::Float64
end

# Override this to allow more active channels to be routed first.
Base.isless(a::CostRoutingChannel, b::CostRoutingChannel) = a.cost < b.cost

function Mapper2.routing_channel(::AbstractKC, start, stop, edge)
    cost = edge.metadata["cost"]
    return CostRoutingChannel(start, stop, cost)
end


text:
using AAABK2018, Test
λ =  1.131874350361868; ψ =  0.037175089899834; ν =  0.205622086535123
α =  0.925520992611854; ϕ =  0.216311889844990; θˡ=  1.390807772162332
θʰ=  1.751156379867318; θᴱ=  0.023600456924424; ε =  2.900000000000000
ρ =  0.020000000000000; γ =  0.500000000000000; γᴱ=  0.500000000000000
σ =  2.000000000000000; Lˢ=  0.165500000000000; ω =  1.0

p = Params(λ,ψ,ν, α, ϕ, θˡ, θʰ, θᴱ, ε, ρ, γ, γᴱ, σ, Lˢ, ω)

# initial guess for equilibirum objects
eqInit = [ 1.3, 0.3, 0.1,
           0.5, 0.6, 1.5]

# solve
eq,res = solveBGP(p,eqInit)

@test res.f_converged | res.x_converged

# solve with 5% incumbent R&D subsidy
pIncSub = Params(p, (sⁱ = 0.05,))
eqIncSub,resIncSub = solveBGP(pIncSub, eqInit)

@test resIncSub.f_converged | resIncSub.x_converged

# solve with 5% operation cost subsidy
pFixSub = Params(p, (sᶠ = 0.05,))
eqFixSub,resFixSub = solveBGP(pFixSub,eqInit)

@test resFixSub.f_converged | resFixSub.x_converged

# solve with 5% entrant R&D subsidy
pEntSub = Params(p, (sᴱ= 0.05,))
eqEntSub,resEntSub = solveBGP(pEntSub,eqInit)

@test resEntSub.f_converged | resEntSub.x_converged

@test all(eqIncSub.x .> eq.x)
@test all(eqFixSub.x .< eq.x)
@test all(eqEntSub.x .< eq.x)

# social planner problem
eqSoc,resSoc = solveSocialPlanner(p)
@test resSoc[:ret] == :FTOL_REACHED
@test abs(resSoc[:cev]/1.0447 - 1.0) < 1.0e-3    # compare with results from the paper
@test abs(eqSoc.x[1]/0.2542 - 1.0)   < 1.0e-3    # compare with results from the paper
@test abs(eqSoc.x[2]/0.4534 - 1.0)   < 1.0e-3    # compare with results from the paper
@test abs(eqSoc.qmin[1]/2.4042 - 1.0)< 1.0e-1    # compare with results from the paper

# optimal policy
# one-tool policies
optPol = policy_opt(p, sⁱ = true)
@test optPol[:ret] == :FTOL_REACHED
@test abs(optPol[:sⁱ] - 0.39) < .01    # compare with results from the paper
@test abs(optPol[:cev]/1.0122 - 1.0) < 1.0e-4    # compare with results from the paper

optPol = policy_opt(p, sᶠ = true)
@test optPol[:ret] == :FTOL_REACHED
@test abs(optPol[:sᶠ] + 0.69) < .01    # compare with results from the paper
@test abs(optPol[:cev]/1.0142 - 1.0) < 1.0e-4    # compare with results from the paper

optPol = policy_opt(p, sᴱ = true)
@test optPol[:ret] == :FTOL_REACHED
@test abs(optPol[:sᴱ] - 0.18) < .01    # compare with results from the paper
@test abs(optPol[:cev]/1.0004 - 1.0) < 1.0e-4    # compare with results from the paper

# two-tool policy
optPol = policy_opt(p, sⁱ = true, sᶠ = true)
@test optPol[:ret] == :FTOL_REACHED
@test abs(optPol[:sⁱ] + 0.03) < .01    # compare with results from the paper
@test abs(optPol[:sᶠ] + 0.74) < .01    # compare with results from the paper
@test abs(optPol[:cev]/1.0142 - 1.0) < 1.0e-4    # compare with results from the paper


text:
include("data_round2.jl")
using Plots

summary_stats("data/full_noincentive.csv")
summary_stats("data/full_bonus50.csv")
summary_stats("data/full_bonus100.csv")

g1n = calculate_gradients("data/step_from_naive.csv", [2.87701, 0.68888], 0.20)
g0n = get_b0_gradient("data/step_from_naive.csv", [2.87701, 0.68888], 0.20)
g1 = calculate_gradients("data/step_from_rrm.csv", [2.83594, 0.547041], 0.10)
g0 = get_b0_gradient("data/step_from_rrm.csv", [2.83594, 0.547041], 0.10)

r2_naive_insample = o_mse_calc(["data/full_noincentive.csv"], [2.87701, 0.68888])
noincentive = ["data/full_noincentive_out.csv"]
r2_naive = o_mse_calc(noincentive, [2.87701, 0.68888])
incentive = ["data/full_bonus100.csv"]
r2_naive_incentive = o_mse_calc(incentive, [2.87701, 0.68888])

b1 = 0.60
b0 = 2.89

r2_ic_ic = o_mse_calc(["data/step_out.csv"], [b0, b1])

plot(["Naive RM (IS)", "Naive RM, No Bonus", "Naive RM, Bonus", "OT-DE, Bonus"],
      [ r2_naive_insample, r2_naive, r2_naive_incentive, r2_ic_ic], seriestype=:bar,
        label=nothing, ylabel="R-Squared")

max_stars_repo_name: evanmunro/personalized-policy

max_issues_repo_name: evanmunro/personalized-policy

max_forks_repo_name: evanmunro/personalized-policy


text:
using ADCME
using PyCall
using LinearAlgebra
using PyPlot
using Random
using AdFem
using SparseArrays
Random.seed!(233)

function fast_solve(rhs)
    fast_solve_ = load_op_and_grad("./build/libFastSolve","fast_solve")
    rhs = convert_to_tensor(rhs, dtype = Float64)
    fast_solve_(rhs)
end

m = 10
n = 5
h = 0.1

K_ = compute_fem_stiffness_matrix1(diagm(0=>ones(2)), m, n, h)
ind = collect(1:m+1)
K_[ind, :] .= 0.0
K_[ind, ind] = spdiagm(0=>ones(length(ind)))
rhs = Float64.(collect((1:(m+1)*(n+1))))
x = K_\rhs
ccall((:factorize_matrix, "./build/libSaver"), Cvoid, (Cint, Cint, Cdouble, Ref{Cint}, Cint),
            Int32(m), Int32(n), h, Int32.(ind .- 1) , Int32(length(ind)))

tf_x = fast_solve(constant(rhs))
sess = Session(); init(sess)
@show norm(run(sess, tf_x)-x)
# uncomment it for testing gradients
# error()

# TODO: change your test parameter to `m`
#       in the case of `multiple=true`, you also need to specify which component you are testings
# gradient check -- v
function scalar_function(m)
    return sum(fast_solve(m)^2)
end

# TODO: change `m_` and `v_` to appropriate values
m_ = constant(rand((m+1)*(n+1)))
v_ = rand((m+1)*(n+1))
y_ = scalar_function(m_)
dy_ = gradients(y_, m_)
ms_ = Array{Any}(undef, 5)
ys_ = Array{Any}(undef, 5)
s_ = Array{Any}(undef, 5)
w_ = Array{Any}(undef, 5)
gs_ =  @. 1 / 10^(1:5)

for i = 1:5
    g_ = gs_[i]
    ms_[i] = m_ + g_*v_
    ys_[i] = scalar_function(ms_[i])
    s_[i] = ys_[i] - y_
    w_[i] = s_[i] - g_*sum(v_.*dy_)
end

sess = Session(); init(sess)
sval_ = run(sess, s_)
wval_ = run(sess, w_)
close("all")
loglog(gs_, abs.(sval_), "*-", label="finite difference")
loglog(gs_, abs.(wval_), "+-", label="automatic differentiation")
loglog(gs_, gs_.^2 * 0.5*abs(wval_[1])/gs_[1]^2, "--",label="\$\\mathcal{O}(\\gamma^2)\$")
loglog(gs_, gs_ * 0.5*abs(sval_[1])/gs_[1], "--",label="\$\\mathcal{O}(\\gamma)\$")

plt.gca().invert_xaxis()
legend()
xlabel("\$\\gamma\$")
ylabel("Error")

max_stars_repo_name: ADCMEMarket/ADCMEImages

max_issues_repo_name: ADCMEMarket/ADCMEImages

max_forks_repo_name: ADCMEMarket/ADCMEImages


text:
using Stan
wd = dirname(@__FILE__)

if !isdefined(Main, :update_R_file)
  function update_R_file{T<:Any}(file::String, dct::Dict{String, T};
      replaceNaNs::Bool=true)
  	isfile(file) && rm(file)
  	strmout = open(file, "w")

  	str = ""
  	for entry in dct
  		str = "\""entry[1]"\""" <- "
  		val = entry[2]
  		if length(val)==1 && length(size(val))==0
  			# Scalar
  			str = str*"$(val)\n"
      elseif length(val)==1 && length(size(val))==1
  			# Single element vector
  			#str = str*"$(val[1])\n"
  		elseif length(val)>1 && length(size(val))==1
  			# Vector
  			str = str*"structure(c("
  			write(strmout, str)
  			str = ""
  			writecsv(strmout, val');
  			str = str*"), .Dim=c($(length(val))))\n"
  		elseif length(val)>1 && length(size(val))>1
  			# Array
  			str = str*"structure(c("
  			write(strmout, str)
  			str = ""
  			writecsv(strmout, val[:]');
  			dimstr = "c"*string(size(val))
  			str = str*"), .Dim=$(dimstr))\n"
  		end
  		write(strmout, str)
  	end
  	close(strmout)
  end
end

dct = Dict{String, Any}[
  Dict("theta" => 0.1),
  Dict("theta" => [0.4]),
  Dict("theta" => [0.1, 0.2, 0.4]),
  Dict("theta" => [0.4 0.5; 0.6 0.7]),
]

for i in 1:length(dct)
  #for j in ["data", "init"]
  for j in ["data"]
    file = joinpath(wd, "test_r_$(i)_$(j).R")
    isfile(file) && rm(file)
    open(file, "w")
    Stan.update_R_file(file, dct[i])
  end
end