text:
__precompile__() # this module is safe to precompile

"""
A Julia wrapper for the code controller for Stanford Pupper, a Raspberry Pi-based quadruped robot that can trot, walk, and jump.
"""
module QuadrupedController
export Configuration, Command
export Robot, run!, behavior_state, behavior_state_string
export toggle_activate, toggle_trot, toggle_hop, turn_left, turn_right, end_turn
export set_height!, set_pitch!, set_roll!, adjust_height!, adjust_roll!

@time using Conda
@time using PyCall

# https://github.com/JuliaPy/PyCall.jl#using-pycall-from-julia-modules
const Configuration = PyNULL()
const four_legs_inverse_kinematics = PyNULL()
const Controller = PyNULL()
const State = PyNULL()

function __init__()
    # Install required packages, if not yet installed
    @static if !Sys.iswindows()
        Conda.add_channel("conda-forge")
    end
    if !haskey(Conda._installed_packages_dict(), "numpy")
        Conda.add("numpy")
    end
    if !haskey(Conda._installed_packages_dict(), "transforms3d")
        Conda.add("transforms3d")
    end

    pushfirst!(PyVector(pyimport("sys")."path"), @__DIR__)    # in order to load a Python module from the current directory

    copy!(Configuration, pyimport("Config").Configuration)
    copy!(four_legs_inverse_kinematics, pyimport("Kinematics").four_legs_inverse_kinematics)
    copy!(Controller, pyimport("Controller").Controller)
    copy!(State, pyimport("State").State)
end

"""
    behavior_state(state)

Returns the behavior_state as an integer

# Examples
```julia-repl
julia> behavior_state(state)
1
```
"""
behavior_state(state) = parse(Int, string(state.behavior_state)[end-2:end-1])

"""
    behavior_state_string(state)

Returns the behavior_state as a string

# Examples
```julia-repl
julia> behavior_state_string(state)
TROT
```
"""
behavior_state_string(state) = string(state.behavior_state)[25:end-1]

mutable struct Command # Stores movement commands
    horizontal_velocity::Array{Float64,1}
    yaw_rate::Float64
    height::Float64
    pitch::Float64
    roll::Float64

    command_interval::Float64

    hop_event::Bool
    trot_event::Bool
    activate_event::Bool

    Command(
        horizontal_velocity::Array{Float64,1} = [0.0, 0.0],
        yaw_rate::Float64 = 0.0,
        height::Float64 = -0.16,
        pitch::Float64 = 0.0,
        roll::Float64 = 0.0,

        hop_event::Bool = false,
        trot_event::Bool = false,
        activate_event::Bool = false;

        command_interval::Float64 = 1.0 / 30
    ) = new(horizontal_velocity, yaw_rate, height, pitch, roll, command_interval, hop_event, trot_event, activate_event)
end

struct Robot
    controller              # Defined in Controller.py
    state                   # Defined in State.py
    command::Command
end

"""
    Robot(config = Configuration(), command::Command = Command())

Robot object contructor
"""
function Robot(config = Configuration(), command::Command = Command())
    controller = Controller(config, four_legs_inverse_kinematics,)
    state = State()   # holds the runtime state of the simulator
    Robot(controller, state, command)
end

# associated functions:
# [sim]step - run
# joystick commands
# reading sensors/environment/state

"""
    run!(robot)

Steps the controller
"""
run!(robot) = robot.controller.run(robot.state, robot.command)

"""
    toggle_activate(robot)

Toggles the robot state from DEACTIVATED to REST, and from any state other than DEACTIVATED to DEACTIVATED

# Examples
```julia-repl
julia> toggle_activate(robot)
Toggling from DEACTIVATED: -1 to REST: 0

julia> toggle_activate(robot)
Toggling from REST: 0 to DEACTIVATED: -1
```
"""
function toggle_activate(robot)
    print("Toggling from ", behavior_state_string(robot.state))
    try
        robot.command.activate_event = 1
    catch
        println("\nError: Illegal state transition in toggle_activate!")
        return
    end
    run!(robot)
    robot.command.activate_event = 0
    println(" to ", behavior_state_string(robot.state))
end

"""
    toggle_trot(robot)

Toggles the robot state between REST or a HOP state to TROT, and from TROT to REST

# Examples
```julia-repl
julia> toggle_activate(robot)
Toggling from DEACTIVATED: -1 to REST: 0

julia> toggle_trot(robot)
Toggling from REST: 0 to TROT: 1
```
"""
function toggle_trot(robot)
    print("Toggling from ", behavior_state_string(robot.state))
    try
        robot.state.behavior_state = robot.controller.trot_transition_mapping[robot.state.behavior_state]
    catch
        println("\nError: Illegal state transition in toggle_trot!")
        return
    end
    run!(robot)
    println(" to ", behavior_state_string(robot.state))
end

"""
    toggle_hop(robot)

Toggles the robot state between REST or TROT to HOP, from HOP to FINISHHOP, and from FINISHHOP to REST

# Examples
```julia-repl
julia> toggle_activate(robot)
Toggling from DEACTIVATED: -1 to REST: 0

julia> toggle_hop(robot)
Toggling from REST: 0 to HOP: 2
```
"""
function toggle_hop(robot)
    print("Toggling from ", behavior_state_string(robot.state))
    try
        robot.state.behavior_state = robot.controller.hop_transition_mapping[robot.state.behavior_state]
    catch
        println("\nError: Illegal state transition in toggle_hop!")
        return
    end
    run!(robot)
    println(" to ", behavior_state_string(robot.state))
end

const turn_yaw_rate = 0.4
const turn_velocity = [0.001, 0]

"""
    turn_left(robot)

Turns the robot roughly 90 degrees to the left by slowing down horizontal velocity
to near 0 and setting a yaw_rate to a positive constant

# Examples
```julia-repl
julia> turn_left(robot)
```
"""
function turn_left(robot)
    robot.command.horizontal_velocity = turn_velocity
    robot.command.yaw_rate = turn_yaw_rate
    #println("Pupper marching in place: ", robot.command.horizontal_velocity)
    #println("Pupper turning in place: ",robot.command.yaw_rate)
end

"""
    turn_right(robot)

Turns the robot roughly 90 degrees to the right by slowing down horizontal velocity
to near 0 and setting a yaw_rate to a negative constant

# Examples
```julia-repl
julia> turn_right(robot)
Pupper marching in place:
```
"""
function turn_right(robot)
    robot.command.horizontal_velocity = turn_velocity
    robot.command.yaw_rate = -turn_yaw_rate
    #println("Pupper marching in place: ", robot.command.horizontal_velocity)
    #println("Pupper turning in place: ",robot.command.yaw_rate)
end

"""
    end_turn(robot)

Ends the turn and proceeds on a forward heading at the pre-turn velocity

# Examples
```julia-repl
julia> end_turn(robot)
```
"""
function end_turn(robot)
    robot.command.horizontal_velocity = [0.2, 0]
    robot.command.yaw_rate = 0.0
end

# Functions ported over from JoystickInterface.py
function deadband(value, band_radius)
    return max(value - band_radius, 0) + min(value + band_radius, 0)
end

function clipped_first_order_filter(input, target, max_rate, tau)
    rate = (target - input) / tau
    return clamp(rate, -max_rate, max_rate)
end

# set commands
set_height!(command::Command, height::Real) = (command.height = clamp(height, -0.25, -0.025))

function set_pitch!(robot::Robot, pitch)
    config  = robot.controller.config

    pitch *= -config.max_pitch
    deadbanded_pitch = deadband(pitch, config.pitch_deadband)
    pitch_rate = clipped_first_order_filter(
        robot.state.pitch,
        deadbanded_pitch,
        config.max_pitch_rate,
        config.pitch_time_constant
    )
    robot.command.pitch = robot.state.pitch + robot.command.command_interval * pitch_rate
end

set_roll!(command::Command, roll::Real) = (command.roll = clamp(roll, -1.0, 1.0))

# adjust commands
function adjust_height!(robot::Robot, adjust_by::Int)
    height = robot.state.height - robot.command.command_interval * robot.controller.config.z_speed * adjust_by
    set_height!(robot.command, height)
end

function adjust_roll!(robot::Robot, adjust_by::Int)
    roll = robot.state.roll + robot.command.command_interval * -robot.controller.config.roll_speed * adjust_by
    set_roll!(robot.command, roll)
end

end


text:
# Estimate cost of first-stage resource allocation decision given scenarios of the demands

function estimateCostOfSoln(z_soln::Array{Float64},demand_scen::Array{Float64,2})

	mod = Model(solver=GurobiSolver(gurobi_env,Presolve=0,OutputFlag=0))

	@variable(mod, v[1:numResources,1:numCustomers] >= 0)
	@variable(mod, w[1:numCustomers] >= 0)

	@objective(mod, Min, sum(recourseCostCoeff[j]*w[j] for j=1:numCustomers))

	@constraint(mod, [i=1:numResources], sum(v[i,j] for j = 1:numCustomers) <= rho[i]*z_soln[i])
	@constraint(mod, con1[j=1:numCustomers], sum(mu[i,j]*v[i,j] for i = 1:numResources) + w[j] >= 0)

	numScenarios::Int64 = size(demand_scen,1)
	solnCost::Float64 = 0.0

	for s = 1:numScenarios
		for j = 1:numCustomers
			JuMP.setRHS(con1[j],demand_scen[s,j])
		end

		status = solve(mod)
		solnCost += getobjectivevalue(mod)
	end

	solnCost /= numScenarios
	solnCost += sum(costVector[i]*z_soln[i] for i=1:numResources)

	return solnCost
end

ext: jl

max_stars_repo_name: rohitkannan/DD-SAA

max_issues_repo_name: rohitkannan/DD-SAA

max_forks_repo_name: rohitkannan/DD-SAA


text:
# ---
# title: 516. Longest Palindromic Subsequence
# id: problem516
# author: Indigo
# date: 2021-06-22
# difficulty: Medium
# categories: Dynamic Programming
# link: <https://leetcode.com/problems/longest-palindromic-subsequence/description/>
# hidden: true
# ---
#
# Given a string s, find the longest palindromic subsequence's length in s. You
# may assume that the maximum length of s is 1000.
#
# **Example 1:**
# Input:
#
#
#
#     "bbbab"
#
#
# Output:
#
#
#
#     4
#
#
# One possible longest palindromic subsequence is "bbbb".
#
#
#
# **Example 2:**
# Input:
#
#
#
#     "cbbd"
#
#
# Output:
#
#
#
#     2
#
#
# One possible longest palindromic subsequence is "bb".
#
#
#
# **Constraints:**
#
#   * `1 <= s.length <= 1000`
#   * `s` consists only of lowercase English letters.
#
#
## @lc code=start
using LeetCode

function longest_palindrome_subseq(s::String)
    len = length(s)
    dp = fill(0, len, len)
    for i in len:-1:1
        dp[i, i] = 1
        for j in i+1:len
            dp[i, j] = (s[i] == s[j]) ? dp[i + 1, j - 1] + 2 : max(dp[i + 1, j], dp[i, j - 1])
        end
    end
    return dp[1, len]
end
## @lc code=end


text:
# Autogenerated wrapper script for bliss_jll for aarch64-linux-gnu
export bliss, libbliss

using GMP_jll
JLLWrappers.@generate_wrapper_header("bliss")
JLLWrappers.@declare_library_product(libbliss, "libbliss.so")
JLLWrappers.@declare_executable_product(bliss)
function __init__()
    JLLWrappers.@generate_init_header(GMP_jll)
    JLLWrappers.@init_library_product(
        libbliss,
        "lib/libbliss.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        bliss,
        "bin/bliss",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# BFS algorithm using Julia stdlib primitives

using SparseArrays

"""
    bfs_primitive(A::SparseMatrixCSC{Bool, U}, source::U) where {U <: Integer}

Compute BFS on the graph represented by `A` starting from `source` using Julia primitives.

# Examples
```julia-repl
julia> dim = 100;

julia> A = sparse_square_matrix(dim, Bool);

julia> source = 1;

julia> bfs_primitive(A,source);
```
"""
function bfs_primitive(
    A::SparseMatrixCSC{Bool, U},
    source::U
    ) where {U <: Integer}

    source > 0 || throw(ArgumentError("source must be positive"))

    dim = size(A, 2)
    # Input matrix must be square
    size(A, 1) == dim || throw(DimensionMismatch())

    q = BitArray(undef, dim)
    q[source] = true

    # Create v, make it dense
    v = zeros(Int64, dim)
    v[source] = 1

    aux = similar(q)
    toupdate = spzeros(Bool, dim)

    @inbounds for level in 1:dim

        aux = broadcast(!iszero, A'*q)
        copy!(q, aux)

        # toupdate[i] = true iff vertex i is
        # in the frontier and was not visited before
        @inbounds for i in eachindex(toupdate)
            toupdate[i] = iszero(v[i]) & aux[i]
        end

        v .+= toupdate .* (level+1)

        if(!reduce(max,q))
            break
        end
    end

    return v
end


text:
# Autogenerated wrapper script for PAPI_jll for armv7l-linux-musleabihf
export libpapi

JLLWrappers.@generate_wrapper_header("PAPI")
JLLWrappers.@declare_library_product(libpapi, "libpapi.so.6.0")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libpapi,
        "lib/libpapi.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using OrdinaryDiffEq
using Trixi

###############################################################################
# semidiscretization of the compressible Euler equations
gamma = 1.4
equations = CompressibleEulerEquations2D(gamma)

"""
    initial_condition_kelvin_helmholtz_instability(x, t, equations::CompressibleEulerEquations2D)

A version of the classical Kelvin-Helmholtz instability based on
- Andrés M. Rueda-Ramírez, Gregor J. Gassner (2021)
  A Subcell Finite Volume Positivity-Preserving Limiter for DGSEM Discretizations
  of the Euler Equations
  [arXiv: 2102.06017](https://arxiv.org/abs/2102.06017)
"""
function initial_condition_kelvin_helmholtz_instability(x, t, equations::CompressibleEulerEquations2D)
  # change discontinuity to tanh
  # typical resolution 128^2, 256^2
  # domain size is [-1,+1]^2
  slope = 15
  amplitude = 0.02
  B = tanh(slope * x[2] + 7.5) - tanh(slope * x[2] - 7.5)
  rho = 0.5 + 0.75 * B
  v1 = 0.5 * (B - 1)
  v2 = 0.1 * sin(2 * pi * x[1])
  p = 1.0
  return prim2cons(SVector(rho, v1, v2, p), equations)
end
initial_condition = initial_condition_kelvin_helmholtz_instability

surface_flux = flux_lax_friedrichs
volume_flux  = flux_ranocha
polydeg = 3
basis = LobattoLegendreBasis(polydeg)
indicator_sc = IndicatorHennemannGassner(equations, basis,
                                         alpha_max=0.002,
                                         alpha_min=0.0001,
                                         alpha_smooth=true,
                                         variable=density_pressure)
volume_integral = VolumeIntegralShockCapturingHG(indicator_sc;
                                                 volume_flux_dg=volume_flux,
                                                 volume_flux_fv=surface_flux)
solver = DGSEM(basis, surface_flux, volume_integral)

coordinates_min = (-1.0, -1.0)
coordinates_max = ( 1.0,  1.0)
mesh = TreeMesh(coordinates_min, coordinates_max,
                initial_refinement_level=5,
                n_cells_max=100_000)
semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)

###############################################################################
# ODE solvers, callbacks etc.

tspan = (0.0, 3.0)
ode = semidiscretize(semi, tspan)

summary_callback = SummaryCallback()

analysis_interval = 100
analysis_callback = AnalysisCallback(semi, interval=analysis_interval)

alive_callback = AliveCallback(analysis_interval=analysis_interval)

save_solution = SaveSolutionCallback(interval=20,
                                     save_initial_solution=true,
                                     save_final_solution=true,
                                     solution_variables=cons2prim)

stepsize_callback = StepsizeCallback(cfl=1.3)

callbacks = CallbackSet(summary_callback,
                        analysis_callback, alive_callback,
                        save_solution,
                        stepsize_callback)

###############################################################################
# run the simulation

sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
            save_everystep=false, callback=callbacks);
summary_callback() # print the timer summary


text:
# Test pointer() functions
let str = to_ascii("this ")
    u8  = utf8(str)
    u16 = utf16(str)
    u32 = utf32(str)
    pa  = pointer(str)
    p8  = pointer(u8)
    p16 = pointer(u16)
    p32 = pointer(u32)
    @test typeof(pa) == Ptr{UInt8}
    @test unsafe_load(pa,1) == 0x74
    @test typeof(p8) == Ptr{UInt8}
    @test unsafe_load(p8,1) == 0x74
    @test typeof(p16) == Ptr{UInt16}
    @test unsafe_load(p16,1) == 0x74
    @test typeof(p32) == Ptr{UInt32}
    @test unsafe_load(p32,1) == 0x74
    pa  = pointer(str, 2)
    p8  = pointer(u8,  2)
    p16 = pointer(u16, 2)
    p32 = pointer(u32, 2)
    @test typeof(pa) == Ptr{UInt8}
    @test unsafe_load(pa,1) == 0x68
    @test typeof(p8) == Ptr{UInt8}
    @test unsafe_load(p8,1) == 0x68
    @test typeof(p16) == Ptr{UInt16}
    @test unsafe_load(p16,1) == 0x68
    @test typeof(p32) == Ptr{UInt32}
    @test unsafe_load(p32,1) == 0x68
    sa  = SubString{ASCIIStr}(str, 3, 5)
    s8  = SubString{UTF8Str}(u8,   3, 5)
    s16 = SubString{UTF16Str}(u16, 3, 5)
    s32 = SubString{UTF32Str}(u32, 3, 5)
    pa  = pointer(sa)
    p8  = pointer(s8)
    p16 = pointer(s16)
    p32 = pointer(s32)
    @test typeof(pa) == Ptr{UInt8}
    @test unsafe_load(pa,1) == 0x69
    @test typeof(p8) == Ptr{UInt8}
    @test unsafe_load(p8,1) == 0x69
    @test typeof(p16) == Ptr{UInt16}
    @test unsafe_load(p16,1) == 0x69
    @test typeof(p32) == Ptr{UInt32}
    @test unsafe_load(p32,1) == 0x69
    pa  = pointer(sa, 2)
    p8  = pointer(s8,  2)
    p16 = pointer(s16, 2)
    p32 = pointer(s32, 2)
    @test typeof(pa) == Ptr{UInt8}
    @test unsafe_load(pa,1) == 0x73
    @test typeof(p8) == Ptr{UInt8}
    @test unsafe_load(p8,1) == 0x73
    @test typeof(p16) == Ptr{UInt16}
    @test unsafe_load(p16,1) == 0x73
    @test typeof(p32) == Ptr{UInt32}
    @test unsafe_load(p32,1) == 0x73
end


text:
struct Constant{T}
    value::T
end

make_zero_derivs(x, ::Val{0}) = ()
make_zero_derivs(x, ::Val{N}) where {N} = (deriv = zero_tangent(x); (deriv, make_zero_derivs(deriv, Val(N - 1))...))

function (traj::Constant{T})(x, ::Val{num_derivs}) where {T, num_derivs}
    traj.value, make_zero_derivs(traj.value, Val(num_derivs))...
end

# Single-argument call overload just returns the value:
(traj::Constant)(x) = first(traj(x, Val(0)))


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule SoXResampler_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("SoXResampler")
JLLWrappers.@generate_main_file("SoXResampler", UUID("fbe68eb6-6641-54c6-99e3-f7c7c4d73a57"))
end  # module SoXResampler_jll


text:
const REWARD_TOLERANCE = 1e-6
const CONSTR_TOLERANCE = 1e-3

_approx_eq(x, y) = abs(x - y) < CONSTR_TOLERANCE
_approx_geq_zero(x) = x > -CONSTR_TOLERANCE
_strictly_geq_zero(x) = x > CONSTR_TOLERANCE
_approx_between(x, l, u) = (x > l - CONSTR_TOLERANCE) && (x < u + CONSTR_TOLERANCE)

function _verify_solution(data::DataSet,
                opt_model::OptModel{Float64},
                only_warn::Bool=true)
    linear_model = opt_model.linear_model
    y_values = opt_model.ys
    z_values = opt_model.zs
    trivially_suboptimal = false
    for (i, sample) in enumerate(data.samples)
        bid_price = sample.bid_price
        pay_price = sample.pay_price
        y_value = y_values[i]
        z_value = z_values[i,:]
        linear_value = dot(linear_model.weights, sample.features) + linear_model.offset
        function _error(description)
            error("Expected: $description.\n Actual: linear_value=$linear_value, y_value=$y_value, pay_price=$pay_price, bid_price=$bid_price, and z_value=$z_value")
        end
        function _warn(description)
            @warn "Expected: $description.\n Actual: linear_value=$linear_value, y_value=$y_value, pay_price=$pay_price, bid_price=$bid_price, and z_value=$z_value"
        end
        reporting_fn = only_warn ? _warn : _error
        if !_approx_geq_zero(y_value - min(0.0, bid_price))
            reporting_fn("y ≥ min(0, bid_price)")
        end
        if _approx_eq(z_value[1], 1.0)
            if !_approx_eq(y_value, pay_price)
                reporting_fn("z[1] = 1 and y = pay_price")
            end
            if linear_value > pay_price + CONSTR_TOLERANCE
                reporting_fn("z[1] = 1 and <β,x> ≤ pay_price")
            end
        elseif _approx_eq(z_value[2], 1.0)
            if !_approx_geq_zero(linear_value - y_value)
                reporting_fn("z[2] = 1 and y ≤ <β,x>")
            end
            if !_approx_between(y_value, pay_price, bid_price)
                reporting_fn("z[2] = 1 and y ∈ [pay_price,bid_price]")
            end
        elseif _approx_eq(z_value[3], 1.0)
            if !_approx_eq(y_value, 0.0)
                reporting_fn("z[3] = 1 and y = 0")
            end
            if !_approx_geq_zero(linear_value - bid_price)
                reporting_fn("z[3] = 1 and <β,x> > bid_price")
            end
        else
            # _error("sum(z) = 1")
        end
    end
end

function _audit_weights(learner::Learner, β_values::Vector{Float64})
    num_boundary = 0
    for i in 1:length(β_values)
        if _approx_eq(β_values[i], learner.param_lb[i]) || _approx_eq(β_values[i], learner.param_ub[i])
            num_boundary += 1
        end
    end
    if num_boundary > 0
        @warn "Solution attained on boundary ($num_boundary of $(length(β_values)) components). Enlarge feasible region and continue."
    end
    return nothing
end

function _unif_random(L::Vector{Float64}, U::Vector{Float64})
    n = length(L)
    @assert n == length(U)
    x = Vector{Float64}(undef, n)
    for i in 1:n
        l = L[i]
        u = U[i]
        if l == u
            x[i] = l
        else
            @assert l < u
            x[i] = rand(Uniform(l, u))
        end
    end
    return x
end

function _reward_function(sample::DataPoint, model::LinearModel)
    linear_value = evaluate(model, sample.features)
    if linear_value <= sample.pay_price
        return sample.pay_price
    elseif linear_value <= sample.bid_price + REWARD_TOLERANCE
        return linear_value
    else
        return 0.0
    end
end

function _grad_reward_function(sample::DataPoint, model::LinearModel)
    num_features = length(sample.features)
    linear_value = evaluate(model, sample.features)
    if linear_value <= sample.pay_price
        return zeros(Float64, num_features), 0.0
    elseif linear_value <= sample.bid_price + REWARD_TOLERANCE
        return vcat(sample.features), 1.0
    else
        return zeros(Float64, num_features), 0.0
    end
end

function _compute_reward(data::DataSet, model::LinearModel)
    val = 0.0
    for sample in data.samples
        val += _reward_function(sample, model)
    end
    return val
end

function _compute_expected_reward(data::DataSet, model::LinearModel)
    return _compute_reward(data, model) / sample_count(data)
end

_compute_expected_reward_if_test_data(data::Nothing, model::LinearModel) = NaN
function _compute_expected_reward_if_test_data(data::DataSet, model::LinearModel)
    return _compute_expected_reward(data, model)
end

function _compute_customer_rewards(data::DataSet, model::LinearModel)
    return [_reward_function(sample, model) for sample in data.samples]
end

function _proportion_sold(data::DataSet, model::LinearModel)
    y_values = _compute_customer_rewards(data, model)
    num_sold = count(_strictly_geq_zero, y_values)
    return num_sold / sample_count(data)
end

_proportion_sold_if_test_data(data::Nothing, model::LinearModel) = NaN
function _proportion_sold_if_test_data(data::DataSet, model::LinearModel)
    return _proportion_sold(data, model)
end

import Distributions

import LinearAlgebra.I

function normal_log_normal_data(num_features, num_samples; ρ=0.9, σ=0.1, α=0.0)

    c = (1 / sqrt(num_features)) * rand(Distributions.Normal(), num_features)
    uncorrelated_c = (1 / sqrt(num_features)) * rand(Distributions.Normal(), num_features)
    correlated_c = ρ * c + sqrt(1 - ρ^2) * uncorrelated_c
    X = [1 / sqrt(num_features) * rand(Distributions.Normal(), num_features) for i in 1:num_samples]

    datapoints = Vector{DataPoint}(undef, num_samples)
    for i in 1:num_samples
        c_dot_x = LinearAlgebra.dot(c, X[i])
        corr_c_dot_x = LinearAlgebra.dot(correlated_c, X[i])
        price_one = rand(Distributions.LogNormal(c_dot_x, σ * abs(c_dot_x)))
        price_two = rand(Distributions.LogNormal(corr_c_dot_x, σ * abs(corr_c_dot_x)))
        datapoints[i] = DataPoint(
            X[i],
            (1 + α) * max(price_one, price_two),
            (1 - α) * min(price_one, price_two)
        )
    end
    mean_bid_price = StatsBase.mean(dp.bid_price for dp in datapoints)
    for i in 1:num_samples
        datapoints[i].bid_price /= mean_bid_price
        datapoints[i].pay_price /= mean_bid_price
    end

    return DataSet(datapoints)
end

max_stars_repo_name: joehuchette/reserve-price-optimization

max_issues_repo_name: joehuchette/reserve-price-optimization

max_forks_repo_name: joehuchette/reserve-price-optimization


text:
using Documenter
using UdpTest

makedocs(
    sitename = "UdpTest",
    modules = [UdpTest]
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
deploydocs(
    repo = "github.com/piebat/UdpTest.jl.git",
    branch = "gh-pages"
)


text:
# module for representing agent
# decide control order to robot
# get action based on policy
# by dynamic programming

mutable struct DpPolicyAgent
  speed
  yaw_rate
  delta_time
  estimator
  prev_spd
  prev_yr
  puddle_coef
  puddle_depth
  total_reward
  in_goal
  final_value
  goal
  pose_min
  pose_max
  widths
  index_nums
  policy_data

  # init
  function DpPolicyAgent(;delta_time::Float64=0.1,
                         estimator=nothing,
                         goal=nothing,
                         puddle_coef=100,
                         widths=[0.2, 0.2, pi/18],
                         lower_left=[-4.0, -4.0],
                         upper_right=[4.0, 4.0])
    self = new()
    self.speed = 0.0
    self.yaw_rate = 0.0
    self.delta_time = delta_time
    self.estimator = estimator
    self.prev_spd = 0.0
    self.prev_yr = 0.0
    self.puddle_coef = puddle_coef
    self.puddle_depth = 0.0
    self.total_reward = 0.0
    self.in_goal = false
    self.final_value = 0.0
    self.goal = goal

    self.pose_min = [lower_left[1], lower_left[2], 0.0]
    self.pose_max = [upper_right[1], upper_right[2], 2*pi]
    self.widths = widths
    self.index_nums = round.(Int64, (self.pose_max - self.pose_min)./self.widths)
    self.policy_data = init_policy(self)
    return self
  end
end

function init_policy(self::DpPolicyAgent)
  tmp = zeros(Tuple([self.index_nums[1], self.index_nums[2], self.index_nums[3], 2]))

  txt_path = "src/decision_making/markov_decision_process/policy.txt"
  open(joinpath(split(@__FILE__, "src")[1], txt_path), "r") do fp
    for line in eachline(fp)
      d = split(line) # [i_x, i_y, i_theta, speed, yaw_rate]
      tmp[parse(Int64, d[1])+1, parse(Int64, d[2])+1, parse(Int64, d[3])+1, :] .= [parse(Float64, d[4]), parse(Float64, d[5])]
    end
  end

  return tmp
end

function reward_per_sec(self::DpPolicyAgent)
  return (-1.0 - self.puddle_depth*self.puddle_coef)
end

function to_index(self::DpPolicyAgent, pose)
  index = Int64.(floor.((pose - self.pose_min)./self.widths))

  # normalize direction index
  index[3] = (index[3] + self.index_nums[3]*1000)%self.index_nums[3]

  for i in 1:2
    if index[i] < 0
      index[i] = 0
    elseif index[i] >= self.index_nums[i]
      index[i] = self.index_nums[i] - 1
    end
  end

  return index
end

function policy(self::DpPolicyAgent, pose)
  index = to_index(self, pose)

  action = self.policy_data[index[1]+1, index[2]+1, index[3]+1, :]

  return action[1], action[2]
end

function draw_decision!(self::DpPolicyAgent, observation)
  # reached goal
  if self.in_goal == true
    return 0.0, 0.0 # stop
  end

  if self.estimator !== nothing
    motion_update(self.estimator, self.prev_spd, self.prev_yr, self.delta_time)
    observation_update(self.estimator, observation)

    self.total_reward += self.delta_time * reward_per_sec(self)

    self.speed, self.yaw_rate = policy(self, self.estimator.estimated_pose)
    self.prev_spd, self.prev_yr = self.speed, self.yaw_rate

    draw!(self.estimator)
  end

  return self.speed, self.yaw_rate
end

max_stars_repo_name: ShisatoYano/JuliaAutonomy

max_issues_repo_name: ShisatoYano/JuliaAutonomy

max_forks_repo_name: ShisatoYano/JuliaAutonomy


text:
#!/usr/bin/env julia

using ArgParse
using JSON
using Word2Vec

s = ArgParseSettings()
@add_arg_table s begin
    "--vwify", "-v"
    help = "convert JSON-GAM input to vw compatible format"
    action = :store_true
    "--importance", "-i"
    help = "give errors a weight proportional to their mapping quality"
    action = :store_true
    "--phred", "-p"
    help = "regressor already outputs phred format, so don't convert"
    action = :store_true
    "--min-mq", "-q"
    help = "take the minimum of the previous and new mapping quality"
    action = :store_true
    "--mq-model", "-m"
    arg_type = String
    help = "set the MQ using the given model"
    default = ""
    "--dna-embedding", "-e"
    arg_type = String
    help = "use the DNA sequence embedding to generate sequence features"
    default = ""
    "--vw-opts", "-o"
    arg_type = String
    help = "additional options to pass to vw"
    default = ""
    "--kmer-size", "-k"
    arg_type = Int
    help = "kmer size to break input DNA sequences into"
    default = 5
    "--kmer-stride", "-s"
    arg_type = Int
    help = "bases between each successive kmer"
    default = 1
    "--cigar", "-c"
    help = "write the cigar base-by-base in a namespace"
    action = :store_true
    "--sequence", "-S"
    help = "write the read sequence base-by-base"
    action = :store_true
end

parsed_args = parse_args(ARGS, s)

embedding_file = parsed_args["dna-embedding"]
embedding = (embedding_file != "" ? wordvectors(embedding_file) : nothing)
kmer_size = parsed_args["kmer-size"]
kmer_stride = parsed_args["kmer-stride"]
write_cigar = parsed_args["cigar"]
write_sequence = parsed_args["sequence"]

global MQ_MAX = 60

function float2phred(prob)
    return min(-10 * log10(prob), MQ_MAX)
end

function phred2float(qual)
    return 10 ^ (qual * -.1)
end

function prob2mq(prob)
    v = 1-prob
    if v == 1
        return 0
    elseif v > 0
        return min(MQ_MAX, Int32(round(float2phred(v))))
    else
        return MQ_MAX
    end
end

function applymodel(alns, model, params)
    (path,io) = mktemp()
    for aln in alns
        println(io, vwify(aln, kmer_size, kmer_stride))
    end
    close(io)
    lines = split(chomp(readstring(pipeline(path, `vw -t -p /dev/stdout -i $model --quiet $(split(params, " "))`))), "\n")
    rm(path)
    assert(length(lines) == length(alns))
    for (line,aln) in zip(lines, alns)
        (result,name) = split(line, " ")
        assert(name == aln["name"])
        result = parsed_args["phred"] ? Int(round(float(result))) : prob2mq(float(result))
        mapping_quality = haskey(aln, "mapping_quality") ? aln["mapping_quality"] : 0
        aln["mapping_quality"] = (parsed_args["min-mq"] ? min(mapping_quality, result) : result)
    end
    return alns
end

function kmers_of(k::Int, j::Int, s::String)
    [s[i:(i+k-1)] for i in range(1,j,Int(floor((length(s)-k)/j))+1)]
end

function seq_embedding(s::String, k::Int, j::Int)
    function in_vocab(kmer::String)
        haskey(embedding.vocab_hash, kmer)
    end
    kmers = filter(in_vocab, kmers_of(k, j, s))
    reduce(+, [get_vector(embedding, kmer) for kmer in kmers])/length(kmers)
end

function vwify(aln, kmer_size::Int, kmer_stride::Int, mqimportance=false)
    name = haskey(aln, "name") ? aln["name"] : ""
    correct = haskey(aln, "correct") ? aln["correct"] : 0
    correct = correct > 0.5 ? 1 : 0
    if haskey(aln, "correct") aln["correct"] = correct end
    score = haskey(aln, "score") ? aln["score"] : 0
    identity = haskey(aln, "identity") ? aln["identity"] : 0
    mapqual = haskey(aln, "mapping_quality") ? aln["mapping_quality"] : 0
    path = haskey(aln, "path") && haskey(aln["path"], "mapping") ? aln["path"]["mapping"] : []
    i = 1
    cigar = []
    if write_cigar
        for mapping in path
            for elem in mapping["edit"]
                #println(elem)
                fromlen = haskey(elem, "from_length") ? elem["from_length"] : 0
                tolen = haskey(elem, "to_length") ? elem["to_length"] : 0
                seq = haskey(elem, "sequence") ? elem["sequence"] : ""
                if fromlen == tolen
                    if seq != ""
                        for j in range(i, tolen) push!(cigar, string(j, "S")) end
                    else
                        for j in range(i, tolen) push!(cigar, string(j, "M")) end
                    end
                elseif fromlen > tolen
                    for j in range(i, fromlen) push!(cigar, string(j, "D")) end
                elseif tolen > fromlen
                    for j in range(i, tolen) push!(cigar, string(j, "I")) end
                end
                i += tolen
            end
        end
    end
    sequence = write_sequence ? string("|bases ", join([string(x, y) for (x,y) in zip(1:length(aln["sequence"]),aln["sequence"])], " ")) : ""
    seqvec = (embedding == nothing ? "" : string("|dnavec ", join(["d$i:$(Float16(j))" for (i,j) in enumerate(seq_embedding(aln["sequence"], kmer_size, kmer_stride))], " ")))
    cigar = write_cigar ? string("|cigar ", join(cigar, " ")) : ""
    secondaries = haskey(aln, "secondary_score") ? aln["secondary_score"] : []
    secondaries = join(["s$i:$score" for (i,score) in zip(1:length(secondaries),secondaries)], " ")
    importance = (mqimportance && correct == 0) ? "$(mapqual+2) " : ((mqimportance && correct == 1 && mapqual < 60) ? "$(60-mapqual+1) " : " ")
    string("$correct $importance '$name |alignment identity:$identity score:$score mapqual:$mapqual |scores s0:$score $secondaries $seqvec $sequence $cigar")
end

alns = []
for line in eachline(STDIN)
    aln = JSON.parse(line)
    if parsed_args["vwify"]
        println(vwify(aln, kmer_size, kmer_stride, parsed_args["importance"]))
    elseif parsed_args["mq-model"] != ""
        # we should run the model on the input and set a new mapping quality
        push!(alns, aln)
        if length(alns) >= 1000
            assert(length(alns) > 0)
            for aln in applymodel(alns, parsed_args["mq-model"], parsed_args["vw-opts"])
                println(JSON.json(aln))
            end
            alns = []
        end
    else
        # default is to do nothing to input
        println(JSON.json(aln))
    end
end

# clean up the buffer
if parsed_args["mq-model"] != "" && length(alns) > 0
    for aln in applymodel(alns, parsed_args["mq-model"], parsed_args["vw-opts"])
        println(JSON.json(aln))
    end
    alns = []
end

max_stars_repo_name: ruolin/vg

max_issues_repo_name: ruolin/vg

max_forks_repo_name: ruolin/vg


text:
#################################################################
#
#    Definition of MCMCModel types
#
#################################################################

export model

######### parameters map info  ############
# These types are used to map scalars in the
#   parameter vector to user facing variables
#
immutable PDims
	pos::Int   # starting position of parameter in the parameter vector
	dims::Tuple    # dimensions of user facing parameter, can be a scalar, vector or matrix
end
typealias PMap Dict{Symbol, PDims}

function ispartition(m::PMap, n::Int)
	c = zeros(n)
	for v in values(m)
		c[v.pos:(v.pos+prod(v.dims)-1)] += 1
	end
	all(c .== 1.)
end

#### misc functions common to all models  ####
hasvectormethod(f::Function) = !isgeneric(f) | length(methods(f, (Vector{Float64},))) == 1
hasgradient{M<:MCMCModel}(m::M) = m.evalg != nothing
hastensor{M<:MCMCModel}(m::M) = m.evalt != nothing
hasdtensor{M<:MCMCModel}(m::M) = m.evaldt != nothing

#### User-facing model creation function  ####

# Currently only a "likelihood" model type makes sense
# Left as is in case other kind of models come up
function model(f::Union(Function, Expr); mtype="likelihood", args...)
	if mtype == "likelihood"
		return MCMCLikelihoodModel(f; args...)
	elseif mtype == "whatever"
	else
	end
end


text:
function load_board(str, go_env::GoEnv)
  reverse_map = Dict{Char, Int}([
      'X' => 1, # BLACK
      'O' => -1,  # WHITE
      '.' => 0, # EMPTY
      '#' => 2, # FILL
      '*' => 3, # Ko
      '?' => 4  # UNKNOWN
  ])
  str = replace(str, r"[^XO\.#]+" => s"")
  @assert length(str) == go_env.N ^ 2 # "Board to load didn't have right dimensions"
  board = zeros(Int8, go_env.N, go_env.N)
  for (i, char) in enumerate(str)
      board[i] = reverse_map[char]
  end
  board = permutedims(board, [2, 1])
  return board
end

function assertEqualArray(array1, array2)
  if !all(array1 .== array2)
    throw(AssertionError("Arrays differed in one or more locations:\n$array1\n%$array2\n"))
  end
  return true
end

function assertEqualLibTracker(lib_tracker1, lib_tracker2)
  # A lib tracker may have differently numbered groups yet still
  # represent the same set of groups.
  # "Sort" the group_ids to ensure they are the same.
  function find_group_mapping(lib_tracker)
    current_gid = 0
    mapping = Dict{Int16, Int}()
    for group_id in lib_tracker.group_index
      if group_id == MISSING_GROUP_ID
        continue
      end
      if group_id ∉ keys(mapping)
        mapping[group_id] = current_gid
        current_gid += 1
      end
    end
    return mapping
  end

  lt1_mapping = find_group_mapping(lib_tracker1)
  lt2_mapping = find_group_mapping(lib_tracker2)

  remapped_group_index1 = [get(lt1_mapping, gid, MISSING_GROUP_ID)
                            for gid in lib_tracker1.group_index]
  remapped_group_index2 = [get(lt2_mapping, gid, MISSING_GROUP_ID)
                            for gid in lib_tracker2.group_index]
  @assert remapped_group_index1 == remapped_group_index2

  remapped_groups1 = Dict(get(lt1_mapping, gid, 0) => group for (gid, group) in lib_tracker1.groups)
  remapped_groups2 = Dict(get(lt2_mapping, gid, 0) => group for (gid, group) in lib_tracker2.groups)
  @assert remapped_groups1 == remapped_groups2

  assertEqualArray(lib_tracker1.liberty_cache, lib_tracker2.liberty_cache)
end

function assertEqualPositions(pos1, pos2)
  assertEqualArray(pos1.board, pos2.board)
  assertEqualLibTracker(pos1.lib_tracker, pos2.lib_tracker)
  @assert pos1.n == pos2.n
  @assert pos1.caps == pos2.caps
  @assert pos1.ko == pos2.ko
  r_len = min(length(pos1.recent), length(pos2.recent))
  if r_len > 0 # if a position has no history, then don't bother testing
    @assert pos1.recent[end - r_len + 1:end] == pos2.recent[end - r_len + 1:end]
  end
  @assert pos1.to_play == pos2.to_play
  return true
end

function assertNoPendingVirtualLosses(root)
  # Raise an error if any node in this subtree has vlosses pending.
  queue = [root]
  while !isempty(queue)
    current = pop!(queue)
    @assert current.losses_applied == 0
    queue = vcat(queue, collect(values(current.children)))
  end
  return true
end


text:
using DataFrames
using Grids
using PrettyTables

struct PolicyValueTest <: TrainingTest
    idx::Integer
    state::AbstractEnvState
    image::AbstractImage
    actions::Vector{<:Integer}
    value::Float32
end

policyvaluetestsvg(::AbstractEnvState, ::PolicyValueTest) = @assert false OVERRIDE_MSG

mutable struct PolicyValueTests <: TrainingTests
    tests::Vector{PolicyValueTest}
    images::AbstractArray{Float32, 4}
    actual_value::AbstractArray{Float32, 2}
    period::Int
    warmup::Int
    function PolicyValueTests(tests; period=1000, warmup=0)
        new(
            tests,
            tensorize(get_image.(tests)),
            Flux.unsqueeze(get_value.(tests), 1),
            period,
            warmup,
        )
    end
end

function runtests(tests::PolicyValueTests, net)
    pvtestfn = gettestfn(tests)
    acc, fails, value_mse = pvtestfn(net)
    acc, fails, value_mse
end

function gettestfn(tests::PolicyValueTests)
    function fn(n::Network)
        # if typeof(n) != UniformPolicy
        #     @show Flux.params(n.representation_network)[1][1]
        # end
        initial_output = initialinferencefn(tests)(n)
        acc, fails = getpolicyacc(tests, initial_output)
        acc, fails, getvaluemse(tests, initial_output)
    end
end

Base.lastindex(x::PolicyValueTests) = x.tests[end]

function df(tests::PolicyValueTests)
    test_num = "Test #"
    optimal_action = "Actual (π*)"
    d = Dict(
        test_num=>[],
        optimal_action=>[],
    )
    for (idx,test) in enumerate(tests.tests)
        push!(d[test_num], idx)
        push!(d[optimal_action], join(test.actions, ","))
    end
    DataFrame(
        test_num=>d[test_num],
        optimal_action=>d[optimal_action],
    )
end

function Base.show(io::IO, m::MIME"text/plain", x::PolicyValueTests)
    len = length(x.tests)
    title_paren = "all $len"
    pretty_table(
        io,
        df(x),
        show_row_number=false,
        nosubheader=true,
        title="Policy Tests ($title_paren)"
    )
end

function Base.show(io::IO, m::MIME"text/plain", x::PolicyValueTest)
    println("\nPOLICY TEST #$(x.idx)")
    println("\nState:\n$(x.state)")
    println("\nTo Play:\n$(x.state.to_play.name)")
    println("\nOptimal Action(s):\n$(join(x.actions, ","))\n")
end

struct PolicyValueTestSVG
    test::PolicyValueTest
end
svg(x::PolicyValueTest) = PolicyValueTestSVG(x)

function Base.show(io::IO, m::MIME"image/svg+xml", x::PolicyValueTestSVG)
    test = x.test
    s = test.state
    player_name = s.to_play.name
    actions = test.actions
    state_grid = getgrid(s; header="to play: $player_name", margin_left=0)
    next_states = Rect[]
    from = 1
    len = length(actions)
    to = min(len, 3)
    for (idx,a) in enumerate(actions[from:to])
        push!(
            next_states,
            getgrid(s, a, 0f0, getnextstate(s, a); header="action: $a", margin_left=0)
        )
    end
    if length(actions) > 3
        push!(next_states, Grid(1,1, (r,c)->"..."; line_thickness=0))
    end
    input_grid = Grid(
        [state_grid],
        :,
        1,
        header="INPUT",
        margin_left=0,
        line_thickness=0
    )
    policy_grid_header = "ACTUAL π*"
    policy_grid = Grid(
        next_states,
        :,
        1,
        header=(text=policy_grid_header, height=40),
        margin_left=0,
        line_thickness=0,
        bg_color=RandColor()
    )
    g = Grid([input_grid,policy_grid]; cell_valign=:top, header="POLICY VALUE TEST #$(test.idx)", line_thickness=0)
    show(io, m, g)
end

struct PolicyValueTestsSVG
    tests::PolicyValueTests
end
svg(x::PolicyValueTests) = PolicyValueTestsSVG(x)


text:
function find_target!(m::Model, p::Pedestrian)
    #  find target based on selected target strategy
    find_target!(m.target_strategy, m, p)

    # checks wheter pedestrian is in exit
    if p.isexit && p.pos == p.target
        p.vel = (0., 0.)
        p.finished = true
        return
    end
end

"""
    Nearest <: TargetStrategy

Basic target strategy that selects nearest (euclidian distance) target.

# Fields

- p::Float64 = 0.5
- ε::Float64 = 0.0

"""
Base.@kwdef struct Nearest <: TargetStrategy
    p::Float64 = 0.5
    ε::Float64 = 0
end

function find_target!(s::Nearest, m::Model, p::Pedestrian)
    dist = Inf
    k = length(m.room.checkpoints)

    # find nearest checkpoint
    for (i, t) in enumerate(vcat(m.room.checkpoints, m.room.exits))
        t_pos = nearest(t, p.pos)
        p.pos[2] > t_pos[2] || continue

        d = distance(p.pos, t_pos)
        if d < dist || abs(d - dist) <= s.ε && rand() > s.p
            dist = d
            p.target = t_pos
            p.isexit = i > k
        end
    end
    return
end


text:
import AlgebraicMultigrid: scale_cols_by_largest_entry!,
            SymmetricStrength, poisson
function symmetric_soc(A::SparseMatrixCSC{T,V}, θ) where {T,V}
    D = abs.(diag(A))
    i,j,v = findnz(A)
    mask = i .!= j
    DD = D[i] .* D[j]
    mask = mask .& (abs.(v.^2) .>= (θ * θ * DD))

    i = i[mask]
    j = j[mask]
    v = v[mask]

    S = sparse(i,j,v, size(A)...) + spdiagm(0=>D)

    for i = 1:size(S.nzval,1)
        S.nzval[i] = abs(S.nzval[i])
    end

    scale_cols_by_largest_entry!(S)

    S
end

# Set up tests
function test_symmetric_soc()

    cases = generate_matrices()

    for matrix in cases
        for θ in (0.0, 0.1, 0.5, 1., 10.)
            ref_matrix = symmetric_soc(matrix, θ)
            calc_matrix, _ = SymmetricStrength(θ)(matrix)

            @test sum(abs2, ref_matrix - calc_matrix) < 1e-6
        end
    end
end

function generate_matrices()

    cases = []

    # Random matrices
    seed!(0)
    for T in (Float32, Float64)

        for s in [2, 3, 5]
            push!(cases, sprand(T, s, s, 1.))
        end

        for s in [2, 3, 5, 7, 10, 11, 19]
            push!(cases, poisson(T, s))
        end
    end

    cases
end

function stand_agg(C)
    n = size(C, 1)

    R = Set(1:n)
    j = 0
    Cpts = Int[]

    aggregates = -ones(Int, n)

    # Pass 1
    for i = 1:n
        Ni = union!(Set(C.rowval[nzrange(C, i)]), Set(i))
        if issubset(Ni, R)
            push!(Cpts, i)
            setdiff!(R, Ni)
            for x in Ni
                aggregates[x] = j
            end
            j += 1
        end
    end

    # Pass 2
    old_R = copy(R)
    for i = 1:n
        if ! (i in R)
            continue
        end

        for x in C.rowval[nzrange(C, i)]
            if !(x in old_R)
                aggregates[i] = aggregates[x]
                setdiff!(R, i)
                break
            end
        end
    end

    # Pass 3
    for i = 1:n
        if !(i in R)
            continue
        end
        Ni = union(Set(C.rowval[nzrange(C,i)]), Set(i))
        push!(Cpts, i)

        for x in Ni
            if x in R
                aggregates[x] = j
            end
            j += 1
        end
    end

    @assert length(R) == 0

    Pj = aggregates .+ 1
    Pp = collect(1:n+1)
    Px = ones(eltype(C), n)

    SparseMatrixCSC(maximum(aggregates .+ 1), n, Pp, Pj, Px)
end

# Standard aggregation tests
function test_standard_aggregation()

    cases = generate_matrices()

    for matrix in cases
        for θ in (0.0, 0.1, 0.5, 1., 10.)
            C = symmetric_soc(matrix, θ)
            calc_matrix = StandardAggregation()(matrix)
            ref_matrix = stand_agg(matrix)
            @test sum(abs2, ref_matrix - calc_matrix) < 1e-6
        end
    end

end

# Test fit_candidates
function test_fit_candidates()

    cases = generate_fit_candidates_cases()

    for (i, (AggOp, fine_candidates)) in enumerate(cases)

        mask_candidates!(AggOp, fine_candidates)

        Q, coarse_candidates = fit_candidates(AggOp, fine_candidates)

        @test isapprox(fine_candidates, Q * coarse_candidates)
        @test isapprox(Q * (Q' * fine_candidates), fine_candidates)
    end
end
function mask_candidates!(A,B)
    B[(diff(A.colptr) .== 0)] .= 0
end

function generate_fit_candidates_cases()
    cases = []

    for T in (Float32, Float64)

        # One candidate
        AggOp = SparseMatrixCSC(2, 5, collect(1:6),
                        [1,1,1,2,2], ones(T,5))
        B =  ones(T,5)
        push!(cases, (AggOp, B))

        AggOp = SparseMatrixCSC(2, 5, collect(1:6),
                        [2,2,1,1,1], ones(T,5))
        B = ones(T, 5)
        push!(cases, (AggOp, B))

        AggOp = SparseMatrixCSC(3, 9, collect(1:10),
                        [1,1,1,2,2,2,3,3,3], ones(T, 9))
        B = ones(T, 9)
        push!(cases, (AggOp, B))

        AggOp = SparseMatrixCSC(3, 9, collect(1:10),
                        [3,2,1,1,2,3,2,1,3], ones(T,9))
        B = T.(collect(1:9))
        push!(cases, (AggOp, B))
    end

    cases
end

# Test approximate spectral radius
function test_approximate_spectral_radius()

    cases = []
    seed!(0)

    push!(cases, [2. 0.
                  0. 1.])

    push!(cases, [-2. 0.
                   0  1])

    push!(cases, [100.   0.  0.
                    0. 101.  0.
                    0.   0. 99.])

    for i in 2:5
        push!(cases, rand(i,i))
    end

    for A in cases
        @static if VERSION < v"0.7-"
            E,V = eig(A)
        else
            E,V = (eigen(A)...,)
        end
        E = abs.(E)
        largest_eig = findall(E .== maximum(E))[1]
        expected_eig = E[largest_eig]

        @test isapprox(approximate_spectral_radius(A), expected_eig)
    end

    # Symmetric matrices
    for A in cases
        A = A + A'
        @static if VERSION < v"0.7-"
            E,V = eig(A)
        else
            E,V = (eigen(A)...,)
        end
        E = abs.(E)
        largest_eig = findall(E .== maximum(E))[1]
        expected_eig = E[largest_eig]

        @test isapprox(approximate_spectral_radius(A), expected_eig)
    end
end

# Test Gauss Seidel
import AlgebraicMultigrid: gs!
function test_gauss_seidel()
    N = 1
    A = spdiagm(0 => 2 * ones(N), -1 => -ones(N-1), 1 => -ones(N-1))
    x = eltype(A).(collect(0:N-1))
    b = zeros(N)
    s = GaussSeidel(ForwardSweep())
    s(A, x, b)
    @test sum(abs2, x - zeros(N)) < 1e-8

    N = 3
    A = spdiagm(0 => 2 * ones(N), -1 => -ones(N-1), 1 => -ones(N-1))
    x = eltype(A).(collect(0:N-1))
    b = zeros(N)
    s = GaussSeidel(ForwardSweep())
    s(A, x, b)
    @test sum(abs2, x - [1.0/2.0, 5.0/4.0, 5.0/8.0]) < 1e-8

    N = 1
    A = spdiagm(0 => 2 * ones(N), -1 => -ones(N-1), 1 => -ones(N-1))
    x = eltype(A).(collect(0:N-1))
    b = zeros(N)
    s = GaussSeidel(BackwardSweep())
    s(A, x, b)
    @test sum(abs2, x - zeros(N)) < 1e-8

    N = 3
    A = spdiagm(0 => 2 * ones(N), -1 => -ones(N-1), 1 => -ones(N-1))
    x = eltype(A).(collect(0:N-1))
    b = zeros(N)
    s = GaussSeidel(BackwardSweep())
    s(A, x, b)
    @test sum(abs2, x - [1.0/8.0, 1.0/4.0, 1.0/2.0]) < 1e-8

    N = 1
    A = spdiagm(0 => 2 * ones(N), -1 => -ones(N-1), 1 => -ones(N-1))
    x = eltype(A).(collect(0:N-1))
    b = eltype(A).([10.])
    s = GaussSeidel(ForwardSweep())
    s(A, x, b)
    @test sum(abs2, x - [5.]) < 1e-8

    N = 3
    A = spdiagm(0 => 2 * ones(N), -1 => -ones(N-1), 1 => -ones(N-1))
    x = eltype(A).(collect(0:N-1))
    b = eltype(A).([10., 20., 30.])
    s = GaussSeidel(ForwardSweep())
    s(A, x, b)
    @test sum(abs2, x - [11.0/2.0, 55.0/4, 175.0/8.0]) < 1e-8

    N = 100
    A = spdiagm(0 => 2 * ones(N), -1 => -ones(N-1), 1 => -ones(N-1))
    x = ones(eltype(A), N)
    b = zeros(eltype(A), N)
    s1 = GaussSeidel(ForwardSweep(), 200)
    s1(A, x, b)
    resid1 = norm(A*x,2)
    x = ones(eltype(A), N)
    s2 = GaussSeidel(BackwardSweep(), 200)
    s2(A, x, b)
    resid2 = norm(A*x,2)
    @test resid1 < 0.01 && resid2 < 0.01
    @test isapprox(resid1, resid2)

end

# Test jacobi smooth prolongator with local weighting
function test_jacobi_prolongator()
    A = poisson(100)
    T = poisson(100)
    x = JacobiProlongation(4/3)(A, T, 1, 1)
    ref = include("ref_R.jl")
    @test sum(abs2, x - ref) < 1e-6
end

# Issue #24
function nodes_not_agg()
    A = include("onetoall.jl")
    ml = smoothed_aggregation(A)
    @test size(ml.levels[2].A) == (11,11)
    @test size(ml.final_A) == (2,2)
end

# Issue 26
function test_symmetric_sweep()
    A = poisson(10)
    s = GaussSeidel(SymmetricSweep(), 4)
    x = ones(size(A,1))
    b = zeros(size(A,1))
    s(A, x, b)
    @test sum(abs2, x - [0.176765; 0.353529; 0.497517; 0.598914;
                            0.653311; 0.659104; 0.615597; 0.52275;
                            0.382787; 0.203251]) < 1e-6
end


text:
module test_Plot1

import ..ModiaMath
using Test
t = range(0.0, stop=10.0, length=100)

result = Dict{AbstractString,Any}()

result["time"] = t
result["phi"]  = sin.(t)

ModiaMath.plot(result, :phi, heading="Sine(time)")

# Test also that warn works for non-existing variables
println("... Next plot should give a warning:")
ModiaMath.plot(result, :SignalNotDefined, heading="Sine(time)", figure=2)

# Print result variables
println("\n... result variables = ", ModiaMath.resultTable(result))

end


text:
abstract type Engine end

get_engine_position(eng::Engine) = [0, 0, 0]
get_engine_orientation(eng::Engine) = [0, 0, 0]
calculate_engine(eng::Engine) = error("abstract method")
get_engine_gyro_effects(eng::Engine) = [0.0, 0.0, 0.0]  # [kg·m²/s]

#Getters
get_pfm(eng::Engine) = eng.pfm
get_cj(eng::Engine) = eng.cj
get_power(eng::Engine) = eng.power
get_efficiency(eng::Engine) = eng.efficiency
get_tanks(eng::Engine) = eng.tanks
get_gyro_effects(eng::Engine) = eng.h  # angular momentum [kg·m²/s]

struct Propulsion
    pfm::PointForcesMoments
    cj::Number
    power::Number
    efficiency::Number
    engines::Array{Engine, 1}
    # Gyroscopic effects
    h::Array{T, 1} where T<:Number  # angular momentum [kg·m²/s]
end

# Getters
get_pfm(prop::Propulsion) = prop.pfm
get_cj(prop::Propulsion) = prop.cj
get_power(prop::Propulsion) = prop.power
get_efficiency(prop::Propulsion) = prop.efficiency
get_engines(prop::Propulsion) = prop.engines
get_gyro_effects(prop::Propulsion) = prop.h

function get_propulsion_position(prop::Propulsion)
    pos = [0.0, 0.0, 0.0]
    engines = get_engines(prop)
    for eng=engines
        pos += get_engine_position(eng)
    end
    return pos / length(engines)
end

function get_tanks(prop::Propulsion)
    tanks = RigidSolid[]
    for eng=get_engines(prop)
        for t=get_tanks(eng)
            push!(tanks, t)
        end
    end
    return tanks
end

get_fuel_mass_props(prop::Propulsion) = sum(get_tanks(prop))

function calculate_propulsion(prop::Propulsion, fcs::FCS, aerostate::AeroState,
                              state::State; consume_fuel=false)

    engines = get_engines(prop)

    power = 0.
    pfm = PointForcesMoments(get_propulsion_position(prop), [0, 0, 0], [0, 0, 0])
    cj = 0
    efficiency = 0
    gyro = [0.0, 0.0, 0.0]

    for (ii, eng)=enumerate(engines)
        engines[ii] = calculate_engine(eng, fcs, aerostate, state;
                                       consume_fuel=consume_fuel)
        # TODO: use get_engine_orientation and rotate every engine pfm to body
        pfm += rotate(get_pfm(engines[ii]), get_engine_orientation(eng)...)
        power += get_power(engines[ii])
        gyro += get_gyro_effects(engines[ii])
        # TODO: cj and efficiency
    end

    Propulsion(pfm, cj, power, efficiency, engines, gyro)
end


text:
struct rt_tuple
    x::Float64
    y::Float64
    z::Float64
    w::Int
end

struct rt_color
    red::Float64
    green::Float64
    blue::Float64
end

function color(r::Real,g::Real,b::Real)
    return rt_color(promote(r,g,b)...)
end

function tuple(x::Real,y::Real,z::Real,w::Int)
    return rt_tuple(promote(x,y,z)...,w)
end

"""
    vector(x,y,z)

Create an `rt_tuple` which represents a vector, where `w=0`.

#Examples
```julia-repl
julia> a = RayTracer.vector(1,2,3)
RayTracer.rt_tuple{Int64}(1, 2, 3, 0)
```
"""
function vector(x,y,z)
    return tuple(x,y,z,0)
end

"""
    point(x,y,z)

Create an `rt_tuple` which represents a point, where `w=1`.

#Examples
```julia-repl
julia> a = RayTracer.point(1,2,3)
RayTracer.rt_tuple{Int64}(1, 2, 3, 1)
```
"""
function point(x,y,z)
    return tuple(x,y,z,1)
end

"Check whether rt_tuple{T} is a vector (w=0) or not."
function is_vector(a::rt_tuple)
    if a.w == 0
        return true
    else
        return false
    end
end

"Check whether rt_tuple{T} is a point (w=1) or not."
function is_point(a::rt_tuple)
    if a.w == 1
        return true
    else
        return false
    end
end

# Tuple Operations

function Base.isapprox(a::rt_tuple,b::rt_tuple,atol::Real=0)
    return (≈(a.x,b.x,atol=atol)) && (≈(a.y,b.y,atol=atol)) && (≈(a.z,b.z,atol=atol)) && (a.w≈b.w)
end

function Base.:+(a::rt_tuple,b::rt_tuple)::rt_tuple
    return tuple(a.x+b.x,a.y+b.y,a.z+b.z,a.w+b.w)
end

function Base.:-(a::rt_tuple,b::rt_tuple)::rt_tuple
    return tuple(a.x-b.x,a.y-b.y,a.z-b.z,a.w-b.w)
end

function Base.:-(a::rt_tuple)::rt_tuple
    return tuple(-a.x,-a.y,-a.z,-a.w)
end

function Base.:*(a::rt_tuple,b::Real)::rt_tuple
    return tuple(b*a.x,b*a.y,b*a.z,a.w)
end

function Base.:*(b::Real,a::rt_tuple)::rt_tuple
    return tuple(b*a.x,b*a.y,b*a.z,a.w)
end

function Base.:/(a::rt_tuple,b::Real)::rt_tuple
    return tuple(a.x/b,a.y/b,a.z/b,a.w)
end

function Base.:*(b::Array{T,2},a::rt_tuple)::rt_tuple where {T<:Real}
    v = [a.x;a.y;a.z;a.w]
    r = b*v
    return tuple(r[1],r[2],r[3],convert(Int, round(r[4], digits=0)))
end

# Tuple Methods

"""
    modulus(v::rt_tuple{T<:Real})

Compute the modulus (also known as magnitude) of the vector `v`.

# Examples
```julia-repl
julia> a = RayTracer.vector(1,2,3); modulus(a)
√14
```
"""
function modulus(v::rt_tuple)
    return sqrt(v.x^2 + v.y^2 + v.z^2)
end

"""
    normalize(v::rt_tuple{T<:Real})

Compute unit vector pointing in the same direction as `v`.

# Examples
```julia-repl
julia> a = RayTracer.vector(1,2,3); normalize(a)
RayTracer.rt_tuple{Float64}(1/√14, 2/√14, 3/√14, 0)
```
"""
function normalize(v::rt_tuple)
    return v/modulus(v)
end

"""
    dot(a::rt_tuple,b::rt_tuple)

Compute dot product between two vectors.

# Examples
```julia-repl
julia> a = RayTracer.vector(1,2,3);
julia> b = RayTracer.vector(2,3,4);
julia> dot(a,b)
20
```
"""
function dot(a::rt_tuple,b::rt_tuple)
    return a.x*b.x + a.y*b.y + a.z*b.z
end

"""
    cross(a::rt_tuple,b::rt_tuple)

Compute cross product between two vectors.

# Examples
```julia-repl
julia> a = RayTracer.vector(1,2,3);
julia> b = RayTracer.vector(2,3,4);
julia> dot(a,b)
20
```
"""
function cross(a::rt_tuple,b::rt_tuple)
    return vector(a.y*b.z-a.z*b.y,a.z*b.x-a.x*b.z,a.x*b.y-a.y*b.x)
end

"""
    reflect(v::rt_tuple,n::rt_tuple)

Reflect a vector 'v' on a surface with a normal vector 'n'.
"""
function reflect(v::rt_tuple,n::rt_tuple)::rt_tuple
    return v - 2*dot(v,n)*n
end

# Color Operations

function Base.isapprox(a::rt_color,b::rt_color,atol::Real=0)
    return (≈(a.red,b.red,atol=atol)) && (≈(a.green,b.green,atol=atol)) && (≈(a.blue,b.blue,atol=atol))
end

function Base.:+(a::rt_color,b::rt_color)::rt_color
    return color(a.red+b.red,a.green+b.green,a.blue+b.blue)
end

function Base.:-(a::rt_color,b::rt_color)::rt_color
    return color(a.red-b.red,a.green-b.green,a.blue-b.blue)
end

function Base.:*(a::rt_color,b::Real)::rt_color
    return color(b*a.red,b*a.green,b*a.blue)
end

function Base.:*(b::Real,a::rt_color)::rt_color
    return color(b*a.red,b*a.green,b*a.blue)
end

function Base.:*(a::rt_color,b::rt_color)::rt_color
    return color(b.red*a.red,b.green*a.green,b.blue*a.blue)
end

# Matrix Operations

function submatrix(a::Array{T,2},r::Int,c::Int)::Array{T,2} where {T<:Real}
    return a[1:end .!= r,1:end .!=c]
end

function minor(a::Array{T,2},r::Int,c::Int)::T where {T<:Real}
    return det(submatrix(a,r,c))
end

function cofactor(a::Array{T,2},r::Int,c::Int)::T where {T<:Real}
    return (1*((r+c)%2==0) - 1*((r+c)%2==1))*minor(a,r,c)
end


text:
using Test, NetCDF, Random

"Function that returns a valid file name of a non-existing file for all OSs"
function tempname2()
    f, h = mktemp()
    close(h)
    rm(f)
    f
end
tlist = [Int8, UInt8, Int16, UInt16, Int32, UInt32, Int64, UInt64, Float32, Float64]

@testset "NetCDF" begin
    include("intermediate.jl")
    include("high.jl")
    include("array-like.jl")
    include("chunks.jl")
end
nothing


text:
using Partial
include("test-partial.jl")


text:
@json struct Identity
    provider::String
    extern_uid::String
end

@json struct User
    id::Int
    username::String
    email::String
    name::String
    state::String
    avatar_url::String
    web_url::String
    created_at::DateTime
    is_admin::Bool
    bio::String
    location::String
    public_email::String
    skype::String
    linkedin::String
    twitter::String
    website_url::String
    organization::String
    last_sign_in_at::DateTime
    confirmed_at::DateTime
    theme_id::Int
    last_activity_on::Date
    color_scheme_id::Int
    projects_limit::Int
    current_sign_in_at::DateTime
    identities::Vector{Identity}
    can_create_group::Bool
    can_create_project::Bool
    two_factor_enabled::Bool
    external::Bool
    private_profile::Bool
    # Undocumented
    shared_runners_minutes_limit::Int
end

endpoint(::GitLabAPI, ::typeof(get_user)) = Endpoint(:GET, "/user")
endpoint(::GitLabAPI, ::typeof(get_user), id::Integer) = Endpoint(:GET, "/users/$id")
endpoint(::GitLabAPI, ::typeof(get_user), name::AStr) =
    Endpoint(:GET, "/users"; query=Dict(:username => name))
postprocessor(::GitLabAPI, ::typeof(get_user)) = DoSomething() do r
    v = JSON2.read(IOBuffer(r.body), Union{User, Vector{User}})
    return if v isa User
        v
    else
        isempty(v) ? nothing : v[1]
    end
end
into(::GitLabAPI, ::typeof(get_user)) = User

endpoint(::GitLabAPI, ::typeof(get_users)) = Endpoint(:GET, "/users")
into(::GitLabAPI, ::typeof(get_users)) = Vector{User}

endpoint(::GitLabAPI, ::typeof(update_user), id::Integer) = Endpoint(:PUT, "/users/$id")
postprocessor(::GitLabAPI, ::typeof(update_user)) = DoNothing()

endpoint(::GitLabAPI, ::typeof(create_user)) = Endpoint(:POST, "/users")
into(::GitLabAPI, ::typeof(create_user)) = User

endpoint(::GitLabAPI, ::typeof(delete_user), id::Integer) = Endpoint(:DELETE, "/users/$id")
postprocessor(::GitLabAPI, ::typeof(delete_user)) = DoNothing()


text:
module OptimalEstimation
using Parameters
using DocStringExtensions
using LinearAlgebra

include("types.jl")
include("linear_least_squares.jl")
include("optimal_estimation.jl")

end # module

max_stars_repo_name: RadiativeTransfer/OptimalEstimation

max_issues_repo_name: RadiativeTransfer/OptimalEstimation

max_forks_repo_name: RadiativeTransfer/OptimalEstimation


text:
# batched routines
@inline indexpos(iAs, i)::Int = findfirst(==(i), iAs)

# can be used in either static or dynamic invoke
function analyse_batched_dim(iAs, iBs, iOuts)
    iABs = iAs ∩ iBs
    pres   = iABs ∩ iOuts
    broad  = setdiff((iAs ∩ iOuts) ∪ (iBs ∩ iOuts), pres)
    summed = setdiff(iABs, pres)

    iAps, iAbs, iAss = pres ∩ iAs, broad ∩ iAs, summed ∩ iAs
    iBps, iBbs, iBss = pres ∩ iBs, broad ∩ iBs, summed ∩ iBs

    pA   = indexpos.(Ref(iAs), vcat(iAbs, iAss, iAps))
    pB   = indexpos.(Ref(iBs), vcat(iBss, iBbs, iBps))
    iABs = vcat(iAbs, iBbs, iAps)
    pOut = indexpos.(Ref(iABs), iOuts)
    return (pA...,), (iAps...,), (iAbs...,), (iAss...,), (pB...,), (iBps...,), (iBbs...,), (iBss...,), pOut
end

function analyse_batched_size(iAs, iAps, iAbs, iAss, sA, iBs, iBps, iBbs, iBss, sB)
    sAbs = getindex.(Ref(sA), indexpos.(Ref(iAs), iAbs))
    sAb = reduce(*, sAbs, init=1)
    sAs = mapreduce(i->sA[indexpos(iAs,i)], *, iAss, init=1)
    sAps = getindex.(Ref(sA), indexpos.(Ref(iAs),iAps))
    sAp = reduce(*, sAps, init=1)

    sBbs = getindex.(Ref(sB), indexpos.(Ref(iBs),iBbs))
    sBb = reduce(*, sBbs, init=1)
    sBs = mapreduce(i->sB[indexpos(iBs, i)], *, iBss, init=1)
    sBp = mapreduce(i->sB[indexpos(iBs, i)], *, iBps, init=1)
    sAB = (sAbs..., sBbs..., sAps...)
    return sAb, sAs, sAp, sBs, sBb, sBp, sAB
end

# batched, dynamic version
@generated function batched_contract(::Val{iAs}, A::AbstractArray, ::Val{iBs}, B::AbstractArray, ::Val{iOuts}) where {iAs, iBs, iOuts, NO,T}
    pA, iAps, iAbs, iAss, pB, iBps, iBbs, iBss, pOut = analyse_batched_dim(iAs, iBs, iOuts)
    quote
        sAb, sAs, sAp, sBs, sBb, sBp, sAB = analyse_batched_size($iAs, $iAps, $iAbs, $iAss, size(A), $iBs, $iBps, $iBbs, $iBss, size(B))

        A, B = align_eltypes(A, B)
        Apr = reshape(conditioned_permutedims(A, $pA, $iAs), sAb, sAs, sAp)
        Bpr = reshape(conditioned_permutedims(B, $pB, $iBs), sBs, sBb, sBp)
        AB = _batched_gemm('N','N', Apr, Bpr)
        AB = conditioned_permutedims(reshape(AB, sAB...), $((pOut...,)), $iOuts)
    end
end

# reload this function for GPU support!
function _batched_gemm(C1::Char, C2::Char, A::StridedArray{T,3}, B::StridedArray{T2,3}) where {T<:BlasFloat, T2<:BlasFloat}
    batched_gemm(C1, C2, A, B)
end

function _batched_gemm(C1::Char, C2::Char, A::AbstractArray{T,3}, B::AbstractArray{T2,3}) where {T<:BlasFloat, T2<:BlasFloat}
    batched_gemm(C1, C2, Array(A), Array(B))
end


text:
using PyPlot, Revise, transD_GP, Random, Statistics
include("skytem_response.jl")
##  modeling parameters
ntimesperdecade = 10
nfreqsperdecade = 5
freqlow = 1e-3
## LM operator
Flm = transD_GP.AEM_VMD_HMD.HFieldDHT(
                      ntimesperdecade = ntimesperdecade,
                      nfreqsperdecade = nfreqsperdecade,
                      lowpassfcs = lowpassfcs,
                      times  = LM_times,
                      ramp   = LM_ramp,
                      zTx    = zTx,
                      rRx    = rRx,
                      rTx    = rTx,
                      zRx    = zRxLM,
                      freqlow = freqlow
                      )
## HM operator
Fhm = transD_GP.AEM_VMD_HMD.HFieldDHT(
                      ntimesperdecade = ntimesperdecade,
                      nfreqsperdecade = nfreqsperdecade,
                      lowpassfcs = lowpassfcs,
                      times  = HM_times,
                      ramp   = HM_ramp,
                      zTx    = zTx,
                      rRx    = rRx,
                      rTx    = rTx,
                      zRx    = zRxHM,
                      freqlow = freqlow
                      )
## get the fields in place
@time transD_GP.AEM_VMD_HMD.getfieldTD!(Flm, zfixed, rho)
@time transD_GP.AEM_VMD_HMD.getfieldTD!(Fhm, zfixed, rho)
## plot
figure()
s1 = subplot(211)
loglog(Flm.times,4*pi*1e-7*abs.(Flm.dBzdt), label="lm_julia")
loglog(Fhm.times,4*pi*1e-7*abs.(Fhm.dBzdt), label="hm_julia")
loglog(Flm.times,dLM, ".", label="lm_skytem")
loglog(Fhm.times,dHM, ".", label="hm_skytem")
legend()
grid(which="both")
ylabel("dBzdt V/Am^4")
s2 = subplot(212, sharex=s1)
semilogx(Flm.times,(4*pi*1e-7*Flm.dBzdt-dLM)./dLM*100, label="lm")
semilogx(Fhm.times,(4*pi*1e-7*Fhm.dBzdt-dHM)./dHM*100, label="hm")
ylabel("% difference")
xlabel("time s")
legend()
grid(which="both")


text:
using FileIO
using Base.Test

@testset "FileIO" begin
    include("query.jl")
    include("loadsave.jl")
    include("error_handling.jl")
end


text:
module LatticeUtilities

using LinearAlgebra
using DocStringExtensions
using Printf

include("UnitCell.jl")
export UnitCell
export loc_to_pos!, loc_to_pos
export displacement_to_vec!, displacement_to_vec

include("Lattice.jl")
export Lattice
export valid_location, pbc!
export unitcell_to_loc!, unitcell_to_loc, loc_to_unitcell

include("Bond.jl")
export Bond

include("utilities.jl")
export get_num_sites, valid_site
export site_to_unitcell, site_to_orbital
export site_to_loc!, site_to_loc
export loc_to_site, site_to_site
export calc_k_point!, calc_k_point
export calc_k_points!, calc_k_points
export bond_to_vec!, bond_to_vec
export build_neighbor_table
export simplify!
export sites_to_displacement!, sites_to_displacement
export sites_to_bond!, sites_to_bond
export map_neighbor_table

end # module


text:
using Revise, InfiniteLinearAlgebra, BlockBandedMatrices, BandedMatrices, BlockArrays, InfiniteArrays, FillArrays, LazyArrays, Test

A = BlockTridiagonal(Vcat([fill(1.0,2,1),Matrix(1.0I,2,2),Matrix(1.0I,2,2),Matrix(1.0I,2,2)],Fill(Matrix(1.0I,2,2), ∞)),
                       Vcat([zeros(1,1)], Fill(zeros(2,2), ∞)),
                       Vcat([fill(1.0,1,2),Matrix(1.0I,2,2)], Fill(Matrix(1.0I,2,2), ∞)))

BlockBandedMatrix(A)
isb


text:
using LinearAlgebra
using MAT
using Plots
using Statistics

file = matread("Data2.mat");
X = file["data"];

# Visualize the first two dimension of the data
scatter(X[:,1], X[:,2], aspect_ratio=:equal, leg=false)

# input: X - data points
# output: E - distance matrix
#
function get_E(X)
    n = size(X,1); # number of points
    E = zeros(n,n); # the diagonal entries are 0

    for ii=1:n
        for jj=ii+1:n
            E[ii,jj] = norm(X[ii,:] - X[jj,:]);
            E[jj,ii] = E[ii,jj];
        end
    end

    return E;
end

# get the distance matrix D
E = get_E(X);

# evaluate sigma
E_sort = sort(E, dims=2);

k = 7;
sigma_loc = E_sort[:, k+1];

# input 1: E - distance matrix
# input 2: sigma - constant
#
# output: K - kernal matrix
#
function get_K(E, sigma_loc)
    n = size(E,1);
    K = ones(n,n);

    for ii = 1:n
        for jj = ii+1:n
            K[ii,jj] = exp(-E[ii,jj]^2/(sigma_loc[ii]*sigma_loc[jj]));
            K[jj,ii] = K[ii,jj]
        end
    end

    return K;
end

# get the kernal matrix K
K = get_K(E, sigma_loc);

# input: K - kernal matrix
#
# output 1: Q
# output 2: d_sq - sqrt{D}
#
function get_Q(K)
    n = size(K,1);
    Q = zeros(n,n);
    d_sq = zeros(n);

    # d_sq = sqrt{D}
    for ii = 1:n
        d_sq[ii] = sqrt(sum(K[ii,:]));
    end

    # get components of Q
    for ii = 1:n
        for jj = 1:n
            Q[ii,jj] = K[ii,jj]/(d_sq[ii]*d_sq[jj]);
        end
    end

    return Q, d_sq;
end

# get Q and d_sq
Q, d_sq = get_Q(K);

# input 1: Q
# input 2: d_sq - sqrt{D}
#
# output 1: v - eigenvectors
# output 2: s - eigenvalues
#
function get_eig(Q, d_sq)

    n = size(Q, 1);

    U,S,V = svd(Q); # U and S contains eigenvectors and eigenvalues repectively
                    # which is arranged in descending power.

    v = zeros(n,n);

     for ii = 1 : n
        V[ii,:] = V[ii,:]/d_sq[ii];
    end

    for ii = 1 : n
        v[:,ii] = V[:,ii]/norm(V[:,ii]);
    end

    return v, S;
end

c = 5 ; # the desired reduced dimension

v, s = get_eig(Q , d_sq);
p1 = scatter(s[1:10], label="eigenvalues 1:10");
p2 = plot(log.(s), label="eigenvalues in log");
plot(p1, p2, layout=2)

function get_Y(v, S, c)

    n = size(v,1);
    Y = zeros(n,c);

    # get components of diffusion map Y
    for ii = 1:c
        Y[:,ii] = v[:,ii+1].*S[ii+1];
    end

    return Y ;
end

Y = get_Y(v, s, c);

# print diffution map
p1 = scatter(Y[:,1], Y[:,3], label="2D", aspect_ratio=:equal)
p2 = scatter(Y[:,1], Y[:,2], Y[:,3], label="3D", aspect_ratio=:equal)
plot(p1, p2, layout=2)


text:
@everywhere include("load_network.jl")

@everywhere using Graphs

@everywhere function create_graph(start_node, end_node)
    @assert Base.length(start_node)==Base.length(end_node)

    no_node = max(maximum(start_node), maximum(end_node))
    no_arc = Base.length(start_node)

    graph = simple_inclist(no_node)
    for i=1:no_arc
        add_edge!(graph, start_node[i], end_node[i])
    end
    return graph
end

@everywhere function get_vector(state, origin, destination, link_dic)
    current = destination
    parent = -1
    x = zeros(Int, maximum(link_dic))

    while parent != origin
        parent = state.parents[current]

        link_idx = link_dic[parent,current]

        if link_idx != 0
            x[link_idx] = 1
        end

        current = parent
    end

    return x
end

ta_data = load_ta_network("Sioux Falls")

# unpacking data from ta_data
network_name = ta_data.network_name

number_of_zones = ta_data.number_of_zones
number_of_nodes = ta_data.number_of_nodes
first_thru_node = ta_data.first_thru_node
number_of_links = ta_data.number_of_links

start_node = ta_data.start_node
end_node = ta_data.end_node
capacity = ta_data.capacity
link_length = ta_data.link_length

free_flow_time = ta_data.free_flow_time
B = ta_data.B
power = ta_data.power
speed_limit = ta_data.speed_limit
toll = ta_data.toll
link_type = ta_data.link_type
number_of_zones = ta_data.number_of_zones
total_od_flow = ta_data.total_od_flow
travel_demand = ta_data.travel_demand
od_pairs = ta_data.od_pairs

toll_factor = ta_data.toll_factor
distance_factor = ta_data.distance_factor

best_objective = ta_data.best_objective

# preparing a graph
graph = create_graph(start_node, end_node)
link_dic = sparse(start_node, end_node, 1:number_of_links)

@everywhere function BPR(x)
    # travel_time = free_flow_time .* ( 1.0 + B .* (x./capacity).^power )
    # generalized_cost = travel_time + toll_factor *toll + distance_factor * link_length
    # return generalized_cost

    bpr = similar(x)
    for i=1:length(bpr)
        bpr[i] = free_flow_time[i] * ( 1.0 + B[i] * (x[i]/capacity[i])^power[i] )
        bpr[i] += toll_factor *toll[i] + distance_factor * link_length[i]
    end
    return bpr
end

@everywhere function objective(x)
    # value = free_flow_time .* ( x + B.* ( x.^(power+1)) ./ (capacity.^power) ./ (power+1))
    # return sum(value)

    sum = 0
    for i=1:length(x)
        sum += free_flow_time[i] * ( x[i] + B[i]* ( x[i]^(power[i]+1)) / (capacity[i]^power[i]) / (power[i]+1))
        sum += toll_factor *toll[i] + distance_factor * link_length[i]
    end
    return sum
end

@everywhere function all_or_nothing_single(travel_time)
    state = []
    path = []
    x = zeros(size(start_node))

    for r=1:size(travel_demand)[1]
        # for each origin node r, find shortest paths to all destination nodes
        state = dijkstra_shortest_paths(graph, travel_time, r)

        for s=1:size(travel_demand)[2]
            # for each destination node s, find the shortest-path vector
            # load travel demand
            x = x + travel_demand[r,s] * get_vector(state, r, s, link_dic)
        end
    end

    return x
end

# # parallel computing version
@everywhere function all_or_nothing_parallel(travel_time)
    state = []
    path = []
    vv = zeros(size(start_node))
    x = zeros(size(start_node))

    x = x + @parallel (+) for r=1:size(travel_demand)[1]
        # for each origin node r, find shortest paths to all destination nodes
        # if there is any travel demand starting from node r.
        vv = zeros(size(start_node))

        if sum(travel_demand, 2)[r] > 0.0
            state = dijkstra_shortest_paths(graph, travel_time, r)

            for s=1:size(travel_demand)[2]
                # for each destination node s, find the shortest-path vector
                # v = get_vector(state, r, s, start_node, end_node)

                if travel_demand[r,s] > 0.0
                    # load travel demand
                    vv = vv + travel_demand[r,s] * get_vector(state, r, s, link_dic)
                end
            end

        end

        vv
    end

    return x
end

@everywhere function all_or_nothing(travel_time)
    if nprocs() > 1 # if multiple CPU processes are available
        all_or_nothing_parallel(travel_time)
    else
        all_or_nothing_single(travel_time)
        # when nprocs()==1, using @parallel just adds unnecessary setup time. I guess.
    end
end

# Finding a starting feasible solution
travel_time = BPR(zeros(number_of_links))
xl = all_or_nothing(travel_time)

max_iter_no = 1e10
l = 1
average_excess_cost = 1
tol = 1e-7

while l < max_iter_no
    # Finding yl
    travel_time = BPR(xl)
    yl = all_or_nothing(travel_time)
    xl = xl + (yl - xl)/ (l + 1)
    l = l + 1

    # Average Excess Cost
    average_excess_cost = ( dot(xl, travel_time) - dot(yl, travel_time) ) / dot(xl, travel_time)

    if average_excess_cost < tol
        break
    end
end

print(xl)

max_stars_repo_name: jingzbu/InverseVITraffic

max_issues_repo_name: jingzbu/InverseVITraffic

max_forks_repo_name: jingzbu/InverseVITraffic


text:
module TweetsController

using App, SearchLight, JSON, Toolbox
@dependencies

function index()
  search_term = haskey(@params, :q) ? lowercase(@params(:q)) : ""
  startswith(search_term, "#") && (search_term = search_term[2:end])

  recent_searches = SearchLight.find(Search, SQLQuery(order = "id DESC", limit = 10))

  tweets, polarity, subjectivity, favorites_count, retweets_count = if isempty(search_term)
    Tweets.EMPTY_DATA
  else
    filters = ["polarity != -100", SQLWhereExpression("text ILIKE ?", "%#$(search_term)%")]
    tweets = SearchLight.find(Tweet, SQLQuery(where = filters, limit = 20, order = "id DESC"))

    search = Search(search = search_term)
    SearchLight.create_or_update_by!!(search, :search)

    isempty(tweets) ? (Tweets.EMPTY_DATA) : (tweets, Tweets.tweets_aggregates(filters)..., Tweets.tweets_stats(tweets)...)
  end

  respond_with_html(:tweets, :index,  search_term = search_term, tweets = tweets, recent_searches = recent_searches,
                                      polarity = JSON.json(polarity), subjectivity = JSON.json(subjectivity),
                                      favorites_count = favorites_count, retweets_count = retweets_count)
end

end

max_stars_repo_name: essenciary/genie-demo-tweet-stats

max_issues_repo_name: essenciary/genie-demo-tweet-stats

max_forks_repo_name: essenciary/genie-demo-tweet-stats


text:
#######
## Winston + Lasem to display LaTeX labels
#######
using Lasem

using Winston
using LaTeXStrings

# cache of file names, which were converted earlier
const fcache = Dict{Uint64,ASCIIString}()

# call pdflatex and pdf2svg to render formula
function latex2svg(str::LaTeXString)

    h = hash(bytestring(str))

    if !haskey(fcache, h)

        filename = tempname()   # get temporary file name

        fs = open("$(filename).tex", "w")
        println(fs, "\\documentclass[12pt]{standalone}")
        println(fs, "\\begin{document}")
        println(fs, bytestring(str))
        println(fs, "\\end{document}")

        close(fs)
        try
            success(`pdflatex -output-directory=$(dirname(filename)) $(filename).tex`)
            success(`pdf2svg $(filename).pdf $(filename).svg`)
            rm("$(filename).tex")
            rm("$(filename).aux")
            rm("$(filename).log")
            rm("$(filename).pdf")

            fcache[h] = filename
        catch
            println("Error creating SVG.")
            rethrow()
        end
    else
        filename = fcache[h]
    end

    return filename*".svg"
end

# Winston text drawing

Winston.textdraw(self::Winston.CairoRenderer, x::Real, y::Real, str::ByteString; kwargs...) =
    Winston.textdraw(self::Winston.CairoRenderer, x::Real, y::Real, latexstring(str); kwargs...)

Winston.textwidth(self::Winston.CairoRenderer, str::ByteString) =
    Winston.textwidth(self::Winston.CairoRenderer, latexstring(str))

Winston.textheight(self::Winston.CairoRenderer, str::ByteString) =
    Winston.textheight(self::Winston.CairoRenderer, latexstring(str))

function Winston.textdraw(self::Winston.CairoRenderer, x::Real, y::Real, str::LaTeXString;
    halign::String = "left", valign::String = "bottom", angle::Real = 0, markup::Bool=false)
# original Winston code
#    return Cairo.text(self.ctx, x, y, set_latex(self.ctx, str, get(self,:fontsize)); markup=true, kwargs...)

    println(str)
    println("align: $halign, $valign")

    # latex string --> svg
    filename = latex2svg(str)

    doc = Lasem.lsm_dom_document_new_from_path(filename)
    doc != C_NULL || error("Processing of $str failed.")

    view = Lasem.lsm_dom_document_create_view(doc)
    view != C_NULL || error("Creating lasem view failed.")

#    Lasem.lsm_dom_view_set_resolution(view, 96.0)
    w, h = Lasem.lsm_dom_view_get_size(view)
    println("view size: $w, $h")
    wv = iceil(w)
    hv = iceil(h)

#    wv, hv = Lasem.lsm_dom_view_get_size_pixels(view)

    # Cairo context
    ctx = self.ctx

    # create Cairo surface and render image into it
    surface = surface_create_similar(ctx.surface, wv, hv)
    Lasem.lsm_dom_view_render(view, CairoContext(surface))
    println("ctx size: $(ctx.surface.width), $(ctx.surface.height)")

    save(ctx)

    # alignment offset
    dxrel = -Cairo.align2offset(halign)
    dyrel = -Cairo.align2offset(valign)

    # create pattern from surface and draw it
    pat = CairoPattern(surface)
    translate(ctx, x, y + hv)

#    set_source_rgba(ctx, 1, 0.1, 0.1, 0.6)
#    arc(ctx, 0., 0., 2.0, 0, 2*pi)
#    fill(ctx)

    rotate(ctx, angle*pi/180.)

#    rectangle(ctx, 0., 0., wv, hv)
#    stroke(ctx)

    translate(ctx, dxrel*wv, -dyrel*hv )
    ctm = get_matrix(ctx)

#    rectangle(ctx, 0., 0., abs(wv), hv)
#    stroke(ctx)

    println("cm: $ctm")
    println("sctx size: $(wv), $(hv)")
    println("pos: $(x + dxrel*wv), $(y -dyrel*hv)")
    println("pos: $(x + 0*dxrel*wv), $(y -0*dyrel*hv)")

    set_source(ctx, pat)
    set_matrix(pat, CairoMatrix(1., 0., 0., -1., 0., -hv))

    paint(ctx)
    restore(ctx)

    Lasem.destroy(view)
    Lasem.destroy(doc)

    ext = Float64[wv,hv]
    #ext = device_to_user_distance!(sctx, ext)
    ext = inv([ctm.xx ctm.xy; ctm.yx ctm.yy])*ext
    w = (ext[1])
    h = abs(ext[2])

    w, h = Base.Graphics.device_to_user(ctx, wv, hv)
    println("size: $w, $h\n")

    BoundingBox(x+dxrel*w, x+(dxrel+1)*w, y-dyrel*h, y+(1-dyrel)*h)
end

function Winston.textwidth(self::Winston.CairoRenderer, str::LaTeXString)
# original Winston code
#    layout_text(self, str)
#    extents = get_layout_size(self.ctx)
#    extents[1]
#
    doc = Lasem.lsm_mathml_document_new_from_itex(bytestring(str))
    doc != C_NULL || error("Processing of $str failed.")

    view = Lasem.lsm_dom_document_create_view(doc)
    view != C_NULL || error("Creating lasem view failed.")

#    Lasem.lsm_dom_view_set_resolution(view, 96.0)
#    w,h = Lasem.lsm_dom_view_get_size_pixels(view)
    w, h = Lasem.lsm_dom_view_get_size(view)

    Lasem.destroy(view)
    Lasem.destroy(doc)

    return w
end

function Winston.textheight(self::Winston.CairoRenderer, str::LaTeXString)
# original Winston code
#    get(self.state, :fontsize) ## XXX: kludge?
#
    doc = Lasem.lsm_mathml_document_new_from_itex(bytestring(str))
    doc != C_NULL || error("Processing of $str failed.")

    view = Lasem.lsm_dom_document_create_view(doc)
    view != C_NULL || error("Creating lasem view failed.")

#    Lasem.lsm_dom_view_set_resolution(view, 96.0)
#    w,h = Lasem.lsm_dom_view_get_size_pixels(view)
    w, h = Lasem.lsm_dom_view_get_size(view)

    Lasem.destroy(view)
    Lasem.destroy(doc)

    return h
end


text:
function arraydeep(a::Array)
    deep = 1
    while eltype(a) <: Array
        deep += 1
        a = a[1]
    end
    deep
end


text:
function main!(r::Robot)
    t = 0
    k = 0
    a = 0
    while !isborder(r,Ost)
        while !isborder(r, Nord)
            move!(r,Nord)
            if ismarker(r)
                k = k + 1
                t = t + temperature(r)
            end
        end
        goback!(r)
        move!(r, Ost)
        if ismarker(r)
            k = k + 1
            t = t + temperature(r)
        end
    end
    while !isborder(r, Nord)
        move!(r, Nord)
        if ismarker(r)
            k = k + 1
            t = t + temperature(r)
        end
    end
    if k > 0
        a = t/k
        println(a)
    end
end


text:
using ModelParameters, InteractModels
using Documenter

makedocs(;
    modules=[ModelParameters, InteractModels],
    authors="Rafael Schouten <rafaelschouten@gmail.com> and contributors",
    repo="https://github.com/rafaqz/ModelParameters.jl/blob/{commit}{path}#L{line}",
    sitename="ModelParameters.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://rafaqz.github.io/ModelParameters.jl",
        assets=String[],
    ),
    pages=[
        "ModelParameters" => "index.md",
        "InteractModels" => "interactmodels.md",
    ],
)

deploydocs(;
    repo="github.com/rafaqz/ModelParameters.jl",
)


text:
#=
---------------------------------------------------------------------
This function computes the Metric Terms for General 3D Hexahedral Grids.

The strategy follows Algorithm 12.4 and 12.5 in F.X. Giraldo's Introduction to Element-based
Galerkin Methods using Tensor-Product Bases: Analysis, Algorithms, and Applications.
Here, we use the cross-product metric terms. As an exercise the student may be asked to build the
Curl-Invariant metric terms.

Written by F.X. Giraldo on 7/2021
           Department of Applied Mathematics
           Naval Postgraduate School
           Monterey; CA 93943-5216
---------------------------------------------------------------------
=#

include("map_deriv.jl")

function compute_metrics(coord,intma,ψ,dψ,Ne,Np,Nq,DFloat)

    #Initialize Global Arrays
    ξ_x=zeros(DFloat,Nq,Nq,Nq,Ne)
    ξ_y=zeros(DFloat,Nq,Nq,Nq,Ne)
    ξ_z=zeros(DFloat,Nq,Nq,Nq,Ne)
    η_x=zeros(DFloat,Nq,Nq,Nq,Ne)
    η_y=zeros(DFloat,Nq,Nq,Nq,Ne)
    η_z=zeros(DFloat,Nq,Nq,Nq,Ne)
    ζ_x=zeros(DFloat,Nq,Nq,Nq,Ne)
    ζ_y=zeros(DFloat,Nq,Nq,Nq,Ne)
    ζ_z=zeros(DFloat,Nq,Nq,Nq,Ne)
    jac=zeros(DFloat,Nq,Nq,Nq,Ne)

    #Initialize Local Arrays
    x=zeros(DFloat,Np,Np,Np)
    y=zeros(DFloat,Np,Np,Np)
    z=zeros(DFloat,Np,Np,Np)

    #loop thru the elements
    for e=1:Ne

        #Store Element Variables
        for k=1:Np, j=1:Np, i=1:Np
            I=intma[i,j,k,e]
            x[i,j,k]=coord[1,I]
            y[i,j,k]=coord[2,I]
            z[i,j,k]=coord[3,I]
        end #j

        #Construct Mapping Derivatives: dx/dξ; dx/dη; dy/dξ; dy/dη
        (x_ξ,x_η,x_ζ)=map_deriv(ψ,dψ,x,Np,Nq,DFloat)
        (y_ξ,y_η,y_ζ)=map_deriv(ψ,dψ,y,Np,Nq,DFloat)
        (z_ξ,z_η,z_ζ)=map_deriv(ψ,dψ,z,Np,Nq,DFloat)

        #Construct Inverse Mapping: dξ/dx; dξ/dy; dη/dx; dη/dy
        for k=1:Nq, j=1:Nq, i=1:Nq
            xj = (x_ξ[i,j,k]*y_η[i,j,k]*z_ζ[i,j,k] - x_ξ[i,j,k]*y_ζ[i,j,k]*z_η[i,j,k]) -
            (y_ξ[i,j,k]*x_η[i,j,k]*z_ζ[i,j,k] - y_ξ[i,j,k]*x_ζ[i,j,k]*z_η[i,j,k]) +
            (z_ξ[i,j,k]*x_η[i,j,k]*y_ζ[i,j,k] - z_ξ[i,j,k]*x_ζ[i,j,k]*y_η[i,j,k])

            ξ_x[i,j,k,e]= (y_η[i,j,k]*z_ζ[i,j,k]-y_ζ[i,j,k]*z_η[i,j,k])/xj
            ξ_y[i,j,k,e]=-(x_η[i,j,k]*z_ζ[i,j,k]-x_ζ[i,j,k]*z_η[i,j,k])/xj
            ξ_z[i,j,k,e]= (x_η[i,j,k]*y_ζ[i,j,k]-x_ζ[i,j,k]*y_η[i,j,k])/xj
            η_x[i,j,k,e]=-(y_ξ[i,j,k]*z_ζ[i,j,k]-y_ζ[i,j,k]*z_ξ[i,j,k])/xj
            η_y[i,j,k,e]= (x_ξ[i,j,k]*z_ζ[i,j,k]-x_ζ[i,j,k]*z_ξ[i,j,k])/xj
            η_z[i,j,k,e]= -(x_ξ[i,j,k]*y_ζ[i,j,k]-x_ζ[i,j,k]*y_ξ[i,j,k])/xj
            ζ_x[i,j,k,e]= (y_ξ[i,j,k]*z_η[i,j,k] -y_η[i,j,k]*z_ξ[i,j,k])/xj
            ζ_y[i,j,k,e]=-(x_ξ[i,j,k]*z_η[i,j,k] -x_η[i,j,k]*z_ξ[i,j,k])/xj
            ζ_z[i,j,k,e]= (x_ξ[i,j,k]*y_η[i,j,k] -x_η[i,j,k]*y_ξ[i,j,k])/xj
            jac[i,j,k,e]=abs(xj)
        end #i, j, k
    end #e

    return (ξ_x,ξ_y,ξ_z,η_x,η_y,η_z,ζ_x,ζ_y,ζ_z,jac)
end


text:
using Plots, CSV, JuMP, Gurobi, Printf, DataFrames, Tables,

### This file reads prices form Balmorel iteraions, driving patterns and heating demands and genrates optimal load scheduling
### the five models are implemented in separate block in order to run and debug stand alone. Results are saved in csv format and
###  read by Pyhton for aggregation as in as in section 7.2.3 - step 1, following project report.

## Change iteration number
iter="5iter"         ### choose #it
Iteration="5"        ### choose #it
scen="Expected"      ### choose between Expected, RapidScenario and SlowScenario
adopt="0.25"         ### choose
path="GitHub/Julia/"

# 4 Weeks optimization
Regions=collect(1:2)
Years=collect(1:4)
CategoriesAA=collect(1:24)
CategoriesSumHouse=collect(1:24)
CategoriesAW=collect(1:12)
Periods=collect(1:672)
yrs=["2020","2025","2030","2040"]
## Prices
DK1=CSV.read(path*"prices/pricesDK1"*scen*adopt*iter*".csv", delim=",")
DK2=CSV.read(path*"prices/pricesDK2"*scen*adopt*iter*".csv", delim=",")
DK1=dropmissing(DK1)
DK2=dropmissing(DK2)
Prices=[DK1,DK2]
λ=zeros(size(Years)[1],size(Regions)[1],size(Periods)[1])
for r in Regions
    for y in Years
        for t in Periods
            λ[y,r,t]=Prices[r][t,y]
        end
    end
end

#plot(DK1[:,2])
#plot(λ[1,1,:])

## HEATING FLEX1  A-A (Air-to-Air)  Flex=1h########
AAnoFlex=CSV.read(path*"DemAA.csv", delim=",")
AAnoFlex=dropmissing(AAnoFlex)
demandAA=vcat(AAnoFlex[1008:1175,:],AAnoFlex[2016:2183,:],AAnoFlex[4032:4199,:],AAnoFlex[6720:6887,:])
PmaxAA=Pmax=ones(size(Periods)[1],size(CategoriesAA)[1],size(Years)[1])*1
COP_AA=2.8
ηgas=0.98
Ps=demandAA[:,2:25]
λgas=100000

OPT=zeros(size(Years)[1],size(Regions)[1])
for y in Years
    for r in Regions
       finalAA=zeros(size(Periods)[1],size(CategoriesAA)[1])
       POWER=zeros(size(Periods)[1],size(CategoriesAA)[1])
       for c in CategoriesAA
          model_flexF1 = Model(solver = GurobiSolver())
          @variable(model_flexF1, PfAA[t=1:672]>=0)   #load consumption at time t for category i
          @variable(model_flexF1, Pgas[t=1:672]>=0)   #load consumption at time t for category i
          @objective(model_flexF1, Min, sum(λ[y,r,t]*PfAA[t]+λgas*Pgas[t] for t in Periods))
          @constraint(model_flexF1, [t in collect(1:668)], PfAW[t]+PfAW[t+1]>=Ps[t,c]+Ps[t+1,c])
          @constraint(model_flexF1, [t in Periods], PfAW[t]<= PmaxAW[t,c,y]) #Pmax
          @constraint(model_flexF1, sum(PfAW[t] for t in collect(1:672))==sum(Ps[t,c] for t in collect(1:672)))
          solve(model_flexF1)
          POWER[:,c]=getvalue.(PfAA[:])
       end
       if y == 3 && r==1
          display(plot(POWER[:,1]))
       end
       FINAL=vcat(DataFrame(repeat(POWER[1:168,:],10)),DataFrame(repeat(POWER[169:336,:],10)),DataFrame(repeat(POWER[337:504,:],10)),DataFrame(repeat(POWER[505:672,:],22)))
       name=string(path*"Iteration"*Iteration*"/FlexEH/FlexEHAA_DK", r ,"_y", yrs[y] ,".csv")
       CSV.write(name, FINAL, writeheader=false)
    end
end

## HEATING FLEX2  A-A (Air-to-Air) SummerHouse  Flex=1h########
SumHousenoFlex=CSV.read(path*"DemSumHouse.csv", delim=",")
SumHousenoFlex=dropmissing(SumHousenoFlex)
demandSH=vcat(SumHousenoFlex[1008:1175,:],SumHousenoFlex[2016:2183,:],SumHousenoFlex[4032:4199,:],SumHousenoFlex[6720:6887,:])
PmaxSumHouse=Pmax=ones(size(Periods)[1],size(CategoriesSumHouse)[1],size(Years)[1])*1
COP_AA=2.8
ηgas=0.98
Ps=demandSH[:,2:25]
λgas=100000
for y in Years
    for r in Regions
       finalSumHouse=zeros(size(Periods)[1],size(CategoriesSumHouse)[1])
       POWER=zeros(size(Periods)[1],size(CategoriesSumHouse)[1])
       for c in CategoriesSumHouse
          model_flexF2 = Model(solver = GurobiSolver())
          @variable(model_flexF2, PfSumHouse[t=1:672]>=0)   #load consumption at time t for category i
          @variable(model_flexF2, Pgas[t=1:672]>=0)   #load consumption at time t for category i
          @objective(model_flexF2, Min, sum(λ[y,r,t]*PfSumHouse[t]+λgas*Pgas[t] for t in Periods))
          @constraint(model_flexF2, [t in collect(1:668)], PfSumHouse[t]+PfSumHouse[t+1]>=Ps[t,c]+Ps[t+1,c])
          @constraint(model_flexF2, [t in Periods], PfSumHouse[t]<= PmaxSumHouse[t,c,y]) #Pmax
          @constraint(model_flexF2, sum(PfSumHouse[t] for t in collect(1:672))==sum(Ps[t,c] for t in collect(1:672)))
          solve(model_flexF2)
          POWER[:,c]=getvalue.(PfSumHouse[:])
       end
       if y == 1 && r==1
          display(plot(POWER[:,1]))
       end
       FINAL=vcat(DataFrame(repeat(POWER[1:168,:],10)),DataFrame(repeat(POWER[169:336,:],10)),DataFrame(repeat(POWER[337:504,:],10)),DataFrame(repeat(POWER[505:672,:],22)))
       name=string(path*"Iteration"*Iteration*"/FlexEH/FlexEHSumHouse_DK", r ,"_y", yrs[y] ,".csv")
       CSV.write(name,  FINAL, writeheader=false)
      end
end

#####
## HEATING FLEX 3 Air-to-Water with storage flex= 3h############
AWnoFlex=CSV.read(path*"DemAW.csv", delim=",")
AWnoFlex=dropmissing(AWnoFlex)
demandAA=vcat(AWnoFlex[1008:1175,:],AWnoFlex[2016:2183,:],AWnoFlex[4032:4199,:],AWnoFlex[6720:6887,:])
PmaxAW=ones(size(Periods)[1],size(CategoriesAW)[1],size(Years)[1])*3
COP_AW=2.8
Ps=demandAA[:,2:13]
final=zeros(size(Periods)[1],size(CategoriesAW)[1])
for y in Years
     for r in Regions
       POWER=zeros(size(Periods)[1],size(CategoriesSumHouse)[1])
       MISSING=zeros(size(Periods)[1],size(CategoriesSumHouse)[1])
       for c in CategoriesAW
       model_flexF3 = Model(solver = GurobiSolver())
       @variable(model_flexF3, PfAW[t=1:672]>=0)   #load consumption at time t for category i
       @objective(model_flexF3, Min, sum(λ[y,r,t]*PfAW[t] for t in Periods))
       @constraint(model_flexF3, [t in collect(1:668)], PfAW[t]+PfAW[t+1]+PfAW[t+2]+PfAW[t+3]>=Ps[t,c]+Ps[t+1,c]+Ps[t+2,c]+Ps[t+3,c])
       @constraint(model_flexF3, [t in Periods], PfAW[t]<= PmaxAW[t,c,y]) #Pmax
       @constraint(model_flexF3, sum(PfAW[t] for t in collect(1:672))==sum(Ps[t,c] for t in collect(1:672)))
       solve(model_flexF3)
       POWER[:,c]=getvalue.(PfAW[:])
       end
       if y == 1 && r==1
          display(plot(POWER[:,2]))
          display(plot!(Ps[:,1]))
          display(plot(POWER[40:88,4],label="Smart",legend=:topleft, right_margin = 15Plots.mm, ylabel = "Power [kW]",ylims=[0,3.6]))
          display(plot!(twinx(),λ[1,1,40:88],label="Price",color=:green, ylabel = "[€/MWh]", right_margin = 15Plots.mm,ylims=[0,140]))
          display(plot!(Ps[40:88,4],label="Passive",legend=:topleft, right_margin = 15Plots.mm, bottom_margin=5Plots.mm))
          xlabel!("Hours")
          plot!(size=(800,280))
          savefig("heatpump.png")
          plot!(label=["Smart","Passive","Price"])
       end
       FINAL=vcat(DataFrame(repeat(POWER[1:168,:],10)),DataFrame(repeat(POWER[169:336,:],10)),DataFrame(repeat(POWER[337:504,:],10)),DataFrame(repeat(POWER[505:672,:],22)))
       name=string(path*"Iteration"*Iteration*"/FlexEH/FlexEHAW_DK", r ,"_y", yrs[y] ,".csv")
       CSV.write(name,  DataFrame(FINAL), writeheader=false)
end
end

## ELECTRIC VEHICLES total
using Plots, CSV, JuMP, Gurobi, Printf, DataFrames, Tables

Regions=collect(1:2)
Years=collect(1:4)
Periods=collect(1:672)
CategoriesAA=collect(1:24)
CategoriesSumHouse=collect(1:24)
CategoriesAW=collect(1:12)
Days=collect(1:2917)*3
WEEKS=[6,12,24,40]*168

## Change iteration number
#iter="2"
path="C:/Users/pietr/Documents/GitHub/Tesi/Julia/"## Prices
DK1=CSV.read(path*"prices/pricesweeklyDK1"*scen*adopt*iter*".csv", delim=",")
DK2=CSV.read(path*"prices/pricesweeklyDK2"*scen*adopt*iter*".csv", delim=",")
DK1=dropmissing(DK1)
DK2=dropmissing(DK2)
Prices=[DK1,DK2]
λ=zeros(size(Years)[1],size(Regions)[1],size(Periods)[1])
for r in Regions
    for y in Years
        for t in Periods
            λ[y,r,t]=Prices[r][t,y]
        end
    end
end

#############################
########### BEV #############
#############################

BEVtypes=collect(1:20)
Availability=CSV.read(path*"Driving/Driving_patterns_avail.csv")
KM_demand=CSV.read(path*"Driving/Driving_patterns.csv")
Ex=CSV.read(path*"Driving/Driving_patterns_exit.csv")
ηcar=6     #Km/KWh
ChMax=3.7
fday=[2,2,1,2,6]*24
ηch=0.95
EM=0.9
BsizeBEV=60
M=50000   #BigM

for y in Years
    AvBEVWD=Availability[1:size(Periods)[1],3:22]#+fday[y],3:22]
    KMBEVWD=KM_demand[1:size(Periods)[1],3:22]#+fday[y],3:22]
    Exit=Ex[1:size(Periods)[1],3:22]#+fday[y],3:22]
    for r in Regions
        SOCC=zeros(size(Periods)[1],size(BEVtypes)[1])
        POWER=zeros(size(Periods)[1],size(BEVtypes)[1])
        for e in BEVtypes
            model_EV = Model(solver = GurobiSolver())
            @variable(model_EV, ChEV[t in Periods]>=0)   #load consumption at time t for type e
            @variable(model_EV, 0<=SOC[t in Periods]<=BsizeBEV)
            @objective(model_EV, Min, sum(λ[y,r,t]*ChEV[t] for t in Periods))
            @constraint(model_EV, SOC[1]==SOC[672]+ChEV[672]) ##t in Periods
            @constraint(model_EV, [t =2:size(Periods)[1]], SOC[t]==SOC[t-1]+ChEV[t]-KMBEVWD[t,e]/ηcar)
            @constraint(model_EV, [t in Periods], ChEV[t]<= AvBEVWD[t,e]*ChMax)
            @constraint(model_EV, [t in Periods], SOC[t]>=Exit[t,e]*BsizeBEV*EM)   #when going out full battery needed
            solve(model_EV)
            POWER[:,e]=getvalue.(ChEV[:])
            SOCC[:,e]=getvalue.(SOC[:])/60
        end
        if y == 1 && r==1
          display(plot(POWER[169:672]))
          display(plot!(λ[y,r,169:672]))
       end
       FINAL=vcat(DataFrame(repeat(POWER[1:168,:],10)),DataFrame(repeat(POWER[169:336,:],10)),DataFrame(repeat(POWER[337:504,:],10)),DataFrame(repeat(POWER[505:672,:],22)))
#       FINAL=vcat(DataFrame(POWER[1:168,:]),DataFrame(POWER[169:336,:]),DataFrame(POWER[337:504,:]),DataFrame(POWER[505:672,:]))
       name=string(path*"Iteration"*Iteration*"/FlexEV/FlexEV_DK", r ,"_y", yrs[y] ,".csv")
       CSV.write(name,  DataFrame(FINAL), writeheader=false)
    end
end

Ps=CSV.read(path*"FLEXNOBEV.csv")
display(plot(POWER[96:144,4],label="Smart",legend=:topleft, right_margin = 15Plots.mm, ylabel = "Power [kW]",ylims=[0,4.9]))
display(plot!(Ps[97:145,4],label="Passive",legend=:topleft, right_margin = 15Plots.mm, bottom_margin=5Plots.mm))
display(plot!(twinx(),λ[1,1,96:144],label="Price",color=:green,legend=:topleft, ylabel = "[€/MWh]", right_margin = 15Plots.mm,ylims=[20,55]))
vline!([18,37],color=:black,label="Arrive")
vline!([10,30],color=:red,label="Departure",legend=:topright)
xlabel!("Hours")
plot!(size=(800,280))
savefig("EV.png")

#############################
########## PHEV #############
#############################
PHEVtypes=collect(1:10)
Availability=CSV.read(path*"Driving/Driving_patterns_avail_PHEV.csv")
KM_demand=CSV.read(path*"Driving/Driving_patterns_PHEV.csv")
Ex=CSV.read(path*"Driving/Driving_patterns_exit_PHEV.csv")
ηcar=5     #Km/KWh
ChMax=3.7
fday=[2,2,1,2,6]*24
ηch=0.95
BsizePHEV=10
M=50000   #BigM

for y in Years
    AvPHEV=Availability[1:size(Periods)[1],3:12]#+fday[y],3:12]
    KMPHEV=KM_demand[1:size(Periods)[1],3:12]#+fday[y],3:12]
    ExitPHEV=Ex[1:size(Periods)[1],3:12]#+fday[y],3:12]
    for r in Regions
        SOCC=zeros(size(Periods)[1],size(PHEVtypes)[1])
        POWERPHEV=zeros(size(Periods)[1],size(PHEVtypes)[1])
        for e in PHEVtypes
            model_PHEV = Model(solver = GurobiSolver())
            @variable(model_PHEV, ChEV[t in Periods]>=0)   #load consumption at time t for type e
            @variable(model_PHEV, 0<=SOC[t in Periods]<=BsizePHEV)
            @variable(model_PHEV, Oil[t in Periods,e in PHEVtypes]>=0)
            @objective(model_PHEV, Min, sum(λ[y,r,t]*ChEV[t]+100000*Oil[t,e] for t in Periods))
            @constraint(model_PHEV, SOC[1]==SOC[672]+ChEV[672]) ##t in Periods
            @constraint(model_PHEV, [t =2:size(Periods)[1]], SOC[t]==SOC[t-1]+Oil[t,e]+ChEV[t]-KMPHEV[t,e]/ηcar)
            @constraint(model_PHEV, [t in Periods], ChEV[t]<= AvPHEV[t,e]*ChMax)
            @constraint(model_PHEV, [t in Periods], SOC[t]>=ExitPHEV[t,e]*BsizePHEV)   #when going out full battery needed
            solve(model_PHEV)
            POWERPHEV[:,e]=getvalue.(ChEV[:])
            SOCC[:,e]=getvalue.(SOC[:])/10
        end
        if r==1
          display(plot(POWERPHEV[:,1:3]))
        end
        FINAL=vcat(DataFrame(repeat(POWERPHEV[1:168,:],10)),DataFrame(repeat(POWERPHEV[169:336,:],10)),DataFrame(repeat(POWERPHEV[337:504,:],10)),DataFrame(repeat(POWERPHEV[505:672,:],22)))
#        FINAL=vcat(DataFrame(POWERPHEV[1:168,:]),DataFrame(POWERPHEV[169:336,:]),DataFrame(POWERPHEV[337:504,:]),DataFrame(POWERPHEV[505:672,:]))
        name=string(path*"Iteration"*Iteration*"/FlexEV/FlexPHEV_DK", r ,"_y", yrs[y] ,".csv")
        CSV.write(name,  DataFrame(FINAL), writeheader=false)
    end
end