text:
#=
  The BUGS Book
  Example 2.5.1, "The how many trick"

  Cf ~/jags/bugs_book_2_5_1.jags
  """
  model {
    for (i in 1:20) {
     Y[i] ~ dgamma(4, 0.04)
    }
    cum[1] <- Y[1]
    for (i in 2:20) {
      cum[i] <- cum[i - 1] + Y[i]
    }
    for (i in 1:20) {
      cum.step[i] <- i*step(1000 - cum[i])
      x[i] <- ifelse(cum[i] < 1001, 1, 0) #
    }
    # number <- ranked(cum.step[], 20) # maximum number in cum.step ## NOTE: ranked is not in JAGS
    number <- sort(cum.step[]) # maximum number in cum.step
    check <- equals(cum.step[20], 0) # always 1 if I=20 big enough
  }

  Output:
           Mean       SD  Naive SE Time-series SE
x[1]  1.000e+00 0.000000 0.000e+00      0.000e+00
x[2]  1.000e+00 0.000000 0.000e+00      0.000e+00
x[3]  1.000e+00 0.000000 0.000e+00      0.000e+00
x[4]  1.000e+00 0.000000 0.000e+00      0.000e+00
x[5]  9.999e-01 0.012247 5.000e-05      5.000e-05
x[6]  9.976e-01 0.049101 2.005e-04      2.043e-04
x[7]  9.810e-01 0.136350 5.566e-04      5.566e-04
x[8]  9.141e-01 0.280293 1.144e-03      1.131e-03
x[9]  7.578e-01 0.428418 1.749e-03      1.741e-03
x[10] 5.222e-01 0.499512 2.039e-03      2.039e-03
x[11] 2.825e-01 0.450235 1.838e-03      1.793e-03
x[12] 1.195e-01 0.324340 1.324e-03      1.307e-03
x[13] 4.015e-02 0.196313 8.014e-04      8.014e-04
x[14] 9.867e-03 0.098841 4.035e-04      4.035e-04
x[15] 2.117e-03 0.045959 1.876e-04      1.897e-04
x[16] 4.500e-04 0.021209 8.658e-05      8.988e-05
x[17] 1.667e-05 0.004082 1.667e-05      1.667e-05
x[18] 0.000e+00 0.000000 0.000e+00      0.000e+00
x[19] 0.000e+00 0.000000 0.000e+00      0.000e+00
x[20] 0.000e+00 0.000000 0.000e+00      0.000e+00
  """

  Cf ~/blog/bugs_book_2_5_1.blog
     ~/webppl/bugs_book_2_5_1.wppl

   Distribution of ix:
12.00000 =>    8898  (0.222450)
11.00000 =>    8239  (0.205975)
13.00000 =>    7365  (0.184125)
10.00000 =>    5209  (0.130225)
14.00000 =>    4494  (0.112350)
9.00000 =>    2373  (0.059325)
15.00000 =>    1829  (0.045725)
8.00000 =>     767  (0.019175)
16.00000 =>     495  (0.012375)
7.00000 =>     157  (0.003925)
17.00000 =>     123  (0.003075)
18.00000 =>      29  (0.000725)
6.00000 =>      18  (0.000450)
19.00000 =>       3  (0.000075)
5.00000 =>       1  (0.000025)

Distribution of x (the change point arra")
[1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]    =>  8968 (0.2242)
[1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]    =>  8020 (0.2005)
[1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]    =>  7477 (0.186925)
[1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]    =>  5173 (0.129325)
[1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]    =>  4491 (0.112275)
[1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]    =>  2387 (0.059675)
[1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0]    =>  1833 (0.045825)
[1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]    =>  753 (0.018825)
[1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0]    =>  549 (0.013725)
[1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]    =>  168 (0.0042)
[1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0]    =>  128 (0.0032)
[1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]    =>  26 (0.00065)
[1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0]    =>  22 (0.00055)
[1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0]    =>  5 (0.000125)

=#

using Turing, StatsPlots, DataFrames
include("jl_utils.jl")

@model function bugs_book_2_5_1()
    n = 20

    y = tzeros(n)
    cum = tzeros(n)
    cum_step = tzeros(n)
    x = tzeros(n)
    for i in 1:n
        # Note:BLOG has second parameter b, webbpl 1/b
        y[i] ~ Gamma(3.0,1/0.04)
        if i == 1
            cum[i] = y[1]+1
        else
            cum[i] = cum[i-1] + y[i] + i
        end
        if cum[i] > 1000.0
            cum_step[i] = i*cum[i]
        else
            cum_step[i] = 0.0
        end

        # A change point
        # x[i] ~ flip(0.5)
        if cum[i] < 1001.0
            x[i] = 1.0
        else
            x[i] = 0.0
        end
    end

    # number ~ Uniform(0,50_000)
    number = maximum(cum_step)

    # Max number with 1
    # ix ~ DiscreteUniform(1,20)
    ix, _ = findmax([i*(x[i]==1) for i in 1:n])

    # return number
    return ix
    # return x
end

model = bugs_book_2_5_1()

num_chains = 4

# chains = sample(model, Prior(), MCMCThreads(), 10_000, num_chains)
chains = sample(model, MH(), MCMCThreads(), 10_000, num_chains)
# chains = sample(model, PG(15), MCMCThreads(), 10_000, num_chains)
# chains = sample(model, SMC(1000), MCMCThreads(), 10_000, num_chains)
# chains = sample(model, IS(), MCMCThreads(), 10_000, num_chains)

display(chains)
# display(plot(chains))

gen = generated_quantities(model, chains)
show_var_dist_pct(gen, 40)

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
##ExoplanetsSysSim/examples/hsu_et_al/process_abc_output.jl
##(c) 2017 Keir A. Ashby
# takes the output from abc_run_christiansen.jl and makes relevant plots
# 1. make_plots sends the output into whatever plotting functions we choose to call
# 2. generation_plot makes a plot for the mean, median and overall distribution of distances or values of theta over each generation of abc. saves each plot

using PyPlot
using ExoplanetsSysSim
include(joinpath(Pkg.dir(), "processabcoutput", "src", "constants.jl"))

#sends the output into two plotting functions
function make_plots(output::ABC.abc_population_type, ss_true::ExoplanetsSysSim.CatalogSummaryStatistics, run::Int64 = 1) #run indicates how many times abc_run_christiansen.jl has generated an output in the event that there are multiple outputs.
  num_thetas = length(output.accept_log.theta[1])
#  gen_start = output.accept_log.generation_starts_at
  run = generation_plot(output.accept_log, output.accept_log.dist, run) 	#plots changing distances over each generation
  for i in 1:num_thetas
    run = generation_plot(output.accept_log, output.accept_log.theta, run, i) 	#plots changing theta values over each generation
  end
end

#takes a log from the output and returns the mean, median, and overall distribution of either distance or a theta parameter over every generation of ABC
function generation_plot(log::ABC.abc_log_type, log_part::Array, run::Int64=1, m::Int64=1) #user can decide which part of theta she wants by specifying m in the input
  if typeof(log_part[1]) == Array{Float64,1} 					#first check whether the user has input an array of distances or of theta
    my_theta = zeros(length(log_part))
    for i in 1:length(my_theta)
      my_theta[i] = log_part[i][m]
    end
    log_part = my_theta
    log_name = "Theta_$m";
    run +=10*m; 								#we keep different figures separated by factors of 10 so as not to overlap
  else
    log_name = "Distance";
  end

  avg_list = zeros(length(log.generation_starts_at)) 				#array holding the average of log_part over each generation of abc
  med_list = zeros(length(log.generation_starts_at)) 				#array holding the median of log_part over each generation of abc
  all_lengths = zeros(length(log.generation_starts_at)) 			#array holding the number of log_part values in each generation of abc
  for i in 1:length(log.generation_starts_at)
    n = 1
    if i == length(log.generation_starts_at)
    gen_list = zeros(length(log_part)+1-log.generation_starts_at[i]) 		#array holding each value of log_part in a given generation
	for j in log.generation_starts_at[i]:length(log_part)
	  gen_list[n] = log_part[j]
	  n +=1
	end
    else
    gen_list = zeros(log.generation_starts_at[i+1]-log.generation_starts_at[i]) #array holding each value of log_part in a given generation
	for j in log.generation_starts_at[i]:log.generation_starts_at[i+1]-1
	  gen_list[n] = log_part[j]
 	  n +=1
	end
    end
    avg_list[i] = mean(gen_list) 						#take the average of log_part over this generation
    med_list[i] = median(gen_list) 						#take the median of log_part over this generation
    all_lengths[i] = length(gen_list) 						#take the number of values of log_part for this generation
  end
  all_gen = ones(1) 								#array holding the generation number corresponding to every value of log_part
  for k in 1:length(log.generation_starts_at)
    this_gen = ones(convert(Int64, all_lengths[k])).*k
    append!(all_gen, this_gen)
  end
  shift!(all_gen)

  figure(run)
  semilogy(1:length(log.generation_starts_at), avg_list, "r") 			#plot of average log_part over each generation of ABC
  title(string(log_name, " vs. Generation"), fontsize=18)
  xlabel("Generation", fontsize=14)
  ylabel(log_name, fontsize=20)
  tick_params(labelsize=13, width=2)

  semilogy(1:length(log.generation_starts_at), med_list, "g") 			#plot of median log_part over each generation of ABC

  xmax= convert(Float64, length(log.generation_starts_at))+0.2 			#just so we can see everything
  semilogy(all_gen, log_part, ".")
  axis([0.8, xmax, 0, maximum(log_part)+0.01])

  savefig(string("plot_", log_name, ".png"))
  return(run)
end

ext: jl


text:
n = rand(1:10)
@test matrixdepot("tridiag", [1,1,1], [1,1,1,1], [1,1,1]) == matrixdepot("tridiag", 4, 1,1,1)
@test matrixdepot("tridiag", Float64, n) == matrixdepot("tridiag", n)

A = matrixdepot("tridiag", n)
@test isposdef(full(A))
@test issym(full(A))
println("'tridiag' passed test...")


text:
module MxArrays

export MxArray, mxarray, jvalue, matlab_homepath, matlab_library_path, open_matlab_library

import Base.eltype, Base.close, Base.size, Base.copy, Base.ndims, Base.convert

include("mxbase.jl")
include("mxarray.jl")

end


text:
import Base: sleep
using Dates: now

"""
    BlockingScheduler(; clock=real_time_clock, delayfunc=_sleep, jobconfig=JobConfig())

`BlockingScheduler` is the simplest scheduler.
It implements `AbstractScheduler`.

This is a monothread implementation of scheduling job.

# Optional arguments
- `clock::AbstractClock`: clock that will be used by scheduler (it's by default `real_time_clock`, which is system UTC time but a `SimClock` struct can also be passed for simulation purpose).
- `delayfunc::DelayFunc`: functor which is responsible (when called) of waiting until next task should be fired (`_sleep` is used by default but a `NoSleep` struct can also be passed for simulation purpose).
- `jobconfig::JobConfig`: job configuration default settings (`misfire_grace_period`...)
"""
mutable struct BlockingScheduler <: AbstractScheduler
    id::String
    clock::AbstractClock
    delayfunc::DelayFunc
    jobstore::AbstractJobStore
    executor::AbstractExecutor
    jobconfig::JobConfig
    running::Bool

    function BlockingScheduler(; clock=real_time_clock, delayfunc=_sleep, jobconfig=JobConfig())
        id = get_scheduler_id()
        jobstore = MemoryJobStore()
        executor = DebugExecutor()
        jobconfig = JobConfig()
        running = true
        now_ = now(clock)
        println("[*] $now_ Init scheduler $id")
        new(id, clock, delayfunc, jobstore, executor, jobconfig, running)
    end
end

"""
    isrunning(sched)

Return `true` when a scheduler named `sched` is running.
"""
isrunning(sched::BlockingScheduler) = sched.running

"""
    isstopped(sched)

Return `true` when a scheduler named `sched` is stopped.
"""
isstopped(sched::BlockingScheduler) = !isrunning(sched)

"""
    JobStore(sched)

Return `jobstore` of a scheduler named `sched`.
"""
JobStore(sched::BlockingScheduler) = sched.jobstore

"""
    Executor(sched)

Return `executor` of a scheduler named `sched`.
"""
Executor(sched::BlockingScheduler) = sched.executor

"""
    sleep(sched, args...; kwargs...)

Block the scheduler for a specified delay.
"""
sleep(sched::BlockingScheduler, args...; kwargs...) = sched.delayfunc(args...; kwargs...)

"""
    run_pending(sched)

Run pending tasks of a scheduler `sched`.

This function should be called instead of `run` when using scheduler in simulation mode.
"""
function run_pending(sched::BlockingScheduler)
    jobstore = JobStore(sched)
    executor = Executor(sched)
    if length(jobstore) > 0
        job = peek(jobstore)
        now_ = now(sched.clock)
        if job.dt_next_fire > now_
            duration = job.dt_next_fire - now_
            println("[ ] $now_ Sleep $duration until job $(job.id) execution")
            sleep(sched, duration)
        else  # job.dt_next_fire <= now_
            job = dequeue!(jobstore)
            job.n_triggered += 1
            if now_ - job.dt_next_fire <= job.config.misfire_grace_period
                println("[*] $now_ RUNNING $(job.id)")
                dump(job)
                run(executor, job)
                sleep(sched, 0)
            else
                println("[x] $now_ Can't execute job $(job.id) (beyond misfire grace period)")
            end
            update(jobstore, job, now_)
        end
    end
end

"""
    run(sched)

Run (in a blocking loop) a scheduler named `sched`.
"""
function run(sched::BlockingScheduler)
    now_ = now(sched.clock)
    println("[*] $now_ Running scheduler $(sched.id)")
    while(isrunning(sched))
        if isempty(JobStore(sched))
            shutdown(sched)
        else
            run_pending(sched)
        end
    end
    now_ = now(sched.clock)
    println("[*] $now_ Scheduler $(sched.id) stopped")
end

"""
    add(sched, action, trigger; name=DEFAULT_JOB_NAME, priority=DEFAULT_PRIORITY)

Schedule when an `Action` named `action` should be triggered (according `trigger`).
"""
function add(sched::AbstractScheduler, action::Action, trigger::AbstractTrigger; name=DEFAULT_JOB_NAME, priority=DEFAULT_PRIORITY)
    if isstopped(sched)
        error("Job can't be add to a scheduler which is stopped")
    end
    jobstore = JobStore(sched)
    job_id = get_job_id(jobstore)
    dt_created = now(sched.clock)
    dt_updated = dt_created
    dt_next_fire = DateTime(0)
    n_triggered = 0
    job = Job(job_id, action, trigger, name, priority, dt_created, dt_updated, dt_next_fire, n_triggered, sched.jobconfig)
    push!(jobstore, job)
end

function add(sched::AbstractScheduler, action::Action, trigger; name=DEFAULT_JOB_NAME, priority=DEFAULT_PRIORITY, n=-1)
    add(sched, action, Trigger(trigger, n=n); name=name, priority=priority)
end

"""
    shutdown(sched)

Shutdown scheduler `sched`.
"""
function shutdown(sched::BlockingScheduler)
    sched.running = false
    now_ = now(sched.clock)
    println("[*] $now_ Stopping scheduler $(sched.id)")
end


text:
# Methods wrapping the api, allocating buffers, converting strings, etc.

### GET ATTRIBUTES METHODS
###################################################################################
#ViStatus niScope_GetAttributeViInt32 (ViSession vi, ViConstString channelList, ViAttr attributeID, ViInt32* value);
function niScope_GetAttributeViInt32(obj::T, ch="", attr=1) where {T<:NIScope}
	chlist = [UInt8.(collect(ch)); UInt8(0)] # terminate with NULL char
	res    = Int32(0)
	pres   = ViPInt32(res)
	status = niScope_GetAttributeViInt32(obj.handle, chlist, Int32(attr), pres)
	return (status, niScope_GetError(obj), pres[])
end

#ViStatus niScope_GetAttributeViInt64 (ViSession vi, ViConstString channelList, ViAttr attributeID, ViInt64* value);
function niScope_GetAttributeViInt64(obj::T, ch="", attr=1) where {T<:NIScope}
	chlist = [UInt8.(collect(ch)); UInt8(0)] # terminate with NULL char
	res    = Int64(0)
	pres   = ViPInt64(res)
	status = niScope_GetAttributeViInt64(obj.handle, chlist, Int32(attr), pres)
	return (status, niScope_GetError(obj), pres[])
end

#ViStatus niScope_GetAttributeViReal64 (ViSession vi, ViConstString channelList, ViAttr attributeID, ViReal64* value);
function niScope_GetAttributeViReal64(obj::T, ch="", attr=1) where {T<:NIScope}
	chlist = [UInt8.(collect(ch)); UInt8(0)] # terminate with NULL char
	res    = Float64(0.0)
	pres   = Ref(res)
	status = niScope_GetAttributeViReal64(obj.handle, chlist, Int32(attr), pres)
	return (status, niScope_GetError(obj), pres[])
end

#ViStatus niScope_GetAttributeViString (ViSession vi, ViConstString channelList, ViAttr attributeID, ViInt32 bufSize, ViChar value[]);
function niScope_GetAttributeViString(obj::T, ch="", attr=1) where {T<:NIScope}
	chlist = [UInt8.(collect(ch)); UInt8(0)] # terminate with NULL char
	buffsz = Int32(256)
	res    = zeros(UInt8,buffsz)
	pres   = Ref(res)
	status = niScope_GetAttributeViString(obj.handle, chlist, Int32(attr), buffsz, res)
	return (status, niScope_GetError(obj), pres[])
end

#ViStatus niScope_GetAttributeViBoolean (ViSession vi, ViConstString channelList, ViAttr attributeID, ViBoolean* value);
function niScope_GetAttributeViBoolean(obj::T, ch="", attr=1) where {T<:NIScope}
	chlist = [UInt8.(collect(ch)); UInt8(0)] # terminate with NULL char
	res    = ViBoolean(false)
	pres   = Ref(res)
	status = niScope_GetAttributeViBoolean(obj.handle, chlist, ViAttr(attr), pres)
	return (status, niScope_GetError(obj), pres[])
end

#ViStatus niScope_GetAttributeViSession (ViSession vi, ViConstString channelList, ViAttr attributeID, ViSession* value);
function niScope_GetAttributeViSession(obj::T, ch="", attr=1) where {T<:NIScope}
	chlist = [UInt8.(collect(ch)); UInt8(0)] # terminate with NULL char
	res    = ViSession(0)
	pres   = Ref(res)
	status = niScope_GetAttributeViSession(obj.handle, chlist, Int32(attr), pres)
	return (status, niScope_GetError(obj), pres[])
end
###################################################################################

### SET ATTRIBUTES METHODS
###################################################################################
#ViStatus niScope_SetAttributeViInt32 (ViSession vi, ViConstString channelList, ViAttr attributeID, ViInt32 value);
function niScope_SetAttributeViInt32(obj::T, ch="", attr=1, val=1) where {T<:NIScope}
	chlist = [UInt8.(collect(ch)); UInt8(0)] # terminate with NULL char
	status = niScope_SetAttributeViInt32(obj.handle, chlist, Int32(attr), Int32(val))
	return (status, niScope_GetError(obj))
end

#ViStatus niScope_SetAttributeViInt64 (ViSession vi, ViConstString channelList, ViAttr attributeID, ViInt64 value);
function niScope_SetAttributeViInt64(obj::T, ch="", attr=1, val=1) where {T<:NIScope}
	chlist = [UInt8.(collect(ch)); UInt8(0)] # terminate with NULL char
	status = niScope_SetAttributeViInt64(obj.handle, chlist, Int32(attr), Int64(val))
	return (status, niScope_GetError(obj))
end

#ViStatus niScope_SetAttributeViReal64 (ViSession vi, ViConstString channelList, ViAttr attributeID, ViReal64 value);
function niScope_SetAttributeViReal64(obj::T, ch="", attr=1, val=1.0) where {T<:NIScope}
	chlist = [UInt8.(collect(ch)); UInt8(0)] # terminate with NULL char
	status = niScope_SetAttributeViReal64(obj.handle, chlist, Int32(attr), Float64(val))
	return (status, niScope_GetError(obj))
end

#ViStatus niScope_SetAttributeViString (ViSession vi, ViConstString channelList, ViAttr attributeID, ViInt32 bufSize, ViChar value[]);
function niScope_SetAttributeViString(obj::T, ch="", attr=1, val=1) where {T<:NIScope}
	chlist = [UInt8.(collect(ch)); UInt8(0)] # terminate with NULL char
	val    = [UInt8.(collect(val)); UInt8(0)]
	status = niScope_SetAttributeViString(obj.handle, chlist, Int32(attr), val)
	return (status, niScope_GetError(obj))
end

#ViStatus niScope_SetAttributeViBoolean (ViSession vi, ViConstString channelList, ViAttr attributeID, ViBoolean value);
function niScope_SetAttributeViBoolean(obj::T, ch="", attr=1, val=false) where {T<:NIScope}
	chlist = [UInt8.(collect(ch)); UInt8(0)] # terminate with NULL char
	status = niScope_SetAttributeViBoolean(obj.handle, chlist, Int32(attr), Bool(val))
	return (status, niScope_GetError(obj))
end

#ViStatus niScope_SetAttributeViSession (ViSession vi, ViConstString channelList, ViAttr attributeID, ViSession value);
function niScope_SetAttributeViSession(obj::T, ch="", attr=1, val=1) where {T<:NIScope}
	chlist = [UInt8.(collect(ch)); UInt8(0)] # terminate with NULL char
	status = niScope_SetAttributeViSession(obj.handle, chlist, Int32(attr), UInt32(val))
	return (status, niScope_GetError(obj))
end
###################################################################################

## GET ERRORS
###################################################################################
#ViStatus _stdcall niScope_GetError ( ViSession vi , ViStatus * errorCode , ViInt32 bufferSize , ViChar description []);
function niScope_GetError(obj)
	sts    = Int32(0)
	psts   = Ref(sts)
	buffsz = Int32(256)
	res    = zeros(UInt8,buffsz)
	status = niScope_GetError(obj.handle, psts, buffsz, res)
	res    = rstrip(String(Char.(res)),['\r', '\n','\0','0'])
	if !isempty(res) # this is brittle
		error(res)
	end
end


text:
@doc raw"""
    function rejection(f, log_l :: Float64, params...)

Samples from the prior with the hard likelihood constraint ``\log L_k > `` `log_l`.

Args:
+ `f` stochastic function. Must have signature `f(t :: Trace, params...)`
+ `log_l`: current log likelihood threshold
+ `params`: any additional arguments to pass to `f`
"""
function rejection(f, log_l :: Float64, params...)
    this_log_l = log_l
    local new_t
    while this_log_l <= log_l
        new_t = trace()
        f(new_t, params...)
        this_log_l = loglikelihood(new_t)
    end
    (new_t, this_log_l)
end

export rejection


text:
# ===========================================================
# Concrete overapproximation with a ball in the infinity norm
# ===========================================================

"""
    ballinf_approximation(S::LazySet)

Overapproximate a set by a tight ball in the infinity norm.

### Input

- `S` -- set

### Output

A tight ball in the infinity norm.

### Algorithm

The center and radius of the box are obtained by evaluating the support function
of the given convex set along the canonical directions.
"""
function ballinf_approximation(S::LazySet{N}) where {N}
    n = dim(S)
    c = Vector{N}(undef, n)
    r = zero(N)
    d = zeros(N, n)

    @inbounds for i in 1:n
        d[i] = one(N)
        htop = ρ(d, S)
        d[i] = -one(N)
        hbottom = -ρ(d, S)
        d[i] = zero(N)
        c[i] = (htop + hbottom) / 2
        rcur = (htop - hbottom) / 2
        if (rcur > r)
            r = rcur
        elseif rcur < 0
            # contradicting bounds => set is empty
            return EmptySet{N}(dim(S))
        end
    end
    return BallInf(c, r)
end

# ===============
# Specializations
# ===============

# empty set specialization
ballinf_approximation(∅::EmptySet) = ∅


text:
using Heuristics, StaticArrays, BenchmarkTools

# Objective function
include("./../test/testProblems.jl")

function main()

    # Size of problem
    N = 2
    M = 200

    # Define problem: N is number of dims
    LB = -100 .* ones(N)
    UB =  100 .* ones(N)
    prob = Problem{N}(rastriginfunc, LB, UB)
    opts = Options(display=true)

    # Initialize Optimizer: M is swarm size!
    optimizer = PSO{M}(prob; initMethod = :LogisticsMap)

    # solve
    #Heuristics.initialize!(optimizer, opts)
    #@btime Heuristics.iterate!($optimizer, $opts)
    optimize!(optimizer, opts)

    return nothing
end

main()


text:
module TestStructView

using UnsafeFields.Demos: StructView
using Test

struct Values2{A,B}
    a::A
    b::B
end

@testset "Tuple{Int, Int}" begin
    A = [(111, 222), (333, 444), (555, 666)]
    mut = StructView(A, 1)
    @test (mut[1], mut[2]) === A[1] === (111, 222)
    mut[1] = 3
    mut[2] = 4
    @test (mut[1], mut[2]) === A[1] === (3, 4)
end

@testset "Values2{String,Symbol}" begin
    A = [Values2("initial string", :initial_symbol)]
    load() = (A[1].a, A[1].b)
    mut = StructView(A, 1)
    @test (mut.a, mut.b) === load() === ("initial string", :initial_symbol)
    mut.a = "second string"
    mut.b = :second_symbol
    @test (mut.a, mut.b) === load() === ("second string", :second_symbol)
end

@testset "Values2{Union{String,Int,Nothing},Union{Int,Nothing,Float64}}" begin
    A = [
        Values2{Union{String,Int,Nothing},Union{Int,Nothing,Float64}}(nothing, nothing)
    ]
    load() = (A[1].a, A[1].b)
    mut = StructView(A, 1)
    @test (mut.a, mut.b) === load() === (nothing, nothing)
    mut.a = "second string"
    mut.b = 222
    @test (mut.a, mut.b) === load() === ("second string", 222)
    mut.a = 111
    mut.b = 333.0
    @test (mut.a, mut.b) === load() === (111, 333.0)
end

@testset "Values2{Any,Union{Integer,Symbol}}" begin
    A = [Values2{Any,Union{Integer,Symbol}}(0, 0)]
    load() = (A[1].a, A[1].b)
    mut = StructView(A, 1)
    @test (mut.a, mut.b) === load() === (0, 0)
    mut.a = "second string"
    mut.b = :second_symbol
    @test (mut.a, mut.b) === load() == ("second string", :second_symbol)
end

end  # module


text:
using Flatten, Unitful, BenchmarkTools, FieldMetadata, Test
import Flatten: flattenable

struct Foo{A,B,C}
    a::A
    b::B
    c::C
end

struct Nest{F,B,C}
    nb::B
    nf::F
    nc::C
end

mutable struct MuFoo{A,B,C}
    a::A
    b::B
    c::C
end

mutable struct MuNest{F,B,C}
    nf::F
    nb::B
    nc::C
end

foo = Foo(1.0, 2.0, 3.0)
nest = Nest(Foo(1,2,3), 4.0, 5.0f0)
nesttuple = Nest((foo, nest), 9, 10)

# Flatten to specific types
@test flatten(nest, Int) == (1, 2, 3)
@test flatten(nest, AbstractFloat) == (4.0, 5.0f0)
@test flatten(nest, Number, Float32) == (1, 2, 3, 4.0)

@test flatten(Foo(1, 2, "3"), String, Nothing) === ("3",)

# Test flattening

@test flatten(Foo(1,2,3)) == (1,2,3)
@test flatten(Foo(1,2,3)) == (1,2,3)
@test flatten(((1,2,3), (4,5))) == (1,2,3,4,5)
@test flatten(Nest(Foo(1,2,3),4,5)) == (1,2,3,4,5)
@test flatten((Nest(Foo(1,2,3),4,5), Nest(Foo(6,7,8), 9, 10))) == (1,2,3,4,5,6,7,8,9,10)
@test flatten(Nest(Foo(1,2,3), (4,5), (6,7))) == (1,2,3,4,5,6,7)
@test flatten((Nest(Foo(1,2,3),4.0,5.0), Nest(Foo(6,7,8), 9, 10))) == (1,2,3,4.0,5.0,6,7,8,9,10)

# Test reconstruction
@test typeof( reconstruct(foo, (1.0, 2.0, 3.0))) == typeof(foo)
@test typeof(reconstruct(nest, flatten(nest))) == typeof(nest)

@test flatten(reconstruct(foo, flatten(foo))) == flatten(foo)
@test flatten(reconstruct(nest, flatten(nest))) == flatten(nest)
@test flatten(reconstruct((nest, nest), flatten((nest, nest)))) == flatten((nest, nest))
@test flatten(reconstruct(nesttuple, flatten(nesttuple))) == flatten(nesttuple)

@test flatten(reconstruct(foo, flatten(foo, Real)), Real) == flatten(foo, Real)
foo2 = Foo(1, "two", :three)
@test flatten(reconstruct(foo2, flatten(foo2, String, Real), String, Real), String, Real) == flatten(foo2, String, Real)

# Test updating mutable structs
mufoo = MuFoo(1.0, 2.0, 3.0)

munest = MuNest(MuFoo(1,2,3), 4.0, 5.0)
@test flatten(update!(mufoo, flatten(mufoo) .* 7)) == (7.0, 14.0, 21.0)
@test flatten(update!(munest, flatten(munest) .* 7)) == (7, 14, 21, 28.0, 35.0)

munesttuple = MuNest((MuFoo(1.0, 2.0, 3.0), MuNest(MuFoo(1,2,3), 4.0, 5.0)), 9, 10)
@test flatten(update!(munesttuple, flatten(munesttuple) .* 7)) == (7.0, 14.0, 21.0, 7, 14, 21, 28.0, 35.0, 63, 70)

# Test retyping
@test typeof(reconstruct(foo, round.(Int, flatten(foo))).a) == Int
@test flatten(reconstruct(nesttuple, round.(Int, flatten(nesttuple)))) == round.(Int, flatten(nesttuple))

# Partial fields with @flattenable

@metadata foobar :nobar

@flattenable @foobar struct Partial{A,B,C}
    " Field a"
    a::A | :foo | true
    " Field b"
    b::B | :bar | true
    " Field c"
    c::C | :foo | false
end
partial = Partial(1.0, 2.0, 3.0)
nestedpartial = Partial(Partial(1.0, 2.0, 3.0), 4, 5)
@test flatten(nestedpartial) === (1.0, 2.0, 4)

@test flatten(reconstruct(nestedpartial, flatten(nestedpartial))) == flatten(nestedpartial)

# Meta flattening
@test metaflatten(partial, foobar) == (:foo, :bar)
@test metaflatten(nestedpartial, foobar) == (:foo, :bar, :bar)
@test metaflatten((nestedpartial, partial), foobar) == (:foo, :bar, :bar, :foo, :bar)

# Helpers
@test fieldnameflatten(nestedpartial) == (:a, :b, :b)
@test fieldtypeflatten((nestedpartial, partial)) == (Float64, Float64, Int, Float64, Float64)
@test parenttypeflatten(nestedpartial) ==
    (Partial{Float64,Float64,Float64}, Partial{Float64,Float64,Float64},
     Partial{Partial{Float64,Float64,Float64},Int,Int})
@test parentnameflatten(nestedpartial) == (:Partial, :Partial, :Partial)

# Partial fields with custom field traits

@metadata flattenable2 true

@reflattenable2 struct Partial
    a::A | false
    b::B | false
    c::C | true
end

@test flatten(nestedpartial, flattenable2) === (5,)
@test flatten(reconstruct(nestedpartial, flatten(nestedpartial, flattenable2), flattenable2), flattenable2) == flatten(nestedpartial, flattenable2)

# Updating metadata updates flattened fields
@reflattenable @refoobar struct Partial
    a | :bar | false
    b | :bar | false
    c | :foo | true
end

@reflattenable mutable struct MuFoo
    a | false
    b | true
    c | false
end

@reflattenable mutable struct MuNest
    nf | true
    nb | true
    nc | false
end

@test flatten(nestedpartial) == (5,)
@test flatten(reconstruct(nestedpartial, flatten(nestedpartial))) == flatten(nestedpartial)
@inferred flatten(reconstruct(nestedpartial, flatten(nestedpartial)))

@test metaflatten(foo, flattenable) == (true, true, true)
@test metaflatten(nest, flattenable) == (true, true, true, true, true)
@test metaflatten(partial, foobar) == (:foo,)
@test metaflatten(nestedpartial, foobar) == (:foo,)

@test fieldnameflatten(foo) == (:a, :b, :c)
@test fieldnameflatten(nest) == (:a, :b, :c, :nf, :nc)
@test fieldnameflatten(nestedpartial) == (:c,)

# Updating mutables
mufoo = MuFoo(1.0, 2.0, 3.0)
munest = MuNest(MuFoo(1,2,3), 4.0, 5.0)
@test flatten(update!(mufoo, flatten(mufoo) .* 7)) == (14.0,)
@test flatten(update!(munest, flatten(munest) .* 7)) == (14, 28.0)

# Test non-parametric types
mutable struct AnyPoint
    x
    y
end
anypoint = AnyPoint(1,2)
@test flatten(anypoint) == (1,2)
@test flatten(reconstruct(anypoint, (1,2))) == (1,2)

# With void
nestvoid = Nest(Foo(1,2,3), nothing, nothing)
munestvoid = MuNest(MuFoo(1,2,3), nothing, nothing)
@test flatten(nestvoid) == (1,2,3)
@test flatten((Nest(Foo(1,2,3), nothing, nothing), Nest(Foo(nothing, nothing, nothing), 9, 10))) == (1,2,3,9,10)
@test flatten(reconstruct(nestvoid, flatten(nestvoid))) == flatten(nestvoid)
@test flatten(update!(munestvoid, flatten(munestvoid))) == flatten(munestvoid)

# Test unit stripping functions

ufoo = Foo(1.0u"m", 2.0u"g", 3.0)
unest = Nest(Foo(1u"m",2u"g",3), 4.0, 5.0u"kPa")
unesttuple = Nest((ufoo, unest), 9, 10u"mol*L^-1")

@test flatten(ufoo, Real) === (1.0,2.0,3.0)
@test flatten(unest, Real) === (1,2,3,4.0,5.0)

@test flatten(reconstruct(ufoo, flatten(ufoo, Real) , Real), Real) == flatten(ufoo, Real)

@test flatten(reconstruct(ufoo, flatten(ufoo), Number), Number) === flatten(ufoo)
@test flatten(reconstruct(ufoo, (1u"m",2u"g",3.0), Number), Number) === (1u"m",2u"g",3.0)

##############################################################################
# Benchmarks

struct BenchFoo{A,B,C}
    a::A
    b::B
    c::C
end

foo = BenchFoo(1.0, 2.0, 3.0)

b = @benchmark flatten($foo, $flattenable, Float64, Nothing)
@test b.allocs == 0

@test reconstruct(foo, (1.0, 2.0, 3.0)) == foo
b = @benchmark reconstruct($foo, (1.0, 2.0, 3.0))
@test b.allocs == 0

v = [1.0, 2.0, 3.0]
@test reconstruct(foo, v) == foo
b = @benchmark reconstruct($foo, $v)
@test b.allocs == 0


text:
module ArtifactInstallation
using Pkg.Artifacts, Test, Libdl
export do_test

function do_test()
    # First, check that `"HelloWorldC"` is installed automatically
    artifacts_toml = joinpath(@__DIR__, "..", "Artifacts.toml")
    @test artifact_exists(artifact_hash("HelloWorldC", artifacts_toml))

    # Test that we can use `artifact""` to get at HelloWorldC
    hello_world_exe = joinpath(artifact"HelloWorldC", "bin", "hello_world")
    if Sys.iswindows()
        hello_world_exe = "$(hello_world_exe).exe"
    end
    @test isfile(hello_world_exe)

    # Test that we can use a variable, not just a literal:
    hello_world = "HelloWorldC"
    hello_world_exe = joinpath(@artifact_str(hello_world), "bin", "hello_world")
    if Sys.iswindows()
        hello_world_exe = "$(hello_world_exe).exe"
    end
    @test isfile(hello_world_exe)
end

end


text:
function topologicalorder(network::AbstractNetwork{T}) where T<:Integer
    nnodes = numnodes(network)
    nodeorder = zeros(T, nnodes)
    deg = indegree.(Ref(network), 1:nnodes)

    S = filter((i) -> deg[i] == 0, 1:nnodes)
    if length(S) == 0
        @warn "No nodes with no incoming links. Cannot define a topological order!"
    elseif length(S) > 1
        @warn "More than one nodes with zero indegree. Topological order not unique!"
    end
    o = 0
    while length(S) > 0
        i = pop!(S)
        o += 1
        nodeorder[i] = o
        for j in outneighbors(network, i)
            deg[j] -= 1
            if deg[j] == 0
                push!(S, j)
            end
        end
    end
    if (o != nnodes)
        @warn "Network is not acyclic. Incomplete topological order!"
    end

    nodeorder
end

orderednodes(network::AbstractNetwork) = invperm(topologicalorder(network))


text:
export LocalCF, LocalCFMsg, LocalController

const LocalCFMsg = NaiveMsg

"""
A controller framework that only performs local delay compensation without
the knowledge of other agent's controller.
"""
struct LocalCF{N,X,Z,U,S,XDy,ZDy} <: ControllerFramework{X,Z,U,NaiveMsg{X,Z},Nothing}
    central::CentralControl{U,S}
    world_model::WorldDynamics{N,XDy,ZDy}
    delay_model::DelayModel
    "Whether to conpensate communication delay by assuming constant future action."
    conpensate_comm::Bool
end

function LocalCF(
    central::CentralControl{U,S}, world_model::WorldDynamics{N,XDy,ZDy}, delay_model;
    X, Z, conpensate_comm
) where {N,U,S,XDy,ZDy}
    LocalCF{N,X,Z,U,S,XDy,ZDy}(central, world_model, delay_model, conpensate_comm)
end

struct LocalController{id,N,X,Z,U,S,XDy,ZDy} <: Controller{X,Z,U,LocalCFMsg,Nothing}
    cf::LocalCF{N,X,Z,U,S,XDy,ZDy}
    "Central control state."
    s_c::Ref{S}
    τ::Ref{𝕋}
    u::Ref{U}                           # t = τ-1
    u_history::TimedQueue{U}            # t ∈ [τ-Tx,τ+Tu-1]
    t0::𝕋
end

function LocalController(
    id,cf::LocalCF{N,X,Z,U,S,XDy,ZDy}, s_c, t, u0, u_history
) where {N,X,Z,U,S,XDy,ZDy}
    LocalController{id,N,X,Z,U,S,XDy,ZDy}(cf, Ref(s_c), Ref(t), Ref(u0), u_history, t)
end

function OneVision.control!(
    ctrl::LocalController{id,N,X,Z,U},
    x::X,
    z::Z,
    msgs::Each{LocalCFMsg{X,Z}},
)::Tuple{U,Each{LocalCFMsg{X,Z}}} where {id,N,X,Z,U}
    @unpack central, world_model, delay_model, conpensate_comm = ctrl.cf
    @unpack Tx, Tu, Tc, ΔT = short_delay_names(delay_model)
    x_dy, z_dy = world_model.dynamics[id], world_model.obs_dynamics[id]
    τ = ctrl.τ[]
    xt, zt = Timed.(τ-Tx, (x,z))

    push_iter(head, tail) = ((i == 0 ? head : tail[i]) for i in 0:length(tail)-1)

    x_history = self_estimate(x_dy,xt,ctrl.u_history.queue)
    z_history = self_z_estimate(z_dy, zt, push_iter(xt, x_history))

    x_self, z_self = isempty(ctrl.u_history) ? (xt, zt) : (x_history[end], z_history[end])
    msg = LocalCFMsg(x_self[τ+Tu], z_self[τ+Tu])
    if mod(τ - ctrl.t0, ΔT) == 0
        msgs[id] = msg
        xs = [m.x for m in msgs]
        zs = [m.z for m in msgs]
        ctrl.u[] = control_one(central, ctrl.s_c[], xs, zs, τ + Tu, id)
    end
    u = ctrl.u[]
    if conpensate_comm
        x_future = self_estimate(x_dy, x_self, Timed(τ + Tu + t - 1, u) for t in 1:Tc)
        z_future = self_z_estimate(z_dy, z_self, push_iter(x_self, x_future))
        t_rec = τ + Tu + Tc
        msg = LocalCFMsg(x_future[end][t_rec], z_future[end][t_rec])
    end
    msgs′ = fill(msg, N)
    pushpop!(ctrl.u_history, Timed(τ + Tu, u))
    ctrl.τ[] += 1
    u, msgs′
end

function OneVision.make_controllers(
    cf::LocalCF{N,X,Z,U,S},
    init_status::Each{Tuple{X,Z,U}},
    t0::𝕋,
)::Tuple where {N,X,Z,U,S}
    @unpack Tx, Tu, ΔT = short_delay_names(cf.delay_model)
    function mk_controller(id)
        x0, z0, u0 = init_status[id]
        u_history = TimedQueue(u0, t0-Tx, t0+Tu-1)
        ideal_s = init_state(cf.central, t0)

        LocalController(id, cf, ideal_s, t0, u0, u_history)
    end
    function mk_q(_)
        receives = [LocalCFMsg(x, z) for (x, z, u) in init_status]
        constant_queue(receives, msg_queue_length(cf.delay_model))
    end

    ctrls = ntuple(mk_controller, N)
    msg_qs = ntuple(mk_q, N)
    ctrls, msg_qs
end

# LocalCF has nothing to log
OneVision.write_logs(::LocalController) = Dict{𝕋,Nothing}()


text:
using Luxor
function dragon(turtle::Turtle, level=4, size=200, direction=45)
    if level != 0
        Turn(turtle, -direction)
        dragon(turtle, level-1, size/sqrt(2), 45)
        Turn(turtle, direction*2)
        dragon(turtle, level-1, size/sqrt(2), -45)
        Turn(turtle, -direction)
    else
        Forward(turtle, size)
    end
end

Drawing(900, 500, "./Dragon.png")
t = Turtle(300, 300, true, 0, (0., 0.0, 0.0));
dragon(t, 10,400)
finish()
preview()

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# mgc2b converts mel generalized cesptrum to MGLSADF filter coefficients.
# TODO: consider `otype` optional argument?
function mgc2b!(mgc::AbstractVector, α, γ)
    b = mgc
    mc2b!(b, α)

    # when gamma = 0, mel-generalized cespstrum corresponds to mel cepstrum
    γ == zero(γ) && return mgc

    gnorm!(b, γ)

    # scale by gamma
    b[1] = log(b[1])
    for i=2:length(b)
        @inbounds b[i] *= γ
    end

    b
end

mgc2b(mgc::AbstractVector, α=0.35, γ=0.0) = mgc2b!(copy(mgc), α, γ)

function mgc2b!{T<:MelGeneralizedCepstrum}(state::SpectralParamState{T})
    assert_not_ready_to_filt(state)

    @assert has_loggain(state)
    γ = glog_gamma(paramdef(state))
    if γ == zero(γ)
        @assert gain_normalized(state)
    else
        @assert !gain_normalized(state)
    end

    def = paramdef(state)
    α = allpass_alpha(def)
    γ = glog_gamma(def)
    data = rawdata(state)
    state.data = mgc2b!(data, α, γ)
    state.has_loggain = true
    state.gain_normalized = true
    state.ready_to_filt = true
    state
end

# mgc is assumed not to be gamma scaled (TODO: handle this explicitly?)
function mgc2b{T<:MelGeneralizedCepstrum}(state::SpectralParamState{T})
    mgc2b!(copy(state))
end


text:
struct WeightSchedulePolicy <: AbstractWeightAdjustmentPolicy
    weights::Vector{Float64}
end

possible_weights(wsp::WeightSchedulePolicy) = wsp.weights

function (wsp::WeightSchedulePolicy)(;num_solutions, kwargs...)
    w = wsp.weights[min(num_solutions + 1, length(wsp.weights))]
    return w
end

function awastar_search_with_scheduled_weights(search_prob::AbstractSearchProblem, weights::AbstractArray{T}, walltime_limit::Real, nodes_budget::Integer) where T <: Real
    wsp = WeightSchedulePolicy(weights)
    return awastar_search_with_dynamic_weights(search_prob, wsp, walltime_limit, nodes_budget)
end


text:
function main()
    coords = Set()
    first_fold = true
    for line in readlines()
        if ',' in line
            x, y = split(line, ',') .|> x -> parse(Int, x)
            push!(coords, (y, x))
        elseif '=' in line
            fold_is_x = occursin("x=", line)
            fold_val = parse(Int, split(line, '=')[2])
            new_coords = Set()
            for (y, x) in coords
                fold_x = fold_is_x ? min(0, fold_val - x) * 2 : 0
                fold_y = fold_is_x ? 0 : min(0, fold_val - y) * 2
                push!(new_coords, (y + fold_y, x + fold_x))
            end
            coords = new_coords
            first_fold && length(coords) |> println
            first_fold = false
        end
    end
    y_max, x_max = maximum(coords) .+ 1
    field = [collect(" " ^ x_max) for _ in 1:y_max]
    for (y, x) in coords
        field[y + 1][x + 1] = '#'
    end
    field .|> join .|> println
end

main()

max_stars_repo_name: LiquidFun/adventofcode

max_issues_repo_name: LiquidFun/adventofcode

max_forks_repo_name: LiquidFun/adventofcode


text:
using LinAlgBigFloat.VectorSpaces
using Random

Random.seed!(1)

l, n, m = 2, 6, 4
A = randn(n, m)
B = randn(l, n)

vz = ZeroSpace(n)
@test rank(vz) == 0 && dim(vz) == n

vf = VectorSpace(n)
@test rank(vf) == n && dim(vf) == n

va = VectorSpace(A)
@test rank(va) == min(m, n) && dim(va) == n
vb = VectorSpace(B)
@test rank(vb) == min(l, n) && dim(vb) == l

@test va == image(A)
@test vb == image(B)

@test vz ⊆ va ⊆ vf

@test kernel(A) == ZeroSpace(m)
@test image(B) == VectorSpace(l)

@test kernel(A) == preimage(A)

max_stars_repo_name: KlausC/LinAlgBigFloat

max_issues_repo_name: KlausC/LinAlgBigFloat

max_forks_repo_name: KlausC/LinAlgBigFloat


text:
mutable struct SpacecraftSim{DEST}
    # Initial Epoch (TAI Seconds)
    t0::Float64

    # Initial Co-States Defining Primer Vector Control
    λ0::Vector{Float64}

    # Time span of trajectory
    ts::Tuple{Float64, Float64}

    # Indirect optimal control spacecraft parameters
    ps::IndirectTrajOpt.CR3BPIndirectParams

    # DifferentialEquations.jl solution object (providing interpolated traj.)
    interp::DEST

    # Constructor
    function SpacecraftSim(t0)
        # Initializing time span, initial co-states, and parameters
        λ0      = [ 5.963713007227692, 12.999481348176722, 1.6584130970663207,
                    -0.04062357612068792, 0.018558730104232632, -0.0004993971529784666, 0.12200533102974148]
        ps      = initCR3BPIndirectParams("Low Thrust 10 CR3BP"); ps.ϵ = 0.0
        ts      = (0.0, 33.1) # [days]
        tsND    = (0.0, ts[2] * 86400 / ps.crp.TU) # [n.d. time units]

        # Boundary conditions and initial state/co-state vector
        μ   = ps.crp.μ; LU = ps.crp.LU; TU = ps.crp.TU
        ics = [-0.019488511458668, -0.016033479812051, 0.0, 8.918881923678198, -4.081793688818725, 0.0, 1.0]
        fcs = [(1.0 - μ) + 6060.483/LU, 19452.284/LU, -34982.968/LU,
                0.082677*(TU/LU), 0.006820*(TU/LU), -0.368434*(TU/LU), 0.0]
        y0  = @SVector [ics[1], ics[2], ics[3], ics[4], ics[5], ics[6], ics[7],
                        λ0[1],  λ0[2],  λ0[3],  λ0[4],  λ0[5],  λ0[6],  λ0[7]]

        # Integrate state and co-state differential equations
        sol = integrate(y0, tsND, ps, CR3BP(), FullSolutionHistory(), MEMF())

        # Create spacecraft sim object
        new{typeof(sol)}(t0, λ0, ts, ps, sol)
    end
end

function GetStateAndControl(ss::SpacecraftSim, t::Real; frame = :inertial, derivs = false)
    # Check that t ∈ t0 <= t <= t0 + ts[2]
    if t < ss.t0 || t > ss.t0 + ss.ts[2]*86400
        throw(ArgumentError("t does not fall within tranjectory time span."))
    end

    # Compute non-dimentioanl time
    tnd     = (t - ss.t0) / (ss.ps.crp.TU)

    # Get state in synodic reference frame
    ysyn    = ss.interp(tnd)

    # Scale physical states to units of km, km/s, and kg
    xsyns   = @SVector [(ysyn[1] + ss.ps.crp.μ)*ss.ps.crp.LU, ysyn[2]*ss.ps.crp.LU, ysyn[3]*ss.ps.crp.LU,
                        ysyn[4]*ss.ps.crp.VU, ysyn[5]*ss.ps.crp.VU, ysyn[6]*ss.ps.crp.VU, ysyn[7]*ss.ps.crp.MU]

    # Rotate to inertial if desired
    xs = zeros(7)
    if frame == :inertial
        # Create rotation matrix to inertial frame
        ω       = 1.0              # [rad / n.d. time]
        ωs      = ω / ss.ps.crp.TU  # [rad / second]
        θ       = -ωs*(ss.t0 + ss.ts[2]*86400 - t)    # [rad]

        if derivs == false
            RsI     = @SMatrix [cos(θ) -sin(θ) 0;
                                sin(θ)  cos(θ) 0;
                                0       0      1]
            dRsI    = @SMatrix [-ωs*sin(θ) -ωs*cos(θ) 0;
                                ωs*cos(θ) -ωs*sin(θ) 0;
                                0         0        0]

            # Rotate to inertial frame
            xs[7] = xsyns[7]
            @views mul!(xs[1:3], RsI,  xsyns[1:3])
            @views mul!(xs[4:6], RsI,  xsyns[4:6])
            @views mul!(xs[4:6], dRsI, xsyns[1:3], 1.0, 1.0)
        else
            AU      = ss.ps.crp.VU / ss.ps.crp.TU
            dysyn   = ss.interp(tnd, Val{1})
            dxsyns  = @SVector [dysyn[1]*ss.ps.crp.VU, dysyn[2]*ss.ps.crp.VU, dysyn[3]*ss.ps.crp.VU,
                                dysyn[4]*AU, dysyn[5]*AU, dysyn[6]*AU, dysyn[7]*ss.ps.crp.MU / ss.ps.crp.TU]

            RsI     = @SMatrix [cos(θ) -sin(θ) 0;
                                sin(θ)  cos(θ) 0;
                                0       0      1]
            dRsI    = @SMatrix [-ωs*sin(θ) -ωs*cos(θ) 0;
                                ωs*cos(θ) -ωs*sin(θ) 0;
                                0         0        0]
            ddRsI   = @SMatrix [-ωs^2*cos(θ)  ωs^2*sin(θ) 0;
                                -ωs^2*sin(θ) -ωs^2*cos(θ) 0;
                                0            0            0]

            # Rotate to inertial frame
            xs[7] = dxsyns[7]
            @views mul!(xs[1:3], RsI,   dxsyns[1:3])
            @views mul!(xs[1:3], dRsI,  xsyns[1:3], 1.0, 1.0)
            @views mul!(xs[4:6], RsI,   dxsyns[4:6])
            @views mul!(xs[4:6], dRsI,  dxsyns[1:3], 2.0, 1.0)
            @views mul!(xs[4:6], ddRsI, xsyns[1:3], 1.0, 1.0)
        end
    else
        xs = xsyns
    end

    # Compute controls
    us      = zeros(4);

    # Scaled exaust velocity
    cs      = ss.ps.sp.c * ss.ps.crp.TU / (ss.ps.crp.LU * 1000.0)

    # Switching function
    λv      = norm(view(ysyn, 11:13))
    S       = IndirectTrajOpt.computeS(ysyn, λv, cs)

    # Throttling factor
    us[1]   = S > 0.0 ? 0.0 : 1.0

    # Thrust direction
    ussyn   = @SVector [-ysyn[11] / λv, -ysyn[12] / λv, -ysyn[13] / λv]

    # Rotate to necisary if desired
    if frame == :inertial
        @views mul!(us[2:4], RsI, ussyn)
    else
        us[2:4] .= ussyn
    end

    # Return state and control
    return (xs, us)
end

function GetStateAndControl(ss::SpacecraftSim, ts::AbstractVector; frame = :inertial)
    xIs = zeros(length(ts), 7)
    us  = zeros(length(ts), 4)
    for i in 1:length(ts)
        # Get time
        t = ts[i]

        # Get state and control
        sc = GetStateAndControl(ss, t; frame = frame)

        # Place in data matricies
        xIs[i,:] .= sc[1]
        us[i,:]  .= sc[2]
    end

    return (xIs, us)
end

function GetControl(ss::SpacecraftSim, t::Real; frame = :inertial)
    # Check that t ∈ t0 <= t <= t0 + ts[2]
    if t < ss.t0 || t > ss.t0 + ss.ts[2]*86400
        throw(ArgumentError("t does not fall within tranjectory time span."))
    end

    # Compute non-dimentioanl time
    tnd     = (t - ss.t0) / ss.ps.crp.TU

    # Get state in synodic reference frame
    ysyn    = ss.interp(tnd)

    # Compute controls
    us      = zeros(4);

    # Scaled exaust velocity
    cs      = ss.ps.sp.c * ss.ps.crp.TU / (ss.ps.crp.LU * 1000.0)

    # Switching function
    λv      = norm(view(ysyn, 11:13))
    S       = IndirectTrajOpt.computeS(ysyn, λv, cs)

    # Throttling factor
    us[1]   = S > 0.0 ? 0.0 : 1.0

    # Thrust direction
    ussyn   = @SVector [-ysyn[11] / λv, -ysyn[12] / λv, -ysyn[13] / λv]

    # Rotate to necisary if desired
    if frame == :inertial
        # Create rotation matrix to inertial frame
        ω       = -1.0              # [rad / n.d. time]
        ωs      = ω / ss.ps.crp.TU  # [rad / second]
        θ       = ωs*(ss.t0 + ss.ts[2]*86400 - t)    # [rad]
        RsI     = @SMatrix [cos(θ) -sin(θ) 0;
                            sin(θ)  cos(θ) 0;
                            0       0      1]

        @views mul!(us[2:4], RsI, ussyn)
    else
        us[2:4] .= ussyn
    end

    # Return state and control
    return us
end

function GetControlAndLunaPos(ss::SpacecraftSim, t::Real; frame = :inertial)
    # Check that t ∈ t0 <= t <= t0 + ts[2]
    if t < ss.t0 || t > ss.t0 + ss.ts[2]*86400
        throw(ArgumentError("t does not fall within tranjectory time span."))
    end

    # Compute non-dimentioanl time
    tnd     = (t - ss.t0) / ss.ps.crp.TU

    # Get state in synodic reference frame
    ysyn    = ss.interp(tnd)

    # Compute controls
    us      = zeros(4);

    # Scaled exaust velocity
    cs      = ss.ps.sp.c * ss.ps.crp.TU / (ss.ps.crp.LU * 1000.0)

    # Switching function
    λv      = norm(view(ysyn, 11:13))
    S       = IndirectTrajOpt.computeS(ysyn, λv, cs)

    # Throttling factor
    us[1]   = S > 0.0 ? 0.0 : 1.0

    # Thrust direction
    ussyn   = @SVector [-ysyn[11] / λv, -ysyn[12] / λv, -ysyn[13] / λv]

    # Luna position
    rLunaSyn    = @SVector [ss.ps.crp.LU, 0.0, 0.0]
    rLuna       = zeros(3)

    # Rotate to necisary if desired
    if frame == :inertial
        # Create rotation matrix to inertial frame
        ω       = -1.0              # [rad / n.d. time]
        ωs      = ω / ss.ps.crp.TU  # [rad / second]
        θ       = ωs*(ss.t0 + ss.ts[2]*86400 - t)    # [rad]
        RsI     = @SMatrix [cos(θ) -sin(θ) 0;
                            sin(θ)  cos(θ) 0;
                            0       0      1]

        @views mul!(us[2:4], RsI, ussyn)
        mul!(rLuna, RsI, rLunaSyn)
    else
        us[2:4] .= ussyn
        rLuna .= rLunaSyn
    end

    # Return state and control
    return (us, rLuna)
end

function PlotTrajectory(ss::SpacecraftSim, ts; frame = :inertial)
    # Get states and controls
    (xs,us) = GetStateAndControl(ss, ts; frame = frame)

    # Get target trajectory
    xt = generateTargetTrajectory()

    # Frame String
    if frame == :inertial
        fs = "i"
    else
        fs = "s"
    end

    # Plot with MATLAB
    mat"""
        % Bulk Data
        ts = $ts;
        xs = $xs;
        us = $us;
        xh = $xt;

        % Split into thrust and coast
        xc = nan(size(xs));
        xt = nan(size(xs));
        ut = nan(size(us));
        lu = us(1,1);
        for i = 1:length(ts)
            if us(i,1) == 1.0 || (us(i,1) == 0.0 && lu == 1.0)
                xt(i,:) = xs(i,:);
                ut(i,:) = us(i,:);
            end
            if us(i,1) == 0.0 || (us(i,1) == 1.0 && lu == 0.0)
                xc(i,:) = xs(i,:);
            end
            lu = us(i,1);
        end

        if $fs == "i"
            figure()
            sp1 = subplot(2,1,1);
            plot(xt(:,1), xt(:,2), 'r')
            hold on
            plot(xc(:,1), xc(:,2), 'b')
            plot(xh(:,1), xh(:,2), '--k')
            ylabel("Y, km", "Interpreter", "latex")

            % Plot Earth
            scatter(0.0, 0.0, "xk")

            % Plot Moon
            scatter($(ss.ps.crp.LU), 0.0, "xk")

            %axis equal
            grid on

            sp2 = subplot(2,1,2);
            plot(xt(:,1), xt(:,3), 'r')
            hold on
            plot(xc(:,1), xc(:,3), 'b')
            plot(xh(:,1), xh(:,3), '--k')
            ylabel("Z, km", "Interpreter", "latex")

            % Plot Earth
            scatter(0.0, 0.0, "xk");

            % Plot Moon
            scatter($(ss.ps.crp.LU), 0.0, "xk")

            %axis equal
            grid on

            %leg1 = legend("Thrusting Arc", "Coasting Arc", "Target NRHO");
            xlabel("X, km", "Interpreter", "latex")

            %linkaxes([sp1,sp2,'x'])
        else
            figure()
            plot(xt(:,1), xt(:,2), 'r')
            hold on
            plot(xc(:,1), xc(:,2), 'b')
            plot(xh(:,1), xh(:,2), '--k')
            ylabel("Y, km", "Interpreter", "latex")
            xlabel("X, km", "Interpreter", "latex")

            % Plot Earth
            scatter(0.0, 0.0, "xk");

            % Plot Moon
            scatter($(ss.ps.crp.LU), 0.0, "xk")

            axis equal
            grid on
        end

        set(gca, "fontname", "Times New Roman", "fontsize", 10)
        %set(leg1, "fontname", "Times New Roman", "fontsize", 10, "Interpreter", "latex")
        set(gcf, "PaperUnits", "inches", "PaperPosition", [0.25, 0.25, 5.0, 4.0])
        set(gcf, "PaperPositionMode", "Manual")
        print("./figures/" + $fs + "Traj.eps", "-depsc", "-r600")
    """
end

function PlotTrajectory(ss::SpacecraftSim; frame = :inertial)
    # Plot full trajectory
    ts = range(ss.t0; stop = ss.t0 + ss.ts[2]*86400, length = 10000)
    PlotTrajectory(ss, ts; frame = frame)
end

function generateTargetTrajectory()
    # Parameters
    ps = initCR3BPIndirectParams("Low Thrust 10 CR3BP")
    μ  = ps.crp.μ
    LU = ps.crp.LU
    TU = ps.crp.TU

    # Target final conditions
    xf = @SVector [(1.0 - μ) + 6060.483/LU, 19452.284/LU, -34982.968/LU, 0.082677*(TU/LU), 0.006820*(TU/LU), -0.368434*(TU/LU)]

    # Integrate target trajectory
    solTarg = integrate(xf, (0.0, 1.6), "Low Thrust 10 CR3BP", CR3BP(), FullSolutionHistoryNoControl())
    rt = zeros(length(solTarg.t), 3)
    for i in 1:length(solTarg.t)
        rt[i,1] = (solTarg.u[i][1] + μ)*LU
        rt[i,2] = solTarg.u[i][2]*LU
        rt[i,3] = solTarg.u[i][3]*LU
    end

    return rt
end

max_stars_repo_name: GrantHecht/OptimalEstimationProject


text:
export FreqInterval, logfreqi, logfreqp, freqi, freqp
export FreqIC, logfreqic, logfreqpc, freqic, freqpc

import Base: +, -, *, ==

# some helper constants
const oct = log(2)
const halfoct = oct/2

# FreqInterval
# ============

"""
    FreqInterval <: Interval

A frequency ratio interval.
Is internally represented by the logarithm of the frequency ratio,
so conversion to Float64 will return a log-frequency ratio!
"""
struct FreqInterval <: Interval
    lfr :: Float64
end

# constructors

"""
    logfreqi(logratio)

Creates a frequency ratio interval from a log-ratio.
"""
logfreqi(lfr) = FreqInterval(lfr)

"""
    logfreqp(logfrequency)

Creates a frequency pitch from a log-frequency.
"""
logfreqp(lf) = Pitch(FreqInterval(lf))

"""
    freqi(ratio)

Creates a frequency ratio interval from a frequency ratio.
"""
freqi(ratio) = FreqInterval(log(ratio))

"""
    freqp(frequency)

Creates a frequency pitch from a frequency.
"""
freqp(freq) = Pitch(FreqInterval(log(freq)))

# Base interface

Base.show(io::IO, i::FreqInterval) = print(io, string("fr", exp(i.lfr)))
Base.show(io::IO, p::Pitch{FreqInterval}) = print(io, string(exp(p.pitch.lfr), "Hz"))

Base.isless(i1::FreqInterval, i2::FreqInterval) =
    i1.lfr < i2.lfr

Base.isequal(i1::FreqInterval, i2::FreqInterval) =
    i1.lfr == i2.lfr

Base.isapprox(a::FreqInterval, b::FreqInterval; kwargs...) =
    Base.isapprox(a.lfr, b.lfr; kwargs...)

Base.hash(i::FreqInterval, x::UInt64) = hash(i.lfr, x)

Base.Float64(i::FreqInterval) = i.lfr
Base.Float64(p::Pitch{FreqInterval}) = p.pitch.lfr

convert(::Type{FreqInterval}, x::N) where {N<:Real} = logfreqi(Float64(x))
convert(::Type{N}, i::FreqInterval) where {N<:Real} = convert(N, i.lfr)
# TODO: conversion to other pitch types, in particular tomidi

# interval interface

+(i1::FreqInterval, i2::FreqInterval) = logfreqi(i1.lfr + i2.lfr)
-(i1::FreqInterval, i2::FreqInterval) = logfreqi(i1.lfr - i2.lfr)
-(i::FreqInterval) = logfreqi(-i.lfr)

*(i::FreqInterval, n::Integer) = logfreqi(n*i.lfr)
*(n::Integer, i::FreqInterval) = logfreqi(n*i.lfr)

Base.sign(i::FreqInterval) = sign(i.lfr)
Base.abs(i::FreqInterval) = logfreqi(abs(i.lfr))

Base.zero(::Type{FreqInterval}) = logfreqi(0.)
Base.zero(::FreqInterval) = logfreqi(0.)
octave(::Type{FreqInterval}) = logfreqi(oct)
chromsemi(::Type{FreqInterval}) = logfreqi(oct/12)
# isstep is left unimplemented

ic(i::FreqInterval) = logfreqic(i.lfr)
embed(i::FreqInterval) = i
intervaltype(::Type{FreqInterval}) = FreqInterval
intervalclasstype(::Type{FreqInterval}) = FreqIC

# FreqIC
# ======

"""
    FreqIC <: Interval

A frequency ratio interval class.
Is internally represented by the logarithm of the frequency ratio,
so conversion to Float64 will return a log-frequency ratio!
Maintains octave equivalence.
"""
struct FreqIC <: IntervalClass
    lfr :: Float64
end

# constructors

"""
    logfreqic(logratio)

Creates a frequency ratio interval class from a log-ratio.
"""
logfreqic(lfr) = FreqIC(mod(lfr, oct))

"""
    logfreqpc(logfrequency)

Creates a frequency pitch class from a log-frequency.
"""
logfreqpc(lf) = Pitch(logfreqic(lf))

"""
    freqic(ratio)

Creates a frequency ratio interval class from a frequency ratio.
"""
freqic(ratio) = logfreqic(log(ratio))

"""
    freqpc(frequency)

Creates a frequency pitch class from a frequency.
"""
freqpc(freq) = Pitch(freqic(freq))

# Base interface

Base.show(io::IO, i::FreqIC) = print(io, string("fr[", exp(i.lfr), "]"))
Base.show(io::IO, p::Pitch{FreqIC}) = print(io, string("[", exp(p.pitch.lfr), "]Hz"))

Base.isless(i1::FreqIC, i2::FreqIC) =
    i1.lfr < i2.lfr

Base.isequal(i1::FreqIC, i2::FreqIC) =
    i1.lfr == i2.lfr

Base.isapprox(a::FreqIC, b::FreqIC; kwargs...) =
    Base.isapprox(a.lfr, b.lfr; kwargs...) ||
    Base.isapprox(a.lfr + log(2), b.lfr; kwargs...) ||
    Base.isapprox(a.lfr, b.lfr + log(2); kwargs...)

Base.hash(i::FreqIC, x::UInt64) = hash(i.lfr, x)

Base.Float64(i::FreqIC) = i.lfr
Base.Float64(p::Pitch{FreqIC}) = p.pitch.lfr

convert(::Type{FreqIC}, x::N) where {N<:Real} = logfreqic(Float64(x))
convert(::Type{N}, i::FreqIC) where {N<:Real} = convert(N, i.lfr)
# TODO: conversion to other pitch types, in particular tomidi

# interval interface

+(i1::FreqIC, i2::FreqIC) = logfreqic(i1.lfr + i2.lfr)
-(i1::FreqIC, i2::FreqIC) = logfreqic(i1.lfr - i2.lfr)
-(i::FreqIC) = logfreqic(-i.lfr)

*(i::FreqIC, n::Integer) = logfreqic(n*i.lfr)
*(n::Integer, i::FreqIC) = logfreqic(n*i.lfr)

Base.sign(i::FreqIC) = i.lfr == 0 ? 0 : -sign(i.lfr - halfoct)
Base.abs(i::FreqIC) = logfreqic(abs(mod(i.lfr + halfoct, oct) - halfoct))

Base.zero(::Type{FreqIC}) = logfreqic(0)
Base.zero(::FreqIC) = logfreqic(0)
octave(::Type{FreqIC}) = logfreqic(0)
chromsemi(::Type{FreqIC}) = logfreqic(oct/12)
# isstep is left unimplemented

ic(i::FreqIC) = i
embed(i::FreqIC) = logfreqi(i.lfr)
intervaltype(::Type{FreqIC}) = FreqInterval
intervalclasstype(::Type{FreqIC}) = FreqIC


text:
function HolmesFeatherstoneAttractionModel(arg0::Frame, arg1::NormalizedSphericalHarmonicsProvider)
    return HolmesFeatherstoneAttractionModel((Frame, NormalizedSphericalHarmonicsProvider), arg0, arg1)
end

function acceleration(obj::HolmesFeatherstoneAttractionModel, arg0::FieldSpacecraftState, arg1::Vector{RealFieldElement})
    return jcall(obj, "acceleration", FieldVector3D, (FieldSpacecraftState, Vector{RealFieldElement}), arg0, arg1)
end

function acceleration(obj::HolmesFeatherstoneAttractionModel, arg0::SpacecraftState, arg1::Vector{jdouble})
    return jcall(obj, "acceleration", Vector3D, (SpacecraftState, Vector{jdouble}), arg0, arg1)
end

function add_contribution(obj::ForceModel, arg0::FieldSpacecraftState, arg1::FieldTimeDerivativesEquations)
    return jcall(obj, "addContribution", void, (FieldSpacecraftState, FieldTimeDerivativesEquations), arg0, arg1)
end

function add_contribution(obj::ForceModel, arg0::SpacecraftState, arg1::TimeDerivativesEquations)
    return jcall(obj, "addContribution", void, (SpacecraftState, TimeDerivativesEquations), arg0, arg1)
end

function depends_on_position_only(obj::HolmesFeatherstoneAttractionModel)
    return jcall(obj, "dependsOnPositionOnly", jboolean, ())
end

function equals(obj::Object, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function get_class(obj::Object)
    return jcall(obj, "getClass", Class, ())
end

function get_events_detectors(obj::HolmesFeatherstoneAttractionModel)
    return jcall(obj, "getEventsDetectors", Stream, ())
end

function get_field_events_detectors(obj::HolmesFeatherstoneAttractionModel, arg0::Field)
    return jcall(obj, "getFieldEventsDetectors", Stream, (Field,), arg0)
end

function get_mu(obj::HolmesFeatherstoneAttractionModel)
    return jcall(obj, "getMu", jdouble, ())
end

function get_parameter_driver(obj::AbstractForceModel, arg0::JString)
    return jcall(obj, "getParameterDriver", ParameterDriver, (JString,), arg0)
end

function get_parameters(obj::ForceModel)
    return jcall(obj, "getParameters", Vector{jdouble}, ())
end

function get_parameters(obj::ForceModel, arg0::Field)
    return jcall(obj, "getParameters", Vector{RealFieldElement}, (Field,), arg0)
end

function get_parameters_drivers(obj::HolmesFeatherstoneAttractionModel)
    return jcall(obj, "getParametersDrivers", Vector{ParameterDriver}, ())
end

function get_tide_system(obj::HolmesFeatherstoneAttractionModel)
    return jcall(obj, "getTideSystem", TideSystem, ())
end

function gradient(obj::HolmesFeatherstoneAttractionModel, arg0::AbsoluteDate, arg1::Vector3D, arg2::jdouble)
    return jcall(obj, "gradient", Vector{jdouble}, (AbsoluteDate, Vector3D, jdouble), arg0, arg1, arg2)
end

function gradient(obj::HolmesFeatherstoneAttractionModel, arg0::FieldAbsoluteDate, arg1::FieldVector3D, arg2::RealFieldElement)
    return jcall(obj, "gradient", Vector{RealFieldElement}, (FieldAbsoluteDate, FieldVector3D, RealFieldElement), arg0, arg1, arg2)
end

function hash_code(obj::Object)
    return jcall(obj, "hashCode", jint, ())
end

function init(obj::ForceModel, arg0::SpacecraftState, arg1::AbsoluteDate)
    return jcall(obj, "init", void, (SpacecraftState, AbsoluteDate), arg0, arg1)
end

function is_supported(obj::AbstractForceModel, arg0::JString)
    return jcall(obj, "isSupported", jboolean, (JString,), arg0)
end

function non_central_part(obj::HolmesFeatherstoneAttractionModel, arg0::AbsoluteDate, arg1::Vector3D, arg2::jdouble)
    return jcall(obj, "nonCentralPart", jdouble, (AbsoluteDate, Vector3D, jdouble), arg0, arg1, arg2)
end

function notify(obj::Object)
    return jcall(obj, "notify", void, ())
end

function notify_all(obj::Object)
    return jcall(obj, "notifyAll", void, ())
end

function to_string(obj::Object)
    return jcall(obj, "toString", JString, ())
end

function value(obj::HolmesFeatherstoneAttractionModel, arg0::AbsoluteDate, arg1::Vector3D, arg2::jdouble)
    return jcall(obj, "value", jdouble, (AbsoluteDate, Vector3D, jdouble), arg0, arg1, arg2)
end

function wait(obj::Object)
    return jcall(obj, "wait", void, ())
end

function wait(obj::Object, arg0::jlong)
    return jcall(obj, "wait", void, (jlong,), arg0)
end

function wait(obj::Object, arg0::jlong, arg1::jint)
    return jcall(obj, "wait", void, (jlong, jint), arg0, arg1)
end


text:
mutable struct IMMCTM
    K::Vector{Int}          # topics
    D::Int                  # documents
    N::Vector{Vector{Int}}  # observations per document modality
    M::Int                  # modalities
    I::Vector{Int}          # features per modality
    J::Vector{Vector{Int}}  # values per modality feature
    V::Vector{Int}          # vocab items per modality

    μ::Vector{Float64}
    Σ::Matrix{Float64}
    invΣ::Matrix{Float64}
    α::Vector{Vector{Float64}}

    ζ::Vector{Vector{Float64}}
    θ::Vector{Vector{Matrix{Float64}}}
    λ::Vector{Vector{Float64}}
    ν::Vector{Vector{Float64}}
    γ::Vector{Vector{Vector{Vector{Float64}}}}
    Elnϕ::Vector{Vector{Vector{Vector{Float64}}}}

    features::Vector{Matrix{Int}}
    X::Vector{Vector{Matrix{Int}}}

    converged::Bool
    elbo::Float64
    ll::Vector{Float64}

    function IMMCTM(k::Vector{Int}, α::Vector{Vector{Float64}},
                    features::Vector{Matrix{Int}},
                    X::Vector{Vector{Matrix{Int}}})
        model = new()

        model.K = copy(k)
        model.α = deepcopy(α)
        model.features = deepcopy(features)
        model.X = X

        model.D = length(X)
        model.M = length(features)
        model.I = [size(features[m])[2] for m in 1:model.M]
        model.J = [vec(maximum(features[m], dims=1)) for m in 1:model.M]
        model.V = [size(features[m])[1] for m in 1:model.M]
        model.N = [[sum(X[d][m][:, 2]) for m in 1:model.M] for d in 1:model.D]

        MK = sum(model.K)

        model.μ = zeros(MK)
        model.Σ = Matrix{Float64}(I, MK, MK)
        model.invΣ = Matrix{Float64}(I, MK, MK)

        model.θ = [
            [
                fill(1.0 / model.K[m], model.K[m], size(model.X[d][m])[1])
                for m in 1:model.M
            ] for d in 1:model.D
        ]

        model.γ = [
            [
                [
                    rand(1:100, model.J[m][i]) for i in 1:model.I[m]
                ] for k in 1:model.K[m]
            ] for m in 1:model.M
        ]
        model.Elnϕ = deepcopy(model.γ)
        update_Elnϕ!(model)

        model.λ = [zeros(MK) for d in 1:model.D]
        model.ν = [ones(MK) for d in 1:model.D]

        model.ζ = [Array{Float64}(undef, model.M) for d in 1:model.D]
        for d in 1:model.D update_ζ!(model, d) end

        model.converged = false

        return model
    end
end

function IMMCTM(k::Vector{Int}, α::Vector{Float64},
                features::Vector{Matrix{Int}},
                X::Vector{Vector{Matrix{Int}}})
    M = length(features)
    I = [size(features[m])[2] for m in 1:M]
    full_α = Vector{Float64}[fill(α[m], I[m]) for m in 1:M]
    return IMMCTM(k, full_α, features, X)
end

function calculate_sumθ(model::IMMCTM, d::Int)
    return vcat(
        [
            vec(sum(model.θ[d][m] .* model.X[d][m][:, 2]', dims=2))
            for m in 1:model.M
        ]...
    )
end

function calculate_Ndivζ(model::IMMCTM, d::Int)
    return vcat(
        [
            fill(model.N[d][m] / model.ζ[d][m], model.K[m]) for m in 1:model.M
        ]...
    )
end

function update_λ!(model::IMMCTM, d::Int)
    opt = Opt(:LD_MMA, sum(model.K))
    xtol_rel!(opt, 1e-4)
    xtol_abs!(opt, 1e-4)

    Ndivζ = calculate_Ndivζ(model, d)
    sumθ = calculate_sumθ(model, d)

    max_objective!(
        opt,
        (λ, ∇λ) -> λ_objective(
            λ, ∇λ, model.ν[d], Ndivζ, sumθ, model.μ, model.invΣ
        )
    )
    (optobj, optλ, ret) = optimize(opt, model.λ[d])
    model.λ[d] .= optλ
end

function update_ν!(model::IMMCTM, d::Int)
    opt = Opt(:LD_MMA, sum(model.K))
    lower_bounds!(opt, 1e-7)
    xtol_rel!(opt, 1e-4)
    xtol_abs!(opt, 1e-4)

    Ndivζ = calculate_Ndivζ(model, d)

    max_objective!(
        opt,
        (ν, ∇ν) -> ν_objective(ν, ∇ν, model.λ[d], Ndivζ, model.μ, model.invΣ)
    )
    (optobj, optν, ret) = optimize(opt, model.ν[d])
    model.ν[d] .= optν
end

function update_ζ!(model::IMMCTM, d::Int)
    start = 1
    for m in 1:model.M
        stop = start + model.K[m] - 1
        model.ζ[d][m] = sum(
            exp.(model.λ[d][start:stop] .+ 0.5 * model.ν[d][start:stop])
        )
        start += model.K[m]
    end
end

function update_θ!(model::IMMCTM, d::Int)
    offset = 0
    for m in 1:model.M
        for w in 1:size(model.X[d][m])[1]
            v = model.X[d][m][w, 1]

            for k in 1:model.K[m]
                model.θ[d][m][k, w] = exp(model.λ[d][offset + k])

                for i in 1:model.I[m]
                    model.θ[d][m][k, w] *= exp(
                        model.Elnϕ[m][k][i][model.features[m][v, i]]
                    )
                end
            end

        end
        model.θ[d][m] ./= sum(model.θ[d][m], dims=1)
        offset += model.K[m]
    end
end

function update_μ!(model::IMMCTM)
    model.μ .= mean(model.λ)
end

function update_Σ!(model::IMMCTM)
    model.Σ .= sum(diagm.(0 .=> model.ν))
    for d in 1:model.D
        diff = model.λ[d] .- model.μ
        model.Σ .+= diff * diff'
    end
    model.Σ ./= model.D
    model.invΣ .= inv(model.Σ)
end

function update_Elnϕ!(model::IMMCTM)
    for m in 1:model.M
        for k in 1:model.K[m]
            for i in 1:model.I[m]
                model.Elnϕ[m][k][i] .= digamma.(model.γ[m][k][i]) .-
                    digamma(sum(model.γ[m][k][i]))
            end
        end
    end
end

function update_γ!(model::IMMCTM)
    for m in 1:model.M
        for k in 1:model.K[m]
            for i in 1:model.I[m]
                for j in 1:model.J[m][i]
                    model.γ[m][k][i][j] = model.α[m][i]
                end
            end
        end
    end
    for d in 1:model.D
        for m in 1:model.M
            Nθ = model.θ[d][m] .* model.X[d][m][:, 2]'
            for w in 1:size(model.X[d][m])[1]
                v = model.X[d][m][w, 1]
                for k in 1:model.K[m]
                    for i in 1:model.I[m]
                        model.γ[m][k][i][model.features[m][v, i]] += Nθ[k, w]
                    end
                end
            end
        end
    end
    update_Elnϕ!(model)
end

function update_α!(model::IMMCTM)
    opt = Opt(:LD_MMA, 1)
    lower_bounds!(opt, 1e-7)
    xtol_rel!(opt, 1e-5)
    xtol_abs!(opt, 1e-5)

    for m in 1:model.M
        for i in 1:model.I[m]
            sum_Elnϕ = sum(sum(model.Elnϕ[m][k][i] for k in 1:model.K[m]))

            max_objective!(
                opt,
                (α, ∇α) -> α_objective(α, ∇α, sum_Elnϕ, model.K[m], model.J[m][i])
            )

            (optobj, optα, ret) = optimize(opt, model.α[m][i:i])
            model.α[m][i] = optα[1]
        end
    end
end

function calculate_ElnPϕ(model::IMMCTM)
    lnp = 0.0

    for m in 1:model.M
        for k in 1:model.K[m]
            for i in 1:model.I[m]
                lnp -= logmvbeta(fill(model.α[m][i], model.J[m][i]))
                for j in 1:model.J[m][i]
                    lnp += (model.α[m][i] - 1) * model.Elnϕ[m][k][i][j]
                end
            end
        end
    end

    return lnp
end

function calculate_ElnPη(model::IMMCTM)
    lnp = 0.0

    for d in 1:model.D
        diff = model.λ[d] .- model.μ
        lnp += 0.5 * (
            logdet(model.invΣ) -
            sum(model.K) * log(2π) -
            tr(diagm(0 => model.ν[d]) * model.invΣ) -
            (diff' * model.invΣ * diff)[1]
        )
    end

    return lnp
end

function calculate_ElnPZ(model::IMMCTM)
    lnp = 0.0

    for d in 1:model.D
        Eeη = exp.(model.λ[d] .+ 0.5model.ν[d])
        sumθ = calculate_sumθ(model, d)
        Ndivζ = calculate_Ndivζ(model, d)

        lnp += sum(model.λ[d] .* sumθ)
        lnp -= sum(Ndivζ .* Eeη) - sum(model.N[d])
        lnp -= sum(model.N[d] .* log.(model.ζ[d]))
    end

    return lnp
end

function calculate_ElnPX(model::IMMCTM)
    lnp = 0.0

    for d in 1:model.D
        for m in 1:model.M
            for w in 1:size(model.X[d][m])[1]
                v = model.X[d][m][w, 1]
                for i in 1:model.I[m]
                    for k in 1:model.K[m]
                        lnp += model.X[d][m][w, 2] * model.θ[d][m][k, w] *
                            model.Elnϕ[m][k][i][model.features[m][v, i]]
                    end
                end
            end
        end
    end

    return lnp
end

function calculate_ElnQϕ(model::IMMCTM)
    lnq = 0.0

    for m in 1:model.M
        for k in 1:model.K[m]
            for i in 1:model.I[m]
                lnq += -logmvbeta(model.γ[m][k][i])
                for j in 1:model.J[m][i]
                    lnq += (model.γ[m][k][i][j] - 1) * model.Elnϕ[m][k][i][j]
                end
            end
        end
    end
    return lnq
end

function calculate_ElnQη(model::IMMCTM)
    lnq = 0.0
    for d in 1:model.D
        lnq += -0.5 * (sum(log.(model.ν[d])) + sum(model.K) * (log(2π) + 1))
    end
    return lnq
end

function calculate_ElnQZ(model::IMMCTM)
    lnq = 0.0
    for d in 1:model.D
        for m in 1:model.M
            lnq += sum(model.X[d][m][:, 2]' .* log.(model.θ[d][m] .^ model.θ[d][m]))
        end
    end
    return lnq
end

function calculate_elbo(model::IMMCTM)
    elbo = 0.0
    elbo += calculate_ElnPϕ(model)
    elbo += calculate_ElnPη(model)
    elbo += calculate_ElnPZ(model)
    elbo += calculate_ElnPX(model)
    elbo -= calculate_ElnQϕ(model)
    elbo -= calculate_ElnQη(model)
    elbo -= calculate_ElnQZ(model)
    return elbo
end

function calculate_docmodality_loglikelihood(X::Matrix{Int},
        η::Vector{Float64}, ϕ::Vector{Vector{Vector{Float64}}},
        features::Matrix{Int})
    props = exp.(η) ./ sum(exp.(η))

    K = length(η)
    I = size(features)[2]

    ll = 0.0
    for w in 1:size(X, 1)
        v = X[w, 1]
        pw = 0.0
        for k in 1:K
            tmp = props[k]
            for i in 1:I
                tmp *= ϕ[k][i][features[v, i]]
            end
            pw += tmp
        end
        ll += X[w, 2] * log(pw)
    end

    return ll / sum(X[:, 2])
end

function calculate_modality_loglikelihood(X::Vector{Matrix{Int}},
        η::Vector{Vector{Float64}}, ϕ::Vector{Vector{Vector{Float64}}},
        features::Matrix{Int})
    D = length(X)

    ll = 0.0
    N = 0
    for d in 1:D
        doc_N = sum(X[d][:, 2])
        if doc_N > 0
            doc_ll = calculate_docmodality_loglikelihood(
                X[d], η[d], ϕ, features
            )
            ll += doc_ll * doc_N
            N += doc_N
        end
    end

    return ll / N
end

function calculate_loglikelihoods(X::Vector{Vector{Matrix{Int}}},
        model::IMMCTM)
    ll = Array{Float64}(undef, model.M)

    offset = 1
    for m in 1:model.M
        mk = offset:(offset + model.K[m] - 1)
        η = [model.λ[d][mk] for d in 1:model.D]
        Xm = [X[d][m] for d in 1:model.D]
        ϕ = [
            [model.γ[m][k][i] ./ sum(model.γ[m][k][i]) for i in 1:model.I[m]]
            for k in 1:model.K[m]
        ]

        ll[m] = calculate_modality_loglikelihood(Xm, η, ϕ, model.features[m])

        offset += model.K[m]
    end

    return ll
end

function fitdoc!(model::IMMCTM, d::Int)
    update_ζ!(model, d)
    update_θ!(model, d)
    update_ν!(model, d)
    update_λ!(model, d)
end

function fit!(model::IMMCTM; maxiter=100, tol=1e-4, verbose=true, autoα=false)
    ll = Vector{Float64}[]

    for iter in 1:maxiter
        for d in 1:model.D
            fitdoc!(model, d)
        end

        update_μ!(model)
        update_Σ!(model)
        update_γ!(model)
        if autoα
            update_α!(model)
        end

        push!(ll, calculate_loglikelihoods(model.X, model))
        if verbose
            println("$iter\tLog-likelihoods: ", join(ll[end], ", "))
        end

        if length(ll) > 10 && check_convergence(ll, tol=tol)
            model.converged = true
            break
        end
    end
    model.elbo = calculate_elbo(model)
    model.ll = ll[end]

    return ll
end

function fit_heldout(Xheldout::Vector{Vector{Matrix{Int}}}, model::IMMCTM;
        maxiter=100, verbose=false)

    heldout_model = IMMCTM(model.K, model.α, model.features, Xheldout)
    heldout_model.μ .= model.μ
    heldout_model.Σ .= model.Σ
    heldout_model.invΣ .= model.invΣ
    heldout_model.γ = deepcopy(model.γ)
    heldout_model.Elnϕ = deepcopy(model.Elnϕ)

    ll = Vector{Float64}[]
    for iter in 1:maxiter
        for d in 1:heldout_model.D
            fitdoc!(heldout_model, d)
        end

        push!(ll, calculate_loglikelihoods(Xheldout, heldout_model))

        if verbose
            println("$iter\tLog-likelihoods: ", join(ll[end], ", "))
        end

        if length(ll) > 10 && check_convergence(ll)
            heldout_model.converged = true
            break
        end
    end

    return heldout_model
end

function predict_modality_η(Xobs::Vector{Vector{Matrix{Int}}}, m::Int,
        model::IMMCTM; maxiter=100, verbose=false)
    obsM = setdiff(1:model.M, m)

    moffset = sum(model.K[1:(m - 1)])
    unobsMK = (moffset + 1):(moffset + model.K[m])
    obsMK = setdiff(1:sum(model.K), unobsMK)

    obsmodel = IMMCTM(model.K[obsM], model.α[obsM], model.features[obsM], Xobs)
    obsmodel.μ .= model.μ[obsMK]
    obsmodel.Σ .= model.Σ[obsMK, obsMK]
    obsmodel.invΣ .= model.invΣ[obsMK, obsMK]
    obsmodel.γ = deepcopy(model.γ[obsM])
    obsmodel.Elnϕ = deepcopy(model.Elnϕ[obsM])

    ll = Vector{Float64}[]
    for iter in 1:maxiter
        for d in 1:obsmodel.D
            fitdoc!(obsmodel, d)
        end

        push!(ll, calculate_loglikelihoods(Xobs, obsmodel))

        if verbose
            println("$iter\tLog-likelihoods: ", join(ll[end], ", "))
        end

        if length(ll) > 10 && check_convergence(ll)
            obsmodel.converged = true
            break
        end
    end

    if !obsmodel.converged
        warn("model not converged.")
    end

    η = [
        (
            model.μ[unobsMK] .+ model.Σ[unobsMK, obsMK] *
            model.invΣ[obsMK, obsMK] * (obsmodel.λ[d] .- model.μ[obsMK])
        )
        for d in 1:obsmodel.D
    ]

    return η
end


text:
#=
Julia EpanetToolkit
=#
module epamodule
using Libdl

lib_path = "/home/augusto/Documents/IC-2020/optimized-calibration/libs-epanet/EPANET/build/lib/libepanet2.so" # absolute path to lib epanet binary file
lib = Libdl.dlopen(lib_path)

#= ----------------------------------------------------------

    Project functions

---------------------------------------------------------- =#

function ENepanet(inpFile::String, rptFile::String="", outFile::String="", vfunc=NaN)

    sym = Libdl.dlsym(lib, :ENepanet)
    err = ccall(sym, Cint, (Cstring, Cstring, Cstring), )
end

# testar
function ENinit(rptFile::String, outFile::String, unitsType::Inf64, headlossType::Int64)
    sym = Libdl.dlsym(lib, :ENinit)
    return ccall(sym, Cint, (Cstring, Cstring, Cint, Cint), rptFile, outFile, unitsType, headlossType)

end

# funciona
function ENopen(inpFile::String, rptFile::String = "", outFile::String = "")
    sym = Libdl.dlsym(lib, :ENopen)
    err = ccall(sym, Cint, (Cstring,Cstring,Cstring), inpFile, rptFile, outFile)
    if err != 0
        return "Error: "*string(err)
    end
end

function ENgettitle()

end

function ENsettitle()
end

function ENgetcomment()
end

function ENsetcomment()
end

function ENgetcount()
end

function ENsaveinpfile()
end

function ENclose()
    sym = Libdl.dlsym(lib, :ENclose)
    err = ccall(sym,Cint,())
    if err != 0
        return "Error: "*string(err)
    end
end

#= ----------------------------------------------------------

    Hydraulic Analysis Functions

---------------------------------------------------------- =#

function ENsolveH()
    sym = Libdl.dlsym(lib, :ENsolveH)
    err = ccall(sym,Cint,())
    if err != 0
        return "Error: "*string(err)
    end
end

function ENsaveH()
end

function ENopenH()
    sym = Libdl.dlsym(lib, :ENopenH)
    err = ccall(sym,Cint,())
    if err != 0
        return "Error: "*string(err)
    end
end

function ENinitH()
end

function ENrunH()
end

function ENnextH()
end

function ENcloseH()
    sym = Libdl.dlsym(lib, :ENcloseH)
    err = ccall(sym,Cint,())
    if err != 0
        return "Error: "*string(err)
    end
end

function ENsavehydfile()
end

function ENusehydfile()
end

#= ----------------------------------------------------------

    Water Quality Analysis Functions

---------------------------------------------------------- =#

#= ----------------------------------------------------------

    Reporting Functions

---------------------------------------------------------- =#

#= ----------------------------------------------------------

    Analysis Options Functions

---------------------------------------------------------- =#

#= ----------------------------------------------------------

    Node Functions

---------------------------------------------------------- =#

#= ----------------------------------------------------------

    Nodal Demand Functions

---------------------------------------------------------- =#

#= ----------------------------------------------------------

    Link Functions

---------------------------------------------------------- =#

#= ----------------------------------------------------------

    Pump Functions

---------------------------------------------------------- =#

#= ----------------------------------------------------------

    Time Pattern Functions

---------------------------------------------------------- =#

#= ----------------------------------------------------------

    Data Curve Functions

---------------------------------------------------------- =#

#= ----------------------------------------------------------

    Simple Controls Functions

---------------------------------------------------------- =#

#= ----------------------------------------------------------

    Rule-Based Controls Functions

---------------------------------------------------------- =#

#= ----------------------------------------------------------

    Analysis Options Functions

---------------------------------------------------------- =#

function ENgetnodeindex(id::String)
    sym = Libdl.dlsym(lib, :ENgetnodeindex)
    index = Ref{Int64}(0)
    err = ccall(sym, Cint, (Cstring, Ref{Int64}), id, index)
    if err != 0
        return "Error: "*string(err)
    end
    return index.x
end

function ENgetnodevalue(index::Int64, paramcode::Int64)
    sym = Libdl.dlsym(lib,:ENgetnodevalue)
    ref = Ref{Float32}(0)
    err = ccall(sym, Cint,(Cint,Cint,Ref{Float32}), index, paramcode, ref)
    if err != 0
        return "Error: "*string(err)
    end
    return ref.x
end

function ENsetnodevalue(index::Int64, paramcode::Int64, value::Float64)
    sym = Libdl.dlsym(lib,:ENsetnodevalue)
    err = ccall(sym, Cint, (Cint, Cint, Cfloat), index, paramcode, value)
    if err != 0
        return "Error: "*string(err)
    end
end

function ENgetlinkindex(id::String)
    sym = Libdl.dlsym(lib, :ENgetlinkindex)
    index = Ref{Int64}(0)
    err = ccall(sym, Cint, (Cstring, Ref{Int64}), id, index)
    if err != 0
        return "Error: "*string(err)
    end
    return index.x
end

function ENgetlinkvalue(index::Int64, paramcode::Int64)
    sym = Libdl.dlsym(lib,:ENgetlinkvalue)
    ref = Ref{Float32}(0)
    err = ccall(sym, Cint,(Cint,Cint,Ref{Float32}), index, paramcode, ref)
    if err != 0
        return "Error: "*string(err)
    end
    return ref.x
end

function ENsetlinkvalue(index::Int64, paramcode::Int64, value::Float64)
    sym = Libdl.dlsym(lib,:ENsetlinkvalue)
    err = ccall(sym, Cint, (Cint, Cint, Cfloat), index, paramcode, value)
    if err != 0
        return "Error: "*string(err)
    end
end

#-------------------------------------------------------------#
#= Node paramcodes =#
EN_ELEVATION     = 0
EN_BASEDEMAND    = 1
EN_PATTERN       = 2
EN_EMITTER       = 3
EN_INITQUAL      = 4
EN_SOURCEQUAL    = 5
EN_SOURCEPAT     = 6
EN_SOURCETYPE    = 7
EN_TANKLEVEL     = 8
EN_DEMAND        = 9
EN_HEAD          = 10
EN_PRESSURE      = 11
EN_QUALITY       = 12
EN_SOURCEMASS    = 13
EN_INITVOLUME    = 14
EN_MIXMODEL      = 15
EN_MIXZONEVOL    = 16

EN_TANKDIAM      = 17
EN_MINVOLUME     = 18
EN_VOLCURVE      = 19
EN_MINLEVEL      = 20
EN_MAXLEVEL      = 21
EN_MIXFRACTION   = 22
EN_TANK_KBULK    = 23

#= Link paramcodes =#
EN_DIAMETER      = 0
EN_LENGTH        = 1
EN_ROUGHNESS     = 2
EN_MINORLOSS     = 3
EN_INITSTATUS    = 4
EN_INITSETTING   = 5
EN_KBULK         = 6
EN_KWALL         = 7
EN_FLOW          = 8
EN_VELOCITY      = 9
EN_HEADLOSS      = 10
EN_STATUS        = 11
EN_SETTING       = 12
EN_ENERGY        = 13

#= Time paramcodes =#
EN_DURATION      = 0      # /* Time parameters */
EN_HYDSTEP       = 1
EN_QUALSTEP      = 2
EN_PATTERNSTEP   = 3
EN_PATTERNSTART  = 4
EN_REPORTSTEP    = 5
EN_REPORTSTART   = 6
EN_RULESTEP      = 7
EN_STATISTIC     = 8
EN_PERIODS       = 9

#= Component counts =#
EN_NODECOUNT     = 0
EN_TANKCOUNT     = 1
EN_LINKCOUNT     = 2
EN_PATCOUNT      = 3
EN_CURVECOUNT    = 4
EN_CONTROLCOUNT  = 5

#= Node types =#
EN_JUNCTION      = 0
EN_RESERVOIR     = 1
EN_TANK          = 2

#= Link types =#
EN_CVPIPE        = 0
EN_PIPE          = 1
EN_PUMP          = 2
EN_PRV           = 3
EN_PSV           = 4
EN_PBV           = 5
EN_FCV           = 6
EN_TCV           = 7
EN_GPV           = 8

#= Quality analysis types =#
EN_NONE          = 0
EN_CHEM          = 1
EN_AGE           = 2
EN_TRACE         = 3

#= Source quality types =#
EN_CONCEN        = 0
EN_MASS          = 1
EN_SETPOINT      = 2
EN_FLOWPACED     = 3

# Flow units types
EN_CFS           = 0
EN_GPM           = 1
EN_MGD           = 2
EN_IMGD          = 3
EN_AFD           = 4
EN_LPS           = 5
EN_LPM           = 6
EN_MLD           = 7
EN_CMH           = 8
EN_CMD           = 9

#= Misc. options =#
EN_TRIALS        = 0
EN_ACCURACY      = 1
EN_TOLERANCE     = 2
EN_EMITEXPON     = 3
EN_DEMANDMULT    = 4

#= Control types =#
EN_LOWLEVEL      = 0
EN_HILEVEL       = 1
EN_TIMER         = 2
EN_TIMEOFDAY     = 3

#= Time statistic types =#
EN_AVERAGE       = 1
EN_MINIMUM       = 2
EN_MAXIMUM       = 3
EN_RANGE         = 4

#= Tank mixing models =#
EN_MIX1          = 0
EN_MIX2          = 1
EN_FIFO          = 2
EN_LIFO          = 3

EN_NOSAVE        = 0      # /* Save-results-to-file flag */
EN_SAVE          = 1
EN_INITFLOW      = 10     # /* Re-initialize flow flag   */

end # Ends module

max_stars_repo_name: augusnunes/epanet-julia

max_issues_repo_name: augusnunes/epanet-julia

max_forks_repo_name: augusnunes/epanet-julia


text:
##### directional funcs

export orthog, orthogs

orthogs(::Val{1}) = (2,3)
orthogs(::Val{2}) = (3,1)
orthogs(::Val{3}) = (1,2)
orthogs(i::Int) = orthogs(Val(i))

orthog(::Val{3}, ::Val{2}) = 1
orthog(::Val{2}, ::Val{3}) = 1
orthog(::Val{3}, ::Val{1}) = 2
orthog(::Val{1}, ::Val{3}) = 2
orthog(::Val{1}, ::Val{2}) = 3
orthog(::Val{2}, ::Val{1}) = 3
orthog(i::Int, j::Int) = orthog(Val(i), Val(j))


text:
# Activating Julia environment
using Pkg
Pkg.activate(".")

# Packages
using JuMP, Gurobi, Random, Distributions, MathOptInterface, Plots
MOI = MathOptInterface

# Generating data
n = 20;
b = 4;
a = rand(20);
c = rand(20);

rhos = collect(0:0.1:2)
costs = []
for ρ in rhos
    m = JuMP.Model(Gurobi.Optimizer)
    @variable(m, x[1:n], Bin)
    # Ellipsoidal uncertainty on the cost of the form...
    # @uncertain(m, u[1:n])
    # @constraint(m, norm(u, 2) <= ρ)
    if ρ != 0
        @constraint(m, [(b - sum(a .* x))/ρ, x...] in MOI.SecondOrderCone(n + 1))
    else
        @constraint(m, sum(a .* x) ≤ b)
    end
    @objective(m, Max, sum(c .* x))
    optimize!(m)
    push!(costs, JuMP.objective_value(m))
end

plot(rhos, costs, title = "Robustness/optimality tradeoff", xlabel="Safety factor", ylabel="Optimal cost", legend=false)


text:
# Autogenerated wrapper script for openPMD_api_jll for x86_64-w64-mingw32-cxx03-julia_version+1.6.0
export libopenPMD, libopenPMD_jl

using ADIOS2_jll
using CompilerSupportLibraries_jll
using MicrosoftMPI_jll
using libcxxwrap_julia_jll
JLLWrappers.@generate_wrapper_header("openPMD_api")
JLLWrappers.@declare_library_product(libopenPMD, "libopenPMD.dll")
JLLWrappers.@declare_library_product(libopenPMD_jl, "libopenPMD.jl.dll")
function __init__()
    JLLWrappers.@generate_init_header(ADIOS2_jll, CompilerSupportLibraries_jll, MicrosoftMPI_jll, libcxxwrap_julia_jll)
    JLLWrappers.@init_library_product(
        libopenPMD,
        "bin\\libopenPMD.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libopenPMD_jl,
        "bin\\libopenPMD.jl.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using Documenter, Checkers

makedocs(modules  = [Checkers],
         doctest  = false,
         clean    = true,
         #format   = :html,
         sitename = "Checkers.jl",
         authors  = "Philip Kalikman, Efim Abrikosov",
         html_prettyurls = !("local" in ARGS),
         #pages    = Any[
         #               "Home" => "index.md",
         #               "Macros" => Any["macros/test-forall.md",
         #                               "macros/test-formany.md",
         #                               "macros/test-exists.md"]
         #           ]
)

deploydocs(deps   = Deps.pip("mkdocs", "mkdocs-material" ,"python-markdown-math", "pygments"),
           repo   = "github.com/pkalikman/Checkers.jl.git",
           target = "site",
           julia  = "0.6",
           latest = "master",
           osname = "linux")


text:
# dependencies
using LFAToolkit
using LinearAlgebra
using CSV
using DataFrames

# setup
dimension = 1
numbercomponents = 1
mesh = []
if dimension == 1
    mesh = Mesh1D(1.0)
elseif dimension == 2
    mesh = Mesh2D(1.0, 1.0)
elseif dimension == 3
    mesh = Mesh3D(1.0, 1.0, 1.0)
end
smoothingfactors = DataFrame()

# test range
for P = 1:4
    println("p = ", 2^P)
    # setup
    # -- bases
    p = 2^P

    # -- diffusion operator
    diffusion = GalleryOperator("diffusion", p + 1, p + 1, mesh)

    # -- Jacobi smoother
    jacobi = Jacobi(diffusion)

    # compute smoothing factor
    # -- setup
    numbersteps = 100
    maxeigenvalue = 0
    θ_min = -π / 2
    θ_min_high = π / 2
    θ_max = 3π / 2
    θ_step = 2π / numbersteps
    θ_range = θ_min:θ_step:(θ_max-θ_step)
    ω_range = 0.50:0.01:1.05
    num_ω = max(size(ω_range)...)

    # -- compute
    maxeigenvalue = zeros(num_ω)
    θ_maxeigenvalue = -1 * ones(num_ω, dimension)
    # -- 1D --
    if dimension == 1
        for i = 1:numbersteps
            θ = [θ_range[i]]
            if abs(θ[1]) > π / 128 && θ[1] > θ_min_high
                A = computesymbols(jacobi, [1.0], θ)
                for w = 1:num_ω
                    eigenvalues = [abs(val) for val in eigvals(I - ω_range[w] * (I - A))]
                    currentmaxeigenvalue = max(eigenvalues...)
                    if (currentmaxeigenvalue > maxeigenvalue[w])
                        maxeigenvalue[w] = currentmaxeigenvalue
                        θ_maxeigenvalue[w, :] = θ / π
                    end
                end
            end
        end
        # -- 2D --
    elseif dimension == 2
        for i = 1:numbersteps, j = 1:numbersteps
            θ = [θ_range[i], θ_range[j]]
            if sqrt(abs(θ[1])^2 + abs(θ[2])^2) > π / 128 &&
               (θ[1] > θ_min_high || θ[2] > θ_min_high)
                A = computesymbols(jacobi, [1.0], θ)
                for w = 1:num_ω
                    eigenvalues = [abs(val) for val in eigvals(I - ω_range[w] * (I - A))]
                    currentmaxeigenvalue = max(eigenvalues...)
                    if (currentmaxeigenvalue > maxeigenvalue[w])
                        maxeigenvalue[w] = currentmaxeigenvalue
                        θ_maxeigenvalue[w, :] = θ / π
                    end
                end
            end
        end
        # -- 3D --
    elseif dimension == 3
        for i = 1:numbersteps, j = 1:numbersteps, k = 1:numbersteps
            θ = [θ_range[i], θ_range[j], θ_range[k]]
            if sqrt(abs(θ[1])^2 + abs(θ[2])^2 + abs(θ[3])^2) > π / 128 &&
               (θ[1] > θ_min_high || θ[2] > θ_min_high || θ[3] > θ_min_high)
                A = computesymbols(jacobi, [1.0], θ)
                for w = 1:num_ω
                    eigenvalues = [abs(val) for val in eigvals(I - ω_range[w] * (I - A))]
                    currentmaxeigenvalue = max(eigenvalues...)
                    if (currentmaxeigenvalue > maxeigenvalue[w])
                        maxeigenvalue[w] = currentmaxeigenvalue
                        θ_maxeigenvalue[w, :] = θ / π
                    end
                end
            end
        end
    end
    (_, ω_index) = findmin(maxeigenvalue)
    append!(
        smoothingfactors,
        DataFrame(
            p = p,
            ω = ω_range[ω_index],
            θ = θ_maxeigenvalue[ω_index, :],
            ρ = maxeigenvalue[ω_index],
        ),
    )
end

CSV.write("jacobi-smoothing-factor.csv", smoothingfactors)


text:
function validar(doc::Doc2011)
	validar_documento(doc)
	for conta in doc.contas
		validar_conta(conta)
	end
	return true
end

function validar_documento(doc::Doc2011)
	assert_contas_repetidas(doc)
	assert_formulas(doc)
end

function validar_conta(conta::Conta)
	assert_codigo_conta(conta)
	assert_soma_detalhes(conta)
	assert_sinal(conta)
	assert_tem_detalhe(conta)
	for detalhe in conta.detalhes
		validar_detalhe(detalhe)
	end
end

function validar_detalhe(detalhe::DetalheConta)
	assert_sinal(detalhe)
end

function assert_contas_repetidas(doc::Doc2011)
	codigos = [conta.codigo for conta in doc.contas]
	@assert length(doc.contas) == length(unique(codigos)) "contas com codigos repetidos"
end

function assert_soma_detalhes(conta::Conta)
	if length(conta.detalhes) > 0
		soma = sum([elem.valor for elem in conta.detalhes])
		@assert conta.valor ≈ soma "conta $(conta.codigo): soma de detalhes ($soma) nao bate com valor ($(conta.valor))"
	end
end

CONTAS_VALIDAS = [
	"111000",
	"121000",
	"121100",
	"131000",
	"132000",
	"141000",
	"151000",
	"161000",
	"171000",
	"181000",
	"210000",
	"220000",
	"230000",
	"240000",
	"310000",
	"310100",
	"310101",
	"310102",
	"310103",
	"310104",
	"310105",
	"410100",
	"410101",
	"410200",
	"410201",
	"410202",
	"410300",
	"410301",
	"410302",
	"410400",
	"410401",
	"410402",
	"410500",
	"410501",
	"410502",
	"410503",
	"410504",
	"410600",
	"410601",
	"410602",
	"410603",
	"410604",
	"410700",
	"410701",
	"410702",
	"410703",
	"410704",
	"410800",
	"410801",
	"410802",
	"410900",
	"410901",
	"410904",
	"410907",
	"410908",
	"410909",
	"410910",
	"501000",
	"502000",
	"503000",
	"504000",
	"505000",
	"506000",
	"507000",
	"508000",
	"610000",
	"620000",
	"620100",
	"620200",
	"620300",
	"620400",
	"620500",
	"620501",
	"620502",
	"620503",
	"620504",
	"620505",
	"630000",
	"630100",
	"630200",
	"630300",
	"630400",
	"630500",
	"630501",
	"630502",
	"630503",
	"630504",
	"630505",
	"710000"
]

function assert_codigo_conta(conta::Conta)
	@assert conta.codigo in CONTAS_VALIDAS "conta invalida: $(conta.codigo)"
end

function assert_sinal(x::Union{Conta, DetalheConta})
	@assert x.valor >= -eps()
end

ELEMENTOS_DETALHES_CONTAS = Dict{String, Vector{DataType}}([
	"111000" => [Moeda; Posicao]
	"121000" => [Moeda; Posicao]
	"121100" => [Moeda; Posicao]
	"131000" => [Moeda; Posicao]
	"132000" => [Moeda; Posicao]
	"141000" => [Moeda; Posicao]
	"151000" => [Moeda; Posicao]
	"161000" => [Moeda; Posicao]
	"171000" => [Moeda; Posicao]
	"181000" => [Moeda; Posicao]
	"210000" => [Pais; Moeda]
	"220000" => [Pais; Moeda]
	"230000" => [Pais; Moeda]
	"240000" => [Pais; Moeda]
	"501000" => [Moeda]
	"502000" => [Moeda]
])

function assert_tem_detalhe(conta::Conta)
	codigo = conta.codigo
	if length(conta.detalhes) == 0
		@assert !haskey(ELEMENTOS_DETALHES_CONTAS, codigo) "conta $codigo deve ter detalhes dos tipos $(ELEMENTOS_DETALHES_CONTAS[codigo])"
	else
		for detalhe in conta.detalhes
			assert_elementos_detalhe(codigo, detalhe)
		end
	end

end

function assert_elementos_detalhe(codigo::String, detalhe::DetalheConta)
	gabarito = ELEMENTOS_DETALHES_CONTAS[codigo]
	tipos = [typeof(el) for el in detalhe.elementos]
	@assert length(tipos) == length(gabarito) "conta $codigo: $detalhe deveria ter $(length(gabarito)) ElementoDetalhe mas tem $(length(tipos))"
	for t in tipos
		@assert t in gabarito "conta $codigo: elemento detalhe invalido ($t)"
	end
end

FORMULAS = Dict{String, Function}([
	"310000" => d -> get_valor(d, "310100") * (get_valor(d, "310105")/100)
	"310100" => d -> get_valor(d, "310101", 0.0) + get_valor(d, "310102", 0.0) + get_valor(d, "310103", 0.0) + get_valor(d, "310104", 0.0)
	"410100" => d -> 100.0
	"410200" => d -> get_valor(d, "410201", 0.0) + get_valor(d, "410202", 0.0)
	"410300" => d -> get_valor(d, "410301", 0.0) + get_valor(d, "410302", 0.0)
	"410400" => d -> get_valor(d, "410401", 0.0) + get_valor(d, "410402", 0.0)
	"410500" => d -> get_valor(d, "410501", 0.0) + get_valor(d, "410502", 0.0) + get_valor(d, "410503", 0.0) + get_valor(d, "410504", 0.0)
	"410600" => d -> get_valor(d, "410601", 0.0) + get_valor(d, "410602", 0.0) + get_valor(d, "410603", 0.0) + get_valor(d, "410604", 0.0)
	"410700" => d -> get_valor(d, "410701", 0.0) + get_valor(d, "410702", 0.0) + get_valor(d, "410703", 0.0) + get_valor(d, "410704", 0.0)
	"410800" => d -> get_valor(d, "410801", 0.0) + get_valor(d, "410802", 0.0)
	"410900" => d -> get_valor(d, "410901", 0.0) + get_valor(d, "410904", 0.0) + get_valor(d, "410907", 0.0) + get_valor(d, "410908", 0.0) + get_valor(d, "410909", 0.0) + get_valor(d, "410910", 0.0)
	"503000" => d -> get_valor(d, "310000", 0.0) + get_valor(d, "410400", 0.0) + get_valor(d, "410500", 0.0) + get_valor(d, "410600", 0.0) + get_valor(d, "410700", 0.0) + get_valor(d, "410800", 0.0) + get_valor(d, "410900", 0.0)
	"620000" => d -> get_valor(d, "620200", 0.0) + get_valor(d, "620300", 0.0) + get_valor(d, "620400", 0.0) + get_valor(d, "620500", 0.0) - get_valor(d, "620100", 0.0)
	"620500" => d -> get_valor(d, "620502", 0.0) + get_valor(d, "620503", 0.0) + get_valor(d, "620504", 0.0) + get_valor(d, "620505", 0.0) - get_valor(d, "620501", 0.0)
	"630000" => d -> get_valor(d, "630200", 0.0) + get_valor(d, "630300", 0.0) + get_valor(d, "630400", 0.0) + get_valor(d, "630500", 0.0) - get_valor(d, "630100", 0.0)
	"630500" => d -> get_valor(d, "630502", 0.0) + get_valor(d, "630503", 0.0) + get_valor(d, "630504", 0.0) + get_valor(d, "630505", 0.0) - get_valor(d, "630501", 0.0)
	"710000" => d -> max(get_valor(d, "610000", 0.0) + get_valor(d, "505000", 0.0),  (get_valor(d, "504000")/100) * get_valor(d, "503000", 0.0))
])

function assert_formulas(doc::Doc2011)

	for (codigo, fun) in FORMULAS
		if has_conta(doc, codigo)
			stored = get_valor(doc, codigo)
			evaluated = fun(doc)
			@assert stored ≈ evaluated "conta $codigo: valor informado ($stored) nao bate com o valor da formula ($evaluated)"
		end
	end

end