text:
function get_ode_integrator(func, u0, dt=nothing, param=nothing, alg=Tsit5(), maxiters=Inf;
    reltol=1e-12, abstol=1e-12)
    """
    Args:
    - func : the update function
    - u0  : Inital vector of value
    - param   : additional param
    - dt  : Time step for integration
    - alg : update algorithm (e.g. Euler(), RK4())
    - maxiters : maxiters of integration. Default is Inf.
    """

    tspan = (0.0, Inf)

    if param == nothing
        prob = ODEProblem(func, u0, tspan)
    else
        prob = ODEProblem(func, u0, tspan, param)
    end

    if dt == nothing
        integrator = init(prob, alg, adaptive=true, maxiters=maxiters)
    else
        integrator = init(prob, alg, dt=dt, adaptive=false, maxiters=maxiters)
    end
    integrator.opts.reltol = reltol
    integrator.opts.abstol = abstol
    return integrator
end

function get_ode_solution(func, u0, tspan, dt=nothing, param=nothing, alg=Tsit5();
    reltol=1e-12, abstol=1e-12, save_everystep=true)
    """
    Args:
    - f   : the update function
    - u0  : Inital vector of value
    - param : additional param
    - dt  : Time step for integration (when adaptive, it is initial dt)
    - alg : update algorithm (e.g. Euler(), RK4())
    """

    if param == nothing
        prob = ODEProblem(func, u0, tspan)
    else
        prob = ODEProblem(func, u0, tspan, param)
    end

    if dt == nothing
        sol = solve(prob, alg, reltol=reltol, abstol=abstol, adaptive=true, save_everystep=save_everystep)
    else
        sol = solve(prob, alg, dt=dt, reltol=reltol, abstol=abstol, adaptive=false, save_everystep=save_everystep)
    end
    return sol
end


text:
using ConcreteInterfaces.ComputedFieldTypes
using Base.Test

# Verify that @computed works on types with no computed fields

@computed struct NoComputed
    x::Int
end

@test typeof(@inferred(NoComputed(1))) === NoComputed
@test fulltype(NoComputed) === NoComputed

@computed struct NoComputedInner
    x::Int
    NoComputedInner(y) = new(y + 1)
end

@test typeof(@inferred(NoComputedInner(2))) === NoComputedInner
@test fulltype(NoComputedInner) === NoComputedInner
@test NoComputedInner(5).x == 6

include("readme_examples.jl")
include("inheritance.jl")

ext: jl


text:
#=
Copyright (c) 2015, Intel Corporation

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

    * Redistributions of source code must retain the above copyright notice,
      this list of conditions and the following disclaimer.
    * Redistributions in binary form must reproduce the above copyright
      notice, this list of conditions and the following disclaimer in the
      documentation and/or other materials provided with the distribution.
    * Neither the name of Intel Corporation nor the names of its contributors
      may be used to endorse or promote products derived from this software
      without specific prior written permission.

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

# $ julia knob-test1.jl ../matrices/bcsstk14.mtx
# assymmetric (2, 8) exists but (8, 2) doesn't
# julia: SymGS.cpp:305: bool SpMP::getSymmetricNnzPattern(const SpMP::CSR *, int **, int **, int **, int **): Assertion `sym.isSymmetric(false, true)' failed.
# Aborted (core dumped)

include("../../src/Sparso.jl")
include("../../src/simple-show.jl")
include("./utils.jl")
using Sparso

function pcg_symgs(x, A, b)
    L, U = Sparso.ilu(A)
    r = b - Sparso.SpMV(1,A,x)
    z = copy(r)
    Sparso.fwdTriSolve!(z,L,r, C_NULL)
end

A = matrix_market_read(ARGS[1], true, true)
m = size(A, 1)
b = ones(Float64, m)
x = zeros(Float64, m)
pcg_symgs(x, A, b)

max_stars_repo_name: IntelLabs/Sparso

max_issues_repo_name: IntelLabs/Sparso

max_forks_repo_name: IntelLabs/Sparso


text:
module OS

export getfqdn
"""
Result of readchomp(`hostname -f`)
"""
getfqdn() = readchomp(`hostname -f`)

export getusername
"""
readchomp(`id -u -n`)
"""
getusername() = readchomp(`id -u -n`)

end

# =============== Pid ====================

module Pid

import OS

debinfofun(on::Bool) = on ? (x...)->info(x...) : (x...)->nothing

statinfofun(on::Bool) = on ? (x...)->info(x...) : (x...)->nothing

export isworks_local
"""
isworks_local(pid::AbstractString; deb::Bool=false) -> Bool
"""
function isworks_local(pid::AbstractString; deb::Bool=false)
    debinfo::Function = debinfofun(deb)
    try cmd = `bash -c "
        ( ps -A | grep -P '\b$pid\b' ) >/dev/null && echo YES || echo NO"`
        cmdrv = readchomp(cmd)
        debinfo("""isworks_local($pid): Command '$cmd'
        - returns '$cmdrv' """)
        cmdrv == "YES"
    catch e
        warn("isworks_local(\"$pid\"): ", e)
        return false
    end

end

export isworks_remote
"""
isworks_remote(pid::AbstractString,
                    host::AbstractString;
                                        deb::Bool=false) -> Bool
"""
function isworks_remote(pid::AbstractString,
            host::AbstractString;
            deb::Bool=false)
    debinfo::Function = debinfofun(deb)
    cmd=`ssh $host "
        ( ps -A | grep $pid ) >/dev/null && echo YES || echo NO"`
    try cmdrv = readchomp(cmd)
        debinfo("""isworks_remote($pid,$host): Command '$cmd'
        -  returns '$cmdrv' """)
        cmdrv == "YES"
    catch e
        warn("isworks_remote($pid,$host): ", e)
        return false
    end

end

export is_local_host
"""
is_local_host(host::AbstractString) -> Bool
host must be a FQDN, like a result of \`hostname -f\` bash command.
"""
is_local_host{S<:AbstractString}(host::S) = host==OS.getfqdn()

export isworks
"""
isworks{S<:AbstractString}(pid::S, host::S; deb::Bool=false) -> Bool
"""
function isworks{S<:AbstractString}( pid::S, host::S; deb::Bool=false )

    if !ismatch(r"^\d+$",pid)
    error("pid can contain only digits (given pid='$pid')")
    end
    debinfo::Function = debinfofun(deb)
    if is_local_host(host)
    debinfo("""isworks:(\"$pid\",$host): - my own host.""")
    return isworks_local(pid, deb=deb)
    else
    debinfo("""isworks(\"$pid\",$host): - remote host.""")
    return isworks_remote(pid, host, deb=deb)
    end
end

end # module

# ======================== Flag ==========================

module Flag

import OS
import Pid

debinfofun(on::Bool) = on ? (x...)->info(x...) : (x...)->nothing

statinfofun(on::Bool) = on ? (x...)->info(x...) : (x...)->nothing

export name
"""
name(\"aa.txt\")      -> \"aa.txt.FLAG\"
name(\"aa.txt.FLAG\") -> \"aa.txt.FLAG\"
"""
function name(f::AbstractString)
 if !ismatch(r"\.FLAG$", f)
    f = string(f,".FLAG")
 end
 f
end

export read
"read(filename) -> Dict"
function read(f::AbstractString)
    f = name(f)
    di=Dict{AbstractString,AbstractString}()
    if isfile(f)
    content = open(readchomp,f)
    if isempty(content)
        return di
    end
    kva::Array = split(content,';')

    for kv in kva
        kv = split(kv,'=')
        if length(kv)<2
        continue
        end
        get!(di,kv[1],kv[2])
    end
    end
    di
end

export isworks
"""
isworks(\"aa.txt\", deb=false)
or the same: isworks(\"aa.txt.FLAG\", deb=false)
"""
function isworks(f, deb::Bool=false, warnings=false)
 debinfo = debinfofun(deb)
 fl = name(f)
 if isfile(fl)
    debinfo("set: File $fl - is exists.")
    if filesize(fl)>1000
    error("Size of flag file $fl > 1000 bytes")
    end
    flag::Dict = Flag.read(fl)
    debinfo("Flag file $fl content: $flag")
    pid = get(flag,"PID","")
    if isempty(pid)
        warnings && warn("empty 'pid' key in flag $fl, unset flag...")
        Flag.unset(fl)
        return false
    end
    if Pid.isworks(pid, flag["HOSTFQDN"], deb=deb)
    return true
    end
 end
 false
end

export newcontent
"Generates and returns Dict with new generated content for flag file."
newcontent()::Dict = Dict("PID"=>getpid(),
            "HOSTFQDN"=>OS.getfqdn(),
            "USER"=>OS.getusername(),
            "STARTED"=>now(),
            "UUID"=>Base.Random.uuid4() )

export set
"""
set(filename::AbstractString; deb::Bool=false) -> Bool
"""
set(f::AbstractString; deb::Bool=false) = set(f, newcontent(), deb=deb)

function set(f::AbstractString, d::Dict; deb::Bool=false)
 debinfo::Function = debinfofun(deb)
 fl = name(f)
 if isworks(fl)
    return false
 end
 open(fl,"w") do io
    str = join(["$k=$v" for (k,v) in d], ";")
    debinfo("Write to $fl: \"$str\"")
    print(io,str)
 end
 true
end

export unset
"""
unset(filename::AbstractString) -> Bool
"""
function unset(f::AbstractString)
 f = name(f)
 if isfile(f)
    try rm(f)
    catch e
    warn("unset: can't unset $f:", e)
    return false
    end
    return true
 else
    false
 end
end

export via
"""
Set flag name, call func() and then unset flag.

via( func::Function, flagfile::String, extrargs... )

via( flagfile::String, extrargs... ) do name, extrargs...
    dosmf(...)
end
"""
function via(fun::Function, flag::AbstractString, extrargs...)
 if !Flag.set(flag)
    return false
 end
 rv = fun(flag, extrargs...)
 Flag.unset(flag)
 rv
end

end # module

# ============================= File ======================

module File

import Flag

export goodsize
"goodsize( filename ) -> Bool"
goodsize(f::AbstractString) = goodsize(f, Dict(".gz"=>20))

"""goodsize( filename, dict ) -> Bool
Define your own dict::Dict{String,Int} to redefine compared filesises.
"""
function goodsize{Str<:AbstractString}(
            f::AbstractString, dict::Dict{Str,Int})
    ext = splitext(f)[2];
    compare_with = get(dict, lowercase(ext), 0)
    filesize(f) > compare_with
end

export iswait
"iswait( filename ) -> Bool"
iswait( filename::AbstractString ) = Flag.isworks( filename )

export istodo
"istodo( filename ) -> Bool"
istodo{S<:AbstractString}(f::S) =
    !goodsize(f) && !iswait(f)

"istodo( filename, dict ) -> Bool"
istodo{S<:AbstractString}(f::S, dict::Dict{S,Int}) =
    !goodsize(f, dict) && !iswait(f)

export tmpname
"""
    tmpname(\"aa.txt\")      -> \"aa.txt.TMP\"
    tmpname(\"aa.txt.TMP\") -> \"aa.txt.TMP\"
"""
function tmpname(f::AbstractString)
 if !ismatch(r"\.TMP$", f)
    f = string(f,".TMP")
 end
 f
end

export files

"""
    files(dir)

    Recursive seach all files in directory dir.
"""
function files(dir::AbstractString)
 f = x->true
 files(dir, @inline(f))
end

"""
    files(dir, filter::Function)

    Recursive search files in directory dir, filtered by filter.
"""
files( dir::AbstractString, only::Function ) =
    @task _files(true, dir, only)

"""
    files(dir, re::Regex, args...)

    Recursive seach files in directory dir where files matched by re.
"""
files( dir::AbstractString, re::Regex, args... ) =
    files(dir, x->ismatch(re,x), args... )

"files(dir, filter::Function, re::Regex, args...)"
files( dir::AbstractString, f1::Function, re::Regex, args...) =
    files(dir, x->f1(x)&&ismatch(re,x), args... )

"files(dir, flt1::Function, flt2::Function, ff::Function...)"
files(dir::AbstractString, f1::Function, f2::Function, ff::Function...) =
 files(dir, x->f1(x)&&f2(x), ff...)

function _files(::Bool, dir::AbstractString, only::Function )
  if !isdir(dir) error("$dir not a directory") end
  for x in readdir(dir)
    pathtox = joinpath(dir,x)
    if only(pathtox)
    produce(pathtox)
    end
    if isdir(pathtox)
        _files(true, pathtox, only)
    end
  end
end

#---------------------- hasnear ------
export hasnear
"""
    Returns true if any existing name
    in the dir of <file> are matches all of [re]
    and optional <test> returns true.

    hasnear(\"a.txt\", [r\".*txt\"])
    hasnear(\"a.txt\", r\".*txt\", f->isfile(f))
"""
function hasnear( file::AbstractString,
                rr::Array{Regex,1},
                tst::Function=goodsize )
    (dir,basename)=splitdir(file)
    if isempty(dir)
    dir="."
    end
    if isempty(basename)
        warn("has_near(\"$file\",$rr): Empty basename in \"$file\"")
        return false
    end
    basenames = readdir(dir)
    for r in rr
    if any( bn->ismatch(r,bn)&&tst(joinpath(dir,bn)) ,basenames )
        return true
    end
    end
    false
end

"""
    Returns true if any existing name
    in the dir of <file> are matches with <re>.

    hasnear( file::String, base_name::Regex, [test::Function] )

    hasnear(\"a.txt\", r\".*txt\")
    hasnear(\"a.txt\", r\".*txt\", f->isfile(f))
    \"a.txt\" |> hasnear(...)

"""
hasnear( file::AbstractString, r::Regex ) =
                hasnear(file, [r])

hasnear( file::AbstractString, r::Regex, tst::Function ) =
                hasnear(file, [r], tst)

"hasnear(basename::Regex) -> file->hasnear(file,[basename])"
hasnear( r::Regex ) = f->hasnear(f,[r])

"hasnear(bn::Regex, tst::Func) -> file->hasnear(file,[bn],tst)"
hasnear( r::Regex, tst::Function ) = f->hasnear(f,[r],tst)

"hasnear(nn::Array{Regex,1}) -> file->hasnear(file,[nn])"
hasnear( rr::Array{Regex,1} ) = f->hasnear(f,rr)

"hasnear(nn::Array{Regex,1}, t::Func) -> file->hasnear(file,[nn],t)"
hasnear( rr::Array{Regex,1}, t::Function ) = f->hasnear(f,rr,t)

# ---------- with string: -----------
"""
    hasnear( file::String, bname::String, [tst::Function] )

    Is there <basename> near the <file> and tst(fullname) returns true ?
"""
function hasnear( file::AbstractString,
           bname::AbstractString,
           tst::Function=goodsize )

    (origdir,origbname) = splitdir(file)
    if bname==origbname
    warn("""hasnear(\"$file\",\"$bname\",...): new basename (\"$bname\") equal old base name.""")
    end

    tst(joinpath(origdir,bname))
end

"""
    hasnear( file::String, change::Function, [tst::Function=goodsize] )

    Gets newname via change(basename) of <file>
    and returns result of tst(newname).
"""
function hasnear( file::AbstractString,
           change::Function,
           tst::Function=goodsize )

 (origdir,origbname) = splitdir(file)
 newbname = change(origbname)
 if newbname==origbname
    warn("""hasnear(\"$file\",...): new basename (\"$newbname\") has not difference whith old base name.""")
 end
 tst( joinpath(origdir,newbname))
end

end # module

max_stars_repo_name: closescreen/File

max_issues_repo_name: closescreen/File

max_forks_repo_name: closescreen/File


text:
module QuasiMonteCarlo

using Reexport
@reexport using Transducers
using Sobol
using StatsFuns

include("quasi.jl")
include("normals.jl")
include("mean.jl")

end # module


text:
int_rules_1_2_4_3 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 1.2.4.3*(d+e*x^(n-q))*(a*x^q+b*x^n+c*x^(2*n-q))^p =#
    @apply_utils Antiderivative((~A + ~(B') * (~x) ^ ~(r')) * (~(a') * (~x) ^ ~(q') + ~(b') * (~x) ^ ~(n') + ~(c') * (~x) ^ ~(j')) ^ ~(p'), ~x) => Antiderivative((~x) ^ (~p * ~q) * (~A + ~B * (~x) ^ (~n - ~q)) * (~a + ~b * (~x) ^ (~n - ~q) + ~c * (~x) ^ (2 * (~n - ~q))) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~A, ~B, ~n, ~q], ~x) && (EqQ(~r, ~n - ~q) && (EqQ(~j, 2 * ~n - ~q) && (IntegerQ(~p) && PosQ(~n - ~q))))
    #= Antiderivative(((~A)+(~B')*(~x)^(~j'))*((~a')*(~x)^(~q')+(~b')*(~x)^(~n')+(~c')*(~x)^(~r'))^(~p),~x)  := sqrt(a*x^q+b*x^n+c*x^(2*n-q))/(x^(q/2)*sqrt(a+b*x^(n-q)+c*x^(2*( n-q))))* Antiderivative(x^(q*p)*(A+B*x^(n-q))*(a+b*x^(n-q)+c*x^(2*(n-q)))^p,x) <--  FreeQ([a,b,c,A,B,n,p,q],x) && EqQ(j,n-q) && EqQ(r,2*n-q) && PosQ(n-q)  && IGtQ(p+1/2,0) =#
    #= Antiderivative(((~A)+(~B')*(~x)^(~j'))*((~a')*(~x)^(~q')+(~b')*(~x)^(~n')+(~c')*(~x)^(~r'))^(~p),~x)  := x^(q/2)*sqrt(a+b*x^(n-q)+c*x^(2*(n-q)))/sqrt(a*x^q+b*x^n+c*x^(2* n-q))* Antiderivative(x^(q*p)*(A+B*x^(n-q))*(a+b*x^(n-q)+c*x^(2*(n-q)))^p,x) <--  FreeQ([a,b,c,A,B,n,p,q],x) && EqQ(j,n-q) && EqQ(r,2*n-q) && PosQ(n-q)  && ILtQ(p-1/2,0) =#
    #= Antiderivative(((~A)+(~B')*(~x)^(~j'))*sqrt((~a')*(~x)^(~q')+(~b')*(~x)^(~n')+(~c')*(~x)^(~r')),(~x)*Symbol) :=  sqrt(a*x^q+b*x^n+c*x^(2*n-q))/(x^(q/2)*sqrt(a+b*x^(n-q)+c*x^(2*(n-q))) )* Antiderivative(x^(q/2)*(A+B*x^(n-q))*sqrt(a+b*x^(n-q)+c*x^(2*(n-q))),x) <--  FreeQ([a,b,c,A,B,n,q],x) && EqQ(j,n-q) && EqQ(r,2*n-q) && PosQ(n-q) =#
    @apply_utils Antiderivative((~A + ~(B') * (~x) ^ ~(j')) / sqrt(~(a') * (~x) ^ ~(q') + ~(b') * (~x) ^ ~(n') + ~(c') * (~x) ^ ~(r')), ~x) => (((~x) ^ (~q / 2) * sqrt(~a + ~b * (~x) ^ (~n - ~q) + ~c * (~x) ^ (2 * (~n - ~q)))) / sqrt(~a * (~x) ^ ~q + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n - ~q))) * Antiderivative((~A + ~B * (~x) ^ (~n - ~q)) / ((~x) ^ (~q / 2) * sqrt(~a + ~b * (~x) ^ (~n - ~q) + ~c * (~x) ^ (2 * (~n - ~q)))), ~x) <-- FreeQ([~a, ~b, ~c, ~A, ~B, ~n, ~q], ~x) && (EqQ(~j, ~n - ~q) && (EqQ(~r, 2 * ~n - ~q) && (PosQ(~n - ~q) && (EqQ(~n, 3) && EqQ(~q, 2)))))
    @apply_utils Antiderivative((~A + ~(B') * (~x) ^ ~(r')) * (~(a') * (~x) ^ ~(q') + ~(b') * (~x) ^ ~(n') + ~(c') * (~x) ^ ~(j')) ^ ~p, ~x) => (~x * (~b * ~B * (~n - ~q) * ~p + ~A * ~c * (~p * ~q + (~n - ~q) * (2 * ~p + 1) + 1) + ~B * ~c * (~p * (2 * ~n - ~q) + 1) * (~x) ^ (~n - ~q)) * (~a * (~x) ^ ~q + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n - ~q)) ^ ~p) / (~c * (~p * (2 * ~n - ~q) + 1) * (~p * ~q + (~n - ~q) * (2 * ~p + 1) + 1)) + (((~n - ~q) * ~p) / (~c * (~p * (2 * ~n - ~q) + 1) * (~p * ~q + (~n - ~q) * (2 * ~p + 1) + 1))) * Antiderivative((~x) ^ ~q * ((2 * ~a * ~A * ~c * (~p * ~q + (~n - ~q) * (2 * ~p + 1) + 1) - ~a * ~b * ~B * (~p * ~q + 1)) + ((2 * ~a * ~B * ~c * (~p * (2 * ~n - ~q) + 1) + ~A * ~b * ~c * (~p * ~q + (~n - ~q) * (2 * ~p + 1) + 1)) - (~b) ^ 2 * ~B * (~p * ~q + (~n - ~q) * ~p + 1)) * (~x) ^ (~n - ~q)) * (~a * (~x) ^ ~q + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n - ~q)) ^ (~p - 1), ~x) <-- FreeQ([~a, ~b, ~c, ~A, ~B, ~n, ~q], ~x) && (EqQ(~r, ~n - ~q) && (EqQ(~j, 2 * ~n - ~q) && (Not(IntegerQ(~p)) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (GtQ(~p, 0) && (NeQ(~p * (2 * ~n - ~q) + 1, 0) && NeQ(~p * ~q + (~n - ~q) * (2 * ~p + 1) + 1, 0)))))))
    @apply_utils Antiderivative((~A + ~(B') * (~x) ^ ~(r')) * (~(a') * (~x) ^ ~(q') + ~(c') * (~x) ^ ~(j')) ^ ~p, ~x) => With([n = ~q + ~r], (~x * (~A * (~p * ~q + (n - ~q) * (2 * ~p + 1) + 1) + ~B * (~p * (2n - ~q) + 1) * (~x) ^ (n - ~q)) * (~a * (~x) ^ ~q + ~c * (~x) ^ (2n - ~q)) ^ ~p) / ((~p * (2n - ~q) + 1) * (~p * ~q + (n - ~q) * (2 * ~p + 1) + 1)) + (((n - ~q) * ~p) / ((~p * (2n - ~q) + 1) * (~p * ~q + (n - ~q) * (2 * ~p + 1) + 1))) * Antiderivative((~x) ^ ~q * (2 * ~a * ~A * (~p * ~q + (n - ~q) * (2 * ~p + 1) + 1) + (2 * ~a * ~B * (~p * (2n - ~q) + 1)) * (~x) ^ (n - ~q)) * (~a * (~x) ^ ~q + ~c * (~x) ^ (2n - ~q)) ^ (~p - 1), ~x) <-- EqQ(~j, 2n - ~q) && (NeQ(~p * (2n - ~q) + 1, 0) && NeQ(~p * ~q + (n - ~q) * (2 * ~p + 1) + 1, 0))) <-- FreeQ([~a, ~c, ~A, ~B, ~q], ~x) && (Not(IntegerQ(~p)) && GtQ(~p, 0))
    @apply_utils Antiderivative((~A + ~(B') * (~x) ^ ~(r')) * (~(a') * (~x) ^ ~(q') + ~(b') * (~x) ^ ~(n') + ~(c') * (~x) ^ ~(j')) ^ ~p, ~x) => (-((~x) ^ (-(~q) + 1)) * (((~A * (~b) ^ 2 - ~a * ~b * ~B) - 2 * ~a * ~A * ~c) + (~A * ~b - 2 * ~a * ~B) * ~c * (~x) ^ (~n - ~q)) * (~a * (~x) ^ ~q + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n - ~q)) ^ (~p + 1)) / (~a * (~n - ~q) * (~p + 1) * ((~b) ^ 2 - 4 * ~a * ~c)) + (1 / (~a * (~n - ~q) * (~p + 1) * ((~b) ^ 2 - 4 * ~a * ~c))) * Antiderivative((~x) ^ -(~q) * ((((~A * (~b) ^ 2 * (~p * ~q + (~n - ~q) * (~p + 1) + 1) - ~a * ~b * ~B * (~p * ~q + 1)) - 2 * ~a * ~A * ~c * (~p * ~q + 2 * (~n - ~q) * (~p + 1) + 1)) + (~p * ~q + (~n - ~q) * (2 * ~p + 3) + 1) * (~A * ~b - 2 * ~a * ~B) * ~c * (~x) ^ (~n - ~q)) * (~a * (~x) ^ ~q + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n - ~q)) ^ (~p + 1)), ~x) <-- FreeQ([~a, ~b, ~c, ~A, ~B, ~n, ~q], ~x) && (EqQ(~r, ~n - ~q) && (EqQ(~j, 2 * ~n - ~q) && (Not(IntegerQ(~p)) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && LtQ(~p, -1)))))
    @apply_utils Antiderivative((~A + ~(B') * (~x) ^ ~(r')) * (~(a') * (~x) ^ ~(q') + ~(c') * (~x) ^ ~(j')) ^ ~p, ~x) => With([n = ~q + ~r], (-((~x) ^ (-(~q) + 1)) * (~a * ~A * ~c + ~a * ~B * ~c * (~x) ^ (n - ~q)) * (~a * (~x) ^ ~q + ~c * (~x) ^ (2n - ~q)) ^ (~p + 1)) / (~a * (n - ~q) * (~p + 1) * (2 * ~a * ~c)) + (1 / (~a * (n - ~q) * (~p + 1) * (2 * ~a * ~c))) * Antiderivative((~x) ^ -(~q) * ((~a * ~A * ~c * (~p * ~q + 2 * (n - ~q) * (~p + 1) + 1) + ~a * ~B * ~c * (~p * ~q + (n - ~q) * (2 * ~p + 3) + 1) * (~x) ^ (n - ~q)) * (~a * (~x) ^ ~q + ~c * (~x) ^ (2n - ~q)) ^ (~p + 1)), ~x) <-- EqQ(~j, 2n - ~q)) <-- FreeQ([~a, ~c, ~A, ~B, ~q], ~x) && (Not(IntegerQ(~p)) && LtQ(~p, -1))
    #= Antiderivative(((~A)+(~B')*(~x)^(~q))*((~a')*(~x)^(~j')+(~b')*(~x)^(~k')+(~c')*(~x)^(~n'))^(~p),~x)  := (a*x^j+b*x^k+c*x^n)^p/(x^(j*p)*(a+b*x^(k-j)+c*x^(2*(k-j)))^p)* Antiderivative(x^(j*p)*(A+B*x^(k-j))*(a+b*x^(k-j)+c*x^(2*(k-j)))^p,x) <--  FreeQ([a,b,c,A,B,j,k,p],x) && EqQ(q,k-j) && EqQ(n,2*k-j) && PosQ(k-j)  && Not(IntegerQ(p)) =#
    @apply_utils Antiderivative((~A + ~(B') * (~x) ^ ~(j')) * (~(a') * (~x) ^ ~(q') + ~(b') * (~x) ^ ~(n') + ~(c') * (~x) ^ ~(r')) ^ ~(p'), ~x) => Unintegrable((~A + ~B * (~x) ^ (~n - ~q)) * (~a * (~x) ^ ~q + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n - ~q)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~A, ~B, ~n, ~p, ~q], ~x) && (EqQ(~j, ~n - ~q) && EqQ(~r, 2 * ~n - ~q))
    @apply_utils Antiderivative((~A + ~(B') * (~u) ^ ~(j')) * (~(a') * (~u) ^ ~(q') + ~(b') * (~u) ^ ~(n') + ~(c') * (~u) ^ ~(r')) ^ ~(p'), ~x) => (1 / Coefficient(~u, ~x, 1)) * Subst(Antiderivative((~A + ~B * (~x) ^ (~n - ~q)) * (~a * (~x) ^ ~q + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n - ~q)) ^ ~p, ~x), ~x, ~u) <-- FreeQ([~a, ~b, ~c, ~A, ~B, ~n, ~p, ~q], ~x) && (EqQ(~j, ~n - ~q) && (EqQ(~r, 2 * ~n - ~q) && (LinearQ(~u, ~x) && NeQ(~u, ~x))))
end

max_stars_repo_path: src/rules/1 Algebraic functions/1.2 Trinomial products/1.2.4 Improper/1.2.4.3 (d+e x^(n-q)) (a x^q+b x^n+c x^(2 n-q))^p.jl


text:
home = joinpath(homedir(),"polar/")
push!(LOAD_PATH, home)
using Light, Spectra
path = "CP"
a = collectspect(Float64,path)
polarspectra(a,"data.csv")
run(`pdflatex plot.tex`)
run(`convert -density 300 plot.pdf -resize 1280x800! plot.png`)
#run(`xfconf-query --channel xfce4-desktop --property /backdrop/screen0/monitor0/image-path --set $home/plot.png`)


text:
module MPC_Module

using JuMP
using OrdinaryDiffEq
using DiffEqBase

include("Base.jl")
using .Base

export
     MPC,
     defineMPC!,
     initOpt!,
     defineIP!,
     mapNames!,
     simIPlant!,
     updateX0!,
     currentIPState,
     goalReached!,
     simMPC!,
     predictX0!

########################################################################################
# MPC structs
########################################################################################

mutable struct IP
 control::Control
 state::State
end

function IP()
 IP(
  Control(),
  State()
  )
end

mutable struct EP
 control::Control
 state::State
end

function EP()
 EP(
  Control(),
  State()
 )
end

mutable struct MPCvariables
 # variables
 t::Float64           # current simulation time (s)
 tp::Any              # prediction time (if finalTimeDV == true -> this is not known before optimization)
 tex::Float64         # execution horizon time
 t0Actual                    # actual initial time TODO ?
 t0::Float64                  # mpc initial time TODO ?
 tf::Float64                  # mpc final time TODO ?
 t0Param::Any        # parameter for mpc t0  TODO ?
 evalNum::Int64       # parameter for keeping track of number of MPC evaluations
 goal                 # goal location w.r.t OCP
 goalTol             # tolerance on goal location
 initOptNum::Int64  # number of initial optimization
 previousSolutionNum::Int64  # number of times the previous solution should be used
end

function MPCvariables()
 MPCvariables(
              0.0,    # t
              Any,    # tp (might be a variable)
              0.5,    # tex
              0.0,
              0.0,
              0.0,
              Any,
              1,
              [],
              [],
              3,
              3
 )
end

mutable struct MPC
 v::MPCvariables
 ip::IP
 ep::EP
end

function MPC()
 MPC(
     MPCvariables(),
     IP(),
     EP()
     )
end

########################################################################################
# MPC functions
########################################################################################
"""
defineMPC!(n)
--------------------------------------------------------------------------------------\n
Author: Huckleberry Febbo, Graduate Student, University of Michigan
Date Create: 4/7/2017, Last Modified: 4/8/2018 \n
--------------------------------------------------------------------------------------\n
"""
function defineMPC!(n;
                   mode::Symbol=:OCP,
                   predictX0::Bool=true,
                   fixedTp::Bool=true,
                   tp::Any=Any,
                   tex::Float64=0.5,
                   IPKnown::Bool=true,
                   saveMode::Symbol=:all,
                   maxSim::Int64=100,
                   goal=n.ocp.XF,
                   goalTol=0.1*abs.(n.ocp.X0 - n.ocp.XF),
                   lastOptimal::Bool=true,
                   printLevel::Int64=2,
                   onlyOptimal::Bool=false)
 n.s.mpc.on = true
 n.mpc::MPC = MPC()
 n.s.mpc.mode = mode
 n.s.mpc.predictX0 = predictX0
 n.s.mpc.fixedTp = fixedTp
 n.mpc.v.tp = tp
 n.mpc.v.tex = tex
 n.s.mpc.IPKnown = IPKnown
 n.s.mpc.saveMode = saveMode
 n.s.mpc.maxSim = maxSim
 n.mpc.v.goal = goal
 n.mpc.v.goalTol = goalTol
 n.s.mpc.lastOptimal = lastOptimal
 n.s.mpc.printLevel = printLevel
 n.s.mpc.onlyOptimal = onlyOptimal

 n.f.mpc.simFailed[1] = false # for some reason it is getting defined as 0.0, not false during initialization
 n.f.mpc.defined = true
 return nothing
end

"""
# TODO consider letting user pass options
--------------------------------------------------------------------------------------\n
Author: Huckleberry Febbo, Graduate Student, University of Michigan
Date Create: 4/08/2018, Last Modified: 12/06/2019 \n
--------------------------------------------------------------------------------------\n
"""
function initOpt!(n;save::Bool=true, evalConstraints::Bool=false)
  if n.s.mpc.on
   error("call initOpt!() before defineMPC!(). initOpt!() will destroy n")
  end
  n.s.ocp.save = false
  n.s.mpc.on = false
  n.s.ocp.evalConstraints = false
  n.s.ocp.cacheOnly = true
  if n.s.ocp.save
   @warn "saving initial optimization results where functions where cached!"
  end
  for k in 1:n.mpc.v.initOptNum # initial optimization (s)
   status = optimize!(n)
   if status==:Optimal; break; end
  end
  # defineSolver!(n,solverConfig(c)) # modifying solver settings NOTE currently not in use
  n.s.ocp.save = save  # set to false if running in parallel to save time
  n.s.ocp.cacheOnly = false
  n.s.ocp.evalConstraints = evalConstraints # set to true to investigate infeasibilities
  return nothing
end

"""
# add a mode that solves as quickly as possible
# consider using the IP always.
defineModel!(n)
--------------------------------------------------------------------------------------\n
Author: Huckleberry Febbo, Graduate Student, University of Michigan
Date Create: 4/12/2018, Last Modified: 4/12/2018 \n
--------------------------------------------------------------------------------------\n
"""
function defineIP!(n,model;stateNames=[],controlNames=[],X0a=[])

   if isequal(n.s.mpc.mode,:OCP) # this function is called automatically for this mode
    if !isempty(stateNames)
     error("stateNames are set automatically for :mode == :OCP and cannot be provided.")
    end
    if !isempty(controlNames)
     error("controlNames are set automatically for :mode == :OCP and cannot be provided.")
    end
    if !isempty(X0a)
     error("X0a is set automatically for :mode == :OCP and cannot be provided.")
    end
    n.r.ip.X0a = copy(n.ocp.X0)  # NEED to append time
    n.mpc.ip.state.model = model
    n.mpc.ip.state.name = n.ocp.state.name
    n.mpc.ip.state.description = n.ocp.state.description
    n.mpc.ip.state.num = n.ocp.state.num
    n.mpc.ip.state.pts = n.ocp.state.pts

    n.mpc.ip.control.name = n.ocp.control.name
    n.mpc.ip.control.description = n.ocp.control.description
    n.mpc.ip.control.num = n.ocp.control.num
    n.mpc.ip.control.pts = n.ocp.control.pts

    # add X0 t0 plant dfs
    n.r.ip.plant[:t] = n.mpc.v.t0
    for st in 1:n.mpc.ip.state.num
      n.r.ip.plant[n.mpc.ip.state.name[st]] = copy(n.ocp.X0)[st]
    end
    for ctr in 1:n.mpc.ip.control.num
      n.r.ip.plant[n.mpc.ip.control.name[ctr]] = 0
    end

   elseif isequal(n.s.mpc.mode,:IP)
    if isempty(stateNames)
     error("unless :mode == :OCP the stateNames must be provided.")
    end
    if isempty(controlNames)
     error("unless :mode == :OCP the controlNames must be provided.")
    end
    if isempty(X0a)
     error("unless :mode == :OCP X0a must be provided.")
    end
    if isempty(model)
     error("A model needs to be passed for the IP mode.")
    else
    if isequal(length(X0a),length(stateNames))
      error(string("\n Length of X0a must match length(stateNames) \n"))
    end

     n.mpc.ip.state::State = State() # reset
     n.mpc.ip.state.num = length(stateNames)
     for i in 1:n.mpc.ip.state.num
       if stateNames[i]==:xxx
         error("xxx is OFF limits for a state name; please choose something else. \n")
       end
       push!(n.mpc.ip.state.name,stateNames[i])
     end

     n.mpc.ip.control::Control = Control() # reset
     n.mpc.ip.control.num = length(controlNames)
     for i in 1:n.mpc.ip.control.num
       if controlNames[i]==:xxx
         error("xxx is OFF limits for a control name; please choose something else. \n")
       end
       push!(n.mpc.ip.control.name,controlNames[i])
     end
     n.mpc.r.ip.X0a = X0a
     n.mpc.ip.state.model = model # TODO validate typeof model
    end
   elseif isequal(n.s.mpc.mode,:EP)
    error("not setup for :EP")
   else
    error("n.mpc.s.mode = ",n.s.mpc.mode," not defined." )
   end

   # consider calling mapNames
   return nothing
end

"""
mapNames!(n)
--------------------------------------------------------------------------------------\n
Author: Huckleberry Febbo, Graduate Student, University of Michigan
Date Create: 4/9/2018, Last Modified: 4/12/2018 \n
--------------------------------------------------------------------------------------\n
"""
function mapNames!(n)
  if isequal(n.s.mpc.mode,:IP)
    s1 = n.ocp.state.name
    c1 = n.ocp.control.name
    s2 = n.mpc.ip.state.name
    c2 = n.mpc.ip.control.name
  elseif isequal(n.s.mpc.mode,:EP)
    error(":EP function not ready")
  else
    error("mode must be either :IP or :EP")
  end

  m = []
  # go through all states in OCP
  idxOCP = 1
  for var in s1
    # go through all states in IP
    idxIP = findall(var.==s2)
    if !isempty(idxIP)
      push!(m, [var; :stOCP; idxOCP; :stIP; idxIP[1]])
    end

    # go through all controls in IP
    idxIP = findall(var.==c2)
    if !isempty(idxIP)
      push!(m, [var; :stOCP; idxOCP; :ctrIP; idxIP[1]])
    end
    idxOCP = idxOCP + 1
  end

  # go through all controls in OCP
  idxOCP = 1
  for var in c1
    # go through all states in IP
    idxIP = findall(var.==s2)
    if !isempty(idxIP)
      push!(m, [var; :ctrOCP; idxOCP; :stIP; idxIP[1]])
    end

    # go through all controls in IP
    idxIP = findall(var.==c2)
    if !isempty(idxIP)
      push!(m, [var; :ctrOCP; idxOCP; :ctrIP; idxIP[1]])
    end
    idxOCP = idxOCP + 1
  end

  if isequal(n.s.mpc.mode,:IP)
    n.mpc.mIP = m
  elseif isequal(n.s.mpc.mode,:EP)
    error(":EP function not ready")
  else
    error("mode must be either :IP or :EP")
  end

  return nothing
end

"""
# TODO fix this so that prediction simulation time is ahead. May not effect results.
# NOTE this may be ok... we are getting X0p for initialization of the OCP
# as long as the OCP pushes the time ahead. which it does then everything is fine!!
# consider making user pass X0, t0, tf
--------------------------------------------------------------------------------------\n
Author: Huckleberry Febbo, Graduate Student, University of Michigan
Date Create: 2/14/2017, Last Modified: 12/06/2019 \n
--------------------------------------------------------------------------------------\n
"""
function simIPlant!(n)
  if isequal(n.mpc.ip.state.pts,0)
   error("isqual(n.mpc.ip.state.pts,0), cannot simulate with zero points.")
  end
  X0 = currentIPState(n)[1]
  t0 = round(n.mpc.v.t, digits = 3) # if rounding is too rough, then the tex will be essentially 0!
  tf = round(n.mpc.v.t + n.mpc.v.tex,digits = 3)

  if isequal(n.s.mpc.mode,:OCP)
   if isequal(n.mpc.v.evalNum,1)
    U = 0*Matrix{Float64}(undef, n.ocp.control.pts,n.ocp.control.num)
    t = Vector(range(t0,tf,length=n.ocp.control.pts))
   elseif n.s.ocp.interpolationOn
    U = n.r.ocp.Upts
    t = n.r.ocp.tpts
   else
    U = n.r.ocp.U
    t = n.r.ocp.tctr
   end
  else
   error("TODO")
  end
  # chop of first control point for bkwEuler as it is typically 0
  if isequal(n.s.ocp.integrationScheme,:bkwEuler)
   U = U[2:end,:]
   t = t[2:end]
  end

  sol, U = n.mpc.ip.state.model(n,X0,t,U,t0,tf)
  return sol, U
end

"""
--------------------------------------------------------------------------------------\n
Author: Huckleberry Febbo, Graduate Student, University of Michigan
Date Create: 4/08/2018, Last Modified: 4/08/2018 \n
--------------------------------------------------------------------------------------\n
"""
function currentIPState(n)
  if isempty(n.r.ip.plant)
    error("there is no data in n.r.ip.plant")
  end

  # even though may have solution for plant ahead of time
  # can only get the state up to n.mpc.v.t
  idx = findall((n.mpc.v.t .- n.r.ip.plant[:t]) .>= 0)
  if isempty(idx)
    error("(n.mpc.v.t - n.r.ip.plant[:t]) .>= 0) is empty.")
  else
    X0 = [zeros(n.mpc.ip.state.num),n.mpc.v.t]
    for st in 1:n.mpc.ip.state.num
      X0[1][st] = n.r.ip.plant[n.mpc.ip.state.name[st]][idx[end]]
    end
  end
  return X0
end

"""
# TODO eventually the "plant" will be different from the "model"
predictX0!(n)
--------------------------------------------------------------------------------------\n
Author: Huckleberry Febbo, Graduate Student, University of Michigan
Date Create: 4/7/2017, Last Modified: 4/25/2018 \n
--------------------------------------------------------------------------------------\n
"""
function predictX0!(n)

  if n.s.mpc.fixedTp
   # NOTE consider passing back (n.mpc.v.t + n.mpc.v.tex) from simIPlant!()
   tp = round(n.mpc.v.t + n.mpc.v.tex, digits = 1)  # TODO add 1 as an MPCparamss
  else
   error("TODO")
  end

  if isequal(n.s.mpc.mode,:OCP)
   sol, U = simIPlant!(n)
   X0p = [sol(sol.t[end])[:],tp]
   push!(n.r.ip.X0p,X0p)
  else
    error("TODO")
  end
 # else
   # with no control signals to follow, X0p is simply the current known location of the plant
 #  X0 = currentIPState(n)
 #  X0p = [X0[1], tp]  # modify X0 to predict the time
  # push!(n.r.ip.X0p,X0p)
 # end
  return nothing
end

"""
--------------------------------------------------------------------------------------\n
Author: Huckleberry Febbo, Graduate Student, University of Michigan
Date Create: 3/06/2018, Last Modified: 4/08/2018 \n
--------------------------------------------------------------------------------------\n
"""
function updateX0!(n,args...)
 # need to map n.r.ip.X0p to n.X0 (states may be different)
 # NOTE for the :OCP mode this is OK

 if !n.s.mpc.predictX0 #  use the current known plant state to update OCP
   push!(n.r.ip.X0p,currentIPState(n))
 else
   predictX0!(n)
 end

 if isequal(n.s.mpc.mode,:OCP)
   if !isequal(length(args),0)
    X0 = args[1]
    if length(X0)!=n.ocp.state.num
      error(string("\n Length of X0 must match number of states \n"));
    end
    n.ocp.X0 = X0
   else
    n.ocp.X0 = n.r.ip.X0p[end][1] # the n.ocp. structure is for running things
   end
   push!(n.r.ocp.X0,n.ocp.X0)    # NOTE this may be for saving data
   setvalue(n.ocp.t0,copy(n.r.ip.X0p[end][2]))
 else
  error("not set up for this mode")
 end

  if n.s.mpc.shiftX0 # TODO consider saving linear shifting occurances
    for st in 1:n.ocp.state.num
      if n.ocp.X0[st] < n.ocp.XL[st]
        n.ocp.X0[st] = n.ocp.XL[st]
      end
      if n.ocp.X0[st] > n.ocp.XU[st]
        n.ocp.X0[st] = n.ocp.XU[st]
      end
    end
  end
  # update states with n.ocp.X0
  for st in 1:n.ocp.state.num
    if n.s.ocp.x0slackVariables
     JuMP.setRHS(n.r.ocp.x0Con[st,1], n.ocp.X0[st])
     JuMP.setRHS(n.r.ocp.x0Con[st,2],-n.ocp.X0[st])
    else
      JuMP.setRHS(n.r.ocp.x0Con[st],n.ocp.X0[st])
    end
  end
  return nothing
end

"""
--------------------------------------------------------------------------------------\n
Author: Huckleberry Febbo, Graduate Student, University of Michigan
Date Create: 4/08/2018, Last Modified: 4/08/2018 \n
--------------------------------------------------------------------------------------\n
"""
function goalReached!(n,args...)
  if isequal(n.s.mpc.mode,:OCP)
    X = currentIPState(n)[1]
  else
    X = args[1]
    error("TODO")
  end
  A = (abs.(X - n.mpc.v.goal) .<= n.mpc.v.goalTol)
  B = isnan.(n.mpc.v.goal)
  C = [A[i]||B[i] for i in 1:length(A)]

  if all(C)
   if isequal(n.s.mpc.printLevel,2)
    println("Goal Attained! \n")
   end
    n.f.mpc.goalReached = true
  elseif n.s.mpc.expandGoal && (getvalue(n.ocp.tf) < n.mpc.v.tex)
    A =( abs.(X - n.mpc.v.goal) .<= n.s.mpc.enlargeGoalTolFactor*n.mpc.v.goalTol)
    C = [A[i]||B[i] for i in 1:length(A)]
    if all(C)
     if isequal(n.s.mpc.printLevel,2)
      println("Expanded Goal Attained! \n")
     end
     n.f.mpc.goalReached = true
    else
     println("Expanded Goal Not Attained! \n
              Stopping Simulation")
     n.f.mpc.simFailed = [true, :expandedGoal]
    end
  end

 return n.f.mpc.goalReached
end
# if the vehicle is very close to the goal sometimes the optimization returns with a small final time
# and it can even be negative (due to tolerances in NLP solver). If this is the case, the goal is slightly
# expanded from the previous check and one final check is performed otherwise the run is failed
#if getvalue(n.ocp.tf) < 0.01
#  if ((n.r.ip.dfplant[end][:x][end]-c["goal"]["x"])^2 + (n.r..ip.dfplant[end][:y][end]-c["goal"]["yVal"])^2)^0.5 < 2*c["goal"]["tol"]
#  println("Expanded Goal Attained! \n"); n.f.mpc.goal_reached=true;
#  break;/
#  else
#  warn("Expanded Goal Not Attained! -> stopping simulation! \n"); break;
#  end
#elseif getvalue(n.ocp.tf) < 0.5 # if the vehicle is near the goal => tf may be less then 0.5 s
#  tf = (n.r.evalNum-1)*n.mpc.v.tex + getvalue(n.ocp.tf)
#else
#  tf = (n.r.evalNum)*n.mpc.v.tex
#end

"""
--------------------------------------------------------------------------------------\n
Author: Huckleberry Febbo, Graduate Student, University of Michigan
Date Create: 2/06/2018, Last Modified: 4/08/2018 \n
--------------------------------------------------------------------------------------\n
"""
function simMPC!(n;updateFunction::Any=[],checkFunction::Any=[])
  for ii = 1:n.s.mpc.maxSim
    if isequal(n.s.mpc.printLevel,2)
     println("Running model for the: ",n.mpc.v.evalNum," time")
    end
    #############################
    # (A) and (B) in "parallel"
    #############################

    # (B) simulate plant
    sol, U = simIPlant!(n) # the plant simulation time will lead the actual time
    plant2dfs!(n,sol,U)

    # check to see if the simulation failed (i.e. plant crashed)
    if !isequal(typeof(checkFunction),Array{Any,1})
     val, sym = checkFunction(n)
     if val
       n.f.mpc.simFailed = [val, sym]
      break
     end
    end

    # check to see if the goal has been reached
    if goalReached!(n); break; end
    if n.f.mpc.simFailed[1]; break; end

    # (A) solve OCP  TODO the time should be ahead here as it runs
    updateX0!(n)  # before updateFunction()
    if !isequal(typeof(updateFunction),Array{Any,1})
      updateFunction(n)
    end

    optimize!(n)
    if n.f.mpc.simFailed[1]; break; end

    # advance time
    n.mpc.v.t = n.mpc.v.t + n.mpc.v.tex
    n.mpc.v.evalNum = n.mpc.v.evalNum + 1
  end
end
# practical concerns/questions
#################################
# 1) predict X0
# 2) shift X0 for NLP feasibility
# 3) ensuring that U passed to the plant is feasible
     # effected by interpolation, demonstrate by varying numPts
     # seems to be a major problem with LGR nodes, possibly due to Runge effect
# 4) fixedTp or variableTp
# 5) usePrevious optimal.
     # at some point will be unable to do this
# 6) infeasibilities, soft constraint on inital conditions

# TODO
# 1) plot the goal, the tolerances on X0p
# 2) calculate the error and plot
end # module


text:
module AbstractOps

import Random
export AbstractOp
export weight, op_index, op_symbols, unsafe_op

abstract type AbstractOp end

"""
    op_index(op::AbstractOp)::Int

Return the `Int` index corresponding to `op`.
"""
function op_index end

function unsafe_op(::Type{T}, ind::Integer) where T
    return T(ind)
end

"""
    op_symbols(::Type{OpT}, ::Type{SymT})

Return a list of objects of type `SymT` that represent operators of type `OpT`.
The list should be in a canonical order.
# Examples
```julia
op_symbols(::Type{<:AbstractPauli}, ::Type{Symbol}) = (:I, :X, :Y, :Z)
```
"""
function op_symbols end

"""
    _AbstractOp(::Type{T}, ind::V)

Create an operator of type `T` from value `ind`.
"""
function _AbstractOp(::Type{T}, ind::V) where {T, V}
    syms = op_symbols(T, V)
    j = 0
    @inbounds for i in eachindex(syms) # 1:length(syms)
        if ind === syms[i]
            j = i
            break
        end
    end
    if j == 0
        throw(ArgumentError("Unrecognized operator symbol"))
    end
    return T(j - 1)  # TODO: abstract this index shift
end

"""
    Vector{T}(opstring::AbstractString) where {T <: AbstractOp}

Initialize a `Vector{T}` by converting each character to type `T`.
"""
Vector{T}(opstring::AbstractString) where {T <: AbstractOp} = [T(s) for s in opstring]

## TODO: If we trust rand_ind_range, we could use unsafe_pauli (or a generic version)
"""
    docstring
"""
function rand_ind_range end

# Using unsafe_op does not seem to be faster
function Random.rand(rng::Random.AbstractRNG, ::Random.SamplerType{OpT}) where {OpT <: AbstractOp}
#    return unsafe_op(OpT, rand(rng, rand_ind_range(OpT)))
    return OpT(rand(rng, rand_ind_range(OpT)))
end

"""
    weight(v::AbstractArray{<:AbstractOp})
    weight(ps::OpTerm{<:AbstractOp})

Count the number of Paulis in the string that are not the identity.
"""
weight(v::AbstractArray{<:AbstractOp}) = count(!isone, v)

"""
    docstring
"""
function accumulate_phase end

"""
    docstring
"""
function accumulate_phase! end

accumulate_phase!(args...) = nothing

"""
    docstring
"""
function compute_phase end

"""
    phase_data(::Type{T}) where T

Return a phase data object for type `T`.
"""
function phase_data end

"""
    docstring
"""
struct DummyPhaseData
end

phase_data(::Type{<:AbstractOp}) = DummyPhaseData()

function _show_op_plain end

function anticommutes end

# Return x and z bits as Tuple{Bool,Bool}
function symplectic end

end # module AbstractOps


text:
mutable struct Soma <: AbstractSoma
   	model::Model.ModelData
   	samples::Model.Samples

   # Axon is the output
   	axon::AbstractAxon

   	dendrite::AbstractDendrite

	# Soma threshold. When exceeded an AP is generated.
   	threshold::Float64
	# Post synaptic potential
   	psp::Float64

	# --------------------------------------------------------
	# Action potential
	# --------------------------------------------------------
	# AP can travel back down the dendrite. The value decays
	# with distance.
   	apFast::Float64 # Fast trace
   	apSlow::Float64 # Slow trace
   	apSlowPrior::Float64 # Slow trace (t-1)

	# The time-mark of the current AP.
   	APt::Float64
	# The previous time-mark
   	preAPt::Float64
   	APMax::Float64

	# --------------------------------------------------------
	# STDP
	# --------------------------------------------------------
	# -----------------------------------
	# AP decay
	# -----------------------------------
   	ntao::Float64 # fast trace
   	ntaoS::Float64 # slow trace

	# Fast Surge
   	nFastSurge::Float64
   	nDynFastSurge::Float64
   	nInitialFastSurge::Float64

	# Slow Surge
   	nSlowSurge::Float64
   	nDynSlowSurge::Float64
   	nInitialSlowSurge::Float64

	# The time-mark at which a spike arrived at a synapse
   	preT::Float64

    refractoryPeriod::Float64
    refractoryCnt::Float64
    refractoryState::Bool

	# -----------------------------------
	# Suppression
	# -----------------------------------
   	ntaoJ::Float64
    efficacyTrace::Float64

    function Soma(axon::AbstractAxon, model::Model.ModelData, samples::Model.Samples)
        o = new()
        o.axon = axon
        o.model = model
        o.samples = samples
        o
    end
end

function initialize!(soma::AbstractSoma)
	# Set properties based on model. These drive the other properties.
   	soma.refractoryPeriod = Model.refractory_period(soma.model)
   	soma.nInitialFastSurge = Model.fast_surge(soma.model)
   	soma.nInitialSlowSurge = Model.slow_surge(soma.model)
   	soma.ntao = Model.tao(soma.model)
   	soma.ntaoS = Model.tao_s(soma.model)
   	soma.ntaoJ = Model.tao_j(soma.model)
   	soma.threshold = Model.threshold(soma.model)
   	soma.APMax = Model.ap_max(soma.model)
    # println(soma.nInitialFastSurge, ", ", soma.nInitialSlowSurge, ", ", soma.ntao, ", ", soma.ntaoS)

   	println("___ Soma properties ___")
    println("| refractoryPeriod: ", soma.refractoryPeriod)
    println("| nInitialFastSurge: ", soma.nInitialFastSurge)
    println("| nInitialSlowSurge: ", soma.nInitialSlowSurge)
    println("| ntao: ", soma.ntao)
    println("| ntaoS: ", soma.ntaoS)
    println("| ntaoJ: ", soma.ntaoJ)
    println("| threshold: ", soma.threshold)
    println("| APMax: ", soma.APMax)
    println("---------------------------")

   	initialize!(soma.dendrite)
end

function reset!(soma::AbstractSoma)
   	# print_trace(stacktrace(), "")
    soma.apFast = 0.0
    soma.apSlow = 0.0
    soma.preT = -1000000000000000.0
    soma.refractoryState = false
    soma.refractoryCnt = 0
    soma.nSlowSurge = 0.0
    soma.nFastSurge = 0.0
    soma.efficacyTrace = 0.0

   	reset!(soma.axon)
    reset!(soma.dendrite)
end

function integrate!(soma::AbstractSoma, span_t::Int64, t::Int64)
   	dt = t - soma.preT

   	soma.efficacyTrace = efficacy(soma, dt)

	# The dendrite will return a value that affects the soma.
    soma.psp = integrate!(soma.dendrite, span_t, t)

	# Default state
    set!(soma.axon, UInt8(0)) # Set output

    if soma.refractoryState
		# this algorithm should be the same as for the synapse or at least very
		# close.

        if soma.refractoryCnt >= soma.refractoryPeriod
            soma.refractoryState = false
            soma.refractoryCnt = 0
			# fmt.Printf("Refractory ended at (%d)\n", int(t))
        else
            soma.refractoryCnt += 1
        end
    else
        if soma.psp > soma.threshold
			# An action potential just occurred.
			# TODO Handle depolarization

            soma.refractoryState = true

			# TODO
			# Generate a back propagating spike that fades spatial/temporally similar to CaDP model.
			# This spike affects forward in time.
			# The value is driven by the time delta of (preAPt - APt)
           	set!(soma.axon, UInt8(1)) # Set output

			# Surge from action potential
            soma.nFastSurge = soma.APMax + soma.apFast * soma.nInitialFastSurge * exp(-soma.apFast / soma.ntao)
            soma.nSlowSurge = soma.APMax + soma.apSlow * soma.nInitialSlowSurge * exp(-soma.apSlow / soma.ntaoS)

			# Reset time deltas
            soma.preT = t
            dt = 0
        end
   	end

	# Prior is for triplet
   	soma.apSlowPrior = soma.apSlow
    # println(soma.nFastSurge, ", ", soma.nSlowSurge, ", ", soma.ntao, ", ", soma.ntaoS)
   	soma.apFast = soma.nFastSurge * exp(-dt / soma.ntao)
   	soma.apSlow = soma.nSlowSurge * exp(-dt / soma.ntaoS)

	# Collecting is centralized in streams.jl for consistency.
   	collect!(soma.samples, soma, span_t)

   	output(soma.axon)
end

function output(soma::AbstractSoma)
   	output(soma.axon)
end

# This is a time based property NOT distance.
# Each spike of the neuron i sets the post spike efficacy j to 0
# whereafter it recovers exponentially to 1 with a time constant toaI.
# In other words, the efficacy of a spike is suppressed by
# the proximity of a previous post spike.
function efficacy_trace(soma::AbstractSoma)
   	soma.efficacyTrace
end

function efficacy(soma::AbstractSoma, dt::Float64)
   	1.0 - exp(-dt / soma.ntaoJ)
end


text:
function f()
	x = 42
	while x >= 0
		println(x)
		x = x - 1
	end
end

f()

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
# -------------------------------------------------------------------------------------------------------------------------------------------------------------
# Script with functions to import and transform raw data
# -------------------------------------------------------------------------------------------------------------------------------------------------------------

# --------------- Raw data ---------------------------------------------------------------------------------------------------------------------------------

"""
    import_data(dir::String, revision::String, year::Integer, N::Integer, S::Integer)

The function loads the Input-Output Table sourced form the OECD or WIOD for a specified year as DataFrame into the environment and
    extracts the raw intermediate, final demand matrices (Z, F) and the inventory matrix IV.

# Arguments
- `dir::String`: directory of raw data (either folder location or entire directory).
- `source::String`: specifies the source of the IO table (OECD, WIOD) as a string.
- `revision::String`: specifies the of the IO table as a string.
- `year::Integer`: specifies the year of the WIOD table which should be imported.
- `N::Integer`: number of origin/destination countries.
- `S::Integer`: number of origin/destination industries.

# Output
- `Z::Matrix{Float64}`: NS×NS, raw origin country-industry destination country-industry intermediate demand matrix Z.
- `F::Matrix{Float64}`: NS×N, raw origin country-industry destination country final demand matrix F.
- `IV::Matrix{Float64}`: NS×N, raw origin country-industry destination country inventory adjustment matrix IV.

# Examples
```julia-repl
julia> Z, F, IV = import_data("C:/Users/u0148308/Desktop/raw/", "2016", 2014, 44, 56)
```
"""

function import_data(dir::String, source::String, revision::String, year::Integer, N::Integer, S::Integer)

    if source == "OECD"

        if revision == "2021"

            # still need to sort out the issue with China and Mexico being split into two!

            file = source * "/" * revision * "/"
            path = ifelse(contains(dir[end-5:end], '.'), dir, dir * file)

            time = 1995:2018
            index = findfirst(time .== year) # get index of specified year

            Z = RData.load(path * "ICIO2021econZ" * ".RData")["ICIO2021econZ"][index,:,:]
            F = RData.load(path * "ICIO2021econFD" * ".RData")["ICIO2021econFD"][index,:,:]
            IV = RData.load(path * "ICIO2021econINVNT" * ".RData")["ICIO2021econINVNT"][index,:,:]

            F = F .- IV # take out inventory from final demand

            println(" ✓ Raw data from $source (rev. $revision) for $year was successfully imported!")

        else
            println(" × The revision - $revision - of $source IO tables is not available!")
        end

    else source == "WIOD"

        if revision == "2016"

            file = source * "/" * revision * "/WIOT" * string(year) * "_October16_ROW.RData"
            path = ifelse(contains(dir[end-5:end], '.'), dir, dir * file)

            df = RData.load(path)["wiot"]
            df = df[1:N*S, 6:end-1] # NS×NS+5N, take out the rows/columns with country/industry names

        elseif revision == "2013"

            name_2013 = ifelse(year >= 2008, "Sep", "Apr")
            file = source * "/" * revision * "/WIOT" * string(year)[end-1:end] * "_ROW_" * name_2013 * "12.xlsx"
            path = ifelse(contains(dir[end-5:end], '.'), dir, dir * file)

            df = DataFrames.DataFrame(XLSX.readxlsx(path)["WIOT_$year"][:], :auto)
            df = df[7:end-8,5:end-1] # NS×NS+5N, take out the rows/columns with country/industry names

        else
            println(" × The revision - $revision - of $source IO tables is not available!")
        end

        IO_table = Matrix(convert.(Float64, df)) # NS×NS+5N
        inventory_columns = N*S+5:5:N*S+5*N

        Z = IO_table[:, 1:N*S] # NS×NS

        F = IO_table[:, Not([1:N*S; inventory_columns])] # NS×4N
        F = [sum(F[i,j:j+3]) for i in 1:N*S, j in 1:4:N*4] # NS×N, sum all types of final demand

        IV = IO_table[:, inventory_columns] # NS×N

        println(" ✓ Raw data from $source (rev. $revision) for $year was successfully imported!")

    end

    return Z, F, IV
end

# --------------- Transformations -----------------------------------------------------------------------------------------------------------------------------

"""
    inventory_adjustment(Z::Matrix, F::Matrix, IV::Matrix, N::Integer, S::Integer)

The function performs the inventory adjustment procedure proposed in Antras et al. (2012).

# Arguments
- `Z::Matrix`: NS×NS, raw origin country-industry destination country-industry intermediate demand matrix Z.
- `F::Matrix`: NS×N, raw origin country-industry destination country final demand matrix F.
- `IV::Matrix`: NS×N, raw origin country-industry destination country inventory adjustment matrix IV.
- `N::Integer`: number of origin/destination countries.
- `S::Integer`: number of origin/destination industries.

# Output
- `Z::Matrix{Float64}`: NS×NS, inventory adjusted origin country-industry destination country-industry intermediate demand matrix Z.
- `F::Matrix{Float64}`: NS×N, inventory adjusted origin country-industry destination country final demand matrix F.

# Examples
```julia-repl
julia> Z, F = inventory_adjustment(Z, F, IV, N, S)
```
"""

function inventory_adjustment(Z::Matrix, F::Matrix, IV::Matrix, N::Integer, S::Integer)

    # remove zero and negative entries (page 37 Antras and Chor, 2018)
    #Z = ifelse.(Z .<= 0.0, 1e-18, Z)
    #F = ifelse.(F .<= 0.0, 1e-18, F)

    Z = ifelse.(Z .< 0.0, 0.0, Z) # no negative values possible in Z
    F = ifelse.(F .< 0.0, 0.0, F) # no negative values possible in F

    # Inventory adjustment (spread IV regardless of origin over all columns)
    IV = [sum(IV[i,:]) for i in 1:N*S] # NS×1, negative values possible
    Y = [sum(Z[i,:]) + sum(F[i,:]) + IV[i] for i in 1:N*S] # NS×1, possibly negative if IV is negative enough
    Y = ifelse.(Y .< 0.0, 0.0, Y) # exclude negative values

    adjustment = Y .- IV # NS×1
    adjustment = ifelse.(adjustment .< 0.0, 0.0, adjustment) # in case negative, as in Antras and Chor (2018)

    # Inventory adjustment
    Z = Z .* ( repeat(Y, 1, N*S) ./ repeat(adjustment, 1, N*S) ) # NS×NS, inventory adjustment
    Z = ifelse.(isnan.(Z), 0.0, Z)
    Z = ifelse.(isinf.(Z), 0.0, Z)
    Z = ifelse.(Z .< 0.0, 0.0, Z) # NS×NS

    F = F .* ( repeat(Y, 1, N) ./ repeat(adjustment, 1, N) ) # NS×N, inventory adjustment
    F = ifelse.(isnan.(F), 0.0, F)
    F = ifelse.(isinf.(F), 0.0, F)
    F = ifelse.(F .< 0.0, 0.0, F) # NS×N

    return Z, F
end

# ---------------

"""
    consistency_adjustment(Z::Matrix, F::Matrix, N::Integer, S::Integer)

The function performs the value added adjustment which borrows from the inventory adjustment procedure proposed in Antras et al. (2012). This makes sure that
    no entries of the value added vector is negative and that the system of Z, F, Y and VA is internally consistent, i.e.
    Z + F = Y_exports = Y_imports = Z + VA.

# Arguments
- `Z::Matrix`: NS×NS, raw origin country-industry destination country-industry intermediate demand matrix Z.
- `F::Matrix`: NS×N, raw origin country-industry destination country final demand matrix F.
- `IV::Matrix`: NS×N, raw origin country-industry destination country inventory adjustment matrix IV.
- `N::Integer`: number of origin/destination countries.
- `S::Integer`: number of origin/destination industries.

# Output
- `Z::Matrix{Float64}`: NS×NS, consistent origin country-industry destination country-industry intermediate demand matrix Z.
- `F::Matrix{Float64}`: NS×N, consistent origin country-industry destination country final demand matrix F.
- `Y::Vector{Float64}`: NS×1, consistent country-industry gross output vector Y.
- `VA::Vector{Float64}`: NS×1, consistent country-industry value added vector VA.

# Examples
```julia-repl
julia> Z, F, Y, VA = consistency_adjustment(Z, F, IV, N, S)
```
"""

function consistency_adjustment(Z::Matrix, F::Matrix, N::Integer, S::Integer)

    # Recompute values for Y and VA to have internal consistency
    Y = [sum(Z[i,:]) + sum(F[i,:]) for i in 1:N*S] # NS×1
    VA = [Y[i] - sum(Z[:,i]) for i in 1:N*S] # NS×1, compute value added (gross output minus sum of inputs)
    #VA = ifelse.(VA .< 0.0, 0.0, VA) # in case one is negative (actually problematic since Z, F, Y internally consistent but not VA!)

    # value added adjustment
    # adjust Z for negative VA values in order to have internal consistency, i.e. Z + F = Y_exports = Y_imports = Z + VA
    adjustment = ifelse.(VA .< 0.0, VA, 0.0) # spread only negative values
    adjustment = Y .- adjustment
    adjustment = ifelse.(adjustment .< 0.0, 0.0, adjustment)
    adjustment = Y ./ adjustment

    Z = Z .* repeat(transpose(adjustment), N*S) # NS×NS, inventory adjustment
    Z = ifelse.(isnan.(Z), 0.0, Z)
    Z = ifelse.(isinf.(Z), 0.0, Z)
    Z = ifelse.(Z .< 0.0, 0.0, Z) # NS×NS

    # rebuild everything
    Y = [sum(Z[i,:]) + sum(F[i,:]) for i in 1:N*S] # NS×1
    VA = [Y[i] - sum(Z[:,i]) for i in 1:N*S] # NS×1, should now be without any negative values

    # calculating gross output via exports/imports+VA should give the same results
    # still not exactly the same! but all differences <1e10
    Y_imports = [sum(Z[:,j]) + VA[j] for j in 1:N*S]

    if Y ≈ Y_imports
        println(" ✓ The matrices are internally consistent!")
    else
        println(" × Failure! The matrices are internally inconsistent!")
    end

    return Z, F, Y, VA
end

# ---------------

"""
    create_expenditure_shares(Z::Matrix, F::Matrix, Y::Vector, N::Integer, S::Integer)

The function computes the origin country-industry destination country-industry intermediate goods expenditure (trade) share matrix π_Z and
    origin country-industry destination country final goods trade share matrix π_F.
For example π_Z[i,j] = 0.2 means that the country-sector pair j sources 20% of its inputs from the country-sector pair i.

# Arguments
- `Z::Matrix`: NS×NS, origin country-industry destination country-industry intermediate demand matrix Z.
- `F::Matrix`: NS×N, origin country-industry destination country final demand matrix F.
- `N::Integer`: number of origin/destination countries.
- `S::Integer`: number of origin/destination industries.

# Output
- `π_Z::Matrix{Float64}`: NS×NS, origin country-industry destination country-industry intermediate goods trade share matrix π_Z.
- `π_F::Matrix{Float64}`: NS×N, origin country-industry destination country foods goods trade share matrix π_Z.

# Examples
```julia-repl
julia> π_Z, π_F = create_expenditure_shares(Z, F, N, S)
```
"""

function create_expenditure_shares(Z::Matrix, F::Matrix, N::Integer, S::Integer)

    Z_agg = [sum(Z[s:S:(N-1)*S+s,j]) for s in 1:S, j in 1:N*S] # S×NS, sum over origin industries
    π_Z = Z ./ repeat(Z_agg, N) # NS×NS
    π_Z = ifelse.(isnan.(π_Z), 0.0, π_Z) # needed in case Z_agg = 0 (highly unlikely  since all entries in Z >= 0)
    !any(0.0 .<= π_Z .<= 1.0) && println(" × Problem! Not all elements in γ are in intervall [0, 1]")

    F_agg = [sum(F[s:S:(N-1)*S+s,j]) for s in 1:S, j in 1:N] # S×N, sum over origin industries
    π_F = F ./ repeat(F_agg, N) # NS×N
    π_F = ifelse.(isnan.(π_F), 0.0, π_F) # needed in case F_agg = 0 (highly unlikely since all entries in F >= 0)
    !any(0.0 .<= π_Z .<= 1.0) && println(" × Problem! Not all elements in γ are in intervall [0, 1]")

    return π_Z, π_F
end

# ---------------

"""
    create_input_shares(Z::Matrix, F::Matrix, Y::Vector, N::Integer, S::Integer)

The function computes the country value added vector VA_ctry respecting inventory adjustments, intermediate goods expenditure share matrix γ
    and final goods import expenditure share matrix α. Additionally, country-industry value added coeffiecient matrix VA_coeff as well as
    country trade balance vector TB_ctry and country final goods import vector F_ctry are produced.
For example if α[s,i] = 0.2 means that country i sources 20% of its final demand from sector s.
For example if VA_coeff[s,i] = 0.2 means that 20% of output in country-industry pair [i,s] was value added.
For example if γ[s,j] = 0.2 means that 20% of inputs in country-industry pair j was sourced from sector s.

# Arguments
- `Z::Matrix`: NS×NS, origin country-industry destination country-industry intermediate demand matrix Z.
- `F::Matrix`: NS×N, origin country-industry destination country final demand matrix F.
- `Y::Vector`: NS×1, origin country-industry gross output vector Y.
- `VA::Vector`: NS×1, country-industry value added vector VA.
- `π_Z::Matrix`: NS×NS, origin country-industry destination country-industry intermediate goods trade share matrix π_Z.
- `π_F::Matrix`: NS×N, origin country-industry destination country foods goods trade share matrix π_Z.
- `N::Integer`: number of origin/destination countries.
- `S::Integer`: number of origin/destination industries.

# Output
- `γ::Matrix{Float64}`: S×NS, country-industry intermediate input expenditure share matrix γ.
- `α::Matrix{Float64}`: S×N, country-industry final expenditure expenditure share matrix α.
- `VA_coeff::Vector{Float64}`: NS×1, country-industry value added coefficients vector VA_coeff.
- `TB_ctry::Vector{Float64}`: NS×1, country trade balance (X-M) vector TB.
- `VA_ctry::Vector{Float64}`: N×1, country value added vector VA_ctry.
- `F_ctry::Vector{Float64}`: N×1, country final import demand vector VA_ctry.

# Examples
```julia-repl
julia> γ, α, VA_coeff, TB_ctry, VA_ctry, F_ctry = create_input_shares(Z, F, Y, π_Z, π_F, VA, N, S)
```
"""

function create_input_shares(Z::Matrix, F::Matrix, Y::Vector, VA::Vector, π_Z::Matrix, π_F::Matrix, N::Integer, S::Integer)

    # Country-industry level intermediate import expenditure shares
    Z_agg = [sum(Z[i:S:(N-1)*S+i, j]) for i in 1:S, j in 1:N*S] # S×NS, sum over origin industries
    Z_agg_VA = [Z_agg; transpose(VA)] # S+1×NS
    γ_VA = [Z_agg_VA[i,j]/sum(Z_agg_VA[:,j]) for i in 1:S+1, j in 1:N*S]

    γ = γ_VA[1:S, 1:N*S] # S×NS
    γ = ifelse.(isnan.(γ), 0.0, γ) # needed in case Z_agg_VA = 0 (unlikely since all entries in Z, VA >= 0, but still necessary)
    !any(0.0 .<= γ .<= 1.0) && println(" × Problem! Not all elements in γ are in intervall [0, 1]")

    # Country-industry level value added coefficients
    VA_coeff = γ_VA[S+1, 1:N*S] # NS×1
    VA_coeff = ifelse.(isnan.(VA_coeff), 1.0, VA_coeff) # needed in case Z_agg_VA = 0 (unlikely since all entries in Z, VA >= 0, but still necessary)
    #VA_coeff = ifelse.(VA_coeff .< 0.0, 0.0, VA_coeff)

    # Country level value added
    VA_ctry = [sum(VA[i:i+S-1]) for i in 1:S:N*S] # N×1

    # Final imports at country level (sum over rows not columns - otherwise exports)
    F_ctry = [sum(F[:,j]) for j in 1:N] # N×1

    # Country-industry final import expenditure shares (columns sum to 1)
    α = [sum(F[i:S:(N-1)*S+i, j])/F_ctry[j] for i in 1:S, j in 1:N] # S×N
    α = ifelse.(isnan.(α), 0.0, α) # needed in case F_ctry = 0 (unlikely since all entries in F >= 0, but still necessary)
    !any(0.0 .<= α .<= 1.0) && println(" × Problem! Not all elements in γ are in intervall [0, 1]")

    # Country level trade balance (exports - imports)
    # need to calculate country level exports/imports (i.e. aggregate and remove intra-country trade)
    E = [sum(Y[i:i+S-1]) - sum(Z[i:i+S-1,i:i+S-1]) - sum(F[i:i+S-1,ceil(Int, i/S)]) for i in 1:S:N*S] # N×1

    M_Z = [sum(Z[:,j:j+S-1]) - sum(Z[j:j+S-1,j:j+S-1]) for j in 1:S:N*S]
    M_F = [sum(F[:,ceil(Int, j/S)]) - sum(F[j:j+S-1,ceil(Int, j/S)]) for j in 1:S:N*S]
    M = M_Z .+ M_F # N×1

    TB_ctry = E .- M # N×1

    #TB_ctry == VA_ctry .- F_ctry # must hold

    return γ, α, VA_coeff, TB_ctry, VA_ctry, F_ctry
end

# -------------------------------------------------------------------------------------------------------------------------------------------------------------
# Summary function
# -------------------------------------------------------------------------------------------------------------------------------------------------------------

"""
    transform_data(dir::String, revision::String, year::Integer, N::Integer, S::Integer)

The function imports and transforms the raw data for further use in the model.

# Arguments
- `dir::String`: directory of raw data (either folder location or entire directory).
- `source::String`: specifies the source of the IO table (OECD, WIOD) as a string.
- `revision::String`: specifies the revision of the IO table as a string.
- `year::Integer`: specifies the year of the WIOD table which should be imported.
- `N::Integer`: number of origin/destination countries.
- `S::Integer`: number of origin/destination industries.

# Output
- `Z::Matrix{Float64}`: NS×NS, origin country-industry destination country-industry intermediate demand matrix Z.
- `F::Matrix{Float64}`: NS×N, origin country-industry destination country final demand matrix F.
- `Y::Vector{Float64}`: NS×1, origin country-industry gross output vector Y.
- `F_ctry::Vector{Float64}`: N×1, country final import demand vector VA_ctry.
- `TB_ctry::Vector{Float64}`: NS×1, country trade balance (X-M) vector TB.
- `VA_ctry::Vector{Float64}`: N×1, country value added vector VA_ctry.
- `VA_coeff::Vector{Float64}`: NS×1, country-industry value added coefficients vector VA_coeff.
- `γ::Matrix{Float64}`: S×NS, country-industry intermediate input expenditure share matrix γ.
- `α::Matrix{Float64}`: S×N, country-industry final expenditure expenditure share matrix α.
- `π_Z::Matrix{Float64}`: NS×NS, origin country-industry destination country-industry intermediate goods trade share matrix π_Z.
- `π_F::Matrix{Float64}`: NS×N, origin country-industry destination country final goods trade share matrix π_F.

# Examples
```julia-repl
julia> Z, F, Y, F_ctry, TB_ctry, VA_ctry, VA_coeff, γ, α, π_Z, π_F = transform_data(dir, "2013", 1995, N, S)
```

"""

function transform_data(dir::String, source::String, revision::String, year::Integer, N::Integer, S::Integer)

    Z, F, IV = import_data(dir, source, revision, year, N, S)
    Z, F = inventory_adjustment(Z, F, IV, N, S)
    Z, F, Y, VA = consistency_adjustment(Z, F, N, S)
    π_Z, π_F = create_expenditure_shares(Z, F, N, S)
    γ, α, VA_coeff, TB_ctry, VA_ctry, F_ctry = create_input_shares(Z, F, Y, VA, π_Z, π_F, N, S)

    println(" ✓ Corresponding transformations from WIOD (rev. $revision) for $year were successfully computed! \n")

    return Z, F, Y, F_ctry, TB_ctry, VA_ctry, VA_coeff, γ, α, π_Z, π_F
end


text:
@testset "sgd" begin
	include("sgdTest.jl")
end

@testset "objFun" begin
	include("dnnObjFctnTest.jl")
end


text: my_f(x, y) = 2x + 3y


text:
module MuFinder

using Parameters
using Statistics
using LinearAlgebra
using Printf
using Statistics

using ..Models: AbstractModel, HolsteinModel, SSHModel, update_model!
using ..GreensFunctions: EstimateGreensFunction, setup!
using ..Measurements: measure_N², measure_density

export MuTuner, update_μ!

mutable struct MuTuner{T<:AbstractFloat}

    active      :: Bool
    μ_traj      :: Vector{T}
    N_traj      :: Vector{T}
    N²_traj     :: Vector{T}
    forgetful_c :: T
    μ           :: T
    N           :: Int
    β           :: T
    Δτ          :: T
    L           :: Int
    target_N    :: T
    μ_bar       :: T
    μ_std       :: T
    κ_bar       :: T
    N_bar       :: T
    N_std       :: T
    N²_bar      :: T
    μ_bar_traj  :: Vector{T}
    κ_bar_traj  :: Vector{T}
    N_bar_traj  :: Vector{T}
    N²_bar_traj :: Vector{T}
    κ_min       :: T
    μ_avg       :: T
    μ_err       :: T
    logfile     :: String

    function MuTuner(active::Bool, init_μ::T, target_N::T, N::Int, β::T, Δτ::T, forgetful_c::T, κ_min::T, logfile::String) where{T<:AbstractFloat}

        L           = round(Int,β/Δτ)
        μ_traj      = [init_μ]
        N_traj      = Vector{T}()
        N²_traj     = Vector{T}()
        μ_bar_traj  = Vector{T}()
        κ_bar_traj  = Vector{T}()
        N_bar_traj  = Vector{T}()
        N²_bar_traj = Vector{T}()

        if !isfile(logfile) && active
            open(logfile,"w") do file
                write(file,"mu_bar kappa_bar n_bar Nsqr_bar mu n Nsqr\n")
            end
        end

        return new{T}(active, μ_traj, N_traj, N²_traj, forgetful_c, init_μ, N, β, Δτ, L, target_N, init_μ, 0.0, κ_min, -1.0, 0.0, -1.0,
                      μ_bar_traj, κ_bar_traj, N_bar_traj, N²_bar_traj, κ_min, init_μ, 0.0, logfile)
    end
end

"""
Update μ values in model.
"""
function update_μ!(model::AbstractModel{T},  tuner::MuTuner{T}, estimator::EstimateGreensFunction{T})::T where {T<:AbstractFloat}

    μ  = model.μ::Vector{T}
    μ₀ = mean(μ)

    if tuner.active

        # ⟨N⟩ and ⟨N²⟩ value
        N  = 0.0
        N² = 0.0

        # iterate over all possible pairs of random vectors
        for i in 1:(estimator.nᵥ-1)
            for j in (i+1):estimator.nᵥ

                # set up estimates
                setup!(estimator,i,j)

                # measure ⟨N⟩
                N  += model.Nsites * measure_density(model,estimator)

                # measure ⟨N²⟩
                N² += real(measure_N²(model,estimator))
            end
        end

        # normalize measurement
        N  /= binomial(estimator.nᵥ,2)
        N² /= binomial(estimator.nᵥ,2)

        # update μ
        μ₁      = update_μ!(tuner,N,N²)
        Δμ      = μ₁-μ₀
        @. μ   += Δμ
        tuner.μ = μ₁

        return μ₁
    else

        tuner.μ = μ₀

        return μ₀
    end
end

"""
Given new measurements for N and N², updates the MuTuner and returns the new value of μ.
"""
function update_μ!(tuner::MuTuner, N::T, N²::T)::T where {T<:AbstractFloat}

    @unpack μ_traj, N_traj, N²_traj, μ_bar_traj, N_bar_traj, N²_bar_traj, κ_bar_traj, forgetful_c, β, target_N, κ_min = tuner

    # record new ⟨N⟩ and ⟨N²⟩ values
    push!(N_traj,  N)
    push!(N²_traj, N²)

    # calculate new averages
    tuner.μ_bar, tuner.μ_std = forgetful_welfords(μ_traj, tuner.μ_bar, tuner.μ_std, forgetful_c)
    tuner.N_bar  = forgetful_mean(N_traj, tuner.N_bar, forgetful_c)
    tuner.N²_bar = forgetful_mean(N²_traj, tuner.N²_bar, forgetful_c)
    push!(μ_bar_traj,tuner.μ_bar)
    push!(N_bar_traj,tuner.N_bar)
    push!(N²_bar_traj,tuner.N²_bar)

    # length of trajectory
    n = length(N_traj)

    # variance of N
    varN = tuner.N²_bar - tuner.N_bar^2

    # calculate lower bound for κ
    κ_lo = κ_min/sqrt(n)

    # calculate upper bound for κ
    if n==1 || varN < 0.0 || tuner.μ_std <= 0.0
        κ_hi = κ_lo
    else
        κ_hi = sqrt(varN)/tuner.μ_std
    end

    # calculate κ
    tuner.κ_bar  = β * varN

    # apply bounds to κ value
    tuner.κ_bar  = min( tuner.κ_bar , κ_hi )
    tuner.κ_bar  = max( tuner.κ_bar , κ_lo )
    push!(κ_bar_traj,tuner.κ_bar)

    # write to log file
    if tuner.active
        open(tuner.logfile,"a") do file
            @printf file "%.8f %.8f %.8f %.8f %.8f %.8f %.8f\n" tuner.μ_bar (tuner.κ_bar/tuner.N) (tuner.N_bar/tuner.N) tuner.N²_bar tuner.μ (N/tuner.N) N²
        end
    end

    # calculate new μ value
    tuner.μ = tuner.μ_bar + (target_N - tuner.N_bar) / tuner.κ_bar
    push!(μ_traj, tuner.μ)

    return tuner.μ
end

"""
Given a MuTuner, returns the best guess for (μ, err_μ) from
previous trajectory.
"""
function estimate_μ(tuner::MuTuner{T}) where {T<:AbstractFloat}

    if tuner.active

        # forgetfulnes parameter
        forgetful_c = tuner.forgetful_c

        # if entire history is used in mu-tuning, then only use
        # half the history to estimate μ
        if forgetful_c == 1.0
            forgetful_c = 0.5
        end

        # calculate best estimate for true μ value
        forgetful_idx = ceil(Int, forgetful_c * length(tuner.μ_traj))
        μ_traj        = @view tuner.μ_traj[forgetful_idx:length(tuner.μ_traj)]
        μ_err         = stdm( μ_traj , median(μ_traj) )
        tuner.μ_avg   = tuner.μ_bar
        tuner.μ_err   = μ_err

        return nothing
    else

        tuner.μ_avg = tuner.μ
        tuner.μ_err = 0.0

        return nothing
    end
end

############################
## PRIVATE MODULE METHODS ##
############################

"""
Update the mean of the vector `x` using only the most recent `c` fraction of the history.
"""
function forgetful_mean(x::AbstractVector{T}, x̄ₙ₋₁::T, c::T)::T where {T<:AbstractFloat}

    N = length(x)
    if N  == 1
        return x[1]
    end
    i′   = 1 + floor(Int, (1.0 - c) * (N - 1))
    n    = N - i′ + 1
    xₙ   = x[N]
    x̄ₙ   = x̄ₙ₋₁ + (xₙ-x̄ₙ₋₁)/n
    i    = 1 + floor(Int, (1.0 - c) * N)
    if i != i′
        n  = N - i + 1
        x₀ = x[i′]
        x̄ₙ = x̄ₙ - (x₀-x̄ₙ)/n
    end
    return x̄ₙ
end

"""
Update the mean and standard deviation of the vector `x` using only the most recent `c` fraction of the history.
Uses the numerically stable Welford's algorithm to update the standard deviation.
"""
function forgetful_welfords(x::AbstractVector{T},x̄ₙ₋₁::T,sₙ₋₁::T,c::T)::Tuple{T,T} where {T<:AbstractFloat}

    N = length(x)
    if N == 1
        return (x[1], 0.0)
    end
    i′   = 1 + floor(Int, (1.0 - c) * (N - 1))
    n    = N - i′ + 1
    xₙ   = x[N]
    Mₙ₋₁ = (n-2)*sₙ₋₁^2
    x̄ₙ   = x̄ₙ₋₁ + (xₙ-x̄ₙ₋₁)/n
    Mₙ   = Mₙ₋₁ + (xₙ-x̄ₙ)*(xₙ-x̄ₙ₋₁)
    i    = 1 + floor(Int, (1.0 - c) * N)
    if i != i′
        n  = N - i + 1
        x₀ = x[i′]
        x̄′ = x̄ₙ
        x̄ₙ = x̄′ - (x₀-x̄′)/n
        Mₙ = Mₙ - (x₀-x̄ₙ)*(x₀-x̄′)
    end
    if n>1
        sₙ = sqrt(Mₙ/(n-1))
    else
        sₙ = 0.0
    end

    return (x̄ₙ, sₙ)
end

end

max_stars_repo_name: cohensbw/Langevin

max_issues_repo_name: cohensbw/ElPhDynamics

max_forks_repo_name: cohensbw/ElPhDynamics


text:
using BinDeps, Compat

@BinDeps.setup

nettle = library_dependency("nettle", aliases = ["libnettle"], runtime = false)
gnutls = library_dependency("gnutls", aliases = ["libgnutls.so.28","libgnutls","libgnutls28", "libgnutls-28", "libgnutls-deb0"], depends = [nettle], validate = function(p,h)
	if !haskey(ENV,"GNUTLS_VERSION")
		return true
	end
	ccall(dlsym(h,:gnutls_check_version),Ptr{UInt8},(Ptr{UInt8},),ENV["GNUTLS_VERSION"]) != C_NULL
end)

provides(Sources, @compat Dict(
	URI("http://www.lysator.liu.se/~nisse/archive/nettle-2.7.tar.gz") => nettle,
	URI("ftp://ftp.gnutls.org/gcrypt/gnutls/v3.2/gnutls-3.2.14.tar.xz") => gnutls))

#provides(Binaries,URI("ftp://ftp.gnutls.org/gcrypt/gnutls/w32/gnutls-3.2.1-w32.zip"),gnutls,os = :Windows)

@windows_only begin
    using WinRPM
    provides(WinRPM.RPM,"nettle",nettle,os = :Windows)
    provides(WinRPM.RPM,"gnutls",gnutls,os = :Windows)
end

if haskey(ENV,"GNUTLS_VERSION")
	requested_version = convert(VersionNumber,ENV["GNUTLS_VERSION"])
	pkgmanager_validate	= (p,dep)->available_version(p) >= requested_version
else
	pkgmanager_validate = true
end

@osx_only begin
	if Pkg.installed("Homebrew") === nothing
		error("Homebrew package not installed, please run Pkg.add(\"Homebrew\")")
	end
	using Homebrew
	provides( Homebrew.HB, "gnutls", gnutls, os = :Darwin )
end

provides(AptGet,"libgnutls28",gnutls,validate = pkgmanager_validate)
provides(AptGet,"libgnutls-deb0-28",gnutls,validate = pkgmanager_validate)
provides(Yum,"libgnutls",gnutls,validate = pkgmanager_validate)
provides(Yum,"gnutls",gnutls,validate = pkgmanager_validate)

julia_usrdir = normpath(JULIA_HOME*"/../") # This is a stopgap, we need a better builtin solution to get the included libraries
libdirs = AbstractString["$(julia_usrdir)/lib"]
includedirs = AbstractString["$(julia_usrdir)/include"]

env = @compat Dict("HOGWEED_LIBS" => "-L$(libdirs[1]) -L$(BinDeps.libdir(nettle)) -lhogweed -lgmp",
		"NETTLE_LIBS" => "-L$(libdirs[1]) -L$(BinDeps.libdir(nettle)) -lnettle -lgmp",
		"LIBS" => "-lgmp ", "LD_LIBRARY_PATH" => join([libdirs[1];BinDeps.libdir(nettle)],":"))

provides(BuildProcess,Autotools(lib_dirs = libdirs, include_dirs = includedirs, env = env),nettle)

# If we're installing gnutls from source we better also installl nettle from source, otherwise we end up with a giant mess.
provides(BuildProcess,Autotools(libtarget = "lib/libgnutls.la", lib_dirs = libdirs, include_dirs = includedirs,
	env = env),gnutls,force_depends = @compat Dict(BuildProcess => nettle))

@BinDeps.install @compat Dict(:gnutls => :gnutls)


text:
using LabView0mqJl
using Test

@testset "LabView0mqJl.jl" begin
    # Write your tests here.

    # Here are the tests:
    @test LabView0mqJl.tmp_test(0) == 0
    @test LabView0mqJl.tmp_test(3) == 3
    @test LabView0mqJl.tmp_test(5) == 5
end


text:
include("shared.jl")
include("../plotting/gnuplot_utils.jl")

BLAS.set_num_threads(1)

# Small neighborhood of the MLE.
const RANGE_GRID  = range(0.25, 5.0,  length=80)
const SMOOTH_GRID = range(1.2,  1.5,  length=70)

function gensurf()
  out = zeros(length(RANGE_GRID), length(SMOOTH_GRID))
  Threads.@threads for j in eachindex(RANGE_GRID)
    @inbounds for k in eachindex(SMOOTH_GRID)
      out[j,k] = profile_nll_replicates((RANGE_GRID[j], SMOOTH_GRID[k]),
                                        PTS_SURF, (SIM_SURF,))
    end
  end
  out .- minimum(out)
end

# Will need to figure out the right color scale here.
if !isinteractive()
  const surf = gensurf()
  gnuplot_save_matrix!("../plotdata/profile_surface.csv", surf,
                       RANGE_GRID, SMOOTH_GRID)
end


text:
#####################################################################################
#                                Pairwse Distance                                   #
#####################################################################################
using NearestNeighbors, StaticArrays, LinearAlgebra
export min_pairwise_distance, orthonormal

# min_pairwise_distance contributed by Kristoffer Carlsson
"""
    min_pairwise_distance(data) -> (min_pair, min_dist)
Calculate the minimum pairwise distance in the data (`Matrix`, `Vector{Vector}` or
`Dataset`). Return the index pair
of the datapoints that have the minimum distance, as well as its value.
"""
function min_pairwise_distance(cts::AbstractMatrix)
    if size(cts, 1) > size(cts, 2)
        error("Points must be close (transpose the Matrix)")
    end
    tree = KDTree(cts)
    min_d = Inf
    min_pair = (0, 0)
    for p in 1:size(cts, 2)
        inds, dists = knn(tree, view(cts, :, p), 1, false, i -> i == p)
        ind, dist = inds[1], dists[1]
        if dist < min_d
            min_d = dist
            min_pair = (p, ind)
        end
    end
    return min_pair, min_d
end

min_pairwise_distance(d::AbstractDataset) = min_pairwise_distance(d.data)

function min_pairwise_distance(
    pts::Vector{SVector{D,T}}) where {D,T<:Real}
    tree = KDTree(pts)
    min_d = eltype(pts[1])(Inf)
    min_pair = (0, 0)
    for p in 1:length(pts)
        inds, dists = knn(tree, pts[p], 1, false, i -> i == p)
        ind, dist = inds[1], dists[1]
        if dist < min_d
            min_d = dist
            min_pair = (p, ind)
        end
    end
    return min_pair, min_d
end

#####################################################################################
#                                Conversions                                        #
#####################################################################################
to_matrix(a::AbstractVector{<:AbstractVector}) = cat(2, a...)
to_matrix(a::AbstractMatrix) = a
function to_Smatrix(m)
    M = to_matrix(m)
    a, b = size(M)
    return SMatrix{a, b}(M)
end
to_vectorSvector(a::AbstractVector{<:SVector}) = a
function to_vectorSvector(a::AbstractMatrix)
    S = eltype(a)
    D, k = size(a)
    ws = Vector{SVector{D, S}}(k)
    for i in 1:k
        ws[i] = SVector{D, S}(a[:, i])
    end
    return ws
end

"""
    orthonormal(D, k) -> ws
Return a matrix `ws` with `k` columns, each being
an `D`-dimensional orthonormal vector.

Always returns `SMatrix` for stability reasons.
"""
function orthonormal end

@inline function orthonormal(D::Int, k::Int)
    k > D && throw(ArgumentError("k must be ≤ D"))
    q = qr(rand(SMatrix{D, k})).Q
end


text:
using LTC
using Plots
gr()
using BenchmarkTools
using DiffEqSensitivity
using OrdinaryDiffEq
using DiffEqFlux
using GalacticOptim
using ModelingToolkit
import Flux: Data.DataLoader

function generate_data()
    in_features = 2
    out_features = 1
    N = 48
    data_x = [sin.(range(0,stop=3π,length=N)), cos.(range(0,stop=3π,length=N))]
    data_x = [reshape([Float32(data_x[1][i]),Float32(data_x[2][i])],2,1) for i in 1:N]# |> f32
    data_y = [reshape([Float32(y)],1,1) for y in sin.(range(0,stop=6π,length=N))]# |> f32
    DataLoader((data_x, data_y), batchsize=N)
end

function data(iter; data_x=nothing, data_y=nothing, short=false, noisy=false)
    ncycle(generate_data(), iter)
end

function train_sine(n, solver=VCABM(), sensealg=InterpolatingAdjoint(autojacvec=ReverseDiffVJP(true)))

  cbg = function (p,l,pred,y;doplot=true)
    display(l)
    if doplot
      fig = plot([ŷ[end,1] for ŷ in pred], label="ŷ")
      plot!(fig, [yi[end,1] for yi in y], label="y")
      #frame(anim)
      display(fig)
    end
    return false
  end

  #train_dl = generate_data()
  batchsize = 1

  wiring = LTC.FWiring(2,1)
  net = LTC.Net(wiring, name=:net)
  sys = ModelingToolkit.structural_simplify(net)

  model = DiffEqFlux.FastChain(LTC.Mapper(wiring.n_in),
                               LTC.RecurMTK(LTC.MTKCell(wiring.n_in, wiring.n_out, net, sys, solver, sensealg)),
                               LTC.Mapper(wiring.n_out),
                               )

  train_dl = data(n)
  opt = GalacticOptim.Flux.Optimiser(ClipValue(0.5), ADAM(0.01))
  # opt = Optim.LBFGS()
  # opt = BBO()
  # opt = ParticleSwarm(;lower=lb, upper=ub)
  # opt = Fminbox(GradientDescent())
  AD = GalacticOptim.AutoZygote()
  # AD = GalacticOptim.AutoModelingToolkit()

  # return model

  LTC.optimize(model, LTC.loss_seq, cbg, opt, AD, train_dl)

end

@time model = train_sine(100)
# @time model = train_sine(100, Tsit5(), InterpolatingAdjoint(autojacvec=DiffEqSensitivity.EnzymeVJP()))
# @time traintest(1000, QNDF())
# @time traintest(1000, TRBDF2())
# @time traintest(1000, AutoTsit5(Rosenbrock23()))


text:
using DSGE
using HDF5, Base.Test, Distributions
include("../../util.jl")

path = dirname(@__FILE__)

### Model
model = Model990()

### Parameters

# Parameters match parameters, bounds, etc. vectors from reference (ϵ = 1e-4)
para = zeros(82)
bounds = zeros(82, 2)
pshape = zeros(82)
pmean = zeros(82)
pstdd = zeros(82)
trspec = zeros(82, 4)

# keys to skip (used to be fixed_parameters)
fixed_parameters = [:δ, :λ_w, :ϵ_p, :ϵ_w, :g_star]

# not all parameters appear in model.parameters
i = 1
for θ in model.parameters
    !isa(θ,AbstractParameter) && error()
    in(θ.key, fixed_parameters) && continue

    para[i] = θ.value

    (left, right) = θ.valuebounds
    bounds[i, 1] = left
    bounds[i, 2] = right

    prior = θ.prior.value

    if isa(prior, DSGE.RootInverseGamma)
        pshape[i] = 4
        (ν, τ) = params(prior)
        pmean[i] = τ
        pstdd[i] = ν
    else
        if isa(prior, Distributions.Beta)
            pshape[i] = 1
        elseif isa(prior, Distributions.Gamma)
            pshape[i] = 2
        elseif isa(prior, Distributions.Normal)
            pshape[i] = 3
        end
        pmean[i] = mean(prior)
        pstdd[i] = std(prior)

    end

    if θ.transform == DSGE.Untransformed()
        trspec[i, 1] = 0
    elseif θ.transform == DSGE.SquareRoot()
        trspec[i, 1] = 1
    elseif  θ.transform == DSGE.Exponential()()
        trspec[i, 1] = 2
    else
       throw(error("This kind of transform not allowed"))
    end

    (left, right) = θ.transform_parameterization
    trspec[i, 2] = left
    trspec[i, 3] = right
    if θ == model[:Iendoα]
        trspec[i, 4] = 0
    else
        trspec[i, 4] = 1
    end

    i += 1
end

### Model indices

# Endogenous states
endo = model.endogenous_states
@test length(endo) == 60
@test endo[:Ez_t] == 60

# Exogenous shocks
exo = model.exogenous_shocks
@test length(exo) == 16
@test exo[:corepce_sh] == 16

# Expectation shocks
ex = model.expected_shocks
@test length(ex) == 13
@test ex[:Erk_f_sh] == 13

# Equations
eq = model.equilibrium_conditions
@test length(eq) == 60
@test eq[:eq_Ez] == 60

# Additional states
endo_new = model.endogenous_states_augmented
@test length(endo_new) == 12
@test endo_new[:y_t1] == 61

# Observables
obs = model.observables
@test length(obs) == 12
@test obs[:obs_tfp] == 12

### Equilibrium conditions
Γ0, Γ1, C, Ψ, Π = eqcond(model)

# Matrices are of expected dimensions
@test size(Γ0) == (60, 60)
@test size(Γ1) == (60, 60)
@test size(C) == (60, 1)
@test size(Ψ) == (60, 16)
@test size(Π) == (60, 13)

# Check output matrices against reference output (ϵ = 1e-4)
h5 = h5open("$path/eqcond.h5")
Γ0_ref = read(h5, "G0")
Γ1_ref = read(h5, "G1")
C_ref  = reshape(read(h5, "C"), 60, 1)
Ψ_ref  = read(h5, "PSI")
Π_ref  = read(h5, "PIE")
close(h5)

@test_matrix_approx_eq Γ0_ref Γ0
@test_matrix_approx_eq Γ1_ref Γ1
@test_matrix_approx_eq C_ref C
@test_matrix_approx_eq Ψ_ref Ψ
@test_matrix_approx_eq Π_ref Π

# ### Measurement equation
expect = Dict{Symbol, Matrix}()
h5 = h5open("$path/measurement.h5")
expect[:ZZ] = read(h5, "ZZ")
expect[:DD] = reshape(read(h5, "DD"), 12, 1)
expect[:QQ] = read(h5, "QQ")
expect[:EE] = read(h5, "EE")
expect[:MM]  = read(h5, "MM")
close(h5)

model = Model990()
TTT, RRR, CCC = solve(model)
actual = measurement(model, TTT, RRR, CCC)
for d in (:ZZ, :DD, :QQ, :EE, :MM)
    @test_matrix_approx_eq expect[d] actual[d]
end

nothing