text:
import NMRDataSetup
import NMRSpectraSimulator
import NMRSpecifyRegions

include("../src/NMRCalibrate.jl")
import .NMRCalibrate

using LinearAlgebra
using FFTW
import PyPlot
import BSON
import JSON

import Statistics
import Random

# include("./helpers/final_helpers.jl")
# include("./helpers/solute_calibrate.jl")
# include("./helpers/loop_entries1.jl")

##### global constants.
#cs_config_path = "/home/roy/Documents/repo/NMRData/src/input/reduced_cs_config.txt"
SH_config_path = "/home/roy/Documents/repo/NMRData/input/SH_configs/select_compounds_SH_configs.json"
surrogate_config_path = "/home/roy/Documents/repo/NMRData/input/surrogate_configs/select_compounds_SH_configs.json"
fit_config_path = "/home/roy/Documents/repo/NMRData/input/fit_configs/calibrate_700MHz_type1_select_compounds.json"

# get mapping from molecule names to their spin system info json files.
H_params_path = "/home/roy/Documents/repo/NMRData/input/coupling_info"
dict_compound_to_filename = JSON.parsefile("/home/roy/Documents/repo/NMRData/input/compound_mapping/select_compounds.json")

# where the experiment bson file is located.
root_folder = "/home/roy/MEGAsync/outputs/NMR/experiments/NRC"
experiment_path = joinpath(root_folder, "NRC-glucose-2018") # second input is subfolder name.
load_path = joinpath(experiment_path, "experiment.bson")

# where the saved calibration project folder is for this experiment.
project_name = "NRC-glucose-2018"

save_base_folder = "/home/roy/MEGAsync/outputs/NMR/calibrate/NRC"
save_folder = joinpath(save_base_folder, project_name)
save_path = joinpath(save_folder, "proxy.bson")
isdir(save_path) || mkpath(save_folder)
#####

##### user input for loading data.

molecule_names = ["D-(+)-Glucose";]
w = [1.0;]

# molecule_names = ["D-(+)-Glucose"; "DSS"; "Ethanol"]
# w = [1.0; 1.0; 1.0]

# molecule_names = ["D-(+)-Glucose"; "DSS"; ]
# w = [1.0; 1.0;]

# note that 0.1% DSS is 0.0046 M = 4.6 mM.

##### end inputs for loading data.

##### user input for surrogate.
dummy_SSFID = NMRSpectraSimulator.SpinSysParamsType1(0.0) # level 2 model.

##### end inputs for surrogate.

println("Now on $(project_name)")

##### load data that was saved in store_exp.jl. See that script for details of the following commands.
load_path = joinpath(load_path)
dic = BSON.load(load_path)
s_t = dic[:s_t]
fs = dic[:fs]
SW = dic[:SW]
ν_0ppm = dic[:ν_0ppm]
α_0ppm = dic[:α_0ppm]
β_0ppm = dic[:β_0ppm]
λ_0ppm = dic[:λ_0ppm]

hz2ppmfunc = uu->(uu - ν_0ppm)*SW/fs
ppm2hzfunc = pp->(ν_0ppm + pp*fs/SW)

# start the first entry at the frequency corresponding to offset_Hz, which we set to 0.3 ppm.
offset_ppm = 0.3
offset_Hz = ν_0ppm - (ppm2hzfunc(offset_ppm)-ppm2hzfunc(0.0))
DFT_s = fft(s_t)
U_DFT, U_y, U_inds = NMRDataSetup.getwraparoundDFTfreqs(length(s_t), fs, offset_Hz)

S_U = DFT_s[U_inds]
P_y = hz2ppmfunc.(U_y)

val, ind = findmin( abs.(U_y .- ν_0ppm) )
Z = abs(S_U[ind])

y = S_U ./ Z

### end loading spectrum.
#@assert 1==2

####### compute surrogate.

# TODO document how to use the configuration files. then delete these default values from tutorial.
tol_coherence = 1e-2 # resonances are pairs of eigenvalues of the Hamiltonian that have quantum numbers that differ by -1. This is the tolerance away from -1 that is allowed.
α_relative_threshold = 0.05 # resonances with relative amplitude less than this factor compared to the maximum resonance in the spin group will be removed. Set to 0.0 to see every single resonance component.
Δc_partition_radius = 0.3 # determines how many resonances get grouped together. Larger number means less groups and thus more resonances per group.
λ0 = λ_0ppm

Δr_default = 1.0 # the samples used to build the surrogate is taken every `Δr` radian on the frequency axis. Decrease for improved accuracy at the expense of computation resources.
Δκ_λ_default = 0.05 # the samples used to build thes urrogate for κ_λ are taken at this sampling spacing. Decrease for improved accuracy at the expense of computation resources.
Δcs_max_scalar_default = 0.2 # In units of ppm. interpolation border that is added to the lowest and highest resonance frequency component of the mixture being simulated.
κ_λ_lb_default = 0.7 # interpolation lower limit for κ_λ.
κ_λ_ub_default = 1.5 # interpolation upper limit for κ_λ.

# get a surrogate where K_{n,i} is encouraged to be no larger than `early_exit_part_size`.
println("Timing: mag equivalence")
@time MEs = NMRSpectraSimulator.getmageqinfomixture(molecule_names,
    H_params_path,
    dict_compound_to_filename;
    unique_cs_atol = 1e-6)

println("Timing: setupmixtureproxies()")
@time mixture_params = NMRSpectraSimulator.setupmixtureSH(molecule_names,
    H_params_path, dict_compound_to_filename, fs, SW,
    ν_0ppm;
    config_path = SH_config_path,
    tol_coherence = tol_coherence,
    α_relative_threshold = α_relative_threshold,
    Δc_partition_radius = Δc_partition_radius)
As = mixture_params

ΩS_ppm = NMRSpectraSimulator.getPsnospininfo(As, hz2ppmfunc)
ΩS_ppm_sorted = sort(NMRSpectraSimulator.combinevectors(ΩS_ppm))

u_offset = 0.5 # in ppm.
u_min = ppm2hzfunc(ΩS_ppm_sorted[1] - u_offset)
u_max = ppm2hzfunc(ΩS_ppm_sorted[end] + u_offset)

println("Timing: fitproxies!()")
dummy_SSFID = NMRSpectraSimulator.SpinSysParamsType1(0.0)
Bs = NMRSpectraSimulator.fitproxies(As, dummy_SSFID, λ0;
    names = molecule_names,
    config_path = surrogate_config_path,
    Δcs_max_scalar_default = Δcs_max_scalar_default,
    κ_λ_lb_default = κ_λ_lb_default,
    κ_λ_ub_default = κ_λ_ub_default,
    u_min = u_min,
    u_max = u_max,
    Δr_default = Δr_default,
    Δκ_λ_default = Δκ_λ_default)

####### end mixture proxy.

### save.
BSON.bson(save_path,
As = As,
Bs = Bs,
surrogate_config_path = surrogate_config_path,
SH_config_path = SH_config_path,
u_min = u_min,
u_max = u_max,
y = y,
U_y = U_y,
SW = SW,
fs = fs,
ν_0ppm  = ν_0ppm,
λ0 = λ0,
molecule_names = molecule_names,
w = w)


text:
# Compile speed tests
# Tests time to get from nothing to certain tasks
using JuMP  # Ensure compiled before running the tests

# Test 1
# Size 100, aff_str into model_str
tic()
run(`$(ENV["_"]) -e """
using JuMP
m = Model()
N = 100
@variable(m, x[1:N])
@constraint(m, sum{i*x[i],i=1:N} >= N)
@time JuMP.aff_str(JuMP.REPLMode, m.linconstr[end].terms)
@time JuMP.model_str(JuMP.REPLMode, m)"""
`)
test1_time = toq()

# Test 2
# Size 100, model_str
tic()
run(`$(ENV["_"]) -e """
using JuMP
m = Model()
N = 100
@variable(m, x[1:N])
@constraint(m, sum{i*x[i],i=1:N} >= N)
@time JuMP.model_str(JuMP.REPLMode, m)"""
`)
test2_time = toq()

# Test 3
# Size 10,000, print model
tic()
run(`$(ENV["_"]) -e """
using JuMP
m = Model()
N = 10000
@variable(m, x[1:N])
@constraint(m, sum{i*x[i],i=1:N} >= N)
@time sprint(print, m)"""
`)
test3_time = toq()

# Test 4
# Size 10, 4 variables, print model
tic()
run(`$(ENV["_"]) -e """
using JuMP
m = Model()
N = 10
@variable(m, x1[1:N])
@variable(m, x2[1:N,f=1:N])
@variable(m, x3[1:N,f=1:2:N])
@variable(m, x4[[:a,:b,:c]])
@variable(m, x5[[:a,:b,:c],[:d,"e",4]])
@constraint(m,
    sum{i*x1[i],i=1:N} +
    sum{i*f*x2[i,f],i=1:N,f=1:N} +
    sum{i*f*x3[i,f],i=1:N,f=1:2:N} +
    sum(x4) >= N)
@time sprint(print, m)"""
`)
test4_time = toq()

println("Julia ", VERSION)
println("Current branch: $(chomp(readall(`git rev-parse --abbrev-ref HEAD`)))")
println("Current commit: $(chomp(readall(`git rev-parse HEAD`)))")
@printf("Test 1: %7.2f s\n", test1_time)
@printf("Test 2: %7.2f s\n", test2_time)
@printf("Test 3: %7.2f s\n", test3_time)
@printf("Test 4: %7.2f s\n", test4_time)

ext: jl


text:
##
cd(@__DIR__())
rot(ϕ) = Complex{Int}(round(Int,cos(deg2rad(ϕ))), round(Int,sin(deg2rad(ϕ))))
moves = Dict(
    'N' => α -> c -> begin (p,d) = c; (p + Complex{Int}(0,α), d) end,
    'S' => α -> c -> begin (p,d) = c; (p + Complex{Int}(0,-α), d) end,
    'E' => α -> c -> begin (p,d) = c; (p + Complex{Int}(α,0), d) end,
    'W' => α -> c -> begin (p,d) = c; (p + Complex{Int}(-α,0), d) end,
    'L' => ϕ -> c -> begin (p,d) = c; (p, d * rot(ϕ)) end,
    'R' => ϕ -> c -> begin (p,d) = c; (p, d * rot(-ϕ)) end,
    'F' => s -> c -> begin (p,d) = c; (p + s * d, d) end,
)
dist(p::Complex) = reim(p) .|> abs |> sum
dist(t::Tuple) = dist(first(t))
dist(t1::Tuple, t2::Tuple) = dist(t2[1] - t1[1]);

##
actions = map(eachline("test.txt")) do line
    moves[line[1]](parse(Int, line[2:end]))
end;

##
start = Complex{Int}(0,0), Complex{Int}(1,0)
finish = foldl(actions; init = start) do c, a
    a(c)
end
@show finish
dist(start, finish)

## PARTE 2
moves = Dict(
    'N' => α -> c -> begin (p,d) = c; (p, d + Complex{Int}(0,α)) end,
    'S' => α -> c -> begin (p,d) = c; (p, d + Complex{Int}(0,-α)) end,
    'E' => α -> c -> begin (p,d) = c; (p, d + Complex{Int}(α,0)) end,
    'W' => α -> c -> begin (p,d) = c; (p, d + Complex{Int}(-α,0)) end,
    'L' => ϕ -> c -> begin (p,d) = c; (p, d * rot(ϕ)) end,
    'R' => ϕ -> c -> begin (p,d) = c; (p, d * rot(-ϕ)) end,
    'F' => s -> c -> begin
        (p,d) = c
        (p + s * d, d)
    end,
)
actions = map(eachline("input.txt")) do line
    moves[line[1]](parse(Int, line[2:end]))
end;

##
start = Complex{Int}(0,0), Complex{Int}(10,1)
finish = foldl(actions; init = start) do c, a
    a(c)
end
@show finish
dist(start, finish)


text:
using IndirectArrays
using Base.Test, Colors

colors = [RGB(1,0,0), RGB(0,1,0), RGB(0,0,1)]
index0 = [1 2;
          3 1]
for index in (index0, map(Int16, index0))
    A = IndirectArray(index, colors)
    @test eltype(A) == RGB{U8}
    @test size(A) == (2,2)
    @test ndims(A) == 2
    @test A[1,1] === A[1] === RGB(1,0,0)
    @test A[2,1] === A[2] === RGB(0,0,1)
    @test A[1,2] === A[3] === RGB(0,1,0)
    @test A[2,2] === A[4] === RGB(1,0,0)
    @test isa(eachindex(A), AbstractUnitRange)
end

# Bounds checking upon construction
index_ob = copy(index0)
index_ob[1] = 5   # out-of-bounds
unsafe_ia(idx, vals) = (@inbounds ret = IndirectArray(idx, vals); ret)
  safe_ia(idx, vals) = (ret = IndirectArray(idx, vals); ret)
@test_throws BoundsError safe_ia(index_ob, colors)
# This requires inlining, which means it fails on Travis since we turn
# off inlining for better coverage stats
# B = unsafe_ia(index_ob, colors)
# @test_throws BoundsError B[1]
# @test B[2] == RGB(0,0,1)


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

mutable struct StorageAccountPropertiesUpdateParameters <: SwaggerModel
    customDomain::Any # spec type: Union{ Nothing, CustomDomain } # spec name: customDomain
    encryption::Any # spec type: Union{ Nothing, Encryption } # spec name: encryption
    accessTier::Any # spec type: Union{ Nothing, String } # spec name: accessTier
    azureFilesIdentityBasedAuthentication::Any # spec type: Union{ Nothing, AzureFilesIdentityBasedAuthentication } # spec name: azureFilesIdentityBasedAuthentication
    supportsHttpsTrafficOnly::Any # spec type: Union{ Nothing, Bool } # spec name: supportsHttpsTrafficOnly
    networkAcls::Any # spec type: Union{ Nothing, NetworkRuleSet } # spec name: networkAcls
    largeFileSharesState::Any # spec type: Union{ Nothing, String } # spec name: largeFileSharesState
    routingPreference::Any # spec type: Union{ Nothing, RoutingPreference } # spec name: routingPreference
    allowBlobPublicAccess::Any # spec type: Union{ Nothing, Bool } # spec name: allowBlobPublicAccess
    minimumTlsVersion::Any # spec type: Union{ Nothing, String } # spec name: minimumTlsVersion

    function StorageAccountPropertiesUpdateParameters(;customDomain=nothing, encryption=nothing, accessTier=nothing, azureFilesIdentityBasedAuthentication=nothing, supportsHttpsTrafficOnly=nothing, networkAcls=nothing, largeFileSharesState=nothing, routingPreference=nothing, allowBlobPublicAccess=nothing, minimumTlsVersion=nothing)
        o = new()
        validate_property(StorageAccountPropertiesUpdateParameters, Symbol("customDomain"), customDomain)
        setfield!(o, Symbol("customDomain"), customDomain)
        validate_property(StorageAccountPropertiesUpdateParameters, Symbol("encryption"), encryption)
        setfield!(o, Symbol("encryption"), encryption)
        validate_property(StorageAccountPropertiesUpdateParameters, Symbol("accessTier"), accessTier)
        setfield!(o, Symbol("accessTier"), accessTier)
        validate_property(StorageAccountPropertiesUpdateParameters, Symbol("azureFilesIdentityBasedAuthentication"), azureFilesIdentityBasedAuthentication)
        setfield!(o, Symbol("azureFilesIdentityBasedAuthentication"), azureFilesIdentityBasedAuthentication)
        validate_property(StorageAccountPropertiesUpdateParameters, Symbol("supportsHttpsTrafficOnly"), supportsHttpsTrafficOnly)
        setfield!(o, Symbol("supportsHttpsTrafficOnly"), supportsHttpsTrafficOnly)
        validate_property(StorageAccountPropertiesUpdateParameters, Symbol("networkAcls"), networkAcls)
        setfield!(o, Symbol("networkAcls"), networkAcls)
        validate_property(StorageAccountPropertiesUpdateParameters, Symbol("largeFileSharesState"), largeFileSharesState)
        setfield!(o, Symbol("largeFileSharesState"), largeFileSharesState)
        validate_property(StorageAccountPropertiesUpdateParameters, Symbol("routingPreference"), routingPreference)
        setfield!(o, Symbol("routingPreference"), routingPreference)
        validate_property(StorageAccountPropertiesUpdateParameters, Symbol("allowBlobPublicAccess"), allowBlobPublicAccess)
        setfield!(o, Symbol("allowBlobPublicAccess"), allowBlobPublicAccess)
        validate_property(StorageAccountPropertiesUpdateParameters, Symbol("minimumTlsVersion"), minimumTlsVersion)
        setfield!(o, Symbol("minimumTlsVersion"), minimumTlsVersion)
        o
    end
end # type StorageAccountPropertiesUpdateParameters

const _property_map_StorageAccountPropertiesUpdateParameters = Dict{Symbol,Symbol}(Symbol("customDomain")=>Symbol("customDomain"), Symbol("encryption")=>Symbol("encryption"), Symbol("accessTier")=>Symbol("accessTier"), Symbol("azureFilesIdentityBasedAuthentication")=>Symbol("azureFilesIdentityBasedAuthentication"), Symbol("supportsHttpsTrafficOnly")=>Symbol("supportsHttpsTrafficOnly"), Symbol("networkAcls")=>Symbol("networkAcls"), Symbol("largeFileSharesState")=>Symbol("largeFileSharesState"), Symbol("routingPreference")=>Symbol("routingPreference"), Symbol("allowBlobPublicAccess")=>Symbol("allowBlobPublicAccess"), Symbol("minimumTlsVersion")=>Symbol("minimumTlsVersion"))
const _property_types_StorageAccountPropertiesUpdateParameters = Dict{Symbol,String}(Symbol("customDomain")=>"CustomDomain", Symbol("encryption")=>"Encryption", Symbol("accessTier")=>"String", Symbol("azureFilesIdentityBasedAuthentication")=>"AzureFilesIdentityBasedAuthentication", Symbol("supportsHttpsTrafficOnly")=>"Bool", Symbol("networkAcls")=>"NetworkRuleSet", Symbol("largeFileSharesState")=>"String", Symbol("routingPreference")=>"RoutingPreference", Symbol("allowBlobPublicAccess")=>"Bool", Symbol("minimumTlsVersion")=>"String")
Base.propertynames(::Type{ StorageAccountPropertiesUpdateParameters }) = collect(keys(_property_map_StorageAccountPropertiesUpdateParameters))
Swagger.property_type(::Type{ StorageAccountPropertiesUpdateParameters }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_StorageAccountPropertiesUpdateParameters[name]))}
Swagger.field_name(::Type{ StorageAccountPropertiesUpdateParameters }, property_name::Symbol) =  _property_map_StorageAccountPropertiesUpdateParameters[property_name]

const _allowed_StorageAccountPropertiesUpdateParameters_accessTier = ["Hot", "Cool"]

const _allowed_StorageAccountPropertiesUpdateParameters_largeFileSharesState = ["Disabled", "Enabled"]

const _allowed_StorageAccountPropertiesUpdateParameters_minimumTlsVersion = ["TLS1_0", "TLS1_1", "TLS1_2"]

function check_required(o::StorageAccountPropertiesUpdateParameters)
    true
end

function validate_property(::Type{ StorageAccountPropertiesUpdateParameters }, name::Symbol, val)
    if name === Symbol("accessTier")
        Swagger.validate_param(name, "StorageAccountPropertiesUpdateParameters", :enum, val, _allowed_StorageAccountPropertiesUpdateParameters_accessTier)
    end
    if name === Symbol("largeFileSharesState")
        Swagger.validate_param(name, "StorageAccountPropertiesUpdateParameters", :enum, val, _allowed_StorageAccountPropertiesUpdateParameters_largeFileSharesState)
    end
    if name === Symbol("minimumTlsVersion")
        Swagger.validate_param(name, "StorageAccountPropertiesUpdateParameters", :enum, val, _allowed_StorageAccountPropertiesUpdateParameters_minimumTlsVersion)
    end
end


text:
# packages -
using Dates
using Optim
using JSON
using DataFrames
using Statistics
using LsqFit
using Reexport
@reexport using PooksoftBase

# include my code -
include("./base/Types.jl")
include("./base/Checks.jl")
include("./base/Intrinsic.jl")
include("./base/Binary.jl")
include("./base/Ternary.jl")
include("./base/Greeks.jl")
include("./base/Compute.jl")
include("./base/Factory.jl")
include("./base/Volatility.jl")
include("./base/Utility.jl")
include("./base/Longstaff.jl")
include("./base/Breakeven.jl")


text:
#=
Illustrates preimages of functions

@author : Spencer Lyon <spencer.lyon@nyu.edu>
          Victoria Gregory <victoria.gregory@nyu.edu>

@date: 07/09/2014
=#
using Plots
pyplot()
using LaTeXStrings

f(x) = 0.6 .* cos(4.0 .* x) .+ 1.3

xmin, xmax = -1.0, 1.0
Nx = 160
x = linspace(xmin, xmax, Nx)
y = f(x)
ya, yb = minimum(y), maximum(y)

p1 = plot(x, y, color=:black, label=[L"$f$" ""], grid=false)
plot!(x, ya*ones(Nx, 1), fill_between=yb*ones(Nx, 1),
     fillalpha=0.1, color=:blue, label="", lw=0)
plot!(zeros(2, 2), [ya ya; yb yb], lw=3, color=:blue, label=[L"range of $f$" ""])
annotate!(0.04, -0.3, L"$0$", ylims=(-0.6, 3.2))
vline!([0], color=:black, label="")
hline!([0], color=:black, label="")
plot!(foreground_color_axis=:white, foreground_color_text=:white,
      foreground_color_border=:white)

ybar = 1.5
plot!(x, x .* 0 .+ ybar, color=:black, linestyle=:dash, label="")
annotate!(0.05, 0.8 * ybar, L"$y$")

x_vals = Array(Float64, 2, 4)
y_vals = Array(Float64, 2, 4)
labels = []
for (i, z) in enumerate([-0.35, 0.35])
  x_vals[:, 2*i-1] = z*ones(2, 1)
  y_vals[2, 2*i-1] = f(z)
  labels = [labels; (z, -0.2, LaTeXString("\$x_$i\$"))]
end
plot!(x_vals, y_vals, color=:black, linestyle=:dash, label="", annotation=labels)

p2 = plot(x, y, color=:black, label=[L"$f$" ""], grid=false)
plot!(x, ya*ones(Nx, 1), fill_between=yb*ones(Nx, 1),
     fillalpha=0.1, color=:blue, label="", lw=0)
plot!(zeros(2, 2), [ya ya; yb yb], lw=3, color=:blue, label=[L"range of $f$" ""])
annotate!(0.04, -0.3, L"$0$", ylims=(-0.6, 3.2))
vline!([0], color=:black, label="")
hline!([0], color=:black, label="")
plot!(foreground_color_axis=:white, foreground_color_text=:white,
      foreground_color_border=:white)

ybar = 2.6
plot!(x, x .* 0 .+ ybar, color=:black, linestyle=:dash, legend=:none)
annotate!(0.04, 0.91 * ybar, L"$y$")

plot(p1, p2, layout=(2, 1), size=(600,700))

max_stars_repo_name: mwhchen/quantecon

max_issues_repo_name: mwhchen/quantecon

max_forks_repo_name: mwhchen/quantecon


text:
@assert (3-im) * simpk' ==  d"(3-im) * (1-3im)< 1 |"
@assert (3-im) * simpk' == (d"(3+im)*(1+3im)| 1 >")'
@assert (-im * projop')[1,1] == 0-10im

@assert length(filter((o,v) -> klabel(o)==StateLabel(3), op')) == 3
@assert b'*b'*b' == (b^3)'
@assert (op'*op') == (op*op)'
@assert convert(OpSum, op') == maplabels(ctranspose, mapcoeffs(ctranspose, op))
@assert op' == map((k,v)->(k', v'), op)
@assert op[2,1] == 8-4im
@assert (k*k)*(b*b) == tensor(k*b,k*b)
@assert op+op == 2 * op
@assert op-op == 0 * op
@assert op-op+op == op
@assert filternz!(xsubspace(tensor(op,op,op),3)) == d"(16 - 88im)| 1,1,1 >< 1,1,1 |"

op_copy = copy(op)
op_copy[3,0] = 32+im
@assert op_copy[3,0] == 32+im
@assert tensor(op_copy, op_copy')[(3,1),(1,2)] == 120

@test_approx_eq norm(qubits) 1
@assert ptrace(belldens, 1) == ptrace(belldens, 2)

@assert d" act_on(< 1 |, | 'a','b','c' >, 2) == < 1 | 'b' >| 'a','c' > "
@assert d" act_on(| 1 >, < 'a','b','c' |, 2) == conj(< 1 | 'b' >)< 'a','c' | "

tranop = d" | 'a','b','c' >< 'd','e','f' | + 2 * | 'i','j','k' >< 'l','m','n' | "
@assert ptranspose(tranop, 2) == d" | 'a','e','c' >< 'd','b','f' | + 2 * | 'i','m','k' >< 'l','j','n' | "


text:
function execute(problem::PotentialFlowProblem{:beam};kwargs...)

  # Parameters
  L = 2*π
  H = 1.0
  xb₀ = 0.8*π
  xb₁ = 1.2*π
  n = 100
  order = 2
  g = 9.81
  ξ = 0.001
  λ = L/2
  k = 2*π/L
  h = L/n
  ω = √(g*k*tanh(k*H))
  t₀ = 0.0
  tf = 2*π
  Δt = 1/(50*ω)#h/(10*λ*ω)
  θ = 0.5

  # Exact solution
  ϕₑ(x,t) = ω/k * ξ * (cosh(k*(x[2]))) / sinh(k*H) * sin(k*x[1] - ω*t)
  ηₑ(x,t) = ξ * cos(k*x[1] - ω*t)
  ϕₑ(t::Real) = x -> ϕₑ(x,t)
  ηₑ(t::Real) = x -> ηₑ(x,t)

  # Interior Domain
  domain = (0,L,0,H)
  partition = (n,n/10)
  model_Ω = CartesianDiscreteModel(domain,partition;isperiodic=(true,false))
  Ω = Triangulation(model_Ω)

  # Surface domains
  function is_beam(coords)
    n = length(coords)
		x = (1/n)*sum(coords)
    (xb₀ <= x[1] <= xb₁ ) * ( x[2] ≈ H )
  end

  # Boundary domain
  labels = get_face_labeling(model_Ω)
  bgface_to_mask = get_face_mask(labels,[3,4,6],1)
  Γface_to_bgface = findall(bgface_to_mask)
  model_Γ = BoundaryDiscreteModel(Polytope{1},model_Ω,Γface_to_bgface)
  Γ = Triangulation(model_Γ)

  # Get beam - free surface triangulations
  Γface_coords = get_cell_coordinates(Γ)
  Γface_mask = lazy_map(is_beam,Γface_coords)
  Γbface_Γface = findall(Γface_mask)
  Γfface_Γface = findall(!,Γface_mask)
  Γb = BoundaryTriangulation(model_Ω,view(Γface_to_bgface,Γbface_Γface))
  Γf = BoundaryTriangulation(model_Ω,view(Γface_to_bgface,Γfface_Γface))
  model_Γb = BoundaryDiscreteModel(Polytope{1},model_Ω,view(Γface_to_bgface,Γbface_Γface))
  Λb = SkeletonTriangulation(model_Γ,Γface_mask)
  nΛb = get_normal_vector(Λb)

  writevtk(model_Γb,"mGb")
  writevtk(Γb,"Gb")
  writevtk(Γf,"Gf")
  writevtk(Λb,"Lb")

  # Quadratures
  dΩ = Measure(Ω,2*order)
  dΓb = Measure(Γb,2*order)
  dΓf = Measure(Γf,2*order)
  dΛb = Measure(Λb,2*order)

  qΛ = get_cell_points(CellQuadrature(Λb,1))
  mean_mask = CellField(mean(CellField(Γface_mask,Γ)),Λb)
  writevtk(Λb,"Lb",cellfields=["mask"=>mean_mask])

  # FE spaces
  reffe = ReferenceFE(lagrangian,Float64,order)
  V_Ω = TestFESpace(model_Ω,reffe,conformity=:H1)
  V_Γ = TestFESpace(model_Γ,reffe,conformity=:H1)
  U_Ω = TransientTrialFESpace(V_Ω)
  U_Γ = TransientTrialFESpace(V_Γ)
  X = TransientMultiFieldFESpace([U_Ω,U_Γ])
  Y = MultiFieldFESpace([V_Ω,V_Γ])

  # Weak form
  γ = 1
  α = 4/Δt^2 + 2/Δt
  m((ϕtt,ηtt),(w,v)) = ∫( v*ηtt )dΓb
  c((ϕt,ηt),(w,v)) = ∫( 0.5*(α/g*w + v)*ϕt - w*ηt )dΓb + ∫( 0.5*(α/g*w + v)*ϕt - w*ηt )dΓf
  a((ϕ,η),(w,v)) = ∫( ∇(ϕ)⋅∇(w) )dΩ + ∫( Δ(v)*Δ(η) + 0.5*(α/g*w + v)*g*η )dΓb + ∫( 0.5*(α/g*w + v)*g*η )dΓf +
                   ∫((mean_mask==1)*( - mean(Δ(η))*jump(∇(v)⋅nΛb) - jump(∇(η)⋅nΛb)*mean(Δ(v)) + γ/h*jump(∇(η)⋅nΛb)*jump(∇(v)⋅nΛb)) )dΛb
  b((w,v)) = ∫( 0.5*(α/g*w + v)*(-0.0) )dΓb
  op = TransientConstantFEOperator(m,c,a,b,X,Y)

  # Solver
  ls = LUSolver()
  odes = Newmark(ls,Δt,0.5,0.25)
  solver = TransientFESolver(odes)

  # Initial solution
  x₀ = interpolate_everywhere([ϕₑ(0.0),ηₑ(0.0),],X(0.0))
  v₀ = interpolate_everywhere([∂t(ϕₑ)(0.0),∂t(ηₑ)(0.0)],X(0.0))
  a₀ = interpolate_everywhere([∂tt(ϕₑ)(0.0),∂tt(ηₑ)(0.0)],X(0.0))

  # Solution
  sol_t = solve(solver,op,(x₀,v₀,a₀),t₀,tf)
  #sol_t = solve(solver,op,x₀,t₀,tf)

  # Post-process
  l2_Ω(w) = √(∑(∫(w*w)dΩ))
  l2_Γb(v) = √(∑(∫(v*v)dΓb))
  l2_Γf(v) = √(∑(∫(v*v)dΓf))
  E_kin(w) = 0.5*∑( ∫(∇(w)⋅∇(w))dΩ )
  E_pot_b(v) = g*0.5*∑( ∫(v*v)dΓb )
  E_pot_f(v) = g*0.5*∑( ∫(v*v)dΓf )
  Eₑ = 0.5*g*ξ^2*L

  folderName = "PFlow-results"
  if !isdir(folderName)
    mkdir(folderName)
  end
  filePath_Ω = folderName * "/fields_O"
  filePath_Γb = folderName * "/fields_Gb"
  filePath_Γf = folderName * "/fields_Gf"
  pvd_Ω = paraview_collection(filePath_Ω, append=false)
  pvd_Γb = paraview_collection(filePath_Γb, append=false)
  pvd_Γf = paraview_collection(filePath_Γf, append=false)
  for ((ϕn,ηn),tn) in sol_t
    # E = E_kin(ϕn) + E_pot_b(ηn) + E_pot_f(ηn)
    # error_ϕ = l2_Ω(ϕn-ϕₑ(tn))
    # error_η = l2_Γb(ηn-ηₑ(tn)) + l2_Γf(ηn-ηₑ(tn))
    # println(E/Eₑ," ", error_ϕ," ",error_η)

    pvd_Ω[tn] = createvtk(Ω,filePath_Ω * "_$tn.vtu",cellfields = ["phi" => ϕn],order=2)
    pvd_Γb[tn] = createvtk(Γb,filePath_Γb * "_$tn.vtu",cellfields = ["eta" => ηn],nsubcells=10)
    pvd_Γf[tn] = createvtk(Γf,filePath_Γf * "_$tn.vtu",cellfields = ["eta" => ηn],nsubcells=10)
  end
  vtk_save(pvd_Ω)
  vtk_save(pvd_Γb)
  vtk_save(pvd_Γf)
end


text:
include("code-2-lib/distributed.jl") # Insert the contents of distributed.jl here.
add_procs() # Calls the fucntin add_procs() from distributed.jl library

@everywhere using Pkg # Using the built-in package manager in julia
@everywhere Pkg.activate(".") # Activating the package manager for required usage

include("solution2-solver.jl") # Insert the contents of solution2-solver.jl here.

function MyJulia2(InFile1::String, InFile2::String, InFile3::String, InFile4::String, TimeLimitInSeconds::Int64, ScoringMethod::Int64, NetworkModel::String) # Calls the function named MyJulia2 with given arguments
    println("running MyJulia2") # Print the statement
    println("  $(InFile1)") # Print the contents of file
    println("  $(InFile2)")
    println("  $(InFile3)")
    println("  $(InFile4)")
    println("  $(TimeLimitInSeconds)")
    println("  $(ScoringMethod)")
    println("  $(NetworkModel)")

    compute_solution2(InFile1, InFile2, InFile3, InFile4, TimeLimitInSeconds, ScoringMethod, NetworkModel) # Calls the functin named compute_solution2 with given arguments
end


text:
# This file is a part of ValueShapes.jl, licensed under the MIT License (MIT).

using ValueShapes
using Test

using Distributions

@testset "distributions_support" begin
    @test @inferred(valshape(Normal())) == ScalarShape{Real}()
    @test @inferred(valshape(MvNormal([2. 1.; 1. 3.]))) == ArrayShape{Real}(2)
    @test @inferred(valshape(MatrixBeta(4, 6, 6))) == ArrayShape{Real}(4, 4)
end


text:
const MAX_SHOW_COMPONENTS = 10
const MAX_SHOW_FORECASTS = 10
const MAX_SHOW_FORECAST_INITIAL_TIMES = 1

function Base.summary(components::Components)
    return "$(typeof(components)): $(get_num_components(components))"
end

function Base.show(io::IO, components::Components)
    i = 1
    for component in iterate_components(components)
        if i <= MAX_SHOW_COMPONENTS
            show(io, component)
            println(io)
        end
        i += 1
    end

    if i > MAX_SHOW_COMPONENTS
        num = i - MAX_SHOW_COMPONENTS
        println(io, "\n***Omitted $num components***\n")
    end
end

function Base.show(io::IO, ::MIME"text/plain", components::Components)
    num_components = get_num_components(components)
    println(io, "Components")
    println(io, "==========")
    println(io, "Num components: $num_components")
    if num_components > 0
        println(io)
        show_components_table(io, components, backend = :auto)
    end
end

function Base.show(io::IO, ::MIME"text/html", components::Components)
    num_components = get_num_components(components)
    println(io, "
Components
")
    println(io, "
Num components: $num_components
")
    if num_components > 0
        show_components_table(io, components, backend = :html)
    end
end

function Base.summary(container::TimeSeriesContainer)
    return "$(typeof(container)): $(length(container))"
end

function Base.show(io::IO, ::MIME"text/plain", container::TimeSeriesContainer)
    println(io, summary(container))
    for key in keys(container.data)
        println(io, "$(key.time_series_type): name=$(key.name)")
    end
end

function Base.summary(time_series::TimeSeriesData)
    return "$(typeof(time_series)) time_series ($(length(time_series)))"
end

function Base.summary(time_series::TimeSeriesMetadata)
    return "$(typeof(time_series)) time_series"
end

function Base.show(io::IO, data::SystemData)
    show(io, data.components)
    println(io, "\n")
    return show(io, data.time_series_params)
end

function Base.show(io::IO, ::MIME"text/plain", data::SystemData)
    show(io, MIME"text/plain"(), data.components)
    println(io, "\n")
    show_time_series_data(io, data, backend = :auto)
    show(io, data.time_series_params)
end

function Base.show(io::IO, ::MIME"text/html", data::SystemData)
    show(io, MIME"text/html"(), data.components)
    println(io, "\n")
    show_time_series_data(io, data, backend = :html)
    show(io, data.time_series_params)
end

function show_time_series_data(io::IO, data::SystemData; kwargs...)
    component_count, ts_count, forecast_count = get_time_series_counts(data)
    res = get_time_series_resolution(data)
    res = res <= Dates.Minute(1) ? Dates.Second(res) : Dates.Minute(res)

    header = ["Property", "Value"]
    table = [
        "Components with time series data" string(component_count)
        "Total StaticTimeSeries" string(ts_count)
        "Total Forecasts" string(forecast_count)
        "Resolution" string(res)
    ]
    if component_count == 0
        return
    end

    if forecast_count > 0
        initial_times = get_forecast_initial_times(data)
        table2 = [
            "First initial time" string(first(initial_times))
            "Last initial time" string(last(initial_times))
            "Horizon" string(get_forecast_horizon(data))
            "Interval" string(Dates.Minute(get_forecast_interval(data)))
            "Forecast window count" string(get_forecast_window_count(data))
        ]
        table = vcat(table, table2)
    end

    PrettyTables.pretty_table(
        io,
        table;
        header = header,
        title = "Time Series Summary",
        alignment = :l,
        kwargs...,
    )
    return
end

function Base.summary(ist::InfrastructureSystemsComponent)
    # All InfrastructureSystemsComponent subtypes are supposed to implement get_name.
    # Some don't.  They need to override this function.
    return "$(get_name(ist)) ($(typeof(ist)))"
end

function Base.show(io::IO, ::MIME"text/plain", system_units::SystemUnitsSettings)
    print(io, summary(system_units), ":")
    for name in fieldnames(typeof(system_units))
        val = getfield(system_units, name)
        print(io, "\n      ", name, ": ", val)
    end
end

function Base.show(io::IO, ::MIME"text/plain", ist::InfrastructureSystemsComponent)
    print(io, summary(ist), ":")
    for (name, field_type) in zip(fieldnames(typeof(ist)), fieldtypes(typeof(ist)))
        obj = getfield(ist, name)
        if obj isa InfrastructureSystemsInternal
            continue
        elseif obj isa TimeSeriesContainer || obj isa InfrastructureSystemsType
            val = summary(getfield(ist, name))
        elseif obj isa Vector{<:InfrastructureSystemsComponent}
            val = summary(getfield(ist, name))
        else
            val = getfield(ist, name)
        end
        # Not allowed to print `nothing`
        if isnothing(val)
            val = "nothing"
        end
        print(io, "\n   ", name, ": ", val)
    end
end

function Base.show(io::IO, ist::InfrastructureSystemsComponent)
    print(io, string(nameof(typeof(ist))), "(")
    is_first = true
    for (name, field_type) in zip(fieldnames(typeof(ist)), fieldtypes(typeof(ist)))
        if field_type <: TimeSeriesContainer || field_type <: InfrastructureSystemsInternal
            continue
        else
            val = getfield(ist, name)
        end
        if is_first
            is_first = false
        else
            print(io, ", ")
        end
        print(io, val)
    end
    print(io, ")")
    return
end

function Base.show(io::IO, ::MIME"text/plain", it::FlattenIteratorWrapper)
    println(io, "$(eltype(it)) Counts: ")
    for (ctype, count) in _get_type_counts(it)
        println(io, "$ctype: $count")
    end
end

function _get_type_counts(it::FlattenIteratorWrapper)
    data = SortedDict()
    for component in it
        ctype = string(typeof(component))
        if !haskey(data, ctype)
            data[ctype] = 1
        else
            data[ctype] += 1
        end
    end

    return data
end

function show_components_table(io::IO, components::Components; kwargs...)
    header = ["Type", "Count", "Has Static Time Series", "Has Forecasts"]
    data = Array{Any, 2}(undef, length(components.data), length(header))

    type_names = [(strip_module_name(string(x)), x) for x in keys(components.data)]
    sort!(type_names, by = x -> x[1])
    for (i, (type_name, type)) in enumerate(type_names)
        vals = components.data[type]
        has_sts = false
        has_forecasts = false
        for val in values(vals)
            if has_time_series(val, StaticTimeSeries)
                has_sts = true
            end
            if has_time_series(val, Forecast)
                has_forecasts = true
            end
            if has_sts && has_forecasts
                break
            end
        end
        data[i, 1] = type_name
        data[i, 2] = length(vals)
        data[i, 3] = has_sts
        data[i, 4] = has_forecasts
    end

    PrettyTables.pretty_table(io, data; header = header, alignment = :l, kwargs...)
    return
end

function show_components(
    io::IO,
    components::Components,
    component_type::Type{<:InfrastructureSystemsComponent},
    additional_columns::Union{Dict, Vector} = [];
    kwargs...,
)
    if !isconcretetype(component_type)
        error("$component_type must be a concrete type")
    end

    title = string(strip_module_name(component_type))
    header = ["name"]
    has_available = false
    if :available in fieldnames(component_type)
        push!(header, "available")
        has_available = true
    end

    if additional_columns isa Dict
        columns = sort!(collect(keys(additional_columns)))
    else
        columns = additional_columns
    end

    for column in columns
        push!(header, string(column))
    end

    comps = get_components(component_type, components)
    data = Array{Any, 2}(undef, length(comps), length(header))
    for (i, component) in enumerate(comps)
        data[i, 1] = get_name(component)
        j = 2
        if has_available
            data[i, 2] = getproperty(component, :available)
            j += 1
        end

        if additional_columns isa Dict
            for column in columns
                data[i, j] = additional_columns[column](component)
                j += 1
            end
        else
            for column in additional_columns
                getter_name = Symbol("get_$column")
                parent = parentmodule(component_type)
                # This logic enables application of system units in PowerSystems through
                # its getter functions.
                val = getproperty(component, column)
                if val isa TimeSeriesContainer ||
                   val isa InfrastructureSystemsType ||
                   val isa Vector{<:InfrastructureSystemsComponent}
                    val = summary(val)
                elseif hasproperty(parent, getter_name)
                    getter_func = getproperty(parent, getter_name)
                    val = getter_func(component)
                end
                data[i, j] = val
                j += 1
            end
        end
    end

    PrettyTables.pretty_table(
        io,
        data;
        header = header,
        title = title,
        alignment = :l,
        kwargs...,
    )
    return
end

## This function takes in a time period or date period and returns a compound period
function convert_compound_period(period::Union{Dates.TimePeriod, Dates.DatePeriod})
    period = time_period_conversion(period)

    milli_weeks = period - (period % Dates.Millisecond(604800000))
    weeks = convert(Dates.Week, milli_weeks)
    period -= milli_weeks

    milli_days = period - (period % Dates.Millisecond(86400000))
    days = convert(Dates.Day, milli_days)
    period -= milli_days

    milli_hours = period - (period % Dates.Millisecond(3600000))
    hours = convert(Dates.Hour, milli_hours)
    period -= milli_hours

    milli_minutes = period - (period % Dates.Millisecond(60000))
    minutes = convert(Dates.Minute, milli_minutes)
    period -= milli_minutes

    seconds = period - (period % Dates.Millisecond(1000)) # finding the seconds
    seconds = convert(Dates.Second, seconds)
    remainder = period % Dates.Millisecond(1000) #finding the remainding milliseconds
    total = weeks + days + hours + minutes + seconds + remainder
    return total
end


text:
module Objectof

export objectof, hasobjectof

objectof(::Type{Any}) = nothing
objectof(::Type{Nothing}) = nothing
objectof(::Type{Missing}) = missing
objectof(::Type{T}) where T<:Number = zero(T)
objectof(::Type{T}) where T<:AbstractString = (T == AbstractString ? String : T)("")
objectof(::Type{Union{}}) = Union{}()
objectof(::Type{<:Ref{T}}) where T = Ref{T}(objectof(T))
objectof(::Type{T}) where T<:AbstractRange = convert(T, 1:0)
objectof(::Type{T}) where T<:Function = T === Function ? identity : T.instance

function objectof(::Type{Array{T,N}}) where {T,N}
    N isa Integer && N >= 0 || throw(ArgumentError("array dimensions must be integers >= 0"))
    A = Array{T,N}(undef, zeros(Int, N)...)
    if N == 0
        A[] = objectof(T)
    end
    A
end

function objectof(::Type{T}) where T
    if typeof(T) == Union
        hasobjectof(T.a) ? objectof(T.a) : objectof(T.b)
    elseif typeof(T) == DataType
        T(objectof.(T.types)...)
    elseif typeof(T) == UnionAll
        objectof(T{T.var.ub})
    else
        throw(ArgumentError("cannot create object of type $T"))
    end
end

hasobjectof(::Type{Any}) = true
hasobjectof(::Type{Nothing}) = true
hasobjectof(::Type{Missing}) = true
hasobjectof(::Type{T}) where T<:Number = true
hasobjectof(::Type{T}) where T<:AbstractString = true
hasobjectof(::Type{Union{}}) = false
hasobjectof(::Type{<:Ref{T}}) where T = hasobjectof(T)
hasobjectof(::Type{T}) where T<:AbstractRange = T in (AbstractRange, OrdinalRange, AbstractStepRange) || typeof(T) != UnionAll
hasobjectof(::Type{T}) where T<:Function = true

hasobjectof(::Type{Array{T,N}}) where {T,N} = true

function hasobjectof(::Type{T}) where T
    if typeof(T) == Union
        hasobjectof(T.a) || hasobjectof(T.b)
    elseif typeof(T) == DataType
        all(hasobjectof.(T.types)) && hasmethod(T, tuple(T.types...))
    elseif typeof(T) == UnionAll
        hasobjectof(T{T.var.ub})
    else
        false
    end
end

end # module


text:
function write_solutions(OPFpbs, variables_jump2, outpath)
    open(joinpath(outpath,"solution1.txt"),"w") do f
      write(f, "--generation dispatch\nbus id,unit id,pg(MW),qg(MVar) \n");
      for (busname, elems) in OPFpbs[basecase_scenario_name()].ds.bus
        for (elemname,element) in elems
                if typeof(element) == GOCGenerator
                  bus =  element.busid
                  gen = element.id
                  Pgen = getvalue(variables_jump2[variable_name("Sgen", busname, elemname, basecase_scenario_name())*"_Re"])
                  Qgen = getvalue(variables_jump2[variable_name("Sgen", busname, elemname, basecase_scenario_name())*"_Im"])
                  write(f, "$bus, $gen, $Pgen, $Qgen\n")
                end
         end
       end
      write(f,"--end of generation dispatch\n");
    end

    Qgen_scen_values = Dict{Tuple{String,Any,Int64,Any}, Float64}()
    volt_values = Dict{Tuple{String,Int64}, Tuple{Float64, Float64}}()
    delta_values = Dict{String,Float64}()
    Slink_values = Dict{Tuple{String,Any,Int64,Int64,Any}, Tuple{Float64, Float64, Float64, Float64}}()

    for (scenario, OPFpb) in OPFpbs
      if scenario==basecase_scenario_name()
        scenario_id = "0"
        ##volt values
         for (busname, elems) in OPFpb.ds.bus
           for (elemid, element) in elems
             if typeof(element) == GOCVolt
               bus = element.busid
               V_re = getvalue(variables_jump2[variable_name("VOLT", busname, "", scenario)*"_Re"])
               V_im = getvalue(variables_jump2[variable_name("VOLT", busname, "", scenario)*"_Im"])
               V_mod = abs(V_re + V_im * im)
               V_theta = angle(V_re + V_im * im)*180/pi
               volt_values[(scenario_id,bus)] = (V_mod, V_theta)
             end
           end
         end
      else
        ##delta values
        scenario_id = String(matchall(r"\d+", scenario)[1])
        delta_values[scenario_id] = getvalue(variables_jump2[get_delta_varname(scenario)])
        ##volt values
         for (busname, elems) in OPFpb.ds.bus
           for (elemid, element) in elems
             if typeof(element) == GOCVolt
               bus = element.busid
               # bus = String(matchall(r"\d+", element.busname)[1])
               V_re = getvalue(variables_jump2[variable_name("VOLT", busname, "", scenario)*"_Re"])
               V_im = getvalue(variables_jump2[variable_name("VOLT", busname, "", scenario)*"_Im"])
               V_mod = abs(V_re + V_im * im)
               V_theta = angle(V_re + V_im * im)*180/pi
               volt_values[(scenario_id,bus)] = (V_mod, V_theta)
             elseif typeof(element) == GOCGenerator
               bus =  element.busid
               gen = element.id
               # Pgen = getvalue(variables_jump2[variable_name("Sgen", busname, elemid, scenario)*"_Re"])
               Qgen = getvalue(variables_jump2[variable_name("Sgen", busname, elemid, scenario)*"_Im"])
               # println("bus ",bus, " : ", element.power_min," << ", Pgen + im * Qgen," << ", element.power_max)
               Qgen_scen_values[(scenario_id,gen,bus, gen)] = Qgen
             end
           end
         end
      end

      for (link, elems) in OPFpb.ds.link
        orig = link.orig
        dest = link.dest
        Vo_re = getvalue(variables_jump2[variable_name("VOLT", orig, "", scenario)*"_Re"])
        Vo_im = getvalue(variables_jump2[variable_name("VOLT", orig, "", scenario)*"_Im"])
        Vd_re = getvalue(variables_jump2[variable_name("VOLT", dest, "", scenario)*"_Re"])
        Vd_im = getvalue(variables_jump2[variable_name("VOLT", dest, "", scenario)*"_Im"])

        link_elems_formulations = OPFpb.mp.link_formulations[link]
        link_vars = OPFpb.mp.link_vars[link]

        pt = Point()
        add_coord!(pt, link_vars["Volt_orig"] , Vo_re + im * Vo_im)
        add_coord!(pt, link_vars["Volt_dest"], Vd_re + im * Vd_im)

        for (elemid, element) in elems
            scenario == basecase_scenario_name() ? scenario_id = "0" : scenario_id = String(matchall(r"\d+", scenario)[1])
            link_id = element.id
            orig_id = element.orig_id
            dest_id = element.dest_id
            elem_formulation = link_elems_formulations[elemid]
            Sor = evaluate(Sorig(element, link, elemid, elem_formulation, link_vars),pt)
            Sde = evaluate(Sdest(element, link, elemid, elem_formulation, link_vars),pt)
            Slink_values[(scenario_id, link_id, orig_id, dest_id, link_id)] = (real(Sor),imag(Sor),real(Sde),imag(Sde))
         end
       end

     end

    open(joinpath(outpath,"solution2.txt"),"w") do f
      write(f, "--contingency generator\nconID,genID,busID,unitID,q(MW) \n");
      for ((sc, genID, busID, unitID),q) in Qgen_scen_values
        write(f,"$sc, $genID, $busID, $unitID, $q\n")
      end

      write(f,"--end of contingency generator\n--bus\ncontingency id,bus id,v(pu),theta(deg) \n");
      for ((sc,bus), (V_mod, V_theta)) in volt_values
        write(f, "$sc, $bus, $V_mod, $V_theta\n")
      end

      write(f,"--end of bus\n--Delta\ncontingency id,Delta(MW) \n");
      for (sc, delta) in delta_values
        write(f, "$sc, $delta \n")
      end

      write(f,"--end of Delta\n--line flow\ncontingency id,line id,origin bus id,destination bus id,circuit id,p_origin(MW),q_origin(MVar),p_destination(MW),q_destination(MVar) \n");
      for ((sc, br_id, orig_id, dest_id, br_id), (porig, qorig, pdest, qdest)) in Slink_values
        write(f, "$sc, $br_id, $orig_id, $dest_id, $br_id, $porig, $qorig, $pdest, $qdest \n")
      end
      write(f,"--end of line flow\n")
    end
  end

  function write_solutions(OPFpbs::OPFProblems, pt::Point, outpath::String)
      open(joinpath(outpath,"solution1.txt"),"w") do f
        write(f, "--generation dispatch\nbus id,unit id,pg(MW),qg(MVar) \n");
        for (busname, elems) in OPFpbs[basecase_scenario_name()].ds.bus
          for (elemname,element) in elems
                  if typeof(element) == GOCGenerator
                    bus =  element.busid
                    gen = element.id
                    Pgen = pt[Variable(variable_name("Sgen", busname, elemname, basecase_scenario_name())*"_Re",Real)]
                    Qgen = pt[Variable(variable_name("Sgen", busname, elemname, basecase_scenario_name())*"_Im",Real)]
                    write(f, "$bus, $gen, $Pgen, $Qgen\n")
                  end
           end
         end
        write(f,"--end of generation dispatch\n");
      end

      Qgen_scen_values = Dict{Tuple{String,Any,Int64,Any}, Float64}()
      volt_values = Dict{Tuple{String,Int64}, Tuple{Float64, Float64}}()
      delta_values = Dict{String,Float64}()
      Slink_values = Dict{Tuple{String,Any,Int64,Int64,Any}, Tuple{Float64, Float64, Float64, Float64}}()

      for (scenario, OPFpb) in OPFpbs
        if scenario==basecase_scenario_name()
          scenario_id = "0"
          ##volt values
           for (busname, elems) in OPFpb.ds.bus
             for (elemid, element) in elems
               if typeof(element) == GOCVolt
                 bus = element.busid
                 V_re = pt[Variable(variable_name("VOLT", busname, "", scenario)*"_Re",Real)]
                 V_im = pt[Variable(variable_name("VOLT", busname, "", scenario)*"_Im",Real)]
                 V_mod = abs(V_re + V_im * im)
                 V_theta = angle(V_re + V_im * im)*180/pi
                 volt_values[(scenario_id,bus)] = (V_mod, V_theta)
               end
             end
           end
        else
          ##delta values
          scenario_id = String(matchall(r"\d+", scenario)[1])
          delta_values[scenario_id] = pt[Variable(get_delta_varname(scenario),Real)]
          ##volt values
           for (busname, elems) in OPFpb.ds.bus
             for (elemid, element) in elems
               if typeof(element) == GOCVolt
                 bus = element.busid
                 # bus = String(matchall(r"\d+", element.busname)[1])
                 V_re = pt[Variable(variable_name("VOLT", busname, "", scenario)*"_Re",Real)]
                 V_im = pt[Variable(variable_name("VOLT", busname, "", scenario)*"_Im",Real)]
                 V_mod = abs(V_re + V_im * im)
                 V_theta = angle(V_re + V_im * im)*180/pi
                 volt_values[(scenario_id,bus)] = (V_mod, V_theta)
               elseif typeof(element) == GOCGenerator
                 bus =  element.busid
                 gen = element.id
                 # Pgen = getvalue(variables_jump2[variable_name("Sgen", busname, elemid, scenario)*"_Re"])
                 Qgen = pt[Variable(variable_name("Sgen", busname, elemid, scenario)*"_Im",Real)]
                 # println("bus ",bus, " : ", element.power_min," << ", Pgen + im * Qgen," << ", element.power_max)
                 Qgen_scen_values[(scenario_id,gen,bus, gen)] = Qgen
               end
             end
           end
        end

        for (link, elems) in OPFpb.ds.link
          orig = link.orig
          dest = link.dest
          Vo_re = pt[Variable(variable_name("VOLT", orig, "", scenario)*"_Re",Real)]
          Vo_im = pt[Variable(variable_name("VOLT", orig, "", scenario)*"_Im",Real)]
          Vd_re = pt[Variable(variable_name("VOLT", dest, "", scenario)*"_Re",Real)]
          Vd_im = pt[Variable(variable_name("VOLT", dest, "", scenario)*"_Im",Real)]

          link_elems_formulations = OPFpb.mp.link_formulations[link]
          link_vars = OPFpb.mp.link_vars[link]

          pt2 = Point()
          add_coord!(pt2, link_vars["Volt_orig"] , Vo_re + im * Vo_im)
          add_coord!(pt2, link_vars["Volt_dest"], Vd_re + im * Vd_im)

          for (elemid, element) in elems
              scenario == basecase_scenario_name() ? scenario_id = "0" : scenario_id = String(matchall(r"\d+", scenario)[1])
              link_id = element.id
              orig_id = element.orig_id
              dest_id = element.dest_id
              elem_formulation = link_elems_formulations[elemid]
              Sor = evaluate(Sorig(element, link, elemid, elem_formulation, link_vars),pt2)
              Sde = evaluate(Sdest(element, link, elemid, elem_formulation, link_vars),pt2)
              Slink_values[(scenario_id, link_id, orig_id, dest_id, link_id)] = (real(Sor),imag(Sor),real(Sde),imag(Sde))
           end
         end

       end

      open(joinpath(outpath,"solution2.txt"),"w") do f
        write(f, "--contingency generator\nconID,genID,busID,unitID,q(MW) \n");
        for ((sc, genID, busID, unitID),q) in Qgen_scen_values
          write(f,"$sc, $genID, $busID, $unitID, $q\n")
        end

        write(f,"--end of contingency generator\n--bus\ncontingency id,bus id,v(pu),theta(deg) \n");
        for ((sc,bus), (V_mod, V_theta)) in volt_values
          write(f, "$sc, $bus, $V_mod, $V_theta\n")
        end

        write(f,"--end of bus\n--Delta\ncontingency id,Delta(MW) \n");
        for (sc, delta) in delta_values
          write(f, "$sc, $delta \n")
        end

        write(f,"--end of Delta\n--line flow\ncontingency id,line id,origin bus id,destination bus id,circuit id,p_origin(MW),q_origin(MVar),p_destination(MW),q_destination(MVar) \n");
        for ((sc, br_id, orig_id, dest_id, br_id), (porig, qorig, pdest, qdest)) in Slink_values
          write(f, "$sc, $br_id, $orig_id, $dest_id, $br_id, $porig, $qorig, $pdest, $qdest \n")
        end
        write(f,"--end of line flow\n")
      end
    end


text:
module GreenVideoExperiments

using CSV
using FilePaths
using FilePathsBase: /
using DataFrames

export load

include("load.jl")

end # module


text:
using SDEs, Base.Test
include("setup.jl")
include("testutils.jl")
#include("test_diffeq.jl")

##
x1, W1 = integrate(prob, x0, t, EM(), 1)
x2, W2 = integrate(prob, x0, t, SDEs.PCEproto(), 1; θ=0., η=0.5)
#x3, W3 = integrate(prob, x0, t, DEQEM(), 1)

println((x2 - x1)[end])

##
using PyPlot
plot(t, x1)
plot(t, x2)

##
figure()
plot(t, W1-W2)

##

max_stars_repo_name: onoderat/SDEs

max_issues_repo_name: onoderat/SDEs

max_forks_repo_name: onoderat/SDEs


text:
module ReactiveMPBenchmarkCI

# Just to make BenchmarkCI.jl happy

end


text:
using TikzGraphs
using Base.Test

@assert success(`lualatex -v`)
using NBInclude
nbinclude(joinpath(dirname(@__FILE__), "..", "doc", "TikzGraphs.ipynb"))


text:
module Range

include("./../../constraints/geometric/constants.jl")
using JuMP
using .Constants

export second_dim_range

function second_dim_range(m)

    #Then, given that x-axis is the same, we now want to check whether the boxes are within the y-range

    #Xj < Xi + Xcoveri
    @variable(m, if2_guillo[i=1:length(boxes)-1, j=1:length(boxes)-1; i!=j], Bin)
    @constraint(m, c1_if2_guillo[i=1:length(boxes)-1, j=1:length(boxes)-1; i!=j],
        L * (1 - m[:if2_guillo][i, j]) >=
        m[:flbx_front_left_bot_x][j] + 0.001 - (m[:flbx_front_left_bot_x][i] + m[:box_x_axis_cover][i])
    )

    @constraint(m, c2_if2_guillo[i=1:length(boxes)-1, j=1:length(boxes)-1; i!=j],
        L * m[:if2_guillo][i, j] >=
        (m[:flbx_front_left_bot_x][i] + m[:box_x_axis_cover][i]) - m[:flbx_front_left_bot_x][j] - 0.001
    )

    # Xj + Xcoverj > Xi
    @variable(m, if3_guillo[i=1:length(boxes)-1, j=1:length(boxes)-1; i!=j], Bin)
    @constraint(m, c1_if3_guillo[i=1:length(boxes)-1, j=1:length(boxes)-1; i!=j],
        L * (1 - m[:if3_guillo][i, j]) >=
        m[:flbx_front_left_bot_x][i] + 0.001 - (m[:flbx_front_left_bot_x][j] + m[:box_x_axis_cover][j])
    )

    @constraint(m, c2_if3_guillo[i=1:length(boxes)-1, j=1:length(boxes)-1; i!=j],
        L * m[:if3_guillo][i, j] >=
        (m[:flbx_front_left_bot_x][j] + m[:box_x_axis_cover][j]) - m[:flbx_front_left_bot_x][i] - 0.001
    )

    @variable(m, x_for_y_guillo[i=1:length(boxes)-1, j=1:length(boxes)-1; i!=j], Bin)
    @constraint(m, c_x_for_y_guillo[i=1:length(boxes)-1, j=1:length(boxes)-1; i!=j],
        m[:if2_guillo][i, j] +
        m[:if3_guillo][i, j] >=
        2 * m[:x_for_y_guillo][i, j]
    )

    # Yj < Yi + Ycoveri
    @variable(m, if4_guillo[i=1:length(boxes)-1, j=1:length(boxes)-1; i!=j], Bin)
    @constraint(m, c1_if4_guillo[i=1:length(boxes)-1, j=1:length(boxes)-1; i!=j],
        W * (1 - m[:if4_guillo][i, j]) >=
        m[:flby_front_left_bot_y][j] + 0.001 - (m[:flby_front_left_bot_y][i] + m[:box_y_axis_cover][i])
    )

    @constraint(m, c2_if4_guillo[i=1:length(boxes)-1, j=1:length(boxes)-1; i!=j],
        W * m[:if4_guillo][i, j] >=
        (m[:flby_front_left_bot_y][i] + m[:box_y_axis_cover][i]) - m[:flby_front_left_bot_y][j] - 0.001
    )

    # Yj + Ycoverj > Yi
    @variable(m, if5_guillo[i=1:length(boxes)-1, j=1:length(boxes)-1; i!=j], Bin)
    @constraint(m, c1_if5_guillo[i=1:length(boxes)-1, j=1:length(boxes)-1; i!=j],
        W * (1 - m[:if5_guillo][i, j]) >=
        m[:flby_front_left_bot_y][i] + 0.001 - (m[:flby_front_left_bot_y][j] + m[:box_y_axis_cover][j])
    )

    @constraint(m, c2_if5_guillo[i=1:length(boxes)-1, j=1:length(boxes)-1; i!=j],
        W * m[:if5_guillo][i, j] >=
        (m[:flby_front_left_bot_y][j] + m[:box_y_axis_cover][j]) - m[:flby_front_left_bot_y][i] - 0.001
    )

    @variable(m, y_for_x_guillo[i=1:length(boxes)-1, j=1:length(boxes)-1; i!=j], Bin)
    @constraint(m, c_y_for_x_guillo[i=1:length(boxes)-1, j=1:length(boxes)-1; i!=j],
        m[:if4_guillo][i, j] +
        m[:if5_guillo][i, j] >=
        2 * m[:y_for_x_guillo][i, j]
    )

    return m
end

end


text:
abstract type AbstractTransport end

"""Transport layer over ZMQ sockets"""
struct ZMQTransport <: AbstractTransport
    ctx::Context
    sock::Socket
    mode::Int # REQ/REP
    bound::Bool

    function ZMQTransport(addr::String, mode::Int, bound::Bool, ctx::Context=Context())
        sock = Socket(ctx, mode)
        if bound
            ZMQ.bind(sock, addr)
        else
            ZMQ.connect(sock, addr)
        end
        new(ctx, sock, mode, bound)
    end
    ZMQTransport(ip, port, mode::Int, bound::Bool, ctx::Context=Context()) = ZMQTransport("tcp://$ip:$port", mode, bound, ctx)
end

function sendrecv(conn::ZMQTransport, msgstr)
    @debug("sending request", msgstr)
    ZMQ.send(conn.sock, ZMQ.Message(msgstr))
    respstr = unsafe_string(ZMQ.recv(conn.sock))
    @debug("received response", respstr)
    respstr
end

function sendresp(conn::ZMQTransport, msgstr)
    @debug("sending response", msgstr)
    ZMQ.send(conn.sock, ZMQ.Message(msgstr))
end

function recvreq(conn::ZMQTransport)
    reqstr = unsafe_string(ZMQ.recv(conn.sock))
    @debug("received request", reqstr)
    reqstr
end

function close(conn::ZMQTransport)
    close(conn.sock)
    # close(conn.ctx)
end

"""Transport layer over in-process Channels"""
struct InProcTransport <: AbstractTransport
    name::Symbol

    function InProcTransport(name::Symbol)
        if !(name in keys(InProcChannels))
            InProcChannels[name] = (Channel{Any}(1), Channel{Any}(1))
        end
        new(name)
    end
end

const InProcChannels = Dict{Symbol,Tuple{Channel{Any},Channel{Any}}}()

function sendrecv(conn::InProcTransport, msg)
    clntq,srvrq = InProcChannels[conn.name]

    @debug("sending request", msg)
    put!(srvrq, msg)
    resp = take!(clntq)
    @debug("received response", resp)
    resp
end

function sendresp(conn::InProcTransport, msg)
    clntq,srvrq = InProcChannels[conn.name]
    @debug("sending response", msg)
    put!(clntq, msg)
    nothing
end

function recvreq(conn::InProcTransport)
    clntq,srvrq = InProcChannels[conn.name]
    req = take!(srvrq)
    @debug("received request", req)
    req
end

function close(conn::InProcTransport)
    if conn.name in keys(InProcChannels)
        delete!(InProcChannels, conn.name)
    end
    nothing
end


text:
@test nv(r1) == 10
@test ne(r1) == 20
@test nv(r2) == 5
@test ne(r2) == 10

@test Graph(10,20,seed=3) == Graph(10,20,seed=3)
@test DiGraph(10,20,seed=3) == DiGraph(10,20,seed=3)
@test Graph(10,20,seed=3) == erdos_renyi(10,20,seed=3)
@test ne(Graph(10,40,seed=3)) == 40
@test ne(DiGraph(10,80,seed=3)) == 80

er = erdos_renyi(10, 0.5)
@test nv(er) == 10
@test is_directed(er) == false
er = erdos_renyi(10, 0.5, is_directed=true)
@test nv(er) == 10
@test is_directed(er) == true

er = erdos_renyi(10, 0.5, seed=17)
@test nv(er) == 10
@test is_directed(er) == false

ws = watts_strogatz(10,4,0.2)
@test nv(ws) == 10
@test ne(ws) == 20
@test is_directed(ws) == false

ws = watts_strogatz(10, 4, 0.2, is_directed=true)
@test nv(ws) == 10
@test ne(ws) == 20
@test is_directed(ws) == true

ba = barabasi_albert(10, 2)
@test nv(ba) == 10
@test ne(ba) == 16
@test is_directed(ba) == false

ba = barabasi_albert(10, 2, 2)
@test nv(ba) == 10
@test ne(ba) == 16
@test is_directed(ba) == false

ba = barabasi_albert(10, 4, 2)
@test nv(ba) == 10
@test ne(ba) == 12
@test is_directed(ba) == false

ba = barabasi_albert(10, 2, complete=true)
@test nv(ba) == 10
@test ne(ba) == 17
@test is_directed(ba) == false

ba = barabasi_albert(10, 2, 2, complete=true)
@test nv(ba) == 10
@test ne(ba) == 17
@test is_directed(ba) == false

ba = barabasi_albert(10, 4, 2, complete=true)
@test nv(ba) == 10
@test ne(ba) == 18
@test is_directed(ba) == false

ba = barabasi_albert(10, 2, is_directed=true)
@test nv(ba) == 10
@test ne(ba) == 16
@test is_directed(ba) == true

ba = barabasi_albert(10, 2, 2, is_directed=true)
@test nv(ba) == 10
@test ne(ba) == 16
@test is_directed(ba) == true

ba = barabasi_albert(10, 4, 2, is_directed=true)
@test nv(ba) == 10
@test ne(ba) == 12
@test is_directed(ba) == true

ba = barabasi_albert(10, 2, is_directed=true, complete=true)
@test nv(ba) == 10
@test ne(ba) == 18
@test is_directed(ba) == true

ba = barabasi_albert(10, 2, 2, is_directed=true, complete=true)
@test nv(ba) == 10
@test ne(ba) == 18
@test is_directed(ba) == true

ba = barabasi_albert(10, 4, 2, is_directed=true, complete=true)
@test nv(ba) == 10
@test ne(ba) == 24
@test is_directed(ba) == true

fm = static_fitness_model(20, rand(10))
@test nv(fm) == 10
@test ne(fm) == 20
@test is_directed(fm) == false

fm = static_fitness_model(20, rand(10), rand(10))
@test nv(fm) == 10
@test ne(fm) == 20
@test is_directed(fm) == true

sf = static_scale_free(10, 20, 2.0)
@test nv(sf) == 10
@test ne(sf) == 20
@test is_directed(sf) == false

sf = static_scale_free(10, 20, 2.0, 2.0)
@test nv(sf) == 10
@test ne(sf) == 20
@test is_directed(sf) == true

rr = random_regular_graph(5, 0)
@test nv(rr) == 5
@test ne(rr) == 0
@test is_directed(rr) == false

rd = random_regular_digraph(10,0)
@test nv(rd) == 10
@test ne(rd) == 0
@test is_directed(rd)

rr = random_regular_graph(6, 3, seed=1)
@test nv(rr) == 6
@test ne(rr) == 9
@test is_directed(rr) == false

rr = random_regular_graph(1000, 50)
@test nv(rr) == 1000
@test ne(rr) == 25000
@test is_directed(rr) == false
for v in vertices(rr)
    @test degree(rr, v) == 50
end

rr = random_configuration_model(10, repmat([2,4] ,5), seed=3)
@test nv(rr) == 10
@test ne(rr) == 15
@test is_directed(rr) == false
num2 = 0; num4 = 0
for v in vertices(rr)
    d = degree(rr, v)
    @test  d == 2 || d == 4
    d == 2 ? num2 += 1 : num4 += 1
end
@test num4 == 5
@test num2 == 5

rr = random_configuration_model(1000, zeros(Int,1000))
@test nv(rr) == 1000
@test ne(rr) == 0
@test is_directed(rr) == false

rr = random_configuration_model(3, [2,2,2], check_graphical=true)
@test nv(rr) == 3
@test ne(rr) == 3
@test is_directed(rr) == false

rd = random_regular_digraph(1000, 4)
@test nv(rd) == 1000
@test ne(rd) == 4000
@test is_directed(rd)
@test std(outdegree(rd)) == 0

rd = random_regular_digraph(1000, 4, dir=:in)
@test nv(rd) == 1000
@test ne(rd) == 4000
@test is_directed(rd)
@test std(indegree(rd)) == 0

rr = random_regular_graph(10, 8, seed=4)
@test nv(rr) == 10
@test ne(rr) == 40
@test is_directed(rr) == false
for v in vertices(rr)
    @test degree(rr, v) == 8
end

rd = random_regular_digraph(10, 8, dir=:out, seed=4)
@test nv(rd) == 10
@test ne(rd) == 80
@test is_directed(rd)

g = stochastic_block_model(2., 3., [100,100])
@test  4.5 < mean(degree(g)) < 5.5
g = stochastic_block_model(3., 4., [100,100,100])
@test  10.5 < mean(degree(g)) < 11.5

function generate_nbp_sbm(numedges, sizes)
    density = 1
    # print(STDERR, "Generating communites with sizes: $sizes\n")
    between = density * 0.90
    intra = density * -0.005
    noise = density * 0.00501
    sbm = nearbipartiteSBM(sizes, between, intra, noise)
    edgestream = @task make_edgestream(sbm)
    g = Graph(sum(sizes), numedges, edgestream)
    return sbm, g
end

numedges = 100
sizes = [10, 10, 10, 10]

n = sum(sizes)
sbm, g = generate_nbp_sbm(numedges, sizes)
bc = blockcounts(sbm, g)
bp = blockfractions(sbm, g) ./ (sizes * sizes')
ratios = bp ./ (sbm.affinities ./ sum(sbm.affinities))
@test norm(ratios) < 0.25

sizes = [200, 200, 100]
internaldeg = 15
externaldeg = 6
internalp = Float64[internaldeg/i for i in sizes]
externalp = externaldeg/sum(sizes)
numedges = internaldeg + externaldeg #+ sum(externaldeg.*sizes[2:end])
numedges *= div(sum(sizes), 2)
sbm = StochasticBlockModel(internalp, externalp, sizes)
g = Graph(sum(sizes), numedges, sbm)
@test ne(g) <= numedges
@test nv(g) == sum(sizes)
bc = blockcounts(sbm, g)
bp = blockfractions(sbm, g) ./ (sizes * sizes')
ratios = bp ./ (sbm.affinities ./ sum(sbm.affinities))
@test norm(ratios) < 0.25

# check that average degree is not too high
# factor of two is cushion for random process
@test mean(degree(g)) <= 4//2*numedges/sum(sizes)
# check that the internal degrees are higher than the external degrees
# 5//4 is cushion for random process.
@test all(sum(bc-diagm(diag(bc)), 1) .<= 5//4 .* diag(bc))

sbm2 = StochasticBlockModel(0.5*ones(4), 0.3, 10*ones(Int,4))
sbm  = StochasticBlockModel(0.5, 0.3, 10, 4)
@test sbm == sbm2
sbm.affinities[1,1] = 0
@test sbm != sbm2


text:
module Systems

export ImpulseSystem,ShockTubeSystem,UUtoEverything

using ..Tools
"""
impulsesystem is a struct containing

    ρ = density
    G = mass flow rate
    P = perimeter of the cross section of the wall
    Ac = Area of the cross section of the wall
    qw = heat flux of the wall when t>0
    Δt = time interval
    Δz = space interval
"""

struct ImpulseSystem
    ρ
    G
    P
    Ac
    qw
    Δz
end

"""
shocktubesystem is a struct containing

    gamma = gamma
    Δx = space interval
"""

struct ShockTubeSystem
    gamma
    Δx
end

"""
UUtoEverything is a struct containing all the thermodyamic values and flow velocity

    u = flow velocity
    ρ = density
    m = mass flow rate
    p = pressure
    e = total energy per mass
    ϵ = internal energy per mass
    h = enthalpy per mass

"""

struct UUtoEverything
    u
    ρ
    m
    p
    e
    ϵ
    h
    gamma
    c
end

function UUtoEverything(uu,pipesystem::ShockTubeSystem)

    gamma = pipesystem.gamma

    u = uu[2,:]./uu[1,:]
    ρ = uu[1,:]
    m = u.*ρ
    p = uutourhop(uu,gamma)[3,:]
    e = uu[3,:]./ρ
    ϵ = p./ρ./(gamma-1)
    h = ϵ.+p./ρ
    c = sqrt.(gamma.*p./ρ)

return UUtoEverything(u,ρ,m,p,e,ϵ,h,gamma,c)

end

end


text:
function logprior(
  rparams::RiskParameters{S},
  rpriors::RiskPriors{S}) where {
  S <: DiseaseStateSequence}
  lprior = 0.0
  for i in 1:length(rpriors.sparks)
    lprior == -Inf && break
    lprior += Distributions.logpdf(rpriors.sparks[i], rparams.sparks[i])
  end
  for i in 1:length(rpriors.susceptibility)
    lprior == -Inf && break
    lprior += Distributions.logpdf(rpriors.susceptibility[i], rparams.susceptibility[i])
  end
  for i in 1:length(rpriors.infectivity)
    lprior == -Inf && break
    lprior += Distributions.logpdf(rpriors.infectivity[i], rparams.infectivity[i])
  end
  for i in 1:length(rpriors.transmissibility)
    lprior == -Inf && break
    lprior += Distributions.logpdf(rpriors.transmissibility[i], rparams.transmissibility[i])
  end
  if State_E ∈ S
    for i in 1:length(rpriors.latency)
      lprior == -Inf && break
      lprior += Distributions.logpdf(rpriors.latency[i], rparams.latency[i])
    end
  end
  if State_R ∈ S
    for i in 1:length(rpriors.removal)
      lprior == -Inf && break
      lprior += Distributions.logpdf(rpriors.removal[i], rparams.removal[i])
    end
  end
  return lprior
end

function logprior(
  sm::NucleicAcidSubstitutionModel,
  sm_priors::Vector{UnivariateDistribution})
  lprior = 0.0
  for i in zip([getproperty(sm, θ) for θ in [propertynames(sm)...]], sm_priors)
    lprior += Distributions.logpdf(i[2], i[1])
    lprior == -Inf && break
  end
  return lprior
end


text:
using MPIMeasurements

startLeft = [-1.0 1.0 0.0]
startRight = [1.0 1.0 0.0]
scaleXYZ = [30.0, 34.0, 1.0]
center=[0.0,-6.0, -5.0]

#scaleXYZ = [1.0, 1.0, 1.0]
shiftX = [scaleXYZ[1]/4, 0.0,0.0]

function rotateInvert(pos)
         temp = pos[1,:]
         pos[1,:] = -pos[2,:]
         pos[2,:] = temp
 return pos
end
# Letter T
posTs = (transpose([startLeft;
          1.0 1.0 0.0;
          0.0 1.0 0.0;
          0.0 -1.0 0.0]).*scaleXYZ.+center)
posTsRotIn=rotateInvert(posTs)
aGT = ArbitraryPositions((posTsRotIn)Unitful.mm)

# Letter H
posHs = (transpose([startLeft;
                  -1.0 -1.0 0.0;
                  -1.0 0.0 0.0;
                  1.0 0.0 0.0
                  1.0 1.0 0.0
                  1.0 -1.0 0.0]).*scaleXYZ.+center)
posHsRotIn=rotateInvert(posHs)
aGH = ArbitraryPositions((posHsRotIn)Unitful.mm)

# Letter E
posEs = (transpose([startRight;
                  startLeft;
                  -1.0 0.0 0.0;
                  1.0 0.0 0.0
                  -1.0 0.0 0.0;
                  -1.0 -1.0 0.0
                  1.0 -1.0 0.0]).*scaleXYZ.+center)
posEsRotIn=rotateInvert(posEs)
aGE = ArbitraryPositions((posEsRotIn)Unitful.mm)

# Letter easy U
posUeasy = (transpose([startLeft;
                  -1.0 -1.0 0.0;
                  1.0 -1.0 0.0;
                  1.0 1.0 0.0]).*scaleXYZ.+center)
posUeasyRotIn=rotateInvert(posUeasy)
aGEesay = ArbitraryPositions((posUeasy)Unitful.mm)

# Letter K
posKs = ((transpose([startLeft;
                  -1.0 -1.0 0.0;
                  -1.0 0.0 0.0;
                  1.0 1.0 0.0;
                  -1.0 0.0 0.0;
                  1.0 -1.0 0.0]).*scaleXYZ).+shiftX.+center)
posKsRotIn=rotateInvert(posKs)
aGK = ArbitraryPositions((posKsRotIn)Unitful.mm)

# Letter U
Udown = -0.0
r = 1.0 - Udown
posUs = ((transpose([startLeft;
                  -1.0 Udown 0.0;
                  -r*cos(0.0) -r*sin(0.0)-Udown 0.0;
                  -r*cos(pi/8) -r*sin(pi/8)-Udown 0.0;
                  -r*cos(pi/4) -r*sin(pi/4)-Udown 0.0;
                  -r*cos(pi*3/8) -r*sin(pi*3/8)-Udown 0.0;
                  -r*cos(pi/2) -r*sin(pi/2)-Udown 0.0;
                  -r*cos(pi*5/8) -r*sin(pi*5/8)-Udown 0.0;
                  -r*cos(pi*6/8) -r*sin(pi*6/8)-Udown 0.0;
                  -r*cos(pi*7/8) -r*sin(pi*7/8)-Udown 0.0;
                  1.0 Udown 0.0;
                  startRight]).*scaleXYZ.+center))
posUsRotIn=rotateInvert(posUs)
aGU = ArbitraryPositions((posUsRotIn)Unitful.mm)

@info "Press Enter to continue"
res=readline()

configFile="HeadScanner.toml"
scanner = MPIScanner(configFile)
robot = getRobot(scanner)
defaultVel = getDefaultVelocity(robot)
setup = getSafety(scanner)
bG = ustrip.(parkPos(robot))
posBG = transpose([bG[1] bG[2] bG[3]])
bGA = ArbitraryPositions((posBG)Unitful.mm)

lettersUKE = [aGU,bGA,aGK,bGA,aGE,bGA]
lettersTUHH = [aGT,bGA,aGU,bGA,aGH,bGA,aGH,bGA]

letters=lettersUKE

for letter in letters
for pos in letter
  isValid = checkCoords(setup, pos, getMinMaxPosX(robot))
end
end
Dates.Time(Dates.now())
PosLog=Any[]
TimeLog=Any[]

setEnabled(robot, true)
#sleep(5)
for letter in letters
  for (index,pos) in enumerate(letter)
   TimeLog=push!(TimeLog,Dates.Time(Dates.now()))
   PosLog=cat(1,PosLog,pos')
   moveAbsUnsafe(robot, pos, defaultVel)
  end
end


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

import Base: peek

macro dotimes(n, body)
    quote
        for i = 1:$(esc(n))
            $(esc(body))
        end
    end
end

const whitespace = " \t\r"

"""
Skip any leading whitespace. Returns io.
"""
function skipwhitespace(io::IO; newlines = true)
    while !eof(io) && (Char(peek(io)) in whitespace || (newlines && peek(io) == UInt8('\n')))
        read(io, Char)
    end
    return io
end

"""
Skip any leading blank lines. Returns the number skipped.
"""
function skipblank(io::IO)
    start = position(io)
    i = 0
    while !eof(io)
        c = read(io, Char)
        c == '\n' && (start = position(io); i+=1; continue)
        c == '\r' && (start = position(io); i+=1; continue)
        c in whitespace || break
    end
    seek(io, start)
    return i
end

"""
Returns true if the line contains only (and, unless allowempty,
at least one of) the characters given.
"""
function linecontains(io::IO, chars; allow_whitespace = true,
                                     eat = true,
                                     allowempty = false)
    start = position(io)
    l = readline(io)
    length(l) == 0 && return allowempty

    result = allowempty
    for c in l
        c in whitespace && (allow_whitespace ? continue : (result = false; break))
        c in chars && (result = true; continue)
        result = false; break
    end
    !(result && eat) && seek(io, start)
    return result
end

blankline(io::IO; eat = true) =
    linecontains(io, "",
                 allow_whitespace = true,
                 allowempty = true,
                 eat = eat)

"""
Test if the stream starts with the given string.
`eat` specifies whether to advance on success (true by default).
`padding` specifies whether leading whitespace should be ignored.
"""
function startswith(stream::IO, s::AbstractString; eat = true, padding = false, newlines = true)
    start = position(stream)
    padding && skipwhitespace(stream, newlines = newlines)
    result = true
    for char in s
        !eof(stream) && read(stream, Char) == char ||
            (result = false; break)
    end
    !(result && eat) && seek(stream, start)
    return result
end

function startswith(stream::IO, c::Char; eat = true)
    if !eof(stream) && peek(stream) == UInt8(c)
        eat && read(stream, Char)
        return true
    else
        return false
    end
end

function startswith{T<:AbstractString}(stream::IO, ss::Vector{T}; kws...)
    any(s->startswith(stream, s; kws...), ss)
end

function startswith(stream::IO, r::Regex; eat = true, padding = false)
    @assert Base.startswith(r.pattern, "^")
    start = position(stream)
    padding && skipwhitespace(stream)
    line = readline(stream)
    seek(stream, start)
    m = match(r, line)
    m === nothing && return ""
    eat && @dotimes length(m.match) read(stream, Char)
    return m.match
end

"""
Executes the block of code, and if the return value is `nothing`,
returns the stream to its initial position.
"""
function withstream(f, stream)
    pos = position(stream)
    result = f()
    (result ≡ nothing || result ≡ false) && seek(stream, pos)
    return result
end

"""
Consume the standard allowed markdown indent of
three spaces. Returns false if there are more than
three present.
"""
function eatindent(io::IO, n = 3)
    withstream(io) do
        m = 0
        while startswith(io, ' ') m += 1 end
        return m <= n
    end
end

"""
Read the stream until startswith(stream, delim)
The delimiter is consumed but not included.
Returns nothing and resets the stream if delim is
not found.
"""
function readuntil(stream::IO, delimiter; newlines = false, match = nothing)
    withstream(stream) do
        buffer = IOBuffer()
        count = 0
        while !eof(stream)
            if startswith(stream, delimiter)
                if count == 0
                    return String(take!(buffer))
                else
                    count -= 1
                    write(buffer, delimiter)
                    continue
                end
            end
            char = read(stream, Char)
            char == match && (count += 1)
            !newlines && char == '\n' && break
            write(buffer, char)
        end
    end
end

# TODO: refactor this. If we're going to assume
# the delimiter is a single character + a minimum
# repeat we may as well just pass that into the
# function.

"""
Parse a symmetrical delimiter which wraps words.
i.e. `*word word*` but not `*word * word`.
`repeat` specifies whether the delimiter can be repeated.
Escaped delimiters are not yet supported.
"""
function parse_inline_wrapper(stream::IO, delimiter::AbstractString; rep = false)
    delimiter, nmin = string(delimiter[1]), length(delimiter)
    withstream(stream) do
        if position(stream) >= 1
            # check the previous byte isn't a delimiter
            skip(stream, -1)
            (read(stream, Char) in delimiter) && return nothing
        end
        n = nmin
        startswith(stream, delimiter^n) || return nothing
        while startswith(stream, delimiter); n += 1; end
        !rep && n > nmin && return nothing
        !eof(stream) && Char(peek(stream)) in whitespace && return nothing

        buffer = IOBuffer()
        while !eof(stream)
            char = read(stream, Char)
            write(buffer, char)
            if !(char in whitespace || char == '\n' || char in delimiter) && startswith(stream, delimiter^n)
                trailing = 0
                while startswith(stream, delimiter); trailing += 1; end
                trailing == 0 && return String(take!(buffer))
                write(buffer, delimiter ^ (n + trailing))
            end
        end
    end
end

function showrest(io::IO)
    start = position(io)
    show(readstring(io))
    println()
    seek(io, start)
end


text:
abstract type DiffusionCache end

mutable struct ExplicitDiffusionCache <: DiffusionCache
    u::Array{Float64,2}
    uprev::Array{Float64,2}
    t::Float64
end
# NOTE: Add types for a matrices
mutable struct CNDiffusionCache <: DiffusionCache
    u::Array{Float64,2}
    uprev::Array{Float64,2}
    rhs::Array{Float64,2}
    a1
    a2
    t::Float64
end
mutable struct ADIDiffusionCache <: DiffusionCache
    u::Array{Float64,2}
    uhalf::Array{Float64,2}
    uprev::Array{Float64,2}
    rhs::Array{Float64,2}
    a1x
    a1y
    a2x
    a2y
    t::Float64
end


text:
# Represents a vector register
immutable VecRegister{N, T}
    data::NTuple{N, VecElement{T}}
end

function Base.getindex(vec_reg::VecRegister, i)
    @inbounds v = vec_reg.data[i].value
    return v
end

function vectupexpr(f,N)
    ex = Expr(:tuple, [f(i) for i=1:N]...)
    return quote
        $ex
    end
end

# Elementwise unary functions
for f in UNARY_FUNCS
    @eval begin
        @generated function Base.$(f){N}(a::VecRegister{N})
            expr = vectupexpr(i -> :(VecElement($($f)(a[$i]))), N)
            return quote
                $(Expr(:meta, :inline))
                VecRegister($expr)
            end
        end
    end
end

# Binary func between VecRegister, VecRegister and
# Number, VecRegister.
for f in BINARY_FUNCS
    @eval begin
        @generated function $(f){N}(a::VecRegister{N}, b::VecRegister{N})
            expr = vectupexpr(i -> :(VecElement( ($($f))(a[$i], b[$i]))), N)
            return quote
            $(Expr(:meta, :inline))
            VecRegister($expr)
            end
        end

        @generated function $(f){N}(a::VecRegister{N}, b::Number)
            expr = vectupexpr(i -> :(VecElement( ($($f))(a[$i], b))), N)
            return quote
            $(Expr(:meta, :inline))
            VecRegister($expr)
            end
        end

        @generated function $(f){N}(b::Number, a::VecRegister{N})
            expr = vectupexpr(i -> :(VecElement( ($($f))(b, a[$i]))), N)
            return quote
            $(Expr(:meta, :inline))
            VecRegister($expr)
            end
        end
    end
end

function recursive_reduction_expr!(expr, i, N, f)
    if i == 2
        return expr
    else
        return Expr(:call, f, recursive_reduction_expr!(expr, i-1, N, f), :(a[$(i)]))
    end
end

for (f, f_scal) in REDUCTION_FUNCS
    @eval begin
        @generated function $(f){N}(a::VecRegister{N})
            if N == 0
                return quote end
            elseif N == 1
                return quote a[1] end
            else
                expr = Expr(:call, $f_scal, :(a[1]), :(a[2]))
                full_red_expr = recursive_reduction_expr!(expr, N, N, $f_scal)
                return quote
                    $(Expr(:meta, :inline))
                    $full_red_expr
                end
            end
        end
    end
end

@generated function Base.zero{N,T}(::Type{VecRegister{N,T}})
    result = vectupexpr(i -> :(VecElement(z)), N)
    return quote
        $(Expr(:meta, :inline))
        z = zero($T)
        VecRegister($result)
    end
end

@generated function Base.one{N,T}(::Type{VecRegister{N,T}})
    result = vectupexpr(i -> :(VecElement(z)), N)
    return quote
        $(Expr(:meta, :inline))
        z = one($T)
        return quote VecRegister($result) end
    end
end

@generated function Base.rand{N, T}(::Type{VecRegister{N,T}})
    expr = vectupexpr(i -> :(VecElement(rand(T))), N)
    return quote
        $(Expr(:meta, :inline))
        VecRegister($result)
    end
end


text:
module TestDatasetsSynthetic

using Test
using MLJBase.DatasetsSynthetic

@testset "make_blobs Tests" begin

    @testset "make_blobs Tests: n=110, p=2, center=3" begin

        n, p, centers = 110, 2, 3
        X, y = DatasetsSynthetic.make_blobs(n; p=p, centers=centers)
        @test (n, p) == size(X)
        @test n == length(y)
        @test centers == length(unique(y))
        end;

    @testset "make_blobs Tests n=90, p=3, center=4" begin
        n, p, centers = 90, 3, 4
        X, y = DatasetsSynthetic.make_blobs(n; p=p, centers=centers)
        @test (n, p) == size(X)
        @test n == length(y)
        @test centers == length(unique(y))
    end;
end;

@testset "make_circles Tests" begin
    n = 55
    X, y = DatasetsSynthetic.make_circles(n)
    @test (n, 2) == size(X)
    @test n == length(y)
    @test 2 == length(unique(y))
end;

@testset "make_moons Tests" begin
    n = 55
    X, y = DatasetsSynthetic.make_moons(n)
    @test (n, 2) == size(X)
    @test n == length(y)
    @test 2 == length(unique(y))
end;

end # module
true


text:
# This file is a part of AstroLib.jl. License is MIT "Expat".
# Copyright (C) 2016 Mosè Giordano.

function _ct2lst{T<:AbstractFloat}(long::T, jd::T)
    t0 = jd - J2000
    t  = t0*inv(JULIANYEAR*100)
    # Compute GST in seconds.
    θ = ct2lst_c[1] + (ct2lst_c[2]*t0) + t*t*(ct2lst_c[3] - t*inv(ct2lst_c[4]))
    return cirrange((θ + long)/15.0, 24.0)
end

"""
    ct2lst(longitude, jd) -> local_sidereal_time
    ct2lst(longitude, tz, date) -> local_sidereal_time

### Purpose ###

Convert from Local Civil Time to Local Mean Sidereal Time.

### Arguments ###

The function can be called in two different ways.  The only argument common to
both methods is `longitude`:

* `longitude`: the longitude in degrees (east of Greenwich) of the place for
  which the local sidereal time is desired, scalar.  The Greenwich mean sidereal
  time (GMST) can be found by setting longitude = `0`.

The civil date to be converted to mean sidereal time can be specified either by
providing the Julian days:

* `jd`: this is number of Julian days for the date to be converted.  It can be a
  scalar or an array.

or the time zone and the date:

* `tz`: the time zone of the site in hours, positive East of the
  Greenwich meridian (ahead of GMT).  Use this parameter to easily
  account for Daylight Savings time (e.g. -4=EDT, -5 = EST/CDT), scalar.
  * `date`: this is the local civil time with type `DateTime`.  It can
  be a scalar or an array.

### Output ###

The local sidereal time for the date/time specified in hours.  This is a scalar
or an array of the same length as `jd` or `date`.

### Method ###

The Julian days of the day and time is question is used to determine the number
of days to have passed since 2000-01-01.  This is used in conjunction with the
GST of that date to extrapolate to the current GST; this is then used to get the
LST.  See Astronomical Algorithms by Jean Meeus, p. 84 (Eq. 11-4) for the
constants used.

### Example ###

Find the Greenwich mean sidereal time (GMST) on 2008-07-30 at 15:53 in
Baltimore, Maryland (longitude=-76.72 degrees).  The timezone is EDT or tz=-4

``` julia
lst = ct2lst(-76.72, -4, DateTime(2008, 7, 30, 15, 53))
# => 11.356505172312609
sixty(lst)
# => 3-element Array{Float64,1}:
#     11.0    # Hours
#     21.0    # Minutes
#     23.4186 # Seconds
```

Find the Greenwich mean sidereal time (GMST) on 2015-11-24 at 13:21 in
Heidelberg, Germany (longitude=08° 43' E).  The timezone is CET or tz=1.
Provide `ct2lst` only with the longitude of the place and the number of
Julian days.

``` julia
# Convert longitude to decimals.
longitude=ten(8, 43);
# Get number of Julian days.  Remember to subtract the time zone in
# order to convert local time to UTC.
jd = jdcnv(DateTime(2015, 11, 24, 13, 21) - Dates.Hour(1));
# Calculate Greenwich Mean Sidereal Time.
lst = ct2lst(longitude, jd)
# => 17.140685171005316
sixty(lst)
# => 3-element Array{Float64,1}:
#     17.0    # Hours
#     8.0    # Minutes
#     26.4666 # Seconds
```

### Notes ###

Code of this function is based on IDL Astronomy User's Library.
"""
ct2lst(long::Real, jd::Real) = _ct2lst(promote(float(long), float(jd))...)

function ct2lst{T<:AbstractFloat}(long::T, tz::T,
                                  date::DateTime)
    # In order to handle time zones, package "TimeZones.jl" is much better, bur
    # here we need only to add the time zone to UTC time.  All time zones I know
    # are either integer, or ±30 or ±45 minutes, so it should be safe enough to
    # convert hours to minutes and subtract minutes from "time".
    date = date - Dates.Minute(round(Int64, tz*60))
    return ct2lst(long, jdcnv(date))
end

ct2lst(long::Real, tz::Real, date::DateTime) =
    ct2lst(promote(float(long), float(tz))..., date)

# Vectorial versions
function ct2lst{D<:DateTime}(long::Real, tz::Real, date::AbstractArray{D})
    lst = similar(date, typeof(float(long)))
    for i in eachindex(date)
        lst[i] = ct2lst(long, tz, date[i])
    end
    return lst
end

function ct2lst{JD<:Real}(long::Real, jd::AbstractArray{JD})
    lst = similar(jd, typeof(float(long)))
    for i in eachindex(jd)
        lst[i] = ct2lst(long, jd[i])
    end
    return lst
end


text:
using ModelingToolkit
using DiffEqJump
using Random
using DataFrames
using Tables
using StatsPlots
using BenchmarkTools

@parameters β c γ
@variables t S(t) I(t) R(t)
N = S+I+R
rate₁   = β*c*S*I/N
affect₁ = [S ~ S - 1, I ~ I + 1]
rate₂   = γ*I
affect₂ = [I ~ I - 1, R ~ R + 1]
j₁      = ConstantRateJump(rate₁,affect₁)
j₂      = ConstantRateJump(rate₂,affect₂)
@named sir_js = JumpSystem([j₁,j₂], t, [S,I,R], [β,c,γ])

tmax = 40.0
tspan = (0.0,tmax)
δt = 1.0;

u0 = [S => 990, I => 10, R => 0];

p = [β => 0.05, c => 10.0, γ => 0.25];

Random.seed!(1234);

sir_dprob = DiscreteProblem(sir_js, u0, tspan,p);

sir_jprob = JumpProblem(sir_js, sir_dprob,Direct());

sol_jump = solve(sir_jprob,SSAStepper());

out_jump = sol_jump(0:δt:tmax);

df_jump = DataFrame(Tables.table(out_jump'))
rename!(df_jump,["S","I","R"])
df_jump[!,:t] = out_jump.t;

@df df_jump plot(:t,
    [:S :I :R],
    xlabel="Time",
    ylabel="Number")

@benchmark solve(sir_jprob,SSAStepper())


text:
using Documenter, LaxZonedDateTimes

makedocs(
    modules = [LaxZonedDateTimes],
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
    pages = [
        "Home" => "index.md",
    ],
    repo = "https://github.ca/JuliaTime/LaxZonedDateTimes.jl/blob/{commit}{path}#L{line}",
    sitename = "LaxZonedDateTimes.jl",
    authors = "Curtis Vogt, Gem Newman",
    assets = [
        "assets/invenia.css"
    ],
    checkdocs = :none,
    strict = true,
)

deploydocs(;
    repo="github.com/JuliaTime/LaxZonedDateTimes.jl"
)


text:
include("../../../src/plotSol.jl")
using LinearAlgebra
using DifferentialEquations
using Revise
# using JLD
# clearconsole()
function SprMassMinReal(m,l,k,c,restLen, tInt, tEnd)
    function xdot!(xdot::Vector{Float64},x::Vector{Float64})
        xdot[1] = x[2]
        xdot[2]  = -k/m*(x[1]-restLen) - c/m*x[2];
    end

    function mainDynMinReal!(dx,x,p,t)
        xdot!(dx,x)
    end

    tEnd = 10.0; tSpan = 0.01;
    x0 = [l;1.0]
    prob = ODEProblem(mainDynMinReal!,x0,(0.0,tEnd))
    sol = solve(prob,Tsit5(),saveat = tSpan,reltol=1e-10,abstol=1e-10)
    tSim = sol.t
    xSol = sol[1,:]
    xDotSol = sol[2,:]
    rSol = [xSol+l/2*ones(length(tSim)) zeros(length(tSim),2)]
    vSol = [xDotSol zeros(length(tSim),2)]
    return rSol, vSol
end

max_stars_repo_name: isrlab/Multibody-Dynamics

max_issues_repo_name: isrlab/Multibody-Dynamics

max_forks_repo_name: isrlab/Multibody-Dynamics


text:
## ------------------------------------------------------------------------------
# Container
mutable struct GDModel{T}
    # current state
    iter::Int
    xi::T
    fi::T
    ϵi::T
    ϵii::T
    Δx::T
    sense::T
    sense_count::T
    damp::T

    # params
    damp_win::Real
    damp_factor::Real
    target::T
    maxΔx::T
    minΔx::T
    maxiter::Int
    gdth::Real
    smooth::Real
    verbose::Bool
end

# constructor
function GDModel(;
        target,
        maxΔx, minΔx,
        maxiter = Int(1e4),
        smooth = 1.0,
        it0 = 1,
        damp = zero(target) .+ one(eltype(target)),
        damp_win = 4,
        damp_factor = 0.9, # damp reduction factor
        gdth,
        verbose = true,
    )

    iter = Int(it0)
    maxiter = Int(maxiter)
    xi = zero(target)
    fi = zero(target)
    ϵi = zero(target)
    ϵii = zero(target)
    Δx = zero(target)
    sense = zero(target)
    sense_count = zero(target)
    GDModel{typeof(target)}(
        # state
        iter, xi, fi, ϵi, ϵii, Δx, sense, sense_count, damp,
        # params
        damp_win, damp_factor, target, maxΔx, minΔx, maxiter, gdth, smooth, verbose
    )
end

gd_value(gdmodel::GDModel) = gdmodel.xi

## ------------------------------------------------------------------------------
# Vectorize version
function _default_Err_fun(gdmodel::GDModel{T}) where T<:AbstractArray
    diff = abs.(gdmodel.target - gdmodel.fi)
    norm = ifelse.(iszero.(gdmodel.target), 1.0, abs.(gdmodel.target))
    diff ./ norm
end

_default_break_cond(gdmodel::GDModel{T}) where T<:AbstractArray =
    maximum(gdmodel.ϵi) < gdmodel.gdth

function _default_get_step(gdmodel::GDModel{T}) where T<:AbstractArray
    max_step = gdmodel.maxΔx
    correction = ifelse.(gdmodel.ϵi .< gdmodel.smooth, min.(1.0, gdmodel.ϵi), 1.0)
    step = max_step .* correction .* gdmodel.damp
    gdmodel.sense .* clamp.(step, gdmodel.minΔx, gdmodel.maxΔx)
end

function _default_get_sense(gdmodel::GDModel{T}) where T<:AbstractArray
    Δϵ = gdmodel.ϵi .- gdmodel.ϵii
    zeros = iszero.(Δϵ)
    sense = gdmodel.sense .* -sign.(gdmodel.ϵi .- gdmodel.ϵii)
    if any(zeros)
        rsense = rand(-1:1, length(Δϵ))
        sense = ifelse.(zeros, rsense, sense)
    end
    return sense
end

function grad_desc_vec!(f::Function, gdmodel::GDModel{T};
        x0::T, x1::T,
        toshow::Vector = [],
        Err::Function = _default_Err_fun,
        break_cond::Function = _default_break_cond,
        oniter!::Function = (dgmodel) -> false,
        get_sense::Function = _default_get_sense,
        get_step::Function = _default_get_step,
    ) where {T<:AbstractArray}

    # initializing
    gdmodel.xi .= x0
    gdmodel.fi .= f(gdmodel)
    gdmodel.ϵi .= Err(gdmodel)
    break_cond(gdmodel) && return gdmodel

    gdmodel.ϵii .= gdmodel.ϵi
    gdmodel.sense .= ones(length(gdmodel.target))
    gdmodel.Δx .= zero(gdmodel.xi)
    gdmodel.xi .= x1

    prog = ProgressThresh(gdmodel.gdth, "Grad desc: ")
    it0 = gdmodel.iter
    while true

        gdmodel.xi += gdmodel.Δx
        gdmodel.fi .= f(gdmodel)
        gdmodel.ϵi .= Err(gdmodel)
        gdmodel.sense .= get_sense(gdmodel)
        gdmodel.Δx .= get_step(gdmodel)

        # callback
        ret = oniter!(gdmodel)
        ret && return gdmodel

        gdmodel.ϵii .= gdmodel.ϵi

        # Damp detection
        if iszero(rem(gdmodel.iter, gdmodel.damp_win))
            isdamping = iszero.(gdmodel.sense_count)
            gdmodel.damp .= ifelse.(isdamping, gdmodel.damp .*= gdmodel.damp_factor, gdmodel.damp)
            gdmodel.sense_count .= zero(gdmodel.sense_count)
        end
        gdmodel.sense_count .+= gdmodel.sense

        # verbose
        if gdmodel.verbose
            maxϵᵢ = maximum(gdmodel.ϵi)
            update!(prog, maxϵᵢ; showvalues = vcat(
                    [
                        ("it", gdmodel.iter),
                        ("maxϵᵢ", maxϵᵢ),
                        ("ϵᵢ", gdmodel.ϵi),
                        ("sense", gdmodel.sense),
                        ("xᵢ", gdmodel.xi),
                        ("Δxᵢ", gdmodel.Δx),
                        ("t", gdmodel.target),
                        ("fi", gdmodel.fi),
                    ], toshow
                )
            )
        end

        # break condition
        break_cond(gdmodel) && break
        gdmodel.iter >= gdmodel.maxiter &&  break

        gdmodel.iter += 1
    end
    gdmodel.verbose && finish!(prog)

    return gdmodel
end

function grad_desc_vec(f::Function;
        x0::T, x1::T,
        target::T,
        maxΔx::T = abs.(x0 - x1),
        minΔx::T = abs.(maxΔx .* 1e-35),
        gdth::Real = 1e-5,
        smooth::Real = 1.0,
        maxiter = Int(1000),
        damp = zero(target) .+ one(eltype(target)),
        damp_win = 4,
        damp_factor = 0.9, # damp reduction factor
        it0 = 1,
        verbose = true,
        kwargs...
    ) where {T<:AbstractArray}

    gdmodel = GDModel(;target, maxΔx, minΔx, maxiter, it0,
        gdth, smooth, damp, damp_win, damp_factor, verbose
    )
    grad_desc_vec!(f, gdmodel; x0, x1, kwargs...)

end

## ------------------------------------------------------------------------------
# Scalar version
function _default_Err_fun(gdmodel::GDModel{T}) where T<:Real
    diff = abs(gdmodel.target - gdmodel.fi)
    norm = ifelse(iszero(gdmodel.target), 1.0, abs(gdmodel.target))
    diff / norm
end

_default_break_cond(gdmodel::GDModel{T}) where T<:Real = gdmodel.ϵi < gdmodel.gdth

function _default_get_step(gdmodel::GDModel{T}) where T<:Real
    max_step = gdmodel.maxΔx
    correction = (gdmodel.ϵi < gdmodel.smooth) ? min.(1.0, gdmodel.ϵi) : 1.0
    step = max_step * correction * gdmodel.damp
    gdmodel.sense * clamp(step, gdmodel.minΔx, gdmodel.maxΔx)
end

function _default_get_sense(gdmodel::GDModel{T}) where T<:Real
    Δϵ = gdmodel.ϵi - gdmodel.ϵii
    iszero(Δϵ) ? rand(-1:1) : gdmodel.sense * -sign(gdmodel.ϵi - gdmodel.ϵii)
end

function grad_desc!(f::Function, gdmodel::GDModel{T};
        x0::T, x1::T,
        toshow::Vector = [],
        Err::Function = _default_Err_fun,
        break_cond::Function = _default_break_cond,
        oniter!::Function = (dgmodel) -> false,
        get_sense::Function = _default_get_sense,
        get_step::Function = _default_get_step,
    ) where {T<:Real}

    # initializing
    gdmodel.xi = x0
    gdmodel.fi = f(gdmodel)
    gdmodel.ϵi = Err(gdmodel)
    break_cond(gdmodel) && return gdmodel

    gdmodel.ϵii = gdmodel.ϵi
    gdmodel.sense = one(gdmodel.target)
    gdmodel.Δx = zero(gdmodel.xi)
    gdmodel.xi = x1

    prog = ProgressThresh(gdmodel.gdth, "Grad desc: ")
    it0 = gdmodel.iter
    while true

        gdmodel.xi += gdmodel.Δx
        gdmodel.fi = f(gdmodel)
        gdmodel.ϵi = Err(gdmodel)
        gdmodel.sense = get_sense(gdmodel)
        gdmodel.Δx = get_step(gdmodel)

        # callback
        ret = oniter!(gdmodel)
        ret && return gdmodel

        if gdmodel.verbose
            update!(prog, gdmodel.ϵi; showvalues = vcat(
                    [
                        ("it", gdmodel.iter),
                        ("ϵᵢ", gdmodel.ϵi),
                        ("sense", gdmodel.sense),
                        ("Δx: min/val/max", (gdmodel.minΔx, gdmodel.Δx, gdmodel.maxΔx)),
                        ("xᵢ", gdmodel.xi),
                        ("t", gdmodel.target),
                        ("fi", gdmodel.fi),
                    ], toshow
                )
            )
        end

        gdmodel.ϵii = gdmodel.ϵi

        # Damp detection
        if iszero(rem(gdmodel.iter, gdmodel.damp_win))
            if iszero(gdmodel.sense_count)
                # Damping detetected
                gdmodel.damp *= gdmodel.damp_factor
            end
            gdmodel.sense_count = 0
        end
        gdmodel.sense_count += gdmodel.sense

        # break condition
        break_cond(gdmodel) && break
        gdmodel.iter >= gdmodel.maxiter &&  break

        gdmodel.iter += 1
    end
    gdmodel.verbose && finish!(prog)

    return gdmodel
end

function grad_desc(f::Function;
        x0::T, x1::T,
        target::T,
        maxΔx::T = abs(x0 - x1),
        minΔx::T = abs(maxΔx * 1e-35),
        smooth = 1.0,
        gdth::Real = 1e-5,
        maxiter = Int(1000),
        it0 = 1,
        damp = 1.0,
        damp_win = 4,
        damp_factor = 0.9, # damp reduction factor
        verbose = true,
        kwargs...
    ) where {T<:Real}

    gdmodel = GDModel(;target, minΔx, maxΔx, maxiter, it0, gdth,
        damp, damp_win, damp_factor, smooth, verbose
    )
    grad_desc!(f, gdmodel; x0, x1, kwargs...)

end


text:
function get_eop_history(obj::ITRFProvider)
    return jcall(obj, "getEOPHistory", EOPHistory, ())
end

function get_non_interpolating_provider(obj::ITRFProvider)
    return jcall(obj, "getNonInterpolatingProvider", ITRFProvider, ())
end

function get_transform(obj::ITRFProvider, arg0::AbsoluteDate)
    return jcall(obj, "getTransform", Transform, (AbsoluteDate,), arg0)
end

function get_transform(obj::ITRFProvider, arg0::FieldAbsoluteDate)
    return jcall(obj, "getTransform", FieldTransform, (FieldAbsoluteDate,), arg0)
end


text:
function mergesort_iter!(a)
    k = 1
    while k < length(a)
        for i in 1:2k:length(a)
            if (i + k - 1 < length(a) < i + 2k -1)
                b = merge(a[i:i+k-1], a[i+k : lastindex(a)])
                a[i:lastindex(a)] = b
            elseif ( i + 2k -1 <= length(a))
                b = merge(a[i:i+k-1], a[i+k : i+2k - 1])
                a[i:i+2k-1] = b
            end
        end
        k *= 2
    end

    return a
end

function merge1!(a, b, c)
    d = merge(a,b)
    for i in 1:length(c)
        c[i] = d[i]
    end
    return c
end

function mergesort_iter2!(a)
    k = 1
    b = similar(a)
    while k < length(a)
        for i in 1:2k:length(a)
            if (i + k - 1 < length(a) < i + 2k -1)
                b = merge(a[i:i+k-1], a[i+k : lastindex(a)])
                a[i:lastindex(a)] = b
            elseif ( i + 2k -1 <= length(a))
                b = merge(a[i:i+k-1], a[i+k : i+2k - 1])
                a[i:i+2k-1] = b
            end
        end
        k *= 2
    end

    a = b
    return a
end


text:
#== # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Tests related to the velocity of the Sun.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
#
#   [1] Vallado, D. A (2013). Fundamentals of Astrodynamics and Applications.
#       Microcosm Press, Hawthorn, CA, USA.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # ==#

# File: ./src/sun/velocity.jl
# ===========================

# Function sun_velocity_i
# -----------------------

################################################################################
#                                 Test Results
################################################################################
#
# The time-derivative of the Sun position vector computed by the function
# `sun_velocity_i` is compared to the numerical differentiation of the Sun
# position vector using the function `sun_position_i`, which was already
# validated.
#
################################################################################

@testset "Function sun_velocity_i" begin
    JD_start = DatetoJD(1950, 1, 1, 0, 0, 0)
    JD_stop  = DatetoJD(2019, 1, 1, 0, 0, 0)

    for trial = 1:100
        JD_UT1 = rand(JD_start:JD_stop)
        Δt     = 0.1
        s_t₁   = sun_position_i(JD_UT1)
        s_t₂   = sun_position_i(JD_UT1 + Δt/86400)
        v_n    = (s_t₂ - s_t₁)/Δt
        v      = sun_velocity_i(JD_UT1)

        @test norm(v-v_n)/norm(v)*100 < 0.055
    end
end


text:
mutable struct drake_signal_t <: LCMType
    dim::Int32
    val::Vector{Float64}
    coord::Vector{String}
    timestamp::Int64
end

@lcmtypesetup(drake_signal_t,
    val => (dim,),
    coord => (dim,),
)


text:
using RandomNumbers
using Random123
import Random: seed!
using Test: @test, @test_throws
using Printf: @printf

function compare_dirs(dir1::AbstractString, dir2::AbstractString)
    files1 = readdir(dir1)
    files2 = readdir(dir2)
    @test files1 == files2

    for file in files1
        file1 = joinpath(dir1, file)
        file2 = joinpath(dir2, file)
        lines1 = readlines(file1)
        lines2 = readlines(file2)
        @test lines1 == lines2
    end
end

strip_cr(line::String) = replace(line, r"\r\n$" => "\n")

@info "Testing Random123"
stdout_ = stdout
pwd_ = pwd()
cd(dirname(@__FILE__))
rm("./actual"; force=true, recursive=true)
mkpath("./actual")

for (rng_name, seed_t, stype, seed, args) in (
    (:Threefry2x, NTuple{2, UInt32}, UInt32, (123, 321), (32,)),
    (:Threefry2x, NTuple{2, UInt64}, UInt64, (123, 321), (32,)),
    (:Threefry4x, NTuple{4, UInt32}, UInt32, (123, 321, 456, 654), (72,)),
    (:Threefry4x, NTuple{4, UInt64}, UInt64, (123, 321, 456, 654), (72,)),
    (:Philox2x,   UInt32, UInt32, 123, (16,)),
    (:Philox2x,   UInt64, UInt64, 123, (16,)),
    (:Philox4x,   NTuple{2, UInt32}, UInt32, (123, 321), (16,)),
    (:Philox4x,   NTuple{2, UInt64}, UInt64, (123, 321), (16,))
)
    outfile = open(string(
        "./actual/check-$(string(lowercase("$rng_name"), sizeof(stype)<<3)).out"
    ), "w")
    redirect_stdout(outfile)

    @eval $rng_name($stype)
    x = @eval $rng_name($stype, $seed, $(args...))
    @test seed_type(x) == seed_t
    @test copyto!(copy(x), x) == x

    x.p = 1
    rand(x, UInt64)
    x.p = 1
    rand(x, UInt128)
    @eval rand($x, NTuple{$(string(rng_name)[end-1]-'0'), $stype})

    set_counter!(x, 0)
    for i in 1:100
        @printf "%.9f\n" rand(x)
    end

    close(outfile)
end
redirect_stdout(stdout_)

compare_dirs("expected", "actual")
cd(pwd_)

include("aesni.jl")
include("ars.jl")


text:
module GridWorlds

using Requires

const GW = GridWorlds
export GW

include("directions.jl")
include("actions.jl")
include("objects.jl")
include("grid_world_base.jl")
include("abstract_grid_world.jl")
include("envs/envs.jl")
include("render_in_terminal.jl")

function __init__()
    @require Makie="ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" include("render_with_Makie.jl")
end

end